LAPACK: SRC/chptrs.f Source File

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