![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zsptrs.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 zsptrs.f CVS log | Up to [local] / rpl / lapack / lapack |
1.1 bertrand 1: SUBROUTINE ZSPTRS( UPLO, N, NRHS, AP, IPIV, B, LDB, INFO )
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: CHARACTER UPLO
10: INTEGER INFO, LDB, N, NRHS
11: * ..
12: * .. Array Arguments ..
13: INTEGER IPIV( * )
14: COMPLEX*16 AP( * ), B( LDB, * )
15: * ..
16: *
17: * Purpose
18: * =======
19: *
20: * ZSPTRS solves a system of linear equations A*X = B with a complex
21: * symmetric matrix A stored in packed format using the factorization
22: * A = U*D*U**T or A = L*D*L**T computed by ZSPTRF.
23: *
24: * Arguments
25: * =========
26: *
27: * UPLO (input) CHARACTER*1
28: * Specifies whether the details of the factorization are stored
29: * as an upper or lower triangular matrix.
30: * = 'U': Upper triangular, form is A = U*D*U**T;
31: * = 'L': Lower triangular, form is A = L*D*L**T.
32: *
33: * N (input) INTEGER
34: * The order of the matrix A. N >= 0.
35: *
36: * NRHS (input) INTEGER
37: * The number of right hand sides, i.e., the number of columns
38: * of the matrix B. NRHS >= 0.
39: *
40: * AP (input) COMPLEX*16 array, dimension (N*(N+1)/2)
41: * The block diagonal matrix D and the multipliers used to
42: * obtain the factor U or L as computed by ZSPTRF, stored as a
43: * packed triangular matrix.
44: *
45: * IPIV (input) INTEGER array, dimension (N)
46: * Details of the interchanges and the block structure of D
47: * as determined by ZSPTRF.
48: *
49: * B (input/output) COMPLEX*16 array, dimension (LDB,NRHS)
50: * On entry, the right hand side matrix B.
51: * On exit, the solution matrix X.
52: *
53: * LDB (input) INTEGER
54: * The leading dimension of the array B. LDB >= max(1,N).
55: *
56: * INFO (output) INTEGER
57: * = 0: successful exit
58: * < 0: if INFO = -i, the i-th argument had an illegal value
59: *
60: * =====================================================================
61: *
62: * .. Parameters ..
63: COMPLEX*16 ONE
64: PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) )
65: * ..
66: * .. Local Scalars ..
67: LOGICAL UPPER
68: INTEGER J, K, KC, KP
69: COMPLEX*16 AK, AKM1, AKM1K, BK, BKM1, DENOM
70: * ..
71: * .. External Functions ..
72: LOGICAL LSAME
73: EXTERNAL LSAME
74: * ..
75: * .. External Subroutines ..
76: EXTERNAL XERBLA, ZGEMV, ZGERU, ZSCAL, ZSWAP
77: * ..
78: * .. Intrinsic Functions ..
79: INTRINSIC MAX
80: * ..
81: * .. Executable Statements ..
82: *
83: INFO = 0
84: UPPER = LSAME( UPLO, 'U' )
85: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
86: INFO = -1
87: ELSE IF( N.LT.0 ) THEN
88: INFO = -2
89: ELSE IF( NRHS.LT.0 ) THEN
90: INFO = -3
91: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
92: INFO = -7
93: END IF
94: IF( INFO.NE.0 ) THEN
95: CALL XERBLA( 'ZSPTRS', -INFO )
96: RETURN
97: END IF
98: *
99: * Quick return if possible
100: *
101: IF( N.EQ.0 .OR. NRHS.EQ.0 )
102: $ RETURN
103: *
104: IF( UPPER ) THEN
105: *
106: * Solve A*X = B, where A = U*D*U'.
107: *
108: * First solve U*D*X = B, overwriting B with X.
109: *
110: * K is the main loop index, decreasing from N to 1 in steps of
111: * 1 or 2, depending on the size of the diagonal blocks.
112: *
113: K = N
114: KC = N*( N+1 ) / 2 + 1
115: 10 CONTINUE
116: *
117: * If K < 1, exit from loop.
118: *
119: IF( K.LT.1 )
120: $ GO TO 30
121: *
122: KC = KC - K
123: IF( IPIV( K ).GT.0 ) THEN
124: *
125: * 1 x 1 diagonal block
126: *
127: * Interchange rows K and IPIV(K).
128: *
129: KP = IPIV( K )
130: IF( KP.NE.K )
131: $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
132: *
133: * Multiply by inv(U(K)), where U(K) is the transformation
134: * stored in column K of A.
135: *
136: CALL ZGERU( K-1, NRHS, -ONE, AP( KC ), 1, B( K, 1 ), LDB,
137: $ B( 1, 1 ), LDB )
138: *
139: * Multiply by the inverse of the diagonal block.
140: *
141: CALL ZSCAL( NRHS, ONE / AP( KC+K-1 ), B( K, 1 ), LDB )
142: K = K - 1
143: ELSE
144: *
145: * 2 x 2 diagonal block
146: *
147: * Interchange rows K-1 and -IPIV(K).
148: *
149: KP = -IPIV( K )
150: IF( KP.NE.K-1 )
151: $ CALL ZSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB )
152: *
153: * Multiply by inv(U(K)), where U(K) is the transformation
154: * stored in columns K-1 and K of A.
155: *
156: CALL ZGERU( K-2, NRHS, -ONE, AP( KC ), 1, B( K, 1 ), LDB,
157: $ B( 1, 1 ), LDB )
158: CALL ZGERU( K-2, NRHS, -ONE, AP( KC-( K-1 ) ), 1,
159: $ B( K-1, 1 ), LDB, B( 1, 1 ), LDB )
160: *
161: * Multiply by the inverse of the diagonal block.
162: *
163: AKM1K = AP( KC+K-2 )
164: AKM1 = AP( KC-1 ) / AKM1K
165: AK = AP( KC+K-1 ) / AKM1K
166: DENOM = AKM1*AK - ONE
167: DO 20 J = 1, NRHS
168: BKM1 = B( K-1, J ) / AKM1K
169: BK = B( K, J ) / AKM1K
170: B( K-1, J ) = ( AK*BKM1-BK ) / DENOM
171: B( K, J ) = ( AKM1*BK-BKM1 ) / DENOM
172: 20 CONTINUE
173: KC = KC - K + 1
174: K = K - 2
175: END IF
176: *
177: GO TO 10
178: 30 CONTINUE
179: *
180: * Next solve U'*X = B, overwriting B with X.
181: *
182: * K is the main loop index, increasing from 1 to N in steps of
183: * 1 or 2, depending on the size of the diagonal blocks.
184: *
185: K = 1
186: KC = 1
187: 40 CONTINUE
188: *
189: * If K > N, exit from loop.
190: *
191: IF( K.GT.N )
192: $ GO TO 50
193: *
194: IF( IPIV( K ).GT.0 ) THEN
195: *
196: * 1 x 1 diagonal block
197: *
198: * Multiply by inv(U'(K)), where U(K) is the transformation
199: * stored in column K of A.
200: *
201: CALL ZGEMV( 'Transpose', K-1, NRHS, -ONE, B, LDB, AP( KC ),
202: $ 1, ONE, B( K, 1 ), LDB )
203: *
204: * Interchange rows K and IPIV(K).
205: *
206: KP = IPIV( K )
207: IF( KP.NE.K )
208: $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
209: KC = KC + K
210: K = K + 1
211: ELSE
212: *
213: * 2 x 2 diagonal block
214: *
215: * Multiply by inv(U'(K+1)), where U(K+1) is the transformation
216: * stored in columns K and K+1 of A.
217: *
218: CALL ZGEMV( 'Transpose', K-1, NRHS, -ONE, B, LDB, AP( KC ),
219: $ 1, ONE, B( K, 1 ), LDB )
220: CALL ZGEMV( 'Transpose', K-1, NRHS, -ONE, B, LDB,
221: $ AP( KC+K ), 1, ONE, B( K+1, 1 ), LDB )
222: *
223: * Interchange rows K and -IPIV(K).
224: *
225: KP = -IPIV( K )
226: IF( KP.NE.K )
227: $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
228: KC = KC + 2*K + 1
229: K = K + 2
230: END IF
231: *
232: GO TO 40
233: 50 CONTINUE
234: *
235: ELSE
236: *
237: * Solve A*X = B, where A = L*D*L'.
238: *
239: * First solve L*D*X = B, overwriting B with X.
240: *
241: * K is the main loop index, increasing from 1 to N in steps of
242: * 1 or 2, depending on the size of the diagonal blocks.
243: *
244: K = 1
245: KC = 1
246: 60 CONTINUE
247: *
248: * If K > N, exit from loop.
249: *
250: IF( K.GT.N )
251: $ GO TO 80
252: *
253: IF( IPIV( K ).GT.0 ) THEN
254: *
255: * 1 x 1 diagonal block
256: *
257: * Interchange rows K and IPIV(K).
258: *
259: KP = IPIV( K )
260: IF( KP.NE.K )
261: $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
262: *
263: * Multiply by inv(L(K)), where L(K) is the transformation
264: * stored in column K of A.
265: *
266: IF( K.LT.N )
267: $ CALL ZGERU( N-K, NRHS, -ONE, AP( KC+1 ), 1, B( K, 1 ),
268: $ LDB, B( K+1, 1 ), LDB )
269: *
270: * Multiply by the inverse of the diagonal block.
271: *
272: CALL ZSCAL( NRHS, ONE / AP( KC ), B( K, 1 ), LDB )
273: KC = KC + N - K + 1
274: K = K + 1
275: ELSE
276: *
277: * 2 x 2 diagonal block
278: *
279: * Interchange rows K+1 and -IPIV(K).
280: *
281: KP = -IPIV( K )
282: IF( KP.NE.K+1 )
283: $ CALL ZSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB )
284: *
285: * Multiply by inv(L(K)), where L(K) is the transformation
286: * stored in columns K and K+1 of A.
287: *
288: IF( K.LT.N-1 ) THEN
289: CALL ZGERU( N-K-1, NRHS, -ONE, AP( KC+2 ), 1, B( K, 1 ),
290: $ LDB, B( K+2, 1 ), LDB )
291: CALL ZGERU( N-K-1, NRHS, -ONE, AP( KC+N-K+2 ), 1,
292: $ B( K+1, 1 ), LDB, B( K+2, 1 ), LDB )
293: END IF
294: *
295: * Multiply by the inverse of the diagonal block.
296: *
297: AKM1K = AP( KC+1 )
298: AKM1 = AP( KC ) / AKM1K
299: AK = AP( KC+N-K+1 ) / AKM1K
300: DENOM = AKM1*AK - ONE
301: DO 70 J = 1, NRHS
302: BKM1 = B( K, J ) / AKM1K
303: BK = B( K+1, J ) / AKM1K
304: B( K, J ) = ( AK*BKM1-BK ) / DENOM
305: B( K+1, J ) = ( AKM1*BK-BKM1 ) / DENOM
306: 70 CONTINUE
307: KC = KC + 2*( N-K ) + 1
308: K = K + 2
309: END IF
310: *
311: GO TO 60
312: 80 CONTINUE
313: *
314: * Next solve L'*X = B, overwriting B with X.
315: *
316: * K is the main loop index, decreasing from N to 1 in steps of
317: * 1 or 2, depending on the size of the diagonal blocks.
318: *
319: K = N
320: KC = N*( N+1 ) / 2 + 1
321: 90 CONTINUE
322: *
323: * If K < 1, exit from loop.
324: *
325: IF( K.LT.1 )
326: $ GO TO 100
327: *
328: KC = KC - ( N-K+1 )
329: IF( IPIV( K ).GT.0 ) THEN
330: *
331: * 1 x 1 diagonal block
332: *
333: * Multiply by inv(L'(K)), where L(K) is the transformation
334: * stored in column K of A.
335: *
336: IF( K.LT.N )
337: $ CALL ZGEMV( 'Transpose', N-K, NRHS, -ONE, B( K+1, 1 ),
338: $ LDB, AP( KC+1 ), 1, ONE, B( K, 1 ), LDB )
339: *
340: * Interchange rows K and IPIV(K).
341: *
342: KP = IPIV( K )
343: IF( KP.NE.K )
344: $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
345: K = K - 1
346: ELSE
347: *
348: * 2 x 2 diagonal block
349: *
350: * Multiply by inv(L'(K-1)), where L(K-1) is the transformation
351: * stored in columns K-1 and K of A.
352: *
353: IF( K.LT.N ) THEN
354: CALL ZGEMV( 'Transpose', N-K, NRHS, -ONE, B( K+1, 1 ),
355: $ LDB, AP( KC+1 ), 1, ONE, B( K, 1 ), LDB )
356: CALL ZGEMV( 'Transpose', N-K, NRHS, -ONE, B( K+1, 1 ),
357: $ LDB, AP( KC-( N-K ) ), 1, ONE, B( K-1, 1 ),
358: $ LDB )
359: END IF
360: *
361: * Interchange rows K and -IPIV(K).
362: *
363: KP = -IPIV( K )
364: IF( KP.NE.K )
365: $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
366: KC = KC - ( N-K+2 )
367: K = K - 2
368: END IF
369: *
370: GO TO 90
371: 100 CONTINUE
372: END IF
373: *
374: RETURN
375: *
376: * End of ZSPTRS
377: *
378: END