Annotation of rpl/lapack/lapack/zhetrs.f, revision 1.18
1.9 bertrand 1: *> \brief \b ZHETRS
2: *
3: * =========== DOCUMENTATION ===========
4: *
1.15 bertrand 5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
1.9 bertrand 7: *
8: *> \htmlonly
1.15 bertrand 9: *> Download ZHETRS + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhetrs.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhetrs.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhetrs.f">
1.9 bertrand 15: *> [TXT]</a>
1.15 bertrand 16: *> \endhtmlonly
1.9 bertrand 17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZHETRS( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO )
1.15 bertrand 22: *
1.9 bertrand 23: * .. Scalar Arguments ..
24: * CHARACTER UPLO
25: * INTEGER INFO, LDA, LDB, N, NRHS
26: * ..
27: * .. Array Arguments ..
28: * INTEGER IPIV( * )
29: * COMPLEX*16 A( LDA, * ), B( LDB, * )
30: * ..
1.15 bertrand 31: *
1.9 bertrand 32: *
33: *> \par Purpose:
34: * =============
35: *>
36: *> \verbatim
37: *>
38: *> ZHETRS solves a system of linear equations A*X = B with a complex
39: *> Hermitian matrix A using the factorization A = U*D*U**H or
40: *> A = L*D*L**H computed by ZHETRF.
41: *> \endverbatim
42: *
43: * Arguments:
44: * ==========
45: *
46: *> \param[in] UPLO
47: *> \verbatim
48: *> UPLO is CHARACTER*1
49: *> Specifies whether the details of the factorization are stored
50: *> as an upper or lower triangular matrix.
51: *> = 'U': Upper triangular, form is A = U*D*U**H;
52: *> = 'L': Lower triangular, form is A = L*D*L**H.
53: *> \endverbatim
54: *>
55: *> \param[in] N
56: *> \verbatim
57: *> N is INTEGER
58: *> The order of the matrix A. N >= 0.
59: *> \endverbatim
60: *>
61: *> \param[in] NRHS
62: *> \verbatim
63: *> NRHS is INTEGER
64: *> The number of right hand sides, i.e., the number of columns
65: *> of the matrix B. NRHS >= 0.
66: *> \endverbatim
67: *>
68: *> \param[in] A
69: *> \verbatim
70: *> A is COMPLEX*16 array, dimension (LDA,N)
71: *> The block diagonal matrix D and the multipliers used to
72: *> obtain the factor U or L as computed by ZHETRF.
73: *> \endverbatim
74: *>
75: *> \param[in] LDA
76: *> \verbatim
77: *> LDA is INTEGER
78: *> The leading dimension of the array A. LDA >= max(1,N).
79: *> \endverbatim
80: *>
81: *> \param[in] IPIV
82: *> \verbatim
83: *> IPIV is INTEGER array, dimension (N)
84: *> Details of the interchanges and the block structure of D
85: *> as determined by ZHETRF.
86: *> \endverbatim
87: *>
88: *> \param[in,out] B
89: *> \verbatim
90: *> B is COMPLEX*16 array, dimension (LDB,NRHS)
91: *> On entry, the right hand side matrix B.
92: *> On exit, the solution matrix X.
93: *> \endverbatim
94: *>
95: *> \param[in] LDB
96: *> \verbatim
97: *> LDB is INTEGER
98: *> The leading dimension of the array B. LDB >= max(1,N).
99: *> \endverbatim
100: *>
101: *> \param[out] INFO
102: *> \verbatim
103: *> INFO is INTEGER
104: *> = 0: successful exit
105: *> < 0: if INFO = -i, the i-th argument had an illegal value
106: *> \endverbatim
107: *
108: * Authors:
109: * ========
110: *
1.15 bertrand 111: *> \author Univ. of Tennessee
112: *> \author Univ. of California Berkeley
113: *> \author Univ. of Colorado Denver
114: *> \author NAG Ltd.
1.9 bertrand 115: *
116: *> \ingroup complex16HEcomputational
117: *
118: * =====================================================================
1.1 bertrand 119: SUBROUTINE ZHETRS( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO )
120: *
1.18 ! bertrand 121: * -- LAPACK computational routine --
1.1 bertrand 122: * -- LAPACK is a software package provided by Univ. of Tennessee, --
123: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
124: *
125: * .. Scalar Arguments ..
126: CHARACTER UPLO
127: INTEGER INFO, LDA, LDB, N, NRHS
128: * ..
129: * .. Array Arguments ..
130: INTEGER IPIV( * )
131: COMPLEX*16 A( LDA, * ), B( LDB, * )
132: * ..
133: *
134: * =====================================================================
135: *
136: * .. Parameters ..
137: COMPLEX*16 ONE
138: PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) )
139: * ..
140: * .. Local Scalars ..
141: LOGICAL UPPER
142: INTEGER J, K, KP
143: DOUBLE PRECISION S
144: COMPLEX*16 AK, AKM1, AKM1K, BK, BKM1, DENOM
145: * ..
146: * .. External Functions ..
147: LOGICAL LSAME
148: EXTERNAL LSAME
149: * ..
150: * .. External Subroutines ..
151: EXTERNAL XERBLA, ZDSCAL, ZGEMV, ZGERU, ZLACGV, ZSWAP
152: * ..
153: * .. Intrinsic Functions ..
154: INTRINSIC DBLE, DCONJG, MAX
155: * ..
156: * .. Executable Statements ..
157: *
158: INFO = 0
159: UPPER = LSAME( UPLO, 'U' )
160: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
161: INFO = -1
162: ELSE IF( N.LT.0 ) THEN
163: INFO = -2
164: ELSE IF( NRHS.LT.0 ) THEN
165: INFO = -3
166: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
167: INFO = -5
168: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
169: INFO = -8
170: END IF
171: IF( INFO.NE.0 ) THEN
172: CALL XERBLA( 'ZHETRS', -INFO )
173: RETURN
174: END IF
175: *
176: * Quick return if possible
177: *
178: IF( N.EQ.0 .OR. NRHS.EQ.0 )
179: $ RETURN
180: *
181: IF( UPPER ) THEN
182: *
1.8 bertrand 183: * Solve A*X = B, where A = U*D*U**H.
1.1 bertrand 184: *
185: * First solve U*D*X = B, overwriting B with X.
186: *
187: * K is the main loop index, decreasing from N to 1 in steps of
188: * 1 or 2, depending on the size of the diagonal blocks.
189: *
190: K = N
191: 10 CONTINUE
192: *
193: * If K < 1, exit from loop.
194: *
195: IF( K.LT.1 )
196: $ GO TO 30
197: *
198: IF( IPIV( K ).GT.0 ) THEN
199: *
200: * 1 x 1 diagonal block
201: *
202: * Interchange rows K and IPIV(K).
203: *
204: KP = IPIV( K )
205: IF( KP.NE.K )
206: $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
207: *
208: * Multiply by inv(U(K)), where U(K) is the transformation
209: * stored in column K of A.
210: *
211: CALL ZGERU( K-1, NRHS, -ONE, A( 1, K ), 1, B( K, 1 ), LDB,
212: $ B( 1, 1 ), LDB )
213: *
214: * Multiply by the inverse of the diagonal block.
215: *
216: S = DBLE( ONE ) / DBLE( A( K, K ) )
217: CALL ZDSCAL( NRHS, S, B( K, 1 ), LDB )
218: K = K - 1
219: ELSE
220: *
221: * 2 x 2 diagonal block
222: *
223: * Interchange rows K-1 and -IPIV(K).
224: *
225: KP = -IPIV( K )
226: IF( KP.NE.K-1 )
227: $ CALL ZSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB )
228: *
229: * Multiply by inv(U(K)), where U(K) is the transformation
230: * stored in columns K-1 and K of A.
231: *
232: CALL ZGERU( K-2, NRHS, -ONE, A( 1, K ), 1, B( K, 1 ), LDB,
233: $ B( 1, 1 ), LDB )
234: CALL ZGERU( K-2, NRHS, -ONE, A( 1, K-1 ), 1, B( K-1, 1 ),
235: $ LDB, B( 1, 1 ), LDB )
236: *
237: * Multiply by the inverse of the diagonal block.
238: *
239: AKM1K = A( K-1, K )
240: AKM1 = A( K-1, K-1 ) / AKM1K
241: AK = A( K, K ) / DCONJG( AKM1K )
242: DENOM = AKM1*AK - ONE
243: DO 20 J = 1, NRHS
244: BKM1 = B( K-1, J ) / AKM1K
245: BK = B( K, J ) / DCONJG( AKM1K )
246: B( K-1, J ) = ( AK*BKM1-BK ) / DENOM
247: B( K, J ) = ( AKM1*BK-BKM1 ) / DENOM
248: 20 CONTINUE
249: K = K - 2
250: END IF
251: *
252: GO TO 10
253: 30 CONTINUE
254: *
1.8 bertrand 255: * Next solve U**H *X = B, overwriting B with X.
1.1 bertrand 256: *
257: * K is the main loop index, increasing from 1 to N in steps of
258: * 1 or 2, depending on the size of the diagonal blocks.
259: *
260: K = 1
261: 40 CONTINUE
262: *
263: * If K > N, exit from loop.
264: *
265: IF( K.GT.N )
266: $ GO TO 50
267: *
268: IF( IPIV( K ).GT.0 ) THEN
269: *
270: * 1 x 1 diagonal block
271: *
1.8 bertrand 272: * Multiply by inv(U**H(K)), where U(K) is the transformation
1.1 bertrand 273: * stored in column K of A.
274: *
275: IF( K.GT.1 ) THEN
276: CALL ZLACGV( NRHS, B( K, 1 ), LDB )
277: CALL ZGEMV( 'Conjugate transpose', K-1, NRHS, -ONE, B,
278: $ LDB, A( 1, K ), 1, ONE, B( K, 1 ), LDB )
279: CALL ZLACGV( NRHS, B( K, 1 ), LDB )
280: END IF
281: *
282: * Interchange rows K and IPIV(K).
283: *
284: KP = IPIV( K )
285: IF( KP.NE.K )
286: $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
287: K = K + 1
288: ELSE
289: *
290: * 2 x 2 diagonal block
291: *
1.8 bertrand 292: * Multiply by inv(U**H(K+1)), where U(K+1) is the transformation
1.1 bertrand 293: * stored in columns K and K+1 of A.
294: *
295: IF( K.GT.1 ) THEN
296: CALL ZLACGV( NRHS, B( K, 1 ), LDB )
297: CALL ZGEMV( 'Conjugate transpose', K-1, NRHS, -ONE, B,
298: $ LDB, A( 1, K ), 1, ONE, B( K, 1 ), LDB )
299: CALL ZLACGV( NRHS, B( K, 1 ), LDB )
300: *
301: CALL ZLACGV( NRHS, B( K+1, 1 ), LDB )
302: CALL ZGEMV( 'Conjugate transpose', K-1, NRHS, -ONE, B,
303: $ LDB, A( 1, K+1 ), 1, ONE, B( K+1, 1 ), LDB )
304: CALL ZLACGV( NRHS, B( K+1, 1 ), LDB )
305: END IF
306: *
307: * Interchange rows K and -IPIV(K).
308: *
309: KP = -IPIV( K )
310: IF( KP.NE.K )
311: $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
312: K = K + 2
313: END IF
314: *
315: GO TO 40
316: 50 CONTINUE
317: *
318: ELSE
319: *
1.8 bertrand 320: * Solve A*X = B, where A = L*D*L**H.
1.1 bertrand 321: *
322: * First solve L*D*X = B, overwriting B with X.
323: *
324: * K is the main loop index, increasing from 1 to N in steps of
325: * 1 or 2, depending on the size of the diagonal blocks.
326: *
327: K = 1
328: 60 CONTINUE
329: *
330: * If K > N, exit from loop.
331: *
332: IF( K.GT.N )
333: $ GO TO 80
334: *
335: IF( IPIV( K ).GT.0 ) THEN
336: *
337: * 1 x 1 diagonal block
338: *
339: * Interchange rows K and IPIV(K).
340: *
341: KP = IPIV( K )
342: IF( KP.NE.K )
343: $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
344: *
345: * Multiply by inv(L(K)), where L(K) is the transformation
346: * stored in column K of A.
347: *
348: IF( K.LT.N )
349: $ CALL ZGERU( N-K, NRHS, -ONE, A( K+1, K ), 1, B( K, 1 ),
350: $ LDB, B( K+1, 1 ), LDB )
351: *
352: * Multiply by the inverse of the diagonal block.
353: *
354: S = DBLE( ONE ) / DBLE( A( K, K ) )
355: CALL ZDSCAL( NRHS, S, B( K, 1 ), LDB )
356: K = K + 1
357: ELSE
358: *
359: * 2 x 2 diagonal block
360: *
361: * Interchange rows K+1 and -IPIV(K).
362: *
363: KP = -IPIV( K )
364: IF( KP.NE.K+1 )
365: $ CALL ZSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB )
366: *
367: * Multiply by inv(L(K)), where L(K) is the transformation
368: * stored in columns K and K+1 of A.
369: *
370: IF( K.LT.N-1 ) THEN
371: CALL ZGERU( N-K-1, NRHS, -ONE, A( K+2, K ), 1, B( K, 1 ),
372: $ LDB, B( K+2, 1 ), LDB )
373: CALL ZGERU( N-K-1, NRHS, -ONE, A( K+2, K+1 ), 1,
374: $ B( K+1, 1 ), LDB, B( K+2, 1 ), LDB )
375: END IF
376: *
377: * Multiply by the inverse of the diagonal block.
378: *
379: AKM1K = A( K+1, K )
380: AKM1 = A( K, K ) / DCONJG( AKM1K )
381: AK = A( K+1, K+1 ) / AKM1K
382: DENOM = AKM1*AK - ONE
383: DO 70 J = 1, NRHS
384: BKM1 = B( K, J ) / DCONJG( AKM1K )
385: BK = B( K+1, J ) / AKM1K
386: B( K, J ) = ( AK*BKM1-BK ) / DENOM
387: B( K+1, J ) = ( AKM1*BK-BKM1 ) / DENOM
388: 70 CONTINUE
389: K = K + 2
390: END IF
391: *
392: GO TO 60
393: 80 CONTINUE
394: *
1.8 bertrand 395: * Next solve L**H *X = B, overwriting B with X.
1.1 bertrand 396: *
397: * K is the main loop index, decreasing from N to 1 in steps of
398: * 1 or 2, depending on the size of the diagonal blocks.
399: *
400: K = N
401: 90 CONTINUE
402: *
403: * If K < 1, exit from loop.
404: *
405: IF( K.LT.1 )
406: $ GO TO 100
407: *
408: IF( IPIV( K ).GT.0 ) THEN
409: *
410: * 1 x 1 diagonal block
411: *
1.8 bertrand 412: * Multiply by inv(L**H(K)), where L(K) is the transformation
1.1 bertrand 413: * stored in column K of A.
414: *
415: IF( K.LT.N ) THEN
416: CALL ZLACGV( NRHS, B( K, 1 ), LDB )
417: CALL ZGEMV( 'Conjugate transpose', N-K, NRHS, -ONE,
418: $ B( K+1, 1 ), LDB, A( K+1, K ), 1, ONE,
419: $ B( K, 1 ), LDB )
420: CALL ZLACGV( NRHS, B( K, 1 ), LDB )
421: END IF
422: *
423: * Interchange rows K and IPIV(K).
424: *
425: KP = IPIV( K )
426: IF( KP.NE.K )
427: $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
428: K = K - 1
429: ELSE
430: *
431: * 2 x 2 diagonal block
432: *
1.8 bertrand 433: * Multiply by inv(L**H(K-1)), where L(K-1) is the transformation
1.1 bertrand 434: * stored in columns K-1 and K of A.
435: *
436: IF( K.LT.N ) THEN
437: CALL ZLACGV( NRHS, B( K, 1 ), LDB )
438: CALL ZGEMV( 'Conjugate transpose', N-K, NRHS, -ONE,
439: $ B( K+1, 1 ), LDB, A( K+1, K ), 1, ONE,
440: $ B( K, 1 ), LDB )
441: CALL ZLACGV( NRHS, B( K, 1 ), LDB )
442: *
443: CALL ZLACGV( NRHS, B( K-1, 1 ), LDB )
444: CALL ZGEMV( 'Conjugate transpose', N-K, NRHS, -ONE,
445: $ B( K+1, 1 ), LDB, A( K+1, K-1 ), 1, ONE,
446: $ B( K-1, 1 ), LDB )
447: CALL ZLACGV( NRHS, B( K-1, 1 ), LDB )
448: END IF
449: *
450: * Interchange rows K and -IPIV(K).
451: *
452: KP = -IPIV( K )
453: IF( KP.NE.K )
454: $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
455: K = K - 2
456: END IF
457: *
458: GO TO 90
459: 100 CONTINUE
460: END IF
461: *
462: RETURN
463: *
464: * End of ZHETRS
465: *
466: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zhetrs.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack