LAPACK: SRC/chetrs2.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CHETRS2( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, 
00002      $                    WORK, INFO )
00003 *
00004 *  -- LAPACK PROTOTYPE routine (version 3.3.0) --
00005 *
00006 *  -- Written by Julie Langou of the Univ. of TN    --
00007 *     November 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 *  CHETRS2 solves a system of linear equations A*X = B with a COMPLEX
00025 *  Hermitian 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**H;
00035 *          = 'L':  Lower triangular, form is A = L*D*L**H.
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 CHETRF.
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 CHETRF.
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       REAL               S
00078       COMPLEX            AK, AKM1, AKM1K, BK, BKM1, DENOM
00079 *     ..
00080 *     .. External Functions ..
00081       LOGICAL            LSAME
00082       EXTERNAL           LSAME
00083 *     ..
00084 *     .. External Subroutines ..
00085       EXTERNAL           CSCAL, CSYCONV, CSWAP, CTRSM, XERBLA
00086 *     ..
00087 *     .. Intrinsic Functions ..
00088       INTRINSIC          CONJG, MAX, REAL
00089 *     ..
00090 *     .. Executable Statements ..
00091 *
00092       INFO = 0
00093       UPPER = LSAME( UPLO, 'U' )
00094       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00095          INFO = -1
00096       ELSE IF( N.LT.0 ) THEN
00097          INFO = -2
00098       ELSE IF( NRHS.LT.0 ) THEN
00099          INFO = -3
00100       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00101          INFO = -5
00102       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00103          INFO = -8
00104       END IF
00105       IF( INFO.NE.0 ) THEN
00106          CALL XERBLA( 'CHETRS2', -INFO )
00107          RETURN
00108       END IF
00109 *
00110 *     Quick return if possible
00111 *
00112       IF( N.EQ.0 .OR. NRHS.EQ.0 )
00113      $   RETURN
00114 *
00115 *     Convert A
00116 *
00117       CALL CSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO )
00118 *
00119       IF( UPPER ) THEN
00120 *
00121 *        Solve A*X = B, where A = U*D*U'.
00122 *
00123 *       P' * B  
00124         K=N
00125         DO WHILE ( K .GE. 1 )
00126          IF( IPIV( K ).GT.0 ) THEN
00127 *           1 x 1 diagonal block
00128 *           Interchange rows K and IPIV(K).
00129             KP = IPIV( K )
00130             IF( KP.NE.K )
00131      $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00132             K=K-1
00133          ELSE
00134 *           2 x 2 diagonal block
00135 *           Interchange rows K-1 and -IPIV(K).
00136             KP = -IPIV( K )
00137             IF( KP.EQ.-IPIV( K-1 ) )
00138      $         CALL CSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB )
00139             K=K-2
00140          END IF
00141         END DO
00142 *
00143 *  Compute (U \P' * B) -> B    [ (U \P' * B) ]
00144 *
00145         CALL CTRSM('L','U','N','U',N,NRHS,ONE,A,N,B,N)
00146 *
00147 *  Compute D \ B -> B   [ D \ (U \P' * B) ]
00148 *       
00149          I=N
00150          DO WHILE ( I .GE. 1 )
00151             IF( IPIV(I) .GT. 0 ) THEN
00152               S = REAL( ONE ) / REAL( A( I, I ) )
00153               CALL CSSCAL( NRHS, S, B( I, 1 ), LDB )
00154             ELSEIF ( I .GT. 1) THEN
00155                IF ( IPIV(I-1) .EQ. IPIV(I) ) THEN
00156                   AKM1K = WORK(I)
00157                   AKM1 = A( I-1, I-1 ) / AKM1K
00158                   AK = A( I, I ) / CONJG( AKM1K )
00159                   DENOM = AKM1*AK - ONE
00160                   DO 15 J = 1, NRHS
00161                      BKM1 = B( I-1, J ) / AKM1K
00162                      BK = B( I, J ) / CONJG( AKM1K )
00163                      B( I-1, J ) = ( AK*BKM1-BK ) / DENOM
00164                      B( I, J ) = ( AKM1*BK-BKM1 ) / DENOM
00165  15              CONTINUE
00166                I = I - 1
00167                ENDIF
00168             ENDIF
00169             I = I - 1
00170          END DO
00171 *
00172 *      Compute (U' \ B) -> B   [ U' \ (D \ (U \P' * B) ) ]
00173 *
00174          CALL CTRSM('L','U','C','U',N,NRHS,ONE,A,N,B,N)
00175 *
00176 *       P * B  [ P * (U' \ (D \ (U \P' * B) )) ]
00177 *
00178         K=1
00179         DO WHILE ( K .LE. N )
00180          IF( IPIV( K ).GT.0 ) THEN
00181 *           1 x 1 diagonal block
00182 *           Interchange rows K and IPIV(K).
00183             KP = IPIV( K )
00184             IF( KP.NE.K )
00185      $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00186             K=K+1
00187          ELSE
00188 *           2 x 2 diagonal block
00189 *           Interchange rows K-1 and -IPIV(K).
00190             KP = -IPIV( K )
00191             IF( K .LT. N .AND. KP.EQ.-IPIV( K+1 ) )
00192      $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00193             K=K+2
00194          ENDIF
00195         END DO
00196 *
00197       ELSE
00198 *
00199 *        Solve A*X = B, where A = L*D*L'.
00200 *
00201 *       P' * B  
00202         K=1
00203         DO WHILE ( K .LE. N )
00204          IF( IPIV( K ).GT.0 ) THEN
00205 *           1 x 1 diagonal block
00206 *           Interchange rows K and IPIV(K).
00207             KP = IPIV( K )
00208             IF( KP.NE.K )
00209      $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00210             K=K+1
00211          ELSE
00212 *           2 x 2 diagonal block
00213 *           Interchange rows K and -IPIV(K+1).
00214             KP = -IPIV( K+1 )
00215             IF( KP.EQ.-IPIV( K ) )
00216      $         CALL CSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB )
00217             K=K+2
00218          ENDIF
00219         END DO
00220 *
00221 *  Compute (L \P' * B) -> B    [ (L \P' * B) ]
00222 *
00223         CALL CTRSM('L','L','N','U',N,NRHS,ONE,A,N,B,N)
00224 *
00225 *  Compute D \ B -> B   [ D \ (L \P' * B) ]
00226 *       
00227          I=1
00228          DO WHILE ( I .LE. N )
00229             IF( IPIV(I) .GT. 0 ) THEN
00230               S = REAL( ONE ) / REAL( A( I, I ) )
00231               CALL CSSCAL( NRHS, S, B( I, 1 ), LDB )
00232             ELSE
00233                   AKM1K = WORK(I)
00234                   AKM1 = A( I, I ) / CONJG( AKM1K )
00235                   AK = A( I+1, I+1 ) / AKM1K
00236                   DENOM = AKM1*AK - ONE
00237                   DO 25 J = 1, NRHS
00238                      BKM1 = B( I, J ) / CONJG( AKM1K )
00239                      BK = B( I+1, J ) / AKM1K
00240                      B( I, J ) = ( AK*BKM1-BK ) / DENOM
00241                      B( I+1, J ) = ( AKM1*BK-BKM1 ) / DENOM
00242  25              CONTINUE
00243                   I = I + 1
00244             ENDIF
00245             I = I + 1
00246          END DO
00247 *
00248 *  Compute (L' \ B) -> B   [ L' \ (D \ (L \P' * B) ) ]
00249 * 
00250         CALL CTRSM('L','L','C','U',N,NRHS,ONE,A,N,B,N)
00251 *
00252 *       P * B  [ P * (L' \ (D \ (L \P' * B) )) ]
00253 *
00254         K=N
00255         DO WHILE ( K .GE. 1 )
00256          IF( IPIV( K ).GT.0 ) THEN
00257 *           1 x 1 diagonal block
00258 *           Interchange rows K and IPIV(K).
00259             KP = IPIV( K )
00260             IF( KP.NE.K )
00261      $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00262             K=K-1
00263          ELSE
00264 *           2 x 2 diagonal block
00265 *           Interchange rows K-1 and -IPIV(K).
00266             KP = -IPIV( K )
00267             IF( K.GT.1 .AND. KP.EQ.-IPIV( K-1 ) )
00268      $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00269             K=K-2
00270          ENDIF
00271         END DO
00272 *
00273       END IF
00274 *
00275 *     Revert A
00276 *
00277       CALL CSYCONV( UPLO, 'R', N, A, LDA, IPIV, WORK, IINFO )
00278 *
00279       RETURN
00280 *
00281 *     End of CHETRS2
00282 *
00283       END