LAPACK: SRC/csptrs.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CSPTRS( 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 *  CSPTRS solves a system of linear equations A*X = B with a complex
00021 *  symmetric matrix A stored in packed format using the factorization
00022 *  A = U*D*U**T or A = L*D*L**T computed by CSPTRF.
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**T;
00031 *          = 'L':  Lower triangular, form is A = L*D*L**T.
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 CSPTRF, 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 CSPTRF.
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       COMPLEX            AK, AKM1, AKM1K, BK, BKM1, DENOM
00070 *     ..
00071 *     .. External Functions ..
00072       LOGICAL            LSAME
00073       EXTERNAL           LSAME
00074 *     ..
00075 *     .. External Subroutines ..
00076       EXTERNAL           CGEMV, CGERU, CSCAL, CSWAP, XERBLA
00077 *     ..
00078 *     .. Intrinsic Functions ..
00079       INTRINSIC          MAX
00080 *     ..
00081 *     .. Executable Statements ..
00082 *
00083       INFO = 0
00084       UPPER = LSAME( UPLO, 'U' )
00085       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00086          INFO = -1
00087       ELSE IF( N.LT.0 ) THEN
00088          INFO = -2
00089       ELSE IF( NRHS.LT.0 ) THEN
00090          INFO = -3
00091       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00092          INFO = -7
00093       END IF
00094       IF( INFO.NE.0 ) THEN
00095          CALL XERBLA( 'CSPTRS', -INFO )
00096          RETURN
00097       END IF
00098 *
00099 *     Quick return if possible
00100 *
00101       IF( N.EQ.0 .OR. NRHS.EQ.0 )
00102      $   RETURN
00103 *
00104       IF( UPPER ) THEN
00105 *
00106 *        Solve A*X = B, where A = U*D*U'.
00107 *
00108 *        First solve U*D*X = B, overwriting B with X.
00109 *
00110 *        K is the main loop index, decreasing from N to 1 in steps of
00111 *        1 or 2, depending on the size of the diagonal blocks.
00112 *
00113          K = N
00114          KC = N*( N+1 ) / 2 + 1
00115    10    CONTINUE
00116 *
00117 *        If K < 1, exit from loop.
00118 *
00119          IF( K.LT.1 )
00120      $      GO TO 30
00121 *
00122          KC = KC - K
00123          IF( IPIV( K ).GT.0 ) THEN
00124 *
00125 *           1 x 1 diagonal block
00126 *
00127 *           Interchange rows K and IPIV(K).
00128 *
00129             KP = IPIV( K )
00130             IF( KP.NE.K )
00131      $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00132 *
00133 *           Multiply by inv(U(K)), where U(K) is the transformation
00134 *           stored in column K of A.
00135 *
00136             CALL CGERU( K-1, NRHS, -ONE, AP( KC ), 1, B( K, 1 ), LDB,
00137      $                  B( 1, 1 ), LDB )
00138 *
00139 *           Multiply by the inverse of the diagonal block.
00140 *
00141             CALL CSCAL( NRHS, ONE / AP( KC+K-1 ), B( K, 1 ), LDB )
00142             K = K - 1
00143          ELSE
00144 *
00145 *           2 x 2 diagonal block
00146 *
00147 *           Interchange rows K-1 and -IPIV(K).
00148 *
00149             KP = -IPIV( K )
00150             IF( KP.NE.K-1 )
00151      $         CALL CSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB )
00152 *
00153 *           Multiply by inv(U(K)), where U(K) is the transformation
00154 *           stored in columns K-1 and K of A.
00155 *
00156             CALL CGERU( K-2, NRHS, -ONE, AP( KC ), 1, B( K, 1 ), LDB,
00157      $                  B( 1, 1 ), LDB )
00158             CALL CGERU( K-2, NRHS, -ONE, AP( KC-( K-1 ) ), 1,
00159      $                  B( K-1, 1 ), LDB, B( 1, 1 ), LDB )
00160 *
00161 *           Multiply by the inverse of the diagonal block.
00162 *
00163             AKM1K = AP( KC+K-2 )
00164             AKM1 = AP( KC-1 ) / AKM1K
00165             AK = AP( KC+K-1 ) / AKM1K
00166             DENOM = AKM1*AK - ONE
00167             DO 20 J = 1, NRHS
00168                BKM1 = B( K-1, J ) / AKM1K
00169                BK = B( K, J ) / AKM1K
00170                B( K-1, J ) = ( AK*BKM1-BK ) / DENOM
00171                B( K, J ) = ( AKM1*BK-BKM1 ) / DENOM
00172    20       CONTINUE
00173             KC = KC - K + 1
00174             K = K - 2
00175          END IF
00176 *
00177          GO TO 10
00178    30    CONTINUE
00179 *
00180 *        Next solve U'*X = B, overwriting B with X.
00181 *
00182 *        K is the main loop index, increasing from 1 to N in steps of
00183 *        1 or 2, depending on the size of the diagonal blocks.
00184 *
00185          K = 1
00186          KC = 1
00187    40    CONTINUE
00188 *
00189 *        If K > N, exit from loop.
00190 *
00191          IF( K.GT.N )
00192      $      GO TO 50
00193 *
00194          IF( IPIV( K ).GT.0 ) THEN
00195 *
00196 *           1 x 1 diagonal block
00197 *
00198 *           Multiply by inv(U'(K)), where U(K) is the transformation
00199 *           stored in column K of A.
00200 *
00201             CALL CGEMV( 'Transpose', K-1, NRHS, -ONE, B, LDB, AP( KC ),
00202      $                  1, ONE, B( K, 1 ), LDB )
00203 *
00204 *           Interchange rows K and IPIV(K).
00205 *
00206             KP = IPIV( K )
00207             IF( KP.NE.K )
00208      $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00209             KC = KC + K
00210             K = K + 1
00211          ELSE
00212 *
00213 *           2 x 2 diagonal block
00214 *
00215 *           Multiply by inv(U'(K+1)), where U(K+1) is the transformation
00216 *           stored in columns K and K+1 of A.
00217 *
00218             CALL CGEMV( 'Transpose', K-1, NRHS, -ONE, B, LDB, AP( KC ),
00219      $                  1, ONE, B( K, 1 ), LDB )
00220             CALL CGEMV( 'Transpose', K-1, NRHS, -ONE, B, LDB,
00221      $                  AP( KC+K ), 1, ONE, B( K+1, 1 ), LDB )
00222 *
00223 *           Interchange rows K and -IPIV(K).
00224 *
00225             KP = -IPIV( K )
00226             IF( KP.NE.K )
00227      $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00228             KC = KC + 2*K + 1
00229             K = K + 2
00230          END IF
00231 *
00232          GO TO 40
00233    50    CONTINUE
00234 *
00235       ELSE
00236 *
00237 *        Solve A*X = B, where A = L*D*L'.
00238 *
00239 *        First solve L*D*X = B, overwriting B with X.
00240 *
00241 *        K is the main loop index, increasing from 1 to N in steps of
00242 *        1 or 2, depending on the size of the diagonal blocks.
00243 *
00244          K = 1
00245          KC = 1
00246    60    CONTINUE
00247 *
00248 *        If K > N, exit from loop.
00249 *
00250          IF( K.GT.N )
00251      $      GO TO 80
00252 *
00253          IF( IPIV( K ).GT.0 ) THEN
00254 *
00255 *           1 x 1 diagonal block
00256 *
00257 *           Interchange rows K and IPIV(K).
00258 *
00259             KP = IPIV( K )
00260             IF( KP.NE.K )
00261      $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00262 *
00263 *           Multiply by inv(L(K)), where L(K) is the transformation
00264 *           stored in column K of A.
00265 *
00266             IF( K.LT.N )
00267      $         CALL CGERU( N-K, NRHS, -ONE, AP( KC+1 ), 1, B( K, 1 ),
00268      $                     LDB, B( K+1, 1 ), LDB )
00269 *
00270 *           Multiply by the inverse of the diagonal block.
00271 *
00272             CALL CSCAL( NRHS, ONE / AP( KC ), B( K, 1 ), LDB )
00273             KC = KC + N - K + 1
00274             K = K + 1
00275          ELSE
00276 *
00277 *           2 x 2 diagonal block
00278 *
00279 *           Interchange rows K+1 and -IPIV(K).
00280 *
00281             KP = -IPIV( K )
00282             IF( KP.NE.K+1 )
00283      $         CALL CSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB )
00284 *
00285 *           Multiply by inv(L(K)), where L(K) is the transformation
00286 *           stored in columns K and K+1 of A.
00287 *
00288             IF( K.LT.N-1 ) THEN
00289                CALL CGERU( N-K-1, NRHS, -ONE, AP( KC+2 ), 1, B( K, 1 ),
00290      $                     LDB, B( K+2, 1 ), LDB )
00291                CALL CGERU( N-K-1, NRHS, -ONE, AP( KC+N-K+2 ), 1,
00292      $                     B( K+1, 1 ), LDB, B( K+2, 1 ), LDB )
00293             END IF
00294 *
00295 *           Multiply by the inverse of the diagonal block.
00296 *
00297             AKM1K = AP( KC+1 )
00298             AKM1 = AP( KC ) / AKM1K
00299             AK = AP( KC+N-K+1 ) / AKM1K
00300             DENOM = AKM1*AK - ONE
00301             DO 70 J = 1, NRHS
00302                BKM1 = B( K, J ) / AKM1K
00303                BK = B( K+1, J ) / AKM1K
00304                B( K, J ) = ( AK*BKM1-BK ) / DENOM
00305                B( K+1, J ) = ( AKM1*BK-BKM1 ) / DENOM
00306    70       CONTINUE
00307             KC = KC + 2*( N-K ) + 1
00308             K = K + 2
00309          END IF
00310 *
00311          GO TO 60
00312    80    CONTINUE
00313 *
00314 *        Next solve L'*X = B, overwriting B with X.
00315 *
00316 *        K is the main loop index, decreasing from N to 1 in steps of
00317 *        1 or 2, depending on the size of the diagonal blocks.
00318 *
00319          K = N
00320          KC = N*( N+1 ) / 2 + 1
00321    90    CONTINUE
00322 *
00323 *        If K < 1, exit from loop.
00324 *
00325          IF( K.LT.1 )
00326      $      GO TO 100
00327 *
00328          KC = KC - ( N-K+1 )
00329          IF( IPIV( K ).GT.0 ) THEN
00330 *
00331 *           1 x 1 diagonal block
00332 *
00333 *           Multiply by inv(L'(K)), where L(K) is the transformation
00334 *           stored in column K of A.
00335 *
00336             IF( K.LT.N )
00337      $         CALL CGEMV( 'Transpose', N-K, NRHS, -ONE, B( K+1, 1 ),
00338      $                     LDB, AP( KC+1 ), 1, ONE, B( K, 1 ), LDB )
00339 *
00340 *           Interchange rows K and IPIV(K).
00341 *
00342             KP = IPIV( K )
00343             IF( KP.NE.K )
00344      $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00345             K = K - 1
00346          ELSE
00347 *
00348 *           2 x 2 diagonal block
00349 *
00350 *           Multiply by inv(L'(K-1)), where L(K-1) is the transformation
00351 *           stored in columns K-1 and K of A.
00352 *
00353             IF( K.LT.N ) THEN
00354                CALL CGEMV( 'Transpose', N-K, NRHS, -ONE, B( K+1, 1 ),
00355      $                     LDB, AP( KC+1 ), 1, ONE, B( K, 1 ), LDB )
00356                CALL CGEMV( 'Transpose', N-K, NRHS, -ONE, B( K+1, 1 ),
00357      $                     LDB, AP( KC-( N-K ) ), 1, ONE, B( K-1, 1 ),
00358      $                     LDB )
00359             END IF
00360 *
00361 *           Interchange rows K and -IPIV(K).
00362 *
00363             KP = -IPIV( K )
00364             IF( KP.NE.K )
00365      $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00366             KC = KC - ( N-K+2 )
00367             K = K - 2
00368          END IF
00369 *
00370          GO TO 90
00371   100    CONTINUE
00372       END IF
00373 *
00374       RETURN
00375 *
00376 *     End of CSPTRS
00377 *
00378       END