Annotation of rpl/lapack/lapack/dsytrs2.f, revision 1.1
1.1 ! bertrand 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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dsytrs2.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack