LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
chetrs.f
Go to the documentation of this file.
00001 *> \brief \b CHETRS
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download CHETRS + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/chetrs.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/chetrs.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/chetrs.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE CHETRS( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO )
00022 * 
00023 *       .. Scalar Arguments ..
00024 *       CHARACTER          UPLO
00025 *       INTEGER            INFO, LDA, LDB, N, NRHS
00026 *       ..
00027 *       .. Array Arguments ..
00028 *       INTEGER            IPIV( * )
00029 *       COMPLEX            A( LDA, * ), B( LDB, * )
00030 *       ..
00031 *  
00032 *
00033 *> \par Purpose:
00034 *  =============
00035 *>
00036 *> \verbatim
00037 *>
00038 *> CHETRS solves a system of linear equations A*X = B with a complex
00039 *> Hermitian matrix A using the factorization A = U*D*U**H or
00040 *> A = L*D*L**H computed by CHETRF.
00041 *> \endverbatim
00042 *
00043 *  Arguments:
00044 *  ==========
00045 *
00046 *> \param[in] UPLO
00047 *> \verbatim
00048 *>          UPLO is CHARACTER*1
00049 *>          Specifies whether the details of the factorization are stored
00050 *>          as an upper or lower triangular matrix.
00051 *>          = 'U':  Upper triangular, form is A = U*D*U**H;
00052 *>          = 'L':  Lower triangular, form is A = L*D*L**H.
00053 *> \endverbatim
00054 *>
00055 *> \param[in] N
00056 *> \verbatim
00057 *>          N is INTEGER
00058 *>          The order of the matrix A.  N >= 0.
00059 *> \endverbatim
00060 *>
00061 *> \param[in] NRHS
00062 *> \verbatim
00063 *>          NRHS is INTEGER
00064 *>          The number of right hand sides, i.e., the number of columns
00065 *>          of the matrix B.  NRHS >= 0.
00066 *> \endverbatim
00067 *>
00068 *> \param[in] A
00069 *> \verbatim
00070 *>          A is COMPLEX array, dimension (LDA,N)
00071 *>          The block diagonal matrix D and the multipliers used to
00072 *>          obtain the factor U or L as computed by CHETRF.
00073 *> \endverbatim
00074 *>
00075 *> \param[in] LDA
00076 *> \verbatim
00077 *>          LDA is INTEGER
00078 *>          The leading dimension of the array A.  LDA >= max(1,N).
00079 *> \endverbatim
00080 *>
00081 *> \param[in] IPIV
00082 *> \verbatim
00083 *>          IPIV is INTEGER array, dimension (N)
00084 *>          Details of the interchanges and the block structure of D
00085 *>          as determined by CHETRF.
00086 *> \endverbatim
00087 *>
00088 *> \param[in,out] B
00089 *> \verbatim
00090 *>          B is COMPLEX array, dimension (LDB,NRHS)
00091 *>          On entry, the right hand side matrix B.
00092 *>          On exit, the solution matrix X.
00093 *> \endverbatim
00094 *>
00095 *> \param[in] LDB
00096 *> \verbatim
00097 *>          LDB is INTEGER
00098 *>          The leading dimension of the array B.  LDB >= max(1,N).
00099 *> \endverbatim
00100 *>
00101 *> \param[out] INFO
00102 *> \verbatim
00103 *>          INFO is INTEGER
00104 *>          = 0:  successful exit
00105 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
00106 *> \endverbatim
00107 *
00108 *  Authors:
00109 *  ========
00110 *
00111 *> \author Univ. of Tennessee 
00112 *> \author Univ. of California Berkeley 
00113 *> \author Univ. of Colorado Denver 
00114 *> \author NAG Ltd. 
00115 *
00116 *> \date November 2011
00117 *
00118 *> \ingroup complexHEcomputational
00119 *
00120 *  =====================================================================
00121       SUBROUTINE CHETRS( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO )
00122 *
00123 *  -- LAPACK computational routine (version 3.4.0) --
00124 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00125 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00126 *     November 2011
00127 *
00128 *     .. Scalar Arguments ..
00129       CHARACTER          UPLO
00130       INTEGER            INFO, LDA, LDB, N, NRHS
00131 *     ..
00132 *     .. Array Arguments ..
00133       INTEGER            IPIV( * )
00134       COMPLEX            A( LDA, * ), B( LDB, * )
00135 *     ..
00136 *
00137 *  =====================================================================
00138 *
00139 *     .. Parameters ..
00140       COMPLEX            ONE
00141       PARAMETER          ( ONE = ( 1.0E+0, 0.0E+0 ) )
00142 *     ..
00143 *     .. Local Scalars ..
00144       LOGICAL            UPPER
00145       INTEGER            J, K, KP
00146       REAL               S
00147       COMPLEX            AK, AKM1, AKM1K, BK, BKM1, DENOM
00148 *     ..
00149 *     .. External Functions ..
00150       LOGICAL            LSAME
00151       EXTERNAL           LSAME
00152 *     ..
00153 *     .. External Subroutines ..
00154       EXTERNAL           CGEMV, CGERU, CLACGV, CSSCAL, CSWAP, XERBLA
00155 *     ..
00156 *     .. Intrinsic Functions ..
00157       INTRINSIC          CONJG, MAX, REAL
00158 *     ..
00159 *     .. Executable Statements ..
00160 *
00161       INFO = 0
00162       UPPER = LSAME( UPLO, 'U' )
00163       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00164          INFO = -1
00165       ELSE IF( N.LT.0 ) THEN
00166          INFO = -2
00167       ELSE IF( NRHS.LT.0 ) THEN
00168          INFO = -3
00169       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00170          INFO = -5
00171       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00172          INFO = -8
00173       END IF
00174       IF( INFO.NE.0 ) THEN
00175          CALL XERBLA( 'CHETRS', -INFO )
00176          RETURN
00177       END IF
00178 *
00179 *     Quick return if possible
00180 *
00181       IF( N.EQ.0 .OR. NRHS.EQ.0 )
00182      $   RETURN
00183 *
00184       IF( UPPER ) THEN
00185 *
00186 *        Solve A*X = B, where A = U*D*U**H.
00187 *
00188 *        First solve U*D*X = B, overwriting B with X.
00189 *
00190 *        K is the main loop index, decreasing from N to 1 in steps of
00191 *        1 or 2, depending on the size of the diagonal blocks.
00192 *
00193          K = N
00194    10    CONTINUE
00195 *
00196 *        If K < 1, exit from loop.
00197 *
00198          IF( K.LT.1 )
00199      $      GO TO 30
00200 *
00201          IF( IPIV( K ).GT.0 ) THEN
00202 *
00203 *           1 x 1 diagonal block
00204 *
00205 *           Interchange rows K and IPIV(K).
00206 *
00207             KP = IPIV( K )
00208             IF( KP.NE.K )
00209      $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00210 *
00211 *           Multiply by inv(U(K)), where U(K) is the transformation
00212 *           stored in column K of A.
00213 *
00214             CALL CGERU( K-1, NRHS, -ONE, A( 1, K ), 1, B( K, 1 ), LDB,
00215      $                  B( 1, 1 ), LDB )
00216 *
00217 *           Multiply by the inverse of the diagonal block.
00218 *
00219             S = REAL( ONE ) / REAL( A( K, K ) )
00220             CALL CSSCAL( NRHS, S, B( K, 1 ), LDB )
00221             K = K - 1
00222          ELSE
00223 *
00224 *           2 x 2 diagonal block
00225 *
00226 *           Interchange rows K-1 and -IPIV(K).
00227 *
00228             KP = -IPIV( K )
00229             IF( KP.NE.K-1 )
00230      $         CALL CSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB )
00231 *
00232 *           Multiply by inv(U(K)), where U(K) is the transformation
00233 *           stored in columns K-1 and K of A.
00234 *
00235             CALL CGERU( K-2, NRHS, -ONE, A( 1, K ), 1, B( K, 1 ), LDB,
00236      $                  B( 1, 1 ), LDB )
00237             CALL CGERU( K-2, NRHS, -ONE, A( 1, K-1 ), 1, B( K-1, 1 ),
00238      $                  LDB, B( 1, 1 ), LDB )
00239 *
00240 *           Multiply by the inverse of the diagonal block.
00241 *
00242             AKM1K = A( K-1, K )
00243             AKM1 = A( K-1, K-1 ) / AKM1K
00244             AK = A( K, K ) / CONJG( AKM1K )
00245             DENOM = AKM1*AK - ONE
00246             DO 20 J = 1, NRHS
00247                BKM1 = B( K-1, J ) / AKM1K
00248                BK = B( K, J ) / CONJG( AKM1K )
00249                B( K-1, J ) = ( AK*BKM1-BK ) / DENOM
00250                B( K, J ) = ( AKM1*BK-BKM1 ) / DENOM
00251    20       CONTINUE
00252             K = K - 2
00253          END IF
00254 *
00255          GO TO 10
00256    30    CONTINUE
00257 *
00258 *        Next solve U**H *X = B, overwriting B with X.
00259 *
00260 *        K is the main loop index, increasing from 1 to N in steps of
00261 *        1 or 2, depending on the size of the diagonal blocks.
00262 *
00263          K = 1
00264    40    CONTINUE
00265 *
00266 *        If K > N, exit from loop.
00267 *
00268          IF( K.GT.N )
00269      $      GO TO 50
00270 *
00271          IF( IPIV( K ).GT.0 ) THEN
00272 *
00273 *           1 x 1 diagonal block
00274 *
00275 *           Multiply by inv(U**H(K)), where U(K) is the transformation
00276 *           stored in column K of A.
00277 *
00278             IF( K.GT.1 ) THEN
00279                CALL CLACGV( NRHS, B( K, 1 ), LDB )
00280                CALL CGEMV( 'Conjugate transpose', K-1, NRHS, -ONE, B,
00281      $                     LDB, A( 1, K ), 1, ONE, B( K, 1 ), LDB )
00282                CALL CLACGV( NRHS, B( K, 1 ), LDB )
00283             END IF
00284 *
00285 *           Interchange rows K and IPIV(K).
00286 *
00287             KP = IPIV( K )
00288             IF( KP.NE.K )
00289      $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00290             K = K + 1
00291          ELSE
00292 *
00293 *           2 x 2 diagonal block
00294 *
00295 *           Multiply by inv(U**H(K+1)), where U(K+1) is the transformation
00296 *           stored in columns K and K+1 of A.
00297 *
00298             IF( K.GT.1 ) THEN
00299                CALL CLACGV( NRHS, B( K, 1 ), LDB )
00300                CALL CGEMV( 'Conjugate transpose', K-1, NRHS, -ONE, B,
00301      $                     LDB, A( 1, K ), 1, ONE, B( K, 1 ), LDB )
00302                CALL CLACGV( NRHS, B( K, 1 ), LDB )
00303 *
00304                CALL CLACGV( NRHS, B( K+1, 1 ), LDB )
00305                CALL CGEMV( 'Conjugate transpose', K-1, NRHS, -ONE, B,
00306      $                     LDB, A( 1, K+1 ), 1, ONE, B( K+1, 1 ), LDB )
00307                CALL CLACGV( NRHS, B( K+1, 1 ), LDB )
00308             END IF
00309 *
00310 *           Interchange rows K and -IPIV(K).
00311 *
00312             KP = -IPIV( K )
00313             IF( KP.NE.K )
00314      $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00315             K = K + 2
00316          END IF
00317 *
00318          GO TO 40
00319    50    CONTINUE
00320 *
00321       ELSE
00322 *
00323 *        Solve A*X = B, where A = L*D*L**H.
00324 *
00325 *        First solve L*D*X = B, overwriting B with X.
00326 *
00327 *        K is the main loop index, increasing from 1 to N in steps of
00328 *        1 or 2, depending on the size of the diagonal blocks.
00329 *
00330          K = 1
00331    60    CONTINUE
00332 *
00333 *        If K > N, exit from loop.
00334 *
00335          IF( K.GT.N )
00336      $      GO TO 80
00337 *
00338          IF( IPIV( K ).GT.0 ) THEN
00339 *
00340 *           1 x 1 diagonal block
00341 *
00342 *           Interchange rows K and IPIV(K).
00343 *
00344             KP = IPIV( K )
00345             IF( KP.NE.K )
00346      $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00347 *
00348 *           Multiply by inv(L(K)), where L(K) is the transformation
00349 *           stored in column K of A.
00350 *
00351             IF( K.LT.N )
00352      $         CALL CGERU( N-K, NRHS, -ONE, A( K+1, K ), 1, B( K, 1 ),
00353      $                     LDB, B( K+1, 1 ), LDB )
00354 *
00355 *           Multiply by the inverse of the diagonal block.
00356 *
00357             S = REAL( ONE ) / REAL( A( K, K ) )
00358             CALL CSSCAL( NRHS, S, B( K, 1 ), LDB )
00359             K = K + 1
00360          ELSE
00361 *
00362 *           2 x 2 diagonal block
00363 *
00364 *           Interchange rows K+1 and -IPIV(K).
00365 *
00366             KP = -IPIV( K )
00367             IF( KP.NE.K+1 )
00368      $         CALL CSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB )
00369 *
00370 *           Multiply by inv(L(K)), where L(K) is the transformation
00371 *           stored in columns K and K+1 of A.
00372 *
00373             IF( K.LT.N-1 ) THEN
00374                CALL CGERU( N-K-1, NRHS, -ONE, A( K+2, K ), 1, B( K, 1 ),
00375      $                     LDB, B( K+2, 1 ), LDB )
00376                CALL CGERU( N-K-1, NRHS, -ONE, A( K+2, K+1 ), 1,
00377      $                     B( K+1, 1 ), LDB, B( K+2, 1 ), LDB )
00378             END IF
00379 *
00380 *           Multiply by the inverse of the diagonal block.
00381 *
00382             AKM1K = A( K+1, K )
00383             AKM1 = A( K, K ) / CONJG( AKM1K )
00384             AK = A( K+1, K+1 ) / AKM1K
00385             DENOM = AKM1*AK - ONE
00386             DO 70 J = 1, NRHS
00387                BKM1 = B( K, J ) / CONJG( AKM1K )
00388                BK = B( K+1, J ) / AKM1K
00389                B( K, J ) = ( AK*BKM1-BK ) / DENOM
00390                B( K+1, J ) = ( AKM1*BK-BKM1 ) / DENOM
00391    70       CONTINUE
00392             K = K + 2
00393          END IF
00394 *
00395          GO TO 60
00396    80    CONTINUE
00397 *
00398 *        Next solve L**H *X = B, overwriting B with X.
00399 *
00400 *        K is the main loop index, decreasing from N to 1 in steps of
00401 *        1 or 2, depending on the size of the diagonal blocks.
00402 *
00403          K = N
00404    90    CONTINUE
00405 *
00406 *        If K < 1, exit from loop.
00407 *
00408          IF( K.LT.1 )
00409      $      GO TO 100
00410 *
00411          IF( IPIV( K ).GT.0 ) THEN
00412 *
00413 *           1 x 1 diagonal block
00414 *
00415 *           Multiply by inv(L**H(K)), where L(K) is the transformation
00416 *           stored in column K of A.
00417 *
00418             IF( K.LT.N ) THEN
00419                CALL CLACGV( NRHS, B( K, 1 ), LDB )
00420                CALL CGEMV( 'Conjugate transpose', N-K, NRHS, -ONE,
00421      $                     B( K+1, 1 ), LDB, A( K+1, K ), 1, ONE,
00422      $                     B( K, 1 ), LDB )
00423                CALL CLACGV( NRHS, B( K, 1 ), LDB )
00424             END IF
00425 *
00426 *           Interchange rows K and IPIV(K).
00427 *
00428             KP = IPIV( K )
00429             IF( KP.NE.K )
00430      $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00431             K = K - 1
00432          ELSE
00433 *
00434 *           2 x 2 diagonal block
00435 *
00436 *           Multiply by inv(L**H(K-1)), where L(K-1) is the transformation
00437 *           stored in columns K-1 and K of A.
00438 *
00439             IF( K.LT.N ) THEN
00440                CALL CLACGV( NRHS, B( K, 1 ), LDB )
00441                CALL CGEMV( 'Conjugate transpose', N-K, NRHS, -ONE,
00442      $                     B( K+1, 1 ), LDB, A( K+1, K ), 1, ONE,
00443      $                     B( K, 1 ), LDB )
00444                CALL CLACGV( NRHS, B( K, 1 ), LDB )
00445 *
00446                CALL CLACGV( NRHS, B( K-1, 1 ), LDB )
00447                CALL CGEMV( 'Conjugate transpose', N-K, NRHS, -ONE,
00448      $                     B( K+1, 1 ), LDB, A( K+1, K-1 ), 1, ONE,
00449      $                     B( K-1, 1 ), LDB )
00450                CALL CLACGV( NRHS, B( K-1, 1 ), LDB )
00451             END IF
00452 *
00453 *           Interchange rows K and -IPIV(K).
00454 *
00455             KP = -IPIV( K )
00456             IF( KP.NE.K )
00457      $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00458             K = K - 2
00459          END IF
00460 *
00461          GO TO 90
00462   100    CONTINUE
00463       END IF
00464 *
00465       RETURN
00466 *
00467 *     End of CHETRS
00468 *
00469       END
 All Files Functions