LAPACK: SRC/csytrs2.f Source File

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