![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE ZHETRS2( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, 2: $ WORK, INFO ) 3: * 4: * -- LAPACK PROTOTYPE routine (version 3.3.0) -- 5: * 6: * -- Written by Julie Langou of the Univ. of TN -- 7: * November 2010 8: * 9: * -- LAPACK is a software package provided by Univ. of Tennessee, -- 10: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 11: * 12: * .. Scalar Arguments .. 13: CHARACTER UPLO 14: INTEGER INFO, LDA, LDB, N, NRHS 15: * .. 16: * .. Array Arguments .. 17: INTEGER IPIV( * ) 18: DOUBLE COMPLEX A( LDA, * ), B( LDB, * ), WORK( * ) 19: * .. 20: * 21: * Purpose 22: * ======= 23: * 24: * ZHETRS2 solves a system of linear equations A*X = B with a real 25: * Hermitian matrix A using the factorization A = U*D*U**T or 26: * A = L*D*L**T computed by ZSYTRF and converted by ZSYCONV. 27: * 28: * Arguments 29: * ========= 30: * 31: * UPLO (input) CHARACTER*1 32: * Specifies whether the details of the factorization are stored 33: * as an upper or lower triangular matrix. 34: * = 'U': Upper triangular, form is A = U*D*U**H; 35: * = 'L': Lower triangular, form is A = L*D*L**H. 36: * 37: * N (input) INTEGER 38: * The order of the matrix A. N >= 0. 39: * 40: * NRHS (input) INTEGER 41: * The number of right hand sides, i.e., the number of columns 42: * of the matrix B. NRHS >= 0. 43: * 44: * A (input) DOUBLE COMPLEX array, dimension (LDA,N) 45: * The block diagonal matrix D and the multipliers used to 46: * obtain the factor U or L as computed by ZHETRF. 47: * 48: * LDA (input) INTEGER 49: * The leading dimension of the array A. LDA >= max(1,N). 50: * 51: * IPIV (input) INTEGER array, dimension (N) 52: * Details of the interchanges and the block structure of D 53: * as determined by ZHETRF. 54: * 55: * B (input/output) DOUBLE COMPLEX array, dimension (LDB,NRHS) 56: * On entry, the right hand side matrix B. 57: * On exit, the solution matrix X. 58: * 59: * LDB (input) INTEGER 60: * The leading dimension of the array B. LDB >= max(1,N). 61: * 62: * WORK (workspace) REAL array, dimension (N) 63: * 64: * INFO (output) INTEGER 65: * = 0: successful exit 66: * < 0: if INFO = -i, the i-th argument had an illegal value 67: * 68: * ===================================================================== 69: * 70: * .. Parameters .. 71: DOUBLE COMPLEX ONE 72: PARAMETER ( ONE = (1.0D+0,0.0D+0) ) 73: * .. 74: * .. Local Scalars .. 75: LOGICAL UPPER 76: INTEGER I, IINFO, J, K, KP 77: DOUBLE PRECISION S 78: DOUBLE COMPLEX AK, AKM1, AKM1K, BK, BKM1, DENOM 79: * .. 80: * .. External Functions .. 81: LOGICAL LSAME 82: EXTERNAL LSAME 83: * .. 84: * .. External Subroutines .. 85: EXTERNAL ZLACGV, ZSCAL, ZSYCONV, ZSWAP, ZTRSM, XERBLA 86: * .. 87: * .. Intrinsic Functions .. 88: INTRINSIC DBLE, DCONJG, MAX 89: * .. 90: * .. Executable Statements .. 91: * 92: INFO = 0 93: UPPER = LSAME( UPLO, 'U' ) 94: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 95: INFO = -1 96: ELSE IF( N.LT.0 ) THEN 97: INFO = -2 98: ELSE IF( NRHS.LT.0 ) THEN 99: INFO = -3 100: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 101: INFO = -5 102: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 103: INFO = -8 104: END IF 105: IF( INFO.NE.0 ) THEN 106: CALL XERBLA( 'ZHETRS2', -INFO ) 107: RETURN 108: END IF 109: * 110: * Quick return if possible 111: * 112: IF( N.EQ.0 .OR. NRHS.EQ.0 ) 113: $ RETURN 114: * 115: * Convert A 116: * 117: CALL ZSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO ) 118: * 119: IF( UPPER ) THEN 120: * 121: * Solve A*X = B, where A = U*D*U'. 122: * 123: * P' * B 124: K=N 125: DO WHILE ( K .GE. 1 ) 126: IF( IPIV( K ).GT.0 ) THEN 127: * 1 x 1 diagonal block 128: * Interchange rows K and IPIV(K). 129: KP = IPIV( K ) 130: IF( KP.NE.K ) 131: $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 132: K=K-1 133: ELSE 134: * 2 x 2 diagonal block 135: * Interchange rows K-1 and -IPIV(K). 136: KP = -IPIV( K ) 137: IF( KP.EQ.-IPIV( K-1 ) ) 138: $ CALL ZSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB ) 139: K=K-2 140: END IF 141: END DO 142: * 143: * Compute (U \P' * B) -> B [ (U \P' * B) ] 144: * 145: CALL ZTRSM('L','U','N','U',N,NRHS,ONE,A,N,B,N) 146: * 147: * Compute D \ B -> B [ D \ (U \P' * B) ] 148: * 149: I=N 150: DO WHILE ( I .GE. 1 ) 151: IF( IPIV(I) .GT. 0 ) THEN 152: S = DBLE( ONE ) / DBLE( A( I, I ) ) 153: CALL ZDSCAL( NRHS, S, B( I, 1 ), LDB ) 154: ELSEIF ( I .GT. 1) THEN 155: IF ( IPIV(I-1) .EQ. IPIV(I) ) THEN 156: AKM1K = WORK(I) 157: AKM1 = A( I-1, I-1 ) / AKM1K 158: AK = A( I, I ) / DCONJG( AKM1K ) 159: DENOM = AKM1*AK - ONE 160: DO 15 J = 1, NRHS 161: BKM1 = B( I-1, J ) / AKM1K 162: BK = B( I, J ) / DCONJG( AKM1K ) 163: B( I-1, J ) = ( AK*BKM1-BK ) / DENOM 164: B( I, J ) = ( AKM1*BK-BKM1 ) / DENOM 165: 15 CONTINUE 166: I = I - 1 167: ENDIF 168: ENDIF 169: I = I - 1 170: END DO 171: * 172: * Compute (U' \ B) -> B [ U' \ (D \ (U \P' * B) ) ] 173: * 174: CALL ZTRSM('L','U','C','U',N,NRHS,ONE,A,N,B,N) 175: * 176: * P * B [ P * (U' \ (D \ (U \P' * B) )) ] 177: * 178: K=1 179: DO WHILE ( K .LE. N ) 180: IF( IPIV( K ).GT.0 ) THEN 181: * 1 x 1 diagonal block 182: * Interchange rows K and IPIV(K). 183: KP = IPIV( K ) 184: IF( KP.NE.K ) 185: $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 186: K=K+1 187: ELSE 188: * 2 x 2 diagonal block 189: * Interchange rows K-1 and -IPIV(K). 190: KP = -IPIV( K ) 191: IF( K .LT. N .AND. KP.EQ.-IPIV( K+1 ) ) 192: $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 193: K=K+2 194: ENDIF 195: END DO 196: * 197: ELSE 198: * 199: * Solve A*X = B, where A = L*D*L'. 200: * 201: * P' * B 202: K=1 203: DO WHILE ( K .LE. N ) 204: IF( IPIV( K ).GT.0 ) THEN 205: * 1 x 1 diagonal block 206: * Interchange rows K and IPIV(K). 207: KP = IPIV( K ) 208: IF( KP.NE.K ) 209: $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 210: K=K+1 211: ELSE 212: * 2 x 2 diagonal block 213: * Interchange rows K and -IPIV(K+1). 214: KP = -IPIV( K+1 ) 215: IF( KP.EQ.-IPIV( K ) ) 216: $ CALL ZSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB ) 217: K=K+2 218: ENDIF 219: END DO 220: * 221: * Compute (L \P' * B) -> B [ (L \P' * B) ] 222: * 223: CALL ZTRSM('L','L','N','U',N,NRHS,ONE,A,N,B,N) 224: * 225: * Compute D \ B -> B [ D \ (L \P' * B) ] 226: * 227: I=1 228: DO WHILE ( I .LE. N ) 229: IF( IPIV(I) .GT. 0 ) THEN 230: S = DBLE( ONE ) / DBLE( A( I, I ) ) 231: CALL ZDSCAL( NRHS, S, B( I, 1 ), LDB ) 232: ELSE 233: AKM1K = WORK(I) 234: AKM1 = A( I, I ) / DCONJG( AKM1K ) 235: AK = A( I+1, I+1 ) / AKM1K 236: DENOM = AKM1*AK - ONE 237: DO 25 J = 1, NRHS 238: BKM1 = B( I, J ) / DCONJG( AKM1K ) 239: BK = B( I+1, J ) / AKM1K 240: B( I, J ) = ( AK*BKM1-BK ) / DENOM 241: B( I+1, J ) = ( AKM1*BK-BKM1 ) / DENOM 242: 25 CONTINUE 243: I = I + 1 244: ENDIF 245: I = I + 1 246: END DO 247: * 248: * Compute (L' \ B) -> B [ L' \ (D \ (L \P' * B) ) ] 249: * 250: CALL ZTRSM('L','L','C','U',N,NRHS,ONE,A,N,B,N) 251: * 252: * P * B [ P * (L' \ (D \ (L \P' * B) )) ] 253: * 254: K=N 255: DO WHILE ( K .GE. 1 ) 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.GT.1 .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: END IF 274: * 275: * Revert A 276: * 277: CALL ZSYCONV( UPLO, 'R', N, A, LDA, IPIV, WORK, IINFO ) 278: * 279: RETURN 280: * 281: * End of ZHETRS2 282: * 283: END