![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zptts2.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack|
Return to zptts2.f CVS log | Up to [local] / rpl / lapack / lapack |
1.1 bertrand 1: SUBROUTINE ZPTTS2( IUPLO, N, NRHS, D, E, B, LDB )
2: *
3: * -- LAPACK routine (version 3.2) --
4: * -- LAPACK is a software package provided by Univ. of Tennessee, --
5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6: * November 2006
7: *
8: * .. Scalar Arguments ..
9: INTEGER IUPLO, LDB, N, NRHS
10: * ..
11: * .. Array Arguments ..
12: DOUBLE PRECISION D( * )
13: COMPLEX*16 B( LDB, * ), E( * )
14: * ..
15: *
16: * Purpose
17: * =======
18: *
19: * ZPTTS2 solves a tridiagonal system of the form
20: * A * X = B
21: * using the factorization A = U'*D*U or A = L*D*L' computed by ZPTTRF.
22: * D is a diagonal matrix specified in the vector D, U (or L) is a unit
23: * bidiagonal matrix whose superdiagonal (subdiagonal) is specified in
24: * the vector E, and X and B are N by NRHS matrices.
25: *
26: * Arguments
27: * =========
28: *
29: * IUPLO (input) INTEGER
30: * Specifies the form of the factorization and whether the
31: * vector E is the superdiagonal of the upper bidiagonal factor
32: * U or the subdiagonal of the lower bidiagonal factor L.
33: * = 1: A = U'*D*U, E is the superdiagonal of U
34: * = 0: A = L*D*L', E is the subdiagonal of L
35: *
36: * N (input) INTEGER
37: * The order of the tridiagonal matrix A. N >= 0.
38: *
39: * NRHS (input) INTEGER
40: * The number of right hand sides, i.e., the number of columns
41: * of the matrix B. NRHS >= 0.
42: *
43: * D (input) DOUBLE PRECISION array, dimension (N)
44: * The n diagonal elements of the diagonal matrix D from the
45: * factorization A = U'*D*U or A = L*D*L'.
46: *
47: * E (input) COMPLEX*16 array, dimension (N-1)
48: * If IUPLO = 1, the (n-1) superdiagonal elements of the unit
49: * bidiagonal factor U from the factorization A = U'*D*U.
50: * If IUPLO = 0, the (n-1) subdiagonal elements of the unit
51: * bidiagonal factor L from the factorization A = L*D*L'.
52: *
53: * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
54: * On entry, the right hand side vectors B for the system of
55: * linear equations.
56: * On exit, the solution vectors, X.
57: *
58: * LDB (input) INTEGER
59: * The leading dimension of the array B. LDB >= max(1,N).
60: *
61: * =====================================================================
62: *
63: * .. Local Scalars ..
64: INTEGER I, J
65: * ..
66: * .. External Subroutines ..
67: EXTERNAL ZDSCAL
68: * ..
69: * .. Intrinsic Functions ..
70: INTRINSIC DCONJG
71: * ..
72: * .. Executable Statements ..
73: *
74: * Quick return if possible
75: *
76: IF( N.LE.1 ) THEN
77: IF( N.EQ.1 )
78: $ CALL ZDSCAL( NRHS, 1.D0 / D( 1 ), B, LDB )
79: RETURN
80: END IF
81: *
82: IF( IUPLO.EQ.1 ) THEN
83: *
84: * Solve A * X = B using the factorization A = U'*D*U,
85: * overwriting each right hand side vector with its solution.
86: *
87: IF( NRHS.LE.2 ) THEN
88: J = 1
89: 10 CONTINUE
90: *
91: * Solve U' * x = b.
92: *
93: DO 20 I = 2, N
94: B( I, J ) = B( I, J ) - B( I-1, J )*DCONJG( E( I-1 ) )
95: 20 CONTINUE
96: *
97: * Solve D * U * x = b.
98: *
99: DO 30 I = 1, N
100: B( I, J ) = B( I, J ) / D( I )
101: 30 CONTINUE
102: DO 40 I = N - 1, 1, -1
103: B( I, J ) = B( I, J ) - B( I+1, J )*E( I )
104: 40 CONTINUE
105: IF( J.LT.NRHS ) THEN
106: J = J + 1
107: GO TO 10
108: END IF
109: ELSE
110: DO 70 J = 1, NRHS
111: *
112: * Solve U' * x = b.
113: *
114: DO 50 I = 2, N
115: B( I, J ) = B( I, J ) - B( I-1, J )*DCONJG( E( I-1 ) )
116: 50 CONTINUE
117: *
118: * Solve D * U * x = b.
119: *
120: B( N, J ) = B( N, J ) / D( N )
121: DO 60 I = N - 1, 1, -1
122: B( I, J ) = B( I, J ) / D( I ) - B( I+1, J )*E( I )
123: 60 CONTINUE
124: 70 CONTINUE
125: END IF
126: ELSE
127: *
128: * Solve A * X = B using the factorization A = L*D*L',
129: * overwriting each right hand side vector with its solution.
130: *
131: IF( NRHS.LE.2 ) THEN
132: J = 1
133: 80 CONTINUE
134: *
135: * Solve L * x = b.
136: *
137: DO 90 I = 2, N
138: B( I, J ) = B( I, J ) - B( I-1, J )*E( I-1 )
139: 90 CONTINUE
140: *
141: * Solve D * L' * x = b.
142: *
143: DO 100 I = 1, N
144: B( I, J ) = B( I, J ) / D( I )
145: 100 CONTINUE
146: DO 110 I = N - 1, 1, -1
147: B( I, J ) = B( I, J ) - B( I+1, J )*DCONJG( E( I ) )
148: 110 CONTINUE
149: IF( J.LT.NRHS ) THEN
150: J = J + 1
151: GO TO 80
152: END IF
153: ELSE
154: DO 140 J = 1, NRHS
155: *
156: * Solve L * x = b.
157: *
158: DO 120 I = 2, N
159: B( I, J ) = B( I, J ) - B( I-1, J )*E( I-1 )
160: 120 CONTINUE
161: *
162: * Solve D * L' * x = b.
163: *
164: B( N, J ) = B( N, J ) / D( N )
165: DO 130 I = N - 1, 1, -1
166: B( I, J ) = B( I, J ) / D( I ) -
167: $ B( I+1, J )*DCONJG( E( I ) )
168: 130 CONTINUE
169: 140 CONTINUE
170: END IF
171: END IF
172: *
173: RETURN
174: *
175: * End of ZPTTS2
176: *
177: END