LAPACK: SRC/dsytrs2.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DSYTRS2( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, 
00002      $                    WORK, INFO )
00003 *
00004 *  -- LAPACK PROTOTYPE routine (version 3.3.0) --
00005 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00006 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00007 *     November 2010
00008 *
00009 *  -- Written by Julie Langou of the Univ. of TN    --
00010 *
00011 *     .. Scalar Arguments ..
00012       CHARACTER          UPLO
00013       INTEGER            INFO, LDA, LDB, N, NRHS
00014 *     ..
00015 *     .. Array Arguments ..
00016       INTEGER            IPIV( * )
00017       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), WORK( * )
00018 *     ..
00019 *
00020 *  Purpose
00021 *  =======
00022 *
00023 *  DSYTRS2 solves a system of linear equations A*X = B with a real
00024 *  symmetric matrix A using the factorization A = U*D*U**T or
00025 *  A = L*D*L**T computed by DSYTRF and converted by DSYCONV.
00026 *
00027 *  Arguments
00028 *  =========
00029 *
00030 *  UPLO    (input) CHARACTER*1
00031 *          Specifies whether the details of the factorization are stored
00032 *          as an upper or lower triangular matrix.
00033 *          = 'U':  Upper triangular, form is A = U*D*U**T;
00034 *          = 'L':  Lower triangular, form is A = L*D*L**T.
00035 *
00036 *  N       (input) INTEGER
00037 *          The order of the matrix A.  N >= 0.
00038 *
00039 *  NRHS    (input) INTEGER
00040 *          The number of right hand sides, i.e., the number of columns
00041 *          of the matrix B.  NRHS >= 0.
00042 *
00043 *  A       (input) DOUBLE PRECISION array, dimension (LDA,N)
00044 *          The block diagonal matrix D and the multipliers used to
00045 *          obtain the factor U or L as computed by DSYTRF.
00046 *
00047 *  LDA     (input) INTEGER
00048 *          The leading dimension of the array A.  LDA >= max(1,N).
00049 *
00050 *  IPIV    (input) INTEGER array, dimension (N)
00051 *          Details of the interchanges and the block structure of D
00052 *          as determined by DSYTRF.
00053 *
00054 *  B       (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
00055 *          On entry, the right hand side matrix B.
00056 *          On exit, the solution matrix X.
00057 *
00058 *  LDB     (input) INTEGER
00059 *          The leading dimension of the array B.  LDB >= max(1,N).
00060 *
00061 *  WORK    (workspace) REAL array, dimension (N)
00062 *
00063 *  INFO    (output) INTEGER
00064 *          = 0:  successful exit
00065 *          < 0:  if INFO = -i, the i-th argument had an illegal value
00066 *
00067 *  =====================================================================
00068 *
00069 *     .. Parameters ..
00070       DOUBLE PRECISION   ONE
00071       PARAMETER          ( ONE = 1.0D+0 )
00072 *     ..
00073 *     .. Local Scalars ..
00074       LOGICAL            UPPER
00075       INTEGER            I, IINFO, J, K, KP
00076       DOUBLE PRECISION   AK, AKM1, AKM1K, BK, BKM1, DENOM
00077 *     ..
00078 *     .. External Functions ..
00079       LOGICAL            LSAME
00080       EXTERNAL           LSAME
00081 *     ..
00082 *     .. External Subroutines ..
00083       EXTERNAL           DSCAL, DSYCONV, DSWAP, DTRSM, XERBLA
00084 *     ..
00085 *     .. Intrinsic Functions ..
00086       INTRINSIC          MAX
00087 *     ..
00088 *     .. Executable Statements ..
00089 *
00090       INFO = 0
00091       UPPER = LSAME( UPLO, 'U' )
00092       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00093          INFO = -1
00094       ELSE IF( N.LT.0 ) THEN
00095          INFO = -2
00096       ELSE IF( NRHS.LT.0 ) THEN
00097          INFO = -3
00098       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00099          INFO = -5
00100       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00101          INFO = -8
00102       END IF
00103       IF( INFO.NE.0 ) THEN
00104          CALL XERBLA( 'DSYTRS2', -INFO )
00105          RETURN
00106       END IF
00107 *
00108 *     Quick return if possible
00109 *
00110       IF( N.EQ.0 .OR. NRHS.EQ.0 )
00111      $   RETURN
00112 *
00113 *     Convert A
00114 *
00115       CALL DSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO )
00116 *
00117       IF( UPPER ) THEN
00118 *
00119 *        Solve A*X = B, where A = U*D*U'.
00120 *
00121 *       P' * B  
00122         K=N
00123         DO WHILE ( K .GE. 1 )
00124          IF( IPIV( K ).GT.0 ) THEN
00125 *           1 x 1 diagonal block
00126 *           Interchange rows K and IPIV(K).
00127             KP = IPIV( K )
00128             IF( KP.NE.K )
00129      $         CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00130             K=K-1
00131          ELSE
00132 *           2 x 2 diagonal block
00133 *           Interchange rows K-1 and -IPIV(K).
00134             KP = -IPIV( K )
00135             IF( KP.EQ.-IPIV( K-1 ) )
00136      $         CALL DSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB )
00137             K=K-2
00138          END IF
00139         END DO
00140 *
00141 *  Compute (U \P' * B) -> B    [ (U \P' * B) ]
00142 *
00143         CALL DTRSM('L','U','N','U',N,NRHS,ONE,A,N,B,N)
00144 *
00145 *  Compute D \ B -> B   [ D \ (U \P' * B) ]
00146 *       
00147          I=N
00148          DO WHILE ( I .GE. 1 )
00149             IF( IPIV(I) .GT. 0 ) THEN
00150               CALL DSCAL( NRHS, ONE / A( I, I ), B( I, 1 ), N )
00151             ELSEIF ( I .GT. 1) THEN
00152                IF ( IPIV(I-1) .EQ. IPIV(I) ) THEN
00153                   AKM1K = WORK(I)
00154                   AKM1 = A( I-1, I-1 ) / AKM1K
00155                   AK = A( I, I ) / AKM1K
00156                   DENOM = AKM1*AK - ONE
00157                   DO 15 J = 1, NRHS
00158                      BKM1 = B( I-1, J ) / AKM1K
00159                      BK = B( I, J ) / AKM1K
00160                      B( I-1, J ) = ( AK*BKM1-BK ) / DENOM
00161                      B( I, J ) = ( AKM1*BK-BKM1 ) / DENOM
00162  15              CONTINUE
00163                I = I - 1
00164                ENDIF
00165             ENDIF
00166             I = I - 1
00167          END DO
00168 *
00169 *      Compute (U' \ B) -> B   [ U' \ (D \ (U \P' * B) ) ]
00170 *
00171          CALL DTRSM('L','U','T','U',N,NRHS,ONE,A,N,B,N)
00172 *
00173 *       P * B  [ P * (U' \ (D \ (U \P' * B) )) ]
00174 *
00175         K=1
00176         DO WHILE ( K .LE. N )
00177          IF( IPIV( K ).GT.0 ) THEN
00178 *           1 x 1 diagonal block
00179 *           Interchange rows K and IPIV(K).
00180             KP = IPIV( K )
00181             IF( KP.NE.K )
00182      $         CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00183             K=K+1
00184          ELSE
00185 *           2 x 2 diagonal block
00186 *           Interchange rows K-1 and -IPIV(K).
00187             KP = -IPIV( K )
00188             IF( K .LT. N .AND. KP.EQ.-IPIV( K+1 ) )
00189      $         CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00190             K=K+2
00191          ENDIF
00192         END DO
00193 *
00194       ELSE
00195 *
00196 *        Solve A*X = B, where A = L*D*L'.
00197 *
00198 *       P' * B  
00199         K=1
00200         DO WHILE ( K .LE. N )
00201          IF( IPIV( K ).GT.0 ) THEN
00202 *           1 x 1 diagonal block
00203 *           Interchange rows K and IPIV(K).
00204             KP = IPIV( K )
00205             IF( KP.NE.K )
00206      $         CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00207             K=K+1
00208          ELSE
00209 *           2 x 2 diagonal block
00210 *           Interchange rows K and -IPIV(K+1).
00211             KP = -IPIV( K+1 )
00212             IF( KP.EQ.-IPIV( K ) )
00213      $         CALL DSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB )
00214             K=K+2
00215          ENDIF
00216         END DO
00217 *
00218 *  Compute (L \P' * B) -> B    [ (L \P' * B) ]
00219 *
00220         CALL DTRSM('L','L','N','U',N,NRHS,ONE,A,N,B,N)
00221 *
00222 *  Compute D \ B -> B   [ D \ (L \P' * B) ]
00223 *       
00224          I=1
00225          DO WHILE ( I .LE. N )
00226             IF( IPIV(I) .GT. 0 ) THEN
00227               CALL DSCAL( NRHS, ONE / A( I, I ), B( I, 1 ), N )
00228             ELSE
00229                   AKM1K = WORK(I)
00230                   AKM1 = A( I, I ) / AKM1K
00231                   AK = A( I+1, I+1 ) / AKM1K
00232                   DENOM = AKM1*AK - ONE
00233                   DO 25 J = 1, NRHS
00234                      BKM1 = B( I, J ) / AKM1K
00235                      BK = B( I+1, J ) / AKM1K
00236                      B( I, J ) = ( AK*BKM1-BK ) / DENOM
00237                      B( I+1, J ) = ( AKM1*BK-BKM1 ) / DENOM
00238  25              CONTINUE
00239                   I = I + 1
00240             ENDIF
00241             I = I + 1
00242          END DO
00243 *
00244 *  Compute (L' \ B) -> B   [ L' \ (D \ (L \P' * B) ) ]
00245 * 
00246         CALL DTRSM('L','L','T','U',N,NRHS,ONE,A,N,B,N)
00247 *
00248 *       P * B  [ P * (L' \ (D \ (L \P' * B) )) ]
00249 *
00250         K=N
00251         DO WHILE ( K .GE. 1 )
00252          IF( IPIV( K ).GT.0 ) THEN
00253 *           1 x 1 diagonal block
00254 *           Interchange rows K and IPIV(K).
00255             KP = IPIV( K )
00256             IF( KP.NE.K )
00257      $         CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00258             K=K-1
00259          ELSE
00260 *           2 x 2 diagonal block
00261 *           Interchange rows K-1 and -IPIV(K).
00262             KP = -IPIV( K )
00263             IF( K.GT.1 .AND. KP.EQ.-IPIV( K-1 ) )
00264      $         CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00265             K=K-2
00266          ENDIF
00267         END DO
00268 *
00269       END IF
00270 *
00271 *     Revert A
00272 *
00273       CALL DSYCONV( UPLO, 'R', N, A, LDA, IPIV, WORK, IINFO )
00274 *
00275       RETURN
00276 *
00277 *     End of DSYTRS2
00278 *
00279       END