![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE ZSYTRS2( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, 2: $ WORK, INFO ) 3: * 4: * -- LAPACK PROTOTYPE routine (version 3.2.2) -- 5: * 6: * -- Written by Julie Langou of the Univ. of TN -- 7: * May 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: * ZSYTRS2 solves a system of linear equations A*X = B with a real 25: * symmetric 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**T; 35: * = 'L': Lower triangular, form is A = L*D*L**T. 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 ZSYTRF. 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 ZSYTRF. 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 COMPLEX AK, AKM1, AKM1K, BK, BKM1, DENOM 78: * .. 79: * .. External Functions .. 80: LOGICAL LSAME 81: EXTERNAL LSAME 82: * .. 83: * .. External Subroutines .. 84: EXTERNAL ZSCAL, ZSYCONV, ZSWAP, ZTRSM, XERBLA 85: * .. 86: * .. Intrinsic Functions .. 87: INTRINSIC MAX 88: * .. 89: * .. Executable Statements .. 90: * 91: INFO = 0 92: UPPER = LSAME( UPLO, 'U' ) 93: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 94: INFO = -1 95: ELSE IF( N.LT.0 ) THEN 96: INFO = -2 97: ELSE IF( NRHS.LT.0 ) THEN 98: INFO = -3 99: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 100: INFO = -5 101: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 102: INFO = -8 103: END IF 104: IF( INFO.NE.0 ) THEN 105: CALL XERBLA( 'ZSYTRS2', -INFO ) 106: RETURN 107: END IF 108: * 109: * Quick return if possible 110: * 111: IF( N.EQ.0 .OR. NRHS.EQ.0 ) 112: $ RETURN 113: * 114: * Convert A 115: * 116: CALL ZSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO ) 117: * 118: IF( UPPER ) THEN 119: * 120: * Solve A*X = B, where A = U*D*U'. 121: * 122: * P' * B 123: K=N 124: DO WHILE ( K .GE. 1 ) 125: IF( IPIV( K ).GT.0 ) THEN 126: * 1 x 1 diagonal block 127: * Interchange rows K and IPIV(K). 128: KP = IPIV( K ) 129: IF( KP.NE.K ) 130: $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 131: K=K-1 132: ELSE 133: * 2 x 2 diagonal block 134: * Interchange rows K-1 and -IPIV(K). 135: KP = -IPIV( K ) 136: IF( KP.EQ.-IPIV( K-1 ) ) 137: $ CALL ZSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB ) 138: K=K-2 139: END IF 140: END DO 141: * 142: * Compute (U \P' * B) -> B [ (U \P' * B) ] 143: * 144: CALL ZTRSM('L','U','N','U',N,NRHS,ONE,A,N,B,N) 145: * 146: * Compute D \ B -> B [ D \ (U \P' * B) ] 147: * 148: I=N 149: DO WHILE ( I .GE. 1 ) 150: IF( IPIV(I) .GT. 0 ) THEN 151: CALL ZSCAL( NRHS, ONE / A( I, I ), B( I, 1 ), N ) 152: ELSEIF ( I .GT. 1) THEN 153: IF ( IPIV(I-1) .EQ. IPIV(I) ) THEN 154: AKM1K = WORK(I) 155: AKM1 = A( I-1, I-1 ) / AKM1K 156: AK = A( I, I ) / AKM1K 157: DENOM = AKM1*AK - ONE 158: DO 15 J = 1, NRHS 159: BKM1 = B( I-1, J ) / AKM1K 160: BK = B( I, J ) / AKM1K 161: B( I-1, J ) = ( AK*BKM1-BK ) / DENOM 162: B( I, J ) = ( AKM1*BK-BKM1 ) / DENOM 163: 15 CONTINUE 164: I = I - 1 165: ENDIF 166: ENDIF 167: I = I - 1 168: END DO 169: * 170: * Compute (U' \ B) -> B [ U' \ (D \ (U \P' * B) ) ] 171: * 172: CALL ZTRSM('L','U','T','U',N,NRHS,ONE,A,N,B,N) 173: * 174: * P * B [ P * (U' \ (D \ (U \P' * B) )) ] 175: * 176: K=1 177: DO WHILE ( K .LE. N ) 178: IF( IPIV( K ).GT.0 ) THEN 179: * 1 x 1 diagonal block 180: * Interchange rows K and IPIV(K). 181: KP = IPIV( K ) 182: IF( KP.NE.K ) 183: $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 184: K=K+1 185: ELSE 186: * 2 x 2 diagonal block 187: * Interchange rows K-1 and -IPIV(K). 188: KP = -IPIV( K ) 189: IF( K .LT. N .AND. KP.EQ.-IPIV( K+1 ) ) 190: $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 191: K=K+2 192: ENDIF 193: END DO 194: * 195: ELSE 196: * 197: * Solve A*X = B, where A = L*D*L'. 198: * 199: * P' * B 200: K=1 201: DO WHILE ( K .LE. N ) 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 and -IPIV(K+1). 212: KP = -IPIV( K+1 ) 213: IF( KP.EQ.-IPIV( K ) ) 214: $ CALL ZSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB ) 215: K=K+2 216: ENDIF 217: END DO 218: * 219: * Compute (L \P' * B) -> B [ (L \P' * B) ] 220: * 221: CALL ZTRSM('L','L','N','U',N,NRHS,ONE,A,N,B,N) 222: * 223: * Compute D \ B -> B [ D \ (L \P' * B) ] 224: * 225: I=1 226: DO WHILE ( I .LE. N ) 227: IF( IPIV(I) .GT. 0 ) THEN 228: CALL ZSCAL( NRHS, ONE / A( I, I ), B( I, 1 ), N ) 229: ELSE 230: AKM1K = WORK(I) 231: AKM1 = A( I, I ) / AKM1K 232: AK = A( I+1, I+1 ) / AKM1K 233: DENOM = AKM1*AK - ONE 234: DO 25 J = 1, NRHS 235: BKM1 = B( I, J ) / AKM1K 236: BK = B( I+1, J ) / AKM1K 237: B( I, J ) = ( AK*BKM1-BK ) / DENOM 238: B( I+1, J ) = ( AKM1*BK-BKM1 ) / DENOM 239: 25 CONTINUE 240: I = I + 1 241: ENDIF 242: I = I + 1 243: END DO 244: * 245: * Compute (L' \ B) -> B [ L' \ (D \ (L \P' * B) ) ] 246: * 247: CALL ZTRSM('L','L','T','U',N,NRHS,ONE,A,N,B,N) 248: * 249: * P * B [ P * (L' \ (D \ (L \P' * B) )) ] 250: * 251: K=N 252: DO WHILE ( K .GE. 1 ) 253: IF( IPIV( K ).GT.0 ) THEN 254: * 1 x 1 diagonal block 255: * Interchange rows K and IPIV(K). 256: KP = IPIV( K ) 257: IF( KP.NE.K ) 258: $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 259: K=K-1 260: ELSE 261: * 2 x 2 diagonal block 262: * Interchange rows K-1 and -IPIV(K). 263: KP = -IPIV( K ) 264: IF( K.GT.1 .AND. KP.EQ.-IPIV( K-1 ) ) 265: $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 266: K=K-2 267: ENDIF 268: END DO 269: * 270: END IF 271: * 272: * Revert A 273: * 274: CALL ZSYCONV( UPLO, 'R', N, A, LDA, IPIV, WORK, IINFO ) 275: * 276: RETURN 277: * 278: * End of ZSYTRS2 279: * 280: END