LAPACK: SRC/csytrs.f Source File

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