Annotation of rpl/lapack/lapack/zhetrs2.f, revision 1.4
1.4 ! bertrand 1: *> \brief \b ZHETRS2
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download ZHETRS2 + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhetrs2.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhetrs2.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhetrs2.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE ZHETRS2( UPLO, N, NRHS, A, LDA, IPIV, B, LDB,
! 22: * WORK, INFO )
! 23: *
! 24: * .. Scalar Arguments ..
! 25: * CHARACTER UPLO
! 26: * INTEGER INFO, LDA, LDB, N, NRHS
! 27: * ..
! 28: * .. Array Arguments ..
! 29: * INTEGER IPIV( * )
! 30: * COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * )
! 31: * ..
! 32: *
! 33: *
! 34: *> \par Purpose:
! 35: * =============
! 36: *>
! 37: *> \verbatim
! 38: *>
! 39: *> ZHETRS2 solves a system of linear equations A*X = B with a complex
! 40: *> Hermitian matrix A using the factorization A = U*D*U**H or
! 41: *> A = L*D*L**H computed by ZHETRF and converted by ZSYCONV.
! 42: *> \endverbatim
! 43: *
! 44: * Arguments:
! 45: * ==========
! 46: *
! 47: *> \param[in] UPLO
! 48: *> \verbatim
! 49: *> UPLO is CHARACTER*1
! 50: *> Specifies whether the details of the factorization are stored
! 51: *> as an upper or lower triangular matrix.
! 52: *> = 'U': Upper triangular, form is A = U*D*U**H;
! 53: *> = 'L': Lower triangular, form is A = L*D*L**H.
! 54: *> \endverbatim
! 55: *>
! 56: *> \param[in] N
! 57: *> \verbatim
! 58: *> N is INTEGER
! 59: *> The order of the matrix A. N >= 0.
! 60: *> \endverbatim
! 61: *>
! 62: *> \param[in] NRHS
! 63: *> \verbatim
! 64: *> NRHS is INTEGER
! 65: *> The number of right hand sides, i.e., the number of columns
! 66: *> of the matrix B. NRHS >= 0.
! 67: *> \endverbatim
! 68: *>
! 69: *> \param[in] A
! 70: *> \verbatim
! 71: *> A is COMPLEX*16 array, dimension (LDA,N)
! 72: *> The block diagonal matrix D and the multipliers used to
! 73: *> obtain the factor U or L as computed by ZHETRF.
! 74: *> \endverbatim
! 75: *>
! 76: *> \param[in] LDA
! 77: *> \verbatim
! 78: *> LDA is INTEGER
! 79: *> The leading dimension of the array A. LDA >= max(1,N).
! 80: *> \endverbatim
! 81: *>
! 82: *> \param[in] IPIV
! 83: *> \verbatim
! 84: *> IPIV is INTEGER array, dimension (N)
! 85: *> Details of the interchanges and the block structure of D
! 86: *> as determined by ZHETRF.
! 87: *> \endverbatim
! 88: *>
! 89: *> \param[in,out] B
! 90: *> \verbatim
! 91: *> B is COMPLEX*16 array, dimension (LDB,NRHS)
! 92: *> On entry, the right hand side matrix B.
! 93: *> On exit, the solution matrix X.
! 94: *> \endverbatim
! 95: *>
! 96: *> \param[in] LDB
! 97: *> \verbatim
! 98: *> LDB is INTEGER
! 99: *> The leading dimension of the array B. LDB >= max(1,N).
! 100: *> \endverbatim
! 101: *>
! 102: *> \param[out] WORK
! 103: *> \verbatim
! 104: *> WORK is REAL array, dimension (N)
! 105: *> \endverbatim
! 106: *>
! 107: *> \param[out] INFO
! 108: *> \verbatim
! 109: *> INFO is INTEGER
! 110: *> = 0: successful exit
! 111: *> < 0: if INFO = -i, the i-th argument had an illegal value
! 112: *> \endverbatim
! 113: *
! 114: * Authors:
! 115: * ========
! 116: *
! 117: *> \author Univ. of Tennessee
! 118: *> \author Univ. of California Berkeley
! 119: *> \author Univ. of Colorado Denver
! 120: *> \author NAG Ltd.
! 121: *
! 122: *> \date November 2011
! 123: *
! 124: *> \ingroup complex16HEcomputational
! 125: *
! 126: * =====================================================================
1.1 bertrand 127: SUBROUTINE ZHETRS2( UPLO, N, NRHS, A, LDA, IPIV, B, LDB,
128: $ WORK, INFO )
129: *
1.4 ! bertrand 130: * -- LAPACK computational routine (version 3.4.0) --
1.1 bertrand 131: * -- LAPACK is a software package provided by Univ. of Tennessee, --
132: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.4 ! bertrand 133: * November 2011
1.1 bertrand 134: *
135: * .. Scalar Arguments ..
136: CHARACTER UPLO
137: INTEGER INFO, LDA, LDB, N, NRHS
138: * ..
139: * .. Array Arguments ..
140: INTEGER IPIV( * )
1.4 ! bertrand 141: COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * )
1.1 bertrand 142: * ..
143: *
144: * =====================================================================
145: *
146: * .. Parameters ..
1.4 ! bertrand 147: COMPLEX*16 ONE
1.1 bertrand 148: PARAMETER ( ONE = (1.0D+0,0.0D+0) )
149: * ..
150: * .. Local Scalars ..
151: LOGICAL UPPER
152: INTEGER I, IINFO, J, K, KP
153: DOUBLE PRECISION S
1.4 ! bertrand 154: COMPLEX*16 AK, AKM1, AKM1K, BK, BKM1, DENOM
1.1 bertrand 155: * ..
156: * .. External Functions ..
157: LOGICAL LSAME
158: EXTERNAL LSAME
159: * ..
160: * .. External Subroutines ..
161: EXTERNAL ZLACGV, ZSCAL, ZSYCONV, ZSWAP, ZTRSM, XERBLA
162: * ..
163: * .. Intrinsic Functions ..
164: INTRINSIC DBLE, DCONJG, MAX
165: * ..
166: * .. Executable Statements ..
167: *
168: INFO = 0
169: UPPER = LSAME( UPLO, 'U' )
170: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
171: INFO = -1
172: ELSE IF( N.LT.0 ) THEN
173: INFO = -2
174: ELSE IF( NRHS.LT.0 ) THEN
175: INFO = -3
176: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
177: INFO = -5
178: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
179: INFO = -8
180: END IF
181: IF( INFO.NE.0 ) THEN
182: CALL XERBLA( 'ZHETRS2', -INFO )
183: RETURN
184: END IF
185: *
186: * Quick return if possible
187: *
188: IF( N.EQ.0 .OR. NRHS.EQ.0 )
189: $ RETURN
190: *
191: * Convert A
192: *
193: CALL ZSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO )
194: *
195: IF( UPPER ) THEN
196: *
1.3 bertrand 197: * Solve A*X = B, where A = U*D*U**H.
1.1 bertrand 198: *
1.3 bertrand 199: * P**T * B
1.1 bertrand 200: K=N
201: DO WHILE ( K .GE. 1 )
202: IF( IPIV( K ).GT.0 ) THEN
203: * 1 x 1 diagonal block
204: * Interchange rows K and IPIV(K).
205: KP = IPIV( K )
206: IF( KP.NE.K )
207: $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
208: K=K-1
209: ELSE
210: * 2 x 2 diagonal block
211: * Interchange rows K-1 and -IPIV(K).
212: KP = -IPIV( K )
213: IF( KP.EQ.-IPIV( K-1 ) )
214: $ CALL ZSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB )
215: K=K-2
216: END IF
217: END DO
218: *
1.3 bertrand 219: * Compute (U \P**T * B) -> B [ (U \P**T * B) ]
1.1 bertrand 220: *
1.3 bertrand 221: CALL ZTRSM('L','U','N','U',N,NRHS,ONE,A,LDA,B,LDB)
1.1 bertrand 222: *
1.3 bertrand 223: * Compute D \ B -> B [ D \ (U \P**T * B) ]
1.1 bertrand 224: *
225: I=N
226: DO WHILE ( I .GE. 1 )
227: IF( IPIV(I) .GT. 0 ) THEN
228: S = DBLE( ONE ) / DBLE( A( I, I ) )
229: CALL ZDSCAL( NRHS, S, B( I, 1 ), LDB )
230: ELSEIF ( I .GT. 1) THEN
231: IF ( IPIV(I-1) .EQ. IPIV(I) ) THEN
232: AKM1K = WORK(I)
233: AKM1 = A( I-1, I-1 ) / AKM1K
234: AK = A( I, I ) / DCONJG( AKM1K )
235: DENOM = AKM1*AK - ONE
236: DO 15 J = 1, NRHS
237: BKM1 = B( I-1, J ) / AKM1K
238: BK = B( I, J ) / DCONJG( AKM1K )
239: B( I-1, J ) = ( AK*BKM1-BK ) / DENOM
240: B( I, J ) = ( AKM1*BK-BKM1 ) / DENOM
241: 15 CONTINUE
242: I = I - 1
243: ENDIF
244: ENDIF
245: I = I - 1
246: END DO
247: *
1.3 bertrand 248: * Compute (U**H \ B) -> B [ U**H \ (D \ (U \P**T * B) ) ]
1.1 bertrand 249: *
1.3 bertrand 250: CALL ZTRSM('L','U','C','U',N,NRHS,ONE,A,LDA,B,LDB)
1.1 bertrand 251: *
1.3 bertrand 252: * P * B [ P * (U**H \ (D \ (U \P**T * B) )) ]
1.1 bertrand 253: *
254: K=1
255: DO WHILE ( K .LE. N )
256: IF( IPIV( K ).GT.0 ) THEN
257: * 1 x 1 diagonal block
258: * Interchange rows K and IPIV(K).
259: KP = IPIV( K )
260: IF( KP.NE.K )
261: $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
262: K=K+1
263: ELSE
264: * 2 x 2 diagonal block
265: * Interchange rows K-1 and -IPIV(K).
266: KP = -IPIV( K )
267: IF( K .LT. N .AND. KP.EQ.-IPIV( K+1 ) )
268: $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
269: K=K+2
270: ENDIF
271: END DO
272: *
273: ELSE
274: *
1.3 bertrand 275: * Solve A*X = B, where A = L*D*L**H.
1.1 bertrand 276: *
1.3 bertrand 277: * P**T * B
1.1 bertrand 278: K=1
279: DO WHILE ( K .LE. N )
280: IF( IPIV( K ).GT.0 ) THEN
281: * 1 x 1 diagonal block
282: * Interchange rows K and IPIV(K).
283: KP = IPIV( K )
284: IF( KP.NE.K )
285: $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
286: K=K+1
287: ELSE
288: * 2 x 2 diagonal block
289: * Interchange rows K and -IPIV(K+1).
290: KP = -IPIV( K+1 )
291: IF( KP.EQ.-IPIV( K ) )
292: $ CALL ZSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB )
293: K=K+2
294: ENDIF
295: END DO
296: *
1.3 bertrand 297: * Compute (L \P**T * B) -> B [ (L \P**T * B) ]
1.1 bertrand 298: *
1.3 bertrand 299: CALL ZTRSM('L','L','N','U',N,NRHS,ONE,A,LDA,B,LDB)
1.1 bertrand 300: *
1.3 bertrand 301: * Compute D \ B -> B [ D \ (L \P**T * B) ]
1.1 bertrand 302: *
303: I=1
304: DO WHILE ( I .LE. N )
305: IF( IPIV(I) .GT. 0 ) THEN
306: S = DBLE( ONE ) / DBLE( A( I, I ) )
307: CALL ZDSCAL( NRHS, S, B( I, 1 ), LDB )
308: ELSE
309: AKM1K = WORK(I)
310: AKM1 = A( I, I ) / DCONJG( AKM1K )
311: AK = A( I+1, I+1 ) / AKM1K
312: DENOM = AKM1*AK - ONE
313: DO 25 J = 1, NRHS
314: BKM1 = B( I, J ) / DCONJG( AKM1K )
315: BK = B( I+1, J ) / AKM1K
316: B( I, J ) = ( AK*BKM1-BK ) / DENOM
317: B( I+1, J ) = ( AKM1*BK-BKM1 ) / DENOM
318: 25 CONTINUE
319: I = I + 1
320: ENDIF
321: I = I + 1
322: END DO
323: *
1.3 bertrand 324: * Compute (L**H \ B) -> B [ L**H \ (D \ (L \P**T * B) ) ]
1.1 bertrand 325: *
1.3 bertrand 326: CALL ZTRSM('L','L','C','U',N,NRHS,ONE,A,LDA,B,LDB)
1.1 bertrand 327: *
1.3 bertrand 328: * P * B [ P * (L**H \ (D \ (L \P**T * B) )) ]
1.1 bertrand 329: *
330: K=N
331: DO WHILE ( K .GE. 1 )
332: IF( IPIV( K ).GT.0 ) THEN
333: * 1 x 1 diagonal block
334: * Interchange rows K and IPIV(K).
335: KP = IPIV( K )
336: IF( KP.NE.K )
337: $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
338: K=K-1
339: ELSE
340: * 2 x 2 diagonal block
341: * Interchange rows K-1 and -IPIV(K).
342: KP = -IPIV( K )
343: IF( K.GT.1 .AND. KP.EQ.-IPIV( K-1 ) )
344: $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
345: K=K-2
346: ENDIF
347: END DO
348: *
349: END IF
350: *
351: * Revert A
352: *
353: CALL ZSYCONV( UPLO, 'R', N, A, LDA, IPIV, WORK, IINFO )
354: *
355: RETURN
356: *
357: * End of ZHETRS2
358: *
359: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zhetrs2.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack