Annotation of rpl/lapack/lapack/dsytrs2.f, revision 1.3

1.1       bertrand    1:       SUBROUTINE DSYTRS2( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, 
                      2:      $                    WORK, INFO )
                      3: *
1.3     ! bertrand    4: *  -- LAPACK PROTOTYPE routine (version 3.3.1) --
1.1       bertrand    5: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                      6: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.3     ! bertrand    7: *  -- April 2011                                                      --
1.1       bertrand    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: *
1.3     ! bertrand  119: *        Solve A*X = B, where A = U*D*U**T.
1.1       bertrand  120: *
1.3     ! bertrand  121: *       P**T * B  
1.1       bertrand  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: *
1.3     ! bertrand  141: *  Compute (U \P**T * B) -> B    [ (U \P**T * B) ]
1.1       bertrand  142: *
1.3     ! bertrand  143:         CALL DTRSM('L','U','N','U',N,NRHS,ONE,A,LDA,B,LDB)
1.1       bertrand  144: *
1.3     ! bertrand  145: *  Compute D \ B -> B   [ D \ (U \P**T * B) ]
1.1       bertrand  146: *       
                    147:          I=N
                    148:          DO WHILE ( I .GE. 1 )
                    149:             IF( IPIV(I) .GT. 0 ) THEN
1.3     ! bertrand  150:               CALL DSCAL( NRHS, ONE / A( I, I ), B( I, 1 ), LDB )
1.1       bertrand  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: *
1.3     ! bertrand  169: *      Compute (U**T \ B) -> B   [ U**T \ (D \ (U \P**T * B) ) ]
1.1       bertrand  170: *
1.3     ! bertrand  171:          CALL DTRSM('L','U','T','U',N,NRHS,ONE,A,LDA,B,LDB)
1.1       bertrand  172: *
1.3     ! bertrand  173: *       P * B  [ P * (U**T \ (D \ (U \P**T * B) )) ]
1.1       bertrand  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: *
1.3     ! bertrand  196: *        Solve A*X = B, where A = L*D*L**T.
1.1       bertrand  197: *
1.3     ! bertrand  198: *       P**T * B  
1.1       bertrand  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: *
1.3     ! bertrand  218: *  Compute (L \P**T * B) -> B    [ (L \P**T * B) ]
1.1       bertrand  219: *
1.3     ! bertrand  220:         CALL DTRSM('L','L','N','U',N,NRHS,ONE,A,LDA,B,LDB)
1.1       bertrand  221: *
1.3     ! bertrand  222: *  Compute D \ B -> B   [ D \ (L \P**T * B) ]
1.1       bertrand  223: *       
                    224:          I=1
                    225:          DO WHILE ( I .LE. N )
                    226:             IF( IPIV(I) .GT. 0 ) THEN
1.3     ! bertrand  227:               CALL DSCAL( NRHS, ONE / A( I, I ), B( I, 1 ), LDB )
1.1       bertrand  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: *
1.3     ! bertrand  244: *  Compute (L**T \ B) -> B   [ L**T \ (D \ (L \P**T * B) ) ]
1.1       bertrand  245: * 
1.3     ! bertrand  246:         CALL DTRSM('L','L','T','U',N,NRHS,ONE,A,LDA,B,LDB)
1.1       bertrand  247: *
1.3     ! bertrand  248: *       P * B  [ P * (L**T \ (D \ (L \P**T * B) )) ]
1.1       bertrand  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>