LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
csytrs2.f
Go to the documentation of this file.
00001 *> \brief \b CSYTRS2
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download CSYTRS2 + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/csytrs2.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/csytrs2.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/csytrs2.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE CSYTRS2( 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 *> CSYTRS2 solves a system of linear equations A*X = B with a COMPLEX
00040 *> symmetric matrix A using the factorization A = U*D*U**T or
00041 *> A = L*D*L**T computed by CSYTRF 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**T;
00053 *>          = 'L':  Lower triangular, form is A = L*D*L**T.
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 CSYTRF.
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 CSYTRF.
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 complexSYcomputational
00125 *
00126 *  =====================================================================
00127       SUBROUTINE CSYTRS2( 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       COMPLEX            AK, AKM1, AKM1K, BK, BKM1, DENOM
00154 *     ..
00155 *     .. External Functions ..
00156       LOGICAL            LSAME
00157       EXTERNAL           LSAME
00158 *     ..
00159 *     .. External Subroutines ..
00160       EXTERNAL           CSCAL, CSYCONV, CSWAP, CTRSM, XERBLA
00161 *     ..
00162 *     .. Intrinsic Functions ..
00163       INTRINSIC          MAX
00164 *     ..
00165 *     .. Executable Statements ..
00166 *
00167       INFO = 0
00168       UPPER = LSAME( UPLO, 'U' )
00169       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00170          INFO = -1
00171       ELSE IF( N.LT.0 ) THEN
00172          INFO = -2
00173       ELSE IF( NRHS.LT.0 ) THEN
00174          INFO = -3
00175       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00176          INFO = -5
00177       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00178          INFO = -8
00179       END IF
00180       IF( INFO.NE.0 ) THEN
00181          CALL XERBLA( 'CSYTRS2', -INFO )
00182          RETURN
00183       END IF
00184 *
00185 *     Quick return if possible
00186 *
00187       IF( N.EQ.0 .OR. NRHS.EQ.0 )
00188      $   RETURN
00189 *
00190 *     Convert A
00191 *
00192       CALL CSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO )
00193 *
00194       IF( UPPER ) THEN
00195 *
00196 *        Solve A*X = B, where A = U*D*U**T.
00197 *
00198 *       P**T * B  
00199         K=N
00200         DO WHILE ( K .GE. 1 )
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 CSWAP( 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-1 and -IPIV(K).
00211             KP = -IPIV( K )
00212             IF( KP.EQ.-IPIV( K-1 ) )
00213      $         CALL CSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB )
00214             K=K-2
00215          END IF
00216         END DO
00217 *
00218 *  Compute (U \P**T * B) -> B    [ (U \P**T * B) ]
00219 *
00220         CALL CTRSM('L','U','N','U',N,NRHS,ONE,A,LDA,B,LDB)
00221 *
00222 *  Compute D \ B -> B   [ D \ (U \P**T * B) ]
00223 *       
00224          I=N
00225          DO WHILE ( I .GE. 1 )
00226             IF( IPIV(I) .GT. 0 ) THEN
00227               CALL CSCAL( NRHS, ONE / A( I, I ), B( I, 1 ), LDB )
00228             ELSEIF ( I .GT. 1) THEN
00229                IF ( IPIV(I-1) .EQ. IPIV(I) ) THEN
00230                   AKM1K = WORK(I)
00231                   AKM1 = A( I-1, I-1 ) / AKM1K
00232                   AK = A( I, I ) / AKM1K
00233                   DENOM = AKM1*AK - ONE
00234                   DO 15 J = 1, NRHS
00235                      BKM1 = B( I-1, J ) / AKM1K
00236                      BK = B( I, J ) / AKM1K
00237                      B( I-1, J ) = ( AK*BKM1-BK ) / DENOM
00238                      B( I, J ) = ( AKM1*BK-BKM1 ) / DENOM
00239  15              CONTINUE
00240                I = I - 1
00241                ENDIF
00242             ENDIF
00243             I = I - 1
00244          END DO
00245 *
00246 *      Compute (U**T \ B) -> B   [ U**T \ (D \ (U \P**T * B) ) ]
00247 *
00248          CALL CTRSM('L','U','T','U',N,NRHS,ONE,A,LDA,B,LDB)
00249 *
00250 *       P * B  [ P * (U**T \ (D \ (U \P**T * B) )) ]
00251 *
00252         K=1
00253         DO WHILE ( K .LE. N )
00254          IF( IPIV( K ).GT.0 ) THEN
00255 *           1 x 1 diagonal block
00256 *           Interchange rows K and IPIV(K).
00257             KP = IPIV( K )
00258             IF( KP.NE.K )
00259      $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00260             K=K+1
00261          ELSE
00262 *           2 x 2 diagonal block
00263 *           Interchange rows K-1 and -IPIV(K).
00264             KP = -IPIV( K )
00265             IF( K .LT. N .AND. KP.EQ.-IPIV( K+1 ) )
00266      $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00267             K=K+2
00268          ENDIF
00269         END DO
00270 *
00271       ELSE
00272 *
00273 *        Solve A*X = B, where A = L*D*L**T.
00274 *
00275 *       P**T * B  
00276         K=1
00277         DO WHILE ( K .LE. N )
00278          IF( IPIV( K ).GT.0 ) THEN
00279 *           1 x 1 diagonal block
00280 *           Interchange rows K and IPIV(K).
00281             KP = IPIV( K )
00282             IF( KP.NE.K )
00283      $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00284             K=K+1
00285          ELSE
00286 *           2 x 2 diagonal block
00287 *           Interchange rows K and -IPIV(K+1).
00288             KP = -IPIV( K+1 )
00289             IF( KP.EQ.-IPIV( K ) )
00290      $         CALL CSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB )
00291             K=K+2
00292          ENDIF
00293         END DO
00294 *
00295 *  Compute (L \P**T * B) -> B    [ (L \P**T * B) ]
00296 *
00297         CALL CTRSM('L','L','N','U',N,NRHS,ONE,A,LDA,B,LDB)
00298 *
00299 *  Compute D \ B -> B   [ D \ (L \P**T * B) ]
00300 *       
00301          I=1
00302          DO WHILE ( I .LE. N )
00303             IF( IPIV(I) .GT. 0 ) THEN
00304               CALL CSCAL( NRHS, ONE / A( I, I ), B( I, 1 ), LDB )
00305             ELSE
00306                   AKM1K = WORK(I)
00307                   AKM1 = A( I, I ) / AKM1K
00308                   AK = A( I+1, I+1 ) / AKM1K
00309                   DENOM = AKM1*AK - ONE
00310                   DO 25 J = 1, NRHS
00311                      BKM1 = B( I, J ) / AKM1K
00312                      BK = B( I+1, J ) / AKM1K
00313                      B( I, J ) = ( AK*BKM1-BK ) / DENOM
00314                      B( I+1, J ) = ( AKM1*BK-BKM1 ) / DENOM
00315  25              CONTINUE
00316                   I = I + 1
00317             ENDIF
00318             I = I + 1
00319          END DO
00320 *
00321 *  Compute (L**T \ B) -> B   [ L**T \ (D \ (L \P**T * B) ) ]
00322 * 
00323         CALL CTRSM('L','L','T','U',N,NRHS,ONE,A,LDA,B,LDB)
00324 *
00325 *       P * B  [ P * (L**T \ (D \ (L \P**T * B) )) ]
00326 *
00327         K=N
00328         DO WHILE ( K .GE. 1 )
00329          IF( IPIV( K ).GT.0 ) THEN
00330 *           1 x 1 diagonal block
00331 *           Interchange rows K and IPIV(K).
00332             KP = IPIV( K )
00333             IF( KP.NE.K )
00334      $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00335             K=K-1
00336          ELSE
00337 *           2 x 2 diagonal block
00338 *           Interchange rows K-1 and -IPIV(K).
00339             KP = -IPIV( K )
00340             IF( K.GT.1 .AND. KP.EQ.-IPIV( K-1 ) )
00341      $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00342             K=K-2
00343          ENDIF
00344         END DO
00345 *
00346       END IF
00347 *
00348 *     Revert A
00349 *
00350       CALL CSYCONV( UPLO, 'R', N, A, LDA, IPIV, WORK, IINFO )
00351 *
00352       RETURN
00353 *
00354 *     End of CSYTRS2
00355 *
00356       END
 All Files Functions