LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
chetrs2.f
Go to the documentation of this file.
00001 *> \brief \b CHETRS2
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download CHETRS2 + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/chetrs2.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/chetrs2.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/chetrs2.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE CHETRS2( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, 
00022 *                           WORK, INFO )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       CHARACTER          UPLO
00026 *       INTEGER            INFO, LDA, LDB, N, NRHS
00027 *       ..
00028 *       .. Array Arguments ..
00029 *       INTEGER            IPIV( * )
00030 *       COMPLEX            A( LDA, * ), B( LDB, * ), WORK( * )
00031 *       ..
00032 *  
00033 *
00034 *> \par Purpose:
00035 *  =============
00036 *>
00037 *> \verbatim
00038 *>
00039 *> CHETRS2 solves a system of linear equations A*X = B with a complex
00040 *> Hermitian matrix A using the factorization A = U*D*U**H or
00041 *> A = L*D*L**H computed by CHETRF and converted by CSYCONV.
00042 *> \endverbatim
00043 *
00044 *  Arguments:
00045 *  ==========
00046 *
00047 *> \param[in] UPLO
00048 *> \verbatim
00049 *>          UPLO is CHARACTER*1
00050 *>          Specifies whether the details of the factorization are stored
00051 *>          as an upper or lower triangular matrix.
00052 *>          = 'U':  Upper triangular, form is A = U*D*U**H;
00053 *>          = 'L':  Lower triangular, form is A = L*D*L**H.
00054 *> \endverbatim
00055 *>
00056 *> \param[in] N
00057 *> \verbatim
00058 *>          N is INTEGER
00059 *>          The order of the matrix A.  N >= 0.
00060 *> \endverbatim
00061 *>
00062 *> \param[in] NRHS
00063 *> \verbatim
00064 *>          NRHS is INTEGER
00065 *>          The number of right hand sides, i.e., the number of columns
00066 *>          of the matrix B.  NRHS >= 0.
00067 *> \endverbatim
00068 *>
00069 *> \param[in] A
00070 *> \verbatim
00071 *>          A is COMPLEX array, dimension (LDA,N)
00072 *>          The block diagonal matrix D and the multipliers used to
00073 *>          obtain the factor U or L as computed by CHETRF.
00074 *> \endverbatim
00075 *>
00076 *> \param[in] LDA
00077 *> \verbatim
00078 *>          LDA is INTEGER
00079 *>          The leading dimension of the array A.  LDA >= max(1,N).
00080 *> \endverbatim
00081 *>
00082 *> \param[in] IPIV
00083 *> \verbatim
00084 *>          IPIV is INTEGER array, dimension (N)
00085 *>          Details of the interchanges and the block structure of D
00086 *>          as determined by CHETRF.
00087 *> \endverbatim
00088 *>
00089 *> \param[in,out] B
00090 *> \verbatim
00091 *>          B is COMPLEX array, dimension (LDB,NRHS)
00092 *>          On entry, the right hand side matrix B.
00093 *>          On exit, the solution matrix X.
00094 *> \endverbatim
00095 *>
00096 *> \param[in] LDB
00097 *> \verbatim
00098 *>          LDB is INTEGER
00099 *>          The leading dimension of the array B.  LDB >= max(1,N).
00100 *> \endverbatim
00101 *>
00102 *> \param[out] WORK
00103 *> \verbatim
00104 *>          WORK is COMPLEX array, dimension (N)
00105 *> \endverbatim
00106 *>
00107 *> \param[out] INFO
00108 *> \verbatim
00109 *>          INFO is INTEGER
00110 *>          = 0:  successful exit
00111 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
00112 *> \endverbatim
00113 *
00114 *  Authors:
00115 *  ========
00116 *
00117 *> \author Univ. of Tennessee 
00118 *> \author Univ. of California Berkeley 
00119 *> \author Univ. of Colorado Denver 
00120 *> \author NAG Ltd. 
00121 *
00122 *> \date November 2011
00123 *
00124 *> \ingroup complexHEcomputational
00125 *
00126 *  =====================================================================
00127       SUBROUTINE CHETRS2( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, 
00128      $                    WORK, INFO )
00129 *
00130 *  -- LAPACK computational routine (version 3.4.0) --
00131 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00132 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00133 *     November 2011
00134 *
00135 *     .. Scalar Arguments ..
00136       CHARACTER          UPLO
00137       INTEGER            INFO, LDA, LDB, N, NRHS
00138 *     ..
00139 *     .. Array Arguments ..
00140       INTEGER            IPIV( * )
00141       COMPLEX            A( LDA, * ), B( LDB, * ), WORK( * )
00142 *     ..
00143 *
00144 *  =====================================================================
00145 *
00146 *     .. Parameters ..
00147       COMPLEX            ONE
00148       PARAMETER          ( ONE = (1.0E+0,0.0E+0) )
00149 *     ..
00150 *     .. Local Scalars ..
00151       LOGICAL            UPPER
00152       INTEGER            I, IINFO, J, K, KP
00153       REAL               S
00154       COMPLEX            AK, AKM1, AKM1K, BK, BKM1, DENOM
00155 *     ..
00156 *     .. External Functions ..
00157       LOGICAL            LSAME
00158       EXTERNAL           LSAME
00159 *     ..
00160 *     .. External Subroutines ..
00161       EXTERNAL           CSCAL, CSYCONV, CSWAP, CTRSM, XERBLA
00162 *     ..
00163 *     .. Intrinsic Functions ..
00164       INTRINSIC          CONJG, MAX, REAL
00165 *     ..
00166 *     .. Executable Statements ..
00167 *
00168       INFO = 0
00169       UPPER = LSAME( UPLO, 'U' )
00170       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00171          INFO = -1
00172       ELSE IF( N.LT.0 ) THEN
00173          INFO = -2
00174       ELSE IF( NRHS.LT.0 ) THEN
00175          INFO = -3
00176       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00177          INFO = -5
00178       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00179          INFO = -8
00180       END IF
00181       IF( INFO.NE.0 ) THEN
00182          CALL XERBLA( 'CHETRS2', -INFO )
00183          RETURN
00184       END IF
00185 *
00186 *     Quick return if possible
00187 *
00188       IF( N.EQ.0 .OR. NRHS.EQ.0 )
00189      $   RETURN
00190 *
00191 *     Convert A
00192 *
00193       CALL CSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO )
00194 *
00195       IF( UPPER ) THEN
00196 *
00197 *        Solve A*X = B, where A = U*D*U**H.
00198 *
00199 *       P**T * B  
00200         K=N
00201         DO WHILE ( K .GE. 1 )
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-1 and -IPIV(K).
00212             KP = -IPIV( K )
00213             IF( KP.EQ.-IPIV( K-1 ) )
00214      $         CALL CSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB )
00215             K=K-2
00216          END IF
00217         END DO
00218 *
00219 *  Compute (U \P**T * B) -> B    [ (U \P**T * B) ]
00220 *
00221         CALL CTRSM('L','U','N','U',N,NRHS,ONE,A,LDA,B,LDB)
00222 *
00223 *  Compute D \ B -> B   [ D \ (U \P**T * B) ]
00224 *       
00225          I=N
00226          DO WHILE ( I .GE. 1 )
00227             IF( IPIV(I) .GT. 0 ) THEN
00228               S = REAL( ONE ) / REAL( A( I, I ) )
00229               CALL CSSCAL( NRHS, S, B( I, 1 ), LDB )
00230             ELSEIF ( I .GT. 1) THEN
00231                IF ( IPIV(I-1) .EQ. IPIV(I) ) THEN
00232                   AKM1K = WORK(I)
00233                   AKM1 = A( I-1, I-1 ) / AKM1K
00234                   AK = A( I, I ) / CONJG( AKM1K )
00235                   DENOM = AKM1*AK - ONE
00236                   DO 15 J = 1, NRHS
00237                      BKM1 = B( I-1, J ) / AKM1K
00238                      BK = B( I, J ) / CONJG( AKM1K )
00239                      B( I-1, J ) = ( AK*BKM1-BK ) / DENOM
00240                      B( I, J ) = ( AKM1*BK-BKM1 ) / DENOM
00241  15              CONTINUE
00242                I = I - 1
00243                ENDIF
00244             ENDIF
00245             I = I - 1
00246          END DO
00247 *
00248 *      Compute (U**H \ B) -> B   [ U**H \ (D \ (U \P**T * B) ) ]
00249 *
00250          CALL CTRSM('L','U','C','U',N,NRHS,ONE,A,LDA,B,LDB)
00251 *
00252 *       P * B  [ P * (U**H \ (D \ (U \P**T * B) )) ]
00253 *
00254         K=1
00255         DO WHILE ( K .LE. N )
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 .LT. N .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       ELSE
00274 *
00275 *        Solve A*X = B, where A = L*D*L**H.
00276 *
00277 *       P**T * B  
00278         K=1
00279         DO WHILE ( K .LE. N )
00280          IF( IPIV( K ).GT.0 ) THEN
00281 *           1 x 1 diagonal block
00282 *           Interchange rows K and IPIV(K).
00283             KP = IPIV( K )
00284             IF( KP.NE.K )
00285      $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00286             K=K+1
00287          ELSE
00288 *           2 x 2 diagonal block
00289 *           Interchange rows K and -IPIV(K+1).
00290             KP = -IPIV( K+1 )
00291             IF( KP.EQ.-IPIV( K ) )
00292      $         CALL CSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB )
00293             K=K+2
00294          ENDIF
00295         END DO
00296 *
00297 *  Compute (L \P**T * B) -> B    [ (L \P**T * B) ]
00298 *
00299         CALL CTRSM('L','L','N','U',N,NRHS,ONE,A,LDA,B,LDB)
00300 *
00301 *  Compute D \ B -> B   [ D \ (L \P**T * B) ]
00302 *       
00303          I=1
00304          DO WHILE ( I .LE. N )
00305             IF( IPIV(I) .GT. 0 ) THEN
00306               S = REAL( ONE ) / REAL( A( I, I ) )
00307               CALL CSSCAL( NRHS, S, B( I, 1 ), LDB )
00308             ELSE
00309                   AKM1K = WORK(I)
00310                   AKM1 = A( I, I ) / CONJG( AKM1K )
00311                   AK = A( I+1, I+1 ) / AKM1K
00312                   DENOM = AKM1*AK - ONE
00313                   DO 25 J = 1, NRHS
00314                      BKM1 = B( I, J ) / CONJG( AKM1K )
00315                      BK = B( I+1, J ) / AKM1K
00316                      B( I, J ) = ( AK*BKM1-BK ) / DENOM
00317                      B( I+1, J ) = ( AKM1*BK-BKM1 ) / DENOM
00318  25              CONTINUE
00319                   I = I + 1
00320             ENDIF
00321             I = I + 1
00322          END DO
00323 *
00324 *  Compute (L**H \ B) -> B   [ L**H \ (D \ (L \P**T * B) ) ]
00325 * 
00326         CALL CTRSM('L','L','C','U',N,NRHS,ONE,A,LDA,B,LDB)
00327 *
00328 *       P * B  [ P * (L**H \ (D \ (L \P**T * B) )) ]
00329 *
00330         K=N
00331         DO WHILE ( K .GE. 1 )
00332          IF( IPIV( K ).GT.0 ) THEN
00333 *           1 x 1 diagonal block
00334 *           Interchange rows K and IPIV(K).
00335             KP = IPIV( K )
00336             IF( KP.NE.K )
00337      $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00338             K=K-1
00339          ELSE
00340 *           2 x 2 diagonal block
00341 *           Interchange rows K-1 and -IPIV(K).
00342             KP = -IPIV( K )
00343             IF( K.GT.1 .AND. KP.EQ.-IPIV( K-1 ) )
00344      $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00345             K=K-2
00346          ENDIF
00347         END DO
00348 *
00349       END IF
00350 *
00351 *     Revert A
00352 *
00353       CALL CSYCONV( UPLO, 'R', N, A, LDA, IPIV, WORK, IINFO )
00354 *
00355       RETURN
00356 *
00357 *     End of CHETRS2
00358 *
00359       END
 All Files Functions