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