Annotation of rpl/lapack/lapack/zhetrs.f, revision 1.9
1.9 ! bertrand 1: *> \brief \b ZHETRS
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 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">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE ZHETRS( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO )
! 22: *
! 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: * ..
! 31: *
! 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: *
! 111: *> \author Univ. of Tennessee
! 112: *> \author Univ. of California Berkeley
! 113: *> \author Univ. of Colorado Denver
! 114: *> \author NAG Ltd.
! 115: *
! 116: *> \date November 2011
! 117: *
! 118: *> \ingroup complex16HEcomputational
! 119: *
! 120: * =====================================================================
1.1 bertrand 121: SUBROUTINE ZHETRS( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO )
122: *
1.9 ! bertrand 123: * -- LAPACK computational routine (version 3.4.0) --
1.1 bertrand 124: * -- LAPACK is a software package provided by Univ. of Tennessee, --
125: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.9 ! bertrand 126: * November 2011
1.1 bertrand 127: *
128: * .. Scalar Arguments ..
129: CHARACTER UPLO
130: INTEGER INFO, LDA, LDB, N, NRHS
131: * ..
132: * .. Array Arguments ..
133: INTEGER IPIV( * )
134: COMPLEX*16 A( LDA, * ), B( LDB, * )
135: * ..
136: *
137: * =====================================================================
138: *
139: * .. Parameters ..
140: COMPLEX*16 ONE
141: PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) )
142: * ..
143: * .. Local Scalars ..
144: LOGICAL UPPER
145: INTEGER J, K, KP
146: DOUBLE PRECISION S
147: COMPLEX*16 AK, AKM1, AKM1K, BK, BKM1, DENOM
148: * ..
149: * .. External Functions ..
150: LOGICAL LSAME
151: EXTERNAL LSAME
152: * ..
153: * .. External Subroutines ..
154: EXTERNAL XERBLA, ZDSCAL, ZGEMV, ZGERU, ZLACGV, ZSWAP
155: * ..
156: * .. Intrinsic Functions ..
157: INTRINSIC DBLE, DCONJG, MAX
158: * ..
159: * .. Executable Statements ..
160: *
161: INFO = 0
162: UPPER = LSAME( UPLO, 'U' )
163: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
164: INFO = -1
165: ELSE IF( N.LT.0 ) THEN
166: INFO = -2
167: ELSE IF( NRHS.LT.0 ) THEN
168: INFO = -3
169: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
170: INFO = -5
171: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
172: INFO = -8
173: END IF
174: IF( INFO.NE.0 ) THEN
175: CALL XERBLA( 'ZHETRS', -INFO )
176: RETURN
177: END IF
178: *
179: * Quick return if possible
180: *
181: IF( N.EQ.0 .OR. NRHS.EQ.0 )
182: $ RETURN
183: *
184: IF( UPPER ) THEN
185: *
1.8 bertrand 186: * Solve A*X = B, where A = U*D*U**H.
1.1 bertrand 187: *
188: * First solve U*D*X = B, overwriting B with X.
189: *
190: * K is the main loop index, decreasing from N to 1 in steps of
191: * 1 or 2, depending on the size of the diagonal blocks.
192: *
193: K = N
194: 10 CONTINUE
195: *
196: * If K < 1, exit from loop.
197: *
198: IF( K.LT.1 )
199: $ GO TO 30
200: *
201: IF( IPIV( K ).GT.0 ) THEN
202: *
203: * 1 x 1 diagonal block
204: *
205: * Interchange rows K and IPIV(K).
206: *
207: KP = IPIV( K )
208: IF( KP.NE.K )
209: $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
210: *
211: * Multiply by inv(U(K)), where U(K) is the transformation
212: * stored in column K of A.
213: *
214: CALL ZGERU( K-1, NRHS, -ONE, A( 1, K ), 1, B( K, 1 ), LDB,
215: $ B( 1, 1 ), LDB )
216: *
217: * Multiply by the inverse of the diagonal block.
218: *
219: S = DBLE( ONE ) / DBLE( A( K, K ) )
220: CALL ZDSCAL( NRHS, S, B( K, 1 ), LDB )
221: K = K - 1
222: ELSE
223: *
224: * 2 x 2 diagonal block
225: *
226: * Interchange rows K-1 and -IPIV(K).
227: *
228: KP = -IPIV( K )
229: IF( KP.NE.K-1 )
230: $ CALL ZSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB )
231: *
232: * Multiply by inv(U(K)), where U(K) is the transformation
233: * stored in columns K-1 and K of A.
234: *
235: CALL ZGERU( K-2, NRHS, -ONE, A( 1, K ), 1, B( K, 1 ), LDB,
236: $ B( 1, 1 ), LDB )
237: CALL ZGERU( K-2, NRHS, -ONE, A( 1, K-1 ), 1, B( K-1, 1 ),
238: $ LDB, B( 1, 1 ), LDB )
239: *
240: * Multiply by the inverse of the diagonal block.
241: *
242: AKM1K = A( K-1, K )
243: AKM1 = A( K-1, K-1 ) / AKM1K
244: AK = A( K, K ) / DCONJG( AKM1K )
245: DENOM = AKM1*AK - ONE
246: DO 20 J = 1, NRHS
247: BKM1 = B( K-1, J ) / AKM1K
248: BK = B( K, J ) / DCONJG( AKM1K )
249: B( K-1, J ) = ( AK*BKM1-BK ) / DENOM
250: B( K, J ) = ( AKM1*BK-BKM1 ) / DENOM
251: 20 CONTINUE
252: K = K - 2
253: END IF
254: *
255: GO TO 10
256: 30 CONTINUE
257: *
1.8 bertrand 258: * Next solve U**H *X = B, overwriting B with X.
1.1 bertrand 259: *
260: * K is the main loop index, increasing from 1 to N in steps of
261: * 1 or 2, depending on the size of the diagonal blocks.
262: *
263: K = 1
264: 40 CONTINUE
265: *
266: * If K > N, exit from loop.
267: *
268: IF( K.GT.N )
269: $ GO TO 50
270: *
271: IF( IPIV( K ).GT.0 ) THEN
272: *
273: * 1 x 1 diagonal block
274: *
1.8 bertrand 275: * Multiply by inv(U**H(K)), where U(K) is the transformation
1.1 bertrand 276: * stored in column K of A.
277: *
278: IF( K.GT.1 ) THEN
279: CALL ZLACGV( NRHS, B( K, 1 ), LDB )
280: CALL ZGEMV( 'Conjugate transpose', K-1, NRHS, -ONE, B,
281: $ LDB, A( 1, K ), 1, ONE, B( K, 1 ), LDB )
282: CALL ZLACGV( NRHS, B( K, 1 ), LDB )
283: END IF
284: *
285: * Interchange rows K and IPIV(K).
286: *
287: KP = IPIV( K )
288: IF( KP.NE.K )
289: $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
290: K = K + 1
291: ELSE
292: *
293: * 2 x 2 diagonal block
294: *
1.8 bertrand 295: * Multiply by inv(U**H(K+1)), where U(K+1) is the transformation
1.1 bertrand 296: * stored in columns K and K+1 of A.
297: *
298: IF( K.GT.1 ) THEN
299: CALL ZLACGV( NRHS, B( K, 1 ), LDB )
300: CALL ZGEMV( 'Conjugate transpose', K-1, NRHS, -ONE, B,
301: $ LDB, A( 1, K ), 1, ONE, B( K, 1 ), LDB )
302: CALL ZLACGV( NRHS, B( K, 1 ), LDB )
303: *
304: CALL ZLACGV( NRHS, B( K+1, 1 ), LDB )
305: CALL ZGEMV( 'Conjugate transpose', K-1, NRHS, -ONE, B,
306: $ LDB, A( 1, K+1 ), 1, ONE, B( K+1, 1 ), LDB )
307: CALL ZLACGV( NRHS, B( K+1, 1 ), LDB )
308: END IF
309: *
310: * Interchange rows K and -IPIV(K).
311: *
312: KP = -IPIV( K )
313: IF( KP.NE.K )
314: $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
315: K = K + 2
316: END IF
317: *
318: GO TO 40
319: 50 CONTINUE
320: *
321: ELSE
322: *
1.8 bertrand 323: * Solve A*X = B, where A = L*D*L**H.
1.1 bertrand 324: *
325: * First solve L*D*X = B, overwriting B with X.
326: *
327: * K is the main loop index, increasing from 1 to N in steps of
328: * 1 or 2, depending on the size of the diagonal blocks.
329: *
330: K = 1
331: 60 CONTINUE
332: *
333: * If K > N, exit from loop.
334: *
335: IF( K.GT.N )
336: $ GO TO 80
337: *
338: IF( IPIV( K ).GT.0 ) THEN
339: *
340: * 1 x 1 diagonal block
341: *
342: * Interchange rows K and IPIV(K).
343: *
344: KP = IPIV( K )
345: IF( KP.NE.K )
346: $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
347: *
348: * Multiply by inv(L(K)), where L(K) is the transformation
349: * stored in column K of A.
350: *
351: IF( K.LT.N )
352: $ CALL ZGERU( N-K, NRHS, -ONE, A( K+1, K ), 1, B( K, 1 ),
353: $ LDB, B( K+1, 1 ), LDB )
354: *
355: * Multiply by the inverse of the diagonal block.
356: *
357: S = DBLE( ONE ) / DBLE( A( K, K ) )
358: CALL ZDSCAL( NRHS, S, B( K, 1 ), LDB )
359: K = K + 1
360: ELSE
361: *
362: * 2 x 2 diagonal block
363: *
364: * Interchange rows K+1 and -IPIV(K).
365: *
366: KP = -IPIV( K )
367: IF( KP.NE.K+1 )
368: $ CALL ZSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB )
369: *
370: * Multiply by inv(L(K)), where L(K) is the transformation
371: * stored in columns K and K+1 of A.
372: *
373: IF( K.LT.N-1 ) THEN
374: CALL ZGERU( N-K-1, NRHS, -ONE, A( K+2, K ), 1, B( K, 1 ),
375: $ LDB, B( K+2, 1 ), LDB )
376: CALL ZGERU( N-K-1, NRHS, -ONE, A( K+2, K+1 ), 1,
377: $ B( K+1, 1 ), LDB, B( K+2, 1 ), LDB )
378: END IF
379: *
380: * Multiply by the inverse of the diagonal block.
381: *
382: AKM1K = A( K+1, K )
383: AKM1 = A( K, K ) / DCONJG( AKM1K )
384: AK = A( K+1, K+1 ) / AKM1K
385: DENOM = AKM1*AK - ONE
386: DO 70 J = 1, NRHS
387: BKM1 = B( K, J ) / DCONJG( AKM1K )
388: BK = B( K+1, J ) / AKM1K
389: B( K, J ) = ( AK*BKM1-BK ) / DENOM
390: B( K+1, J ) = ( AKM1*BK-BKM1 ) / DENOM
391: 70 CONTINUE
392: K = K + 2
393: END IF
394: *
395: GO TO 60
396: 80 CONTINUE
397: *
1.8 bertrand 398: * Next solve L**H *X = B, overwriting B with X.
1.1 bertrand 399: *
400: * K is the main loop index, decreasing from N to 1 in steps of
401: * 1 or 2, depending on the size of the diagonal blocks.
402: *
403: K = N
404: 90 CONTINUE
405: *
406: * If K < 1, exit from loop.
407: *
408: IF( K.LT.1 )
409: $ GO TO 100
410: *
411: IF( IPIV( K ).GT.0 ) THEN
412: *
413: * 1 x 1 diagonal block
414: *
1.8 bertrand 415: * Multiply by inv(L**H(K)), where L(K) is the transformation
1.1 bertrand 416: * stored in column K of A.
417: *
418: IF( K.LT.N ) THEN
419: CALL ZLACGV( NRHS, B( K, 1 ), LDB )
420: CALL ZGEMV( 'Conjugate transpose', N-K, NRHS, -ONE,
421: $ B( K+1, 1 ), LDB, A( K+1, K ), 1, ONE,
422: $ B( K, 1 ), LDB )
423: CALL ZLACGV( NRHS, B( K, 1 ), LDB )
424: END IF
425: *
426: * Interchange rows K and IPIV(K).
427: *
428: KP = IPIV( K )
429: IF( KP.NE.K )
430: $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
431: K = K - 1
432: ELSE
433: *
434: * 2 x 2 diagonal block
435: *
1.8 bertrand 436: * Multiply by inv(L**H(K-1)), where L(K-1) is the transformation
1.1 bertrand 437: * stored in columns K-1 and K of A.
438: *
439: IF( K.LT.N ) THEN
440: CALL ZLACGV( NRHS, B( K, 1 ), LDB )
441: CALL ZGEMV( 'Conjugate transpose', N-K, NRHS, -ONE,
442: $ B( K+1, 1 ), LDB, A( K+1, K ), 1, ONE,
443: $ B( K, 1 ), LDB )
444: CALL ZLACGV( NRHS, B( K, 1 ), LDB )
445: *
446: CALL ZLACGV( NRHS, B( K-1, 1 ), LDB )
447: CALL ZGEMV( 'Conjugate transpose', N-K, NRHS, -ONE,
448: $ B( K+1, 1 ), LDB, A( K+1, K-1 ), 1, ONE,
449: $ B( K-1, 1 ), LDB )
450: CALL ZLACGV( NRHS, B( K-1, 1 ), LDB )
451: END IF
452: *
453: * Interchange rows K and -IPIV(K).
454: *
455: KP = -IPIV( K )
456: IF( KP.NE.K )
457: $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
458: K = K - 2
459: END IF
460: *
461: GO TO 90
462: 100 CONTINUE
463: END IF
464: *
465: RETURN
466: *
467: * End of ZHETRS
468: *
469: 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