LAPACK: SRC/chetrs.f Source File

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