![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zhptrs.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 zhptrs.f CVS log | Up to [local] / rpl / lapack / lapack |
1.1 bertrand 1: SUBROUTINE ZHPTRS( 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: * ZHPTRS solves a system of linear equations A*X = B with a complex
21: * Hermitian matrix A stored in packed format using the factorization
22: * A = U*D*U**H or A = L*D*L**H computed by ZHPTRF.
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**H;
31: * = 'L': Lower triangular, form is A = L*D*L**H.
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 ZHPTRF, 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 ZHPTRF.
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: DOUBLE PRECISION S
70: COMPLEX*16 AK, AKM1, AKM1K, BK, BKM1, DENOM
71: * ..
72: * .. External Functions ..
73: LOGICAL LSAME
74: EXTERNAL LSAME
75: * ..
76: * .. External Subroutines ..
77: EXTERNAL XERBLA, ZDSCAL, ZGEMV, ZGERU, ZLACGV, ZSWAP
78: * ..
79: * .. Intrinsic Functions ..
80: INTRINSIC DBLE, DCONJG, MAX
81: * ..
82: * .. Executable Statements ..
83: *
84: INFO = 0
85: UPPER = LSAME( UPLO, 'U' )
86: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
87: INFO = -1
88: ELSE IF( N.LT.0 ) THEN
89: INFO = -2
90: ELSE IF( NRHS.LT.0 ) THEN
91: INFO = -3
92: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
93: INFO = -7
94: END IF
95: IF( INFO.NE.0 ) THEN
96: CALL XERBLA( 'ZHPTRS', -INFO )
97: RETURN
98: END IF
99: *
100: * Quick return if possible
101: *
102: IF( N.EQ.0 .OR. NRHS.EQ.0 )
103: $ RETURN
104: *
105: IF( UPPER ) THEN
106: *
107: * Solve A*X = B, where A = U*D*U'.
108: *
109: * First solve U*D*X = B, overwriting B with X.
110: *
111: * K is the main loop index, decreasing from N to 1 in steps of
112: * 1 or 2, depending on the size of the diagonal blocks.
113: *
114: K = N
115: KC = N*( N+1 ) / 2 + 1
116: 10 CONTINUE
117: *
118: * If K < 1, exit from loop.
119: *
120: IF( K.LT.1 )
121: $ GO TO 30
122: *
123: KC = KC - K
124: IF( IPIV( K ).GT.0 ) THEN
125: *
126: * 1 x 1 diagonal block
127: *
128: * Interchange rows K and IPIV(K).
129: *
130: KP = IPIV( K )
131: IF( KP.NE.K )
132: $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
133: *
134: * Multiply by inv(U(K)), where U(K) is the transformation
135: * stored in column K of A.
136: *
137: CALL ZGERU( K-1, NRHS, -ONE, AP( KC ), 1, B( K, 1 ), LDB,
138: $ B( 1, 1 ), LDB )
139: *
140: * Multiply by the inverse of the diagonal block.
141: *
142: S = DBLE( ONE ) / DBLE( AP( KC+K-1 ) )
143: CALL ZDSCAL( NRHS, S, B( K, 1 ), LDB )
144: K = K - 1
145: ELSE
146: *
147: * 2 x 2 diagonal block
148: *
149: * Interchange rows K-1 and -IPIV(K).
150: *
151: KP = -IPIV( K )
152: IF( KP.NE.K-1 )
153: $ CALL ZSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB )
154: *
155: * Multiply by inv(U(K)), where U(K) is the transformation
156: * stored in columns K-1 and K of A.
157: *
158: CALL ZGERU( K-2, NRHS, -ONE, AP( KC ), 1, B( K, 1 ), LDB,
159: $ B( 1, 1 ), LDB )
160: CALL ZGERU( K-2, NRHS, -ONE, AP( KC-( K-1 ) ), 1,
161: $ B( K-1, 1 ), LDB, B( 1, 1 ), LDB )
162: *
163: * Multiply by the inverse of the diagonal block.
164: *
165: AKM1K = AP( KC+K-2 )
166: AKM1 = AP( KC-1 ) / AKM1K
167: AK = AP( KC+K-1 ) / DCONJG( AKM1K )
168: DENOM = AKM1*AK - ONE
169: DO 20 J = 1, NRHS
170: BKM1 = B( K-1, J ) / AKM1K
171: BK = B( K, J ) / DCONJG( AKM1K )
172: B( K-1, J ) = ( AK*BKM1-BK ) / DENOM
173: B( K, J ) = ( AKM1*BK-BKM1 ) / DENOM
174: 20 CONTINUE
175: KC = KC - K + 1
176: K = K - 2
177: END IF
178: *
179: GO TO 10
180: 30 CONTINUE
181: *
182: * Next solve U'*X = B, overwriting B with X.
183: *
184: * K is the main loop index, increasing from 1 to N in steps of
185: * 1 or 2, depending on the size of the diagonal blocks.
186: *
187: K = 1
188: KC = 1
189: 40 CONTINUE
190: *
191: * If K > N, exit from loop.
192: *
193: IF( K.GT.N )
194: $ GO TO 50
195: *
196: IF( IPIV( K ).GT.0 ) THEN
197: *
198: * 1 x 1 diagonal block
199: *
200: * Multiply by inv(U'(K)), where U(K) is the transformation
201: * stored in column K of A.
202: *
203: IF( K.GT.1 ) THEN
204: CALL ZLACGV( NRHS, B( K, 1 ), LDB )
205: CALL ZGEMV( 'Conjugate transpose', K-1, NRHS, -ONE, B,
206: $ LDB, AP( KC ), 1, ONE, B( K, 1 ), LDB )
207: CALL ZLACGV( NRHS, B( K, 1 ), LDB )
208: END IF
209: *
210: * Interchange rows K and IPIV(K).
211: *
212: KP = IPIV( K )
213: IF( KP.NE.K )
214: $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
215: KC = KC + K
216: K = K + 1
217: ELSE
218: *
219: * 2 x 2 diagonal block
220: *
221: * Multiply by inv(U'(K+1)), where U(K+1) is the transformation
222: * stored in columns K and K+1 of A.
223: *
224: IF( K.GT.1 ) THEN
225: CALL ZLACGV( NRHS, B( K, 1 ), LDB )
226: CALL ZGEMV( 'Conjugate transpose', K-1, NRHS, -ONE, B,
227: $ LDB, AP( KC ), 1, ONE, B( K, 1 ), LDB )
228: CALL ZLACGV( NRHS, B( K, 1 ), LDB )
229: *
230: CALL ZLACGV( NRHS, B( K+1, 1 ), LDB )
231: CALL ZGEMV( 'Conjugate transpose', K-1, NRHS, -ONE, B,
232: $ LDB, AP( KC+K ), 1, ONE, B( K+1, 1 ), LDB )
233: CALL ZLACGV( NRHS, B( K+1, 1 ), LDB )
234: END IF
235: *
236: * Interchange rows K and -IPIV(K).
237: *
238: KP = -IPIV( K )
239: IF( KP.NE.K )
240: $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
241: KC = KC + 2*K + 1
242: K = K + 2
243: END IF
244: *
245: GO TO 40
246: 50 CONTINUE
247: *
248: ELSE
249: *
250: * Solve A*X = B, where A = L*D*L'.
251: *
252: * First solve L*D*X = B, overwriting B with X.
253: *
254: * K is the main loop index, increasing from 1 to N in steps of
255: * 1 or 2, depending on the size of the diagonal blocks.
256: *
257: K = 1
258: KC = 1
259: 60 CONTINUE
260: *
261: * If K > N, exit from loop.
262: *
263: IF( K.GT.N )
264: $ GO TO 80
265: *
266: IF( IPIV( K ).GT.0 ) THEN
267: *
268: * 1 x 1 diagonal block
269: *
270: * Interchange rows K and IPIV(K).
271: *
272: KP = IPIV( K )
273: IF( KP.NE.K )
274: $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
275: *
276: * Multiply by inv(L(K)), where L(K) is the transformation
277: * stored in column K of A.
278: *
279: IF( K.LT.N )
280: $ CALL ZGERU( N-K, NRHS, -ONE, AP( KC+1 ), 1, B( K, 1 ),
281: $ LDB, B( K+1, 1 ), LDB )
282: *
283: * Multiply by the inverse of the diagonal block.
284: *
285: S = DBLE( ONE ) / DBLE( AP( KC ) )
286: CALL ZDSCAL( NRHS, S, B( K, 1 ), LDB )
287: KC = KC + N - K + 1
288: K = K + 1
289: ELSE
290: *
291: * 2 x 2 diagonal block
292: *
293: * Interchange rows K+1 and -IPIV(K).
294: *
295: KP = -IPIV( K )
296: IF( KP.NE.K+1 )
297: $ CALL ZSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB )
298: *
299: * Multiply by inv(L(K)), where L(K) is the transformation
300: * stored in columns K and K+1 of A.
301: *
302: IF( K.LT.N-1 ) THEN
303: CALL ZGERU( N-K-1, NRHS, -ONE, AP( KC+2 ), 1, B( K, 1 ),
304: $ LDB, B( K+2, 1 ), LDB )
305: CALL ZGERU( N-K-1, NRHS, -ONE, AP( KC+N-K+2 ), 1,
306: $ B( K+1, 1 ), LDB, B( K+2, 1 ), LDB )
307: END IF
308: *
309: * Multiply by the inverse of the diagonal block.
310: *
311: AKM1K = AP( KC+1 )
312: AKM1 = AP( KC ) / DCONJG( AKM1K )
313: AK = AP( KC+N-K+1 ) / AKM1K
314: DENOM = AKM1*AK - ONE
315: DO 70 J = 1, NRHS
316: BKM1 = B( K, J ) / DCONJG( AKM1K )
317: BK = B( K+1, J ) / AKM1K
318: B( K, J ) = ( AK*BKM1-BK ) / DENOM
319: B( K+1, J ) = ( AKM1*BK-BKM1 ) / DENOM
320: 70 CONTINUE
321: KC = KC + 2*( N-K ) + 1
322: K = K + 2
323: END IF
324: *
325: GO TO 60
326: 80 CONTINUE
327: *
328: * Next solve L'*X = B, overwriting B with X.
329: *
330: * K is the main loop index, decreasing from N to 1 in steps of
331: * 1 or 2, depending on the size of the diagonal blocks.
332: *
333: K = N
334: KC = N*( N+1 ) / 2 + 1
335: 90 CONTINUE
336: *
337: * If K < 1, exit from loop.
338: *
339: IF( K.LT.1 )
340: $ GO TO 100
341: *
342: KC = KC - ( N-K+1 )
343: IF( IPIV( K ).GT.0 ) THEN
344: *
345: * 1 x 1 diagonal block
346: *
347: * Multiply by inv(L'(K)), where L(K) is the transformation
348: * stored in column K of A.
349: *
350: IF( K.LT.N ) THEN
351: CALL ZLACGV( NRHS, B( K, 1 ), LDB )
352: CALL ZGEMV( 'Conjugate transpose', N-K, NRHS, -ONE,
353: $ B( K+1, 1 ), LDB, AP( KC+1 ), 1, ONE,
354: $ B( K, 1 ), LDB )
355: CALL ZLACGV( NRHS, B( K, 1 ), LDB )
356: END IF
357: *
358: * Interchange rows K and IPIV(K).
359: *
360: KP = IPIV( K )
361: IF( KP.NE.K )
362: $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
363: K = K - 1
364: ELSE
365: *
366: * 2 x 2 diagonal block
367: *
368: * Multiply by inv(L'(K-1)), where L(K-1) is the transformation
369: * stored in columns K-1 and K of A.
370: *
371: IF( K.LT.N ) THEN
372: CALL ZLACGV( NRHS, B( K, 1 ), LDB )
373: CALL ZGEMV( 'Conjugate transpose', N-K, NRHS, -ONE,
374: $ B( K+1, 1 ), LDB, AP( KC+1 ), 1, ONE,
375: $ B( K, 1 ), LDB )
376: CALL ZLACGV( NRHS, B( K, 1 ), LDB )
377: *
378: CALL ZLACGV( NRHS, B( K-1, 1 ), LDB )
379: CALL ZGEMV( 'Conjugate transpose', N-K, NRHS, -ONE,
380: $ B( K+1, 1 ), LDB, AP( KC-( N-K ) ), 1, ONE,
381: $ B( K-1, 1 ), LDB )
382: CALL ZLACGV( NRHS, B( K-1, 1 ), LDB )
383: END IF
384: *
385: * Interchange rows K and -IPIV(K).
386: *
387: KP = -IPIV( K )
388: IF( KP.NE.K )
389: $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
390: KC = KC - ( N-K+2 )
391: K = K - 2
392: END IF
393: *
394: GO TO 90
395: 100 CONTINUE
396: END IF
397: *
398: RETURN
399: *
400: * End of ZHPTRS
401: *
402: END