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