LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
csysvx.f
Go to the documentation of this file.
00001 *> \brief <b> CSYSVX computes the solution to system of linear equations A * X = B for SY matrices</b>
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download CSYSVX + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/csysvx.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/csysvx.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/csysvx.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE CSYSVX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B,
00022 *                          LDB, X, LDX, RCOND, FERR, BERR, WORK, LWORK,
00023 *                          RWORK, INFO )
00024 * 
00025 *       .. Scalar Arguments ..
00026 *       CHARACTER          FACT, UPLO
00027 *       INTEGER            INFO, LDA, LDAF, LDB, LDX, LWORK, N, NRHS
00028 *       REAL               RCOND
00029 *       ..
00030 *       .. Array Arguments ..
00031 *       INTEGER            IPIV( * )
00032 *       REAL               BERR( * ), FERR( * ), RWORK( * )
00033 *       COMPLEX            A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
00034 *      $                   WORK( * ), X( LDX, * )
00035 *       ..
00036 *  
00037 *
00038 *> \par Purpose:
00039 *  =============
00040 *>
00041 *> \verbatim
00042 *>
00043 *> CSYSVX uses the diagonal pivoting factorization to compute the
00044 *> solution to a complex system of linear equations A * X = B,
00045 *> where A is an N-by-N symmetric matrix and X and B are N-by-NRHS
00046 *> matrices.
00047 *>
00048 *> Error bounds on the solution and a condition estimate are also
00049 *> provided.
00050 *> \endverbatim
00051 *
00052 *> \par Description:
00053 *  =================
00054 *>
00055 *> \verbatim
00056 *>
00057 *> The following steps are performed:
00058 *>
00059 *> 1. If FACT = 'N', the diagonal pivoting method is used to factor A.
00060 *>    The form of the factorization is
00061 *>       A = U * D * U**T,  if UPLO = 'U', or
00062 *>       A = L * D * L**T,  if UPLO = 'L',
00063 *>    where U (or L) is a product of permutation and unit upper (lower)
00064 *>    triangular matrices, and D is symmetric and block diagonal with
00065 *>    1-by-1 and 2-by-2 diagonal blocks.
00066 *>
00067 *> 2. If some D(i,i)=0, so that D is exactly singular, then the routine
00068 *>    returns with INFO = i. Otherwise, the factored form of A is used
00069 *>    to estimate the condition number of the matrix A.  If the
00070 *>    reciprocal of the condition number is less than machine precision,
00071 *>    INFO = N+1 is returned as a warning, but the routine still goes on
00072 *>    to solve for X and compute error bounds as described below.
00073 *>
00074 *> 3. The system of equations is solved for X using the factored form
00075 *>    of A.
00076 *>
00077 *> 4. Iterative refinement is applied to improve the computed solution
00078 *>    matrix and calculate error bounds and backward error estimates
00079 *>    for it.
00080 *> \endverbatim
00081 *
00082 *  Arguments:
00083 *  ==========
00084 *
00085 *> \param[in] FACT
00086 *> \verbatim
00087 *>          FACT is CHARACTER*1
00088 *>          Specifies whether or not the factored form of A has been
00089 *>          supplied on entry.
00090 *>          = 'F':  On entry, AF and IPIV contain the factored form
00091 *>                  of A.  A, AF and IPIV will not be modified.
00092 *>          = 'N':  The matrix A will be copied to AF and factored.
00093 *> \endverbatim
00094 *>
00095 *> \param[in] UPLO
00096 *> \verbatim
00097 *>          UPLO is CHARACTER*1
00098 *>          = 'U':  Upper triangle of A is stored;
00099 *>          = 'L':  Lower triangle of A is stored.
00100 *> \endverbatim
00101 *>
00102 *> \param[in] N
00103 *> \verbatim
00104 *>          N is INTEGER
00105 *>          The number of linear equations, i.e., the order of the
00106 *>          matrix A.  N >= 0.
00107 *> \endverbatim
00108 *>
00109 *> \param[in] NRHS
00110 *> \verbatim
00111 *>          NRHS is INTEGER
00112 *>          The number of right hand sides, i.e., the number of columns
00113 *>          of the matrices B and X.  NRHS >= 0.
00114 *> \endverbatim
00115 *>
00116 *> \param[in] A
00117 *> \verbatim
00118 *>          A is COMPLEX array, dimension (LDA,N)
00119 *>          The symmetric matrix A.  If UPLO = 'U', the leading N-by-N
00120 *>          upper triangular part of A contains the upper triangular part
00121 *>          of the matrix A, and the strictly lower triangular part of A
00122 *>          is not referenced.  If UPLO = 'L', the leading N-by-N lower
00123 *>          triangular part of A contains the lower triangular part of
00124 *>          the matrix A, and the strictly upper triangular part of A is
00125 *>          not referenced.
00126 *> \endverbatim
00127 *>
00128 *> \param[in] LDA
00129 *> \verbatim
00130 *>          LDA is INTEGER
00131 *>          The leading dimension of the array A.  LDA >= max(1,N).
00132 *> \endverbatim
00133 *>
00134 *> \param[in,out] AF
00135 *> \verbatim
00136 *>          AF is COMPLEX array, dimension (LDAF,N)
00137 *>          If FACT = 'F', then AF is an input argument and on entry
00138 *>          contains the block diagonal matrix D and the multipliers used
00139 *>          to obtain the factor U or L from the factorization
00140 *>          A = U*D*U**T or A = L*D*L**T as computed by CSYTRF.
00141 *>
00142 *>          If FACT = 'N', then AF is an output argument and on exit
00143 *>          returns the block diagonal matrix D and the multipliers used
00144 *>          to obtain the factor U or L from the factorization
00145 *>          A = U*D*U**T or A = L*D*L**T.
00146 *> \endverbatim
00147 *>
00148 *> \param[in] LDAF
00149 *> \verbatim
00150 *>          LDAF is INTEGER
00151 *>          The leading dimension of the array AF.  LDAF >= max(1,N).
00152 *> \endverbatim
00153 *>
00154 *> \param[in,out] IPIV
00155 *> \verbatim
00156 *>          IPIV is INTEGER array, dimension (N)
00157 *>          If FACT = 'F', then IPIV is an input argument and on entry
00158 *>          contains details of the interchanges and the block structure
00159 *>          of D, as determined by CSYTRF.
00160 *>          If IPIV(k) > 0, then rows and columns k and IPIV(k) were
00161 *>          interchanged and D(k,k) is a 1-by-1 diagonal block.
00162 *>          If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and
00163 *>          columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
00164 *>          is a 2-by-2 diagonal block.  If UPLO = 'L' and IPIV(k) =
00165 *>          IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were
00166 *>          interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
00167 *>
00168 *>          If FACT = 'N', then IPIV is an output argument and on exit
00169 *>          contains details of the interchanges and the block structure
00170 *>          of D, as determined by CSYTRF.
00171 *> \endverbatim
00172 *>
00173 *> \param[in] B
00174 *> \verbatim
00175 *>          B is COMPLEX array, dimension (LDB,NRHS)
00176 *>          The N-by-NRHS right hand side matrix B.
00177 *> \endverbatim
00178 *>
00179 *> \param[in] LDB
00180 *> \verbatim
00181 *>          LDB is INTEGER
00182 *>          The leading dimension of the array B.  LDB >= max(1,N).
00183 *> \endverbatim
00184 *>
00185 *> \param[out] X
00186 *> \verbatim
00187 *>          X is COMPLEX array, dimension (LDX,NRHS)
00188 *>          If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X.
00189 *> \endverbatim
00190 *>
00191 *> \param[in] LDX
00192 *> \verbatim
00193 *>          LDX is INTEGER
00194 *>          The leading dimension of the array X.  LDX >= max(1,N).
00195 *> \endverbatim
00196 *>
00197 *> \param[out] RCOND
00198 *> \verbatim
00199 *>          RCOND is REAL
00200 *>          The estimate of the reciprocal condition number of the matrix
00201 *>          A.  If RCOND is less than the machine precision (in
00202 *>          particular, if RCOND = 0), the matrix is singular to working
00203 *>          precision.  This condition is indicated by a return code of
00204 *>          INFO > 0.
00205 *> \endverbatim
00206 *>
00207 *> \param[out] FERR
00208 *> \verbatim
00209 *>          FERR is REAL array, dimension (NRHS)
00210 *>          The estimated forward error bound for each solution vector
00211 *>          X(j) (the j-th column of the solution matrix X).
00212 *>          If XTRUE is the true solution corresponding to X(j), FERR(j)
00213 *>          is an estimated upper bound for the magnitude of the largest
00214 *>          element in (X(j) - XTRUE) divided by the magnitude of the
00215 *>          largest element in X(j).  The estimate is as reliable as
00216 *>          the estimate for RCOND, and is almost always a slight
00217 *>          overestimate of the true error.
00218 *> \endverbatim
00219 *>
00220 *> \param[out] BERR
00221 *> \verbatim
00222 *>          BERR is REAL array, dimension (NRHS)
00223 *>          The componentwise relative backward error of each solution
00224 *>          vector X(j) (i.e., the smallest relative change in
00225 *>          any element of A or B that makes X(j) an exact solution).
00226 *> \endverbatim
00227 *>
00228 *> \param[out] WORK
00229 *> \verbatim
00230 *>          WORK is COMPLEX array, dimension (MAX(1,LWORK))
00231 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00232 *> \endverbatim
00233 *>
00234 *> \param[in] LWORK
00235 *> \verbatim
00236 *>          LWORK is INTEGER
00237 *>          The length of WORK.  LWORK >= max(1,2*N), and for best
00238 *>          performance, when FACT = 'N', LWORK >= max(1,2*N,N*NB), where
00239 *>          NB is the optimal blocksize for CSYTRF.
00240 *>
00241 *>          If LWORK = -1, then a workspace query is assumed; the routine
00242 *>          only calculates the optimal size of the WORK array, returns
00243 *>          this value as the first entry of the WORK array, and no error
00244 *>          message related to LWORK is issued by XERBLA.
00245 *> \endverbatim
00246 *>
00247 *> \param[out] RWORK
00248 *> \verbatim
00249 *>          RWORK is REAL array, dimension (N)
00250 *> \endverbatim
00251 *>
00252 *> \param[out] INFO
00253 *> \verbatim
00254 *>          INFO is INTEGER
00255 *>          = 0: successful exit
00256 *>          < 0: if INFO = -i, the i-th argument had an illegal value
00257 *>          > 0: if INFO = i, and i is
00258 *>                <= N:  D(i,i) is exactly zero.  The factorization
00259 *>                       has been completed but the factor D is exactly
00260 *>                       singular, so the solution and error bounds could
00261 *>                       not be computed. RCOND = 0 is returned.
00262 *>                = N+1: D is nonsingular, but RCOND is less than machine
00263 *>                       precision, meaning that the matrix is singular
00264 *>                       to working precision.  Nevertheless, the
00265 *>                       solution and error bounds are computed because
00266 *>                       there are a number of situations where the
00267 *>                       computed solution can be more accurate than the
00268 *>                       value of RCOND would suggest.
00269 *> \endverbatim
00270 *
00271 *  Authors:
00272 *  ========
00273 *
00274 *> \author Univ. of Tennessee 
00275 *> \author Univ. of California Berkeley 
00276 *> \author Univ. of Colorado Denver 
00277 *> \author NAG Ltd. 
00278 *
00279 *> \date April 2012
00280 *
00281 *> \ingroup complexSYsolve
00282 *
00283 *  =====================================================================
00284       SUBROUTINE CSYSVX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B,
00285      $                   LDB, X, LDX, RCOND, FERR, BERR, WORK, LWORK,
00286      $                   RWORK, INFO )
00287 *
00288 *  -- LAPACK driver routine (version 3.4.1) --
00289 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00290 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00291 *     April 2012
00292 *
00293 *     .. Scalar Arguments ..
00294       CHARACTER          FACT, UPLO
00295       INTEGER            INFO, LDA, LDAF, LDB, LDX, LWORK, N, NRHS
00296       REAL               RCOND
00297 *     ..
00298 *     .. Array Arguments ..
00299       INTEGER            IPIV( * )
00300       REAL               BERR( * ), FERR( * ), RWORK( * )
00301       COMPLEX            A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
00302      $                   WORK( * ), X( LDX, * )
00303 *     ..
00304 *
00305 *  =====================================================================
00306 *
00307 *     .. Parameters ..
00308       REAL               ZERO
00309       PARAMETER          ( ZERO = 0.0E+0 )
00310 *     ..
00311 *     .. Local Scalars ..
00312       LOGICAL            LQUERY, NOFACT
00313       INTEGER            LWKOPT, NB
00314       REAL               ANORM
00315 *     ..
00316 *     .. External Functions ..
00317       LOGICAL            LSAME
00318       INTEGER            ILAENV
00319       REAL               CLANSY, SLAMCH
00320       EXTERNAL           ILAENV, LSAME, CLANSY, SLAMCH
00321 *     ..
00322 *     .. External Subroutines ..
00323       EXTERNAL           CLACPY, CSYCON, CSYRFS, CSYTRF, CSYTRS, XERBLA
00324 *     ..
00325 *     .. Intrinsic Functions ..
00326       INTRINSIC          MAX
00327 *     ..
00328 *     .. Executable Statements ..
00329 *
00330 *     Test the input parameters.
00331 *
00332       INFO = 0
00333       NOFACT = LSAME( FACT, 'N' )
00334       LQUERY = ( LWORK.EQ.-1 )
00335       IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) ) THEN
00336          INFO = -1
00337       ELSE IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) )
00338      $          THEN
00339          INFO = -2
00340       ELSE IF( N.LT.0 ) THEN
00341          INFO = -3
00342       ELSE IF( NRHS.LT.0 ) THEN
00343          INFO = -4
00344       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00345          INFO = -6
00346       ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN
00347          INFO = -8
00348       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00349          INFO = -11
00350       ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
00351          INFO = -13
00352       ELSE IF( LWORK.LT.MAX( 1, 2*N ) .AND. .NOT.LQUERY ) THEN
00353          INFO = -18
00354       END IF
00355 *
00356       IF( INFO.EQ.0 ) THEN
00357          LWKOPT = MAX( 1, 2*N )
00358          IF( NOFACT ) THEN
00359             NB = ILAENV( 1, 'CSYTRF', UPLO, N, -1, -1, -1 )
00360             LWKOPT = MAX( LWKOPT, N*NB )
00361          END IF
00362          WORK( 1 ) = LWKOPT
00363       END IF
00364 *
00365       IF( INFO.NE.0 ) THEN
00366          CALL XERBLA( 'CSYSVX', -INFO )
00367          RETURN
00368       ELSE IF( LQUERY ) THEN
00369          RETURN
00370       END IF
00371 *
00372       IF( NOFACT ) THEN
00373 *
00374 *        Compute the factorization A = U*D*U**T or A = L*D*L**T.
00375 *
00376          CALL CLACPY( UPLO, N, N, A, LDA, AF, LDAF )
00377          CALL CSYTRF( UPLO, N, AF, LDAF, IPIV, WORK, LWORK, INFO )
00378 *
00379 *        Return if INFO is non-zero.
00380 *
00381          IF( INFO.GT.0 )THEN
00382             RCOND = ZERO
00383             RETURN
00384          END IF
00385       END IF
00386 *
00387 *     Compute the norm of the matrix A.
00388 *
00389       ANORM = CLANSY( 'I', UPLO, N, A, LDA, RWORK )
00390 *
00391 *     Compute the reciprocal of the condition number of A.
00392 *
00393       CALL CSYCON( UPLO, N, AF, LDAF, IPIV, ANORM, RCOND, WORK, INFO )
00394 *
00395 *     Compute the solution vectors X.
00396 *
00397       CALL CLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
00398       CALL CSYTRS( UPLO, N, NRHS, AF, LDAF, IPIV, X, LDX, INFO )
00399 *
00400 *     Use iterative refinement to improve the computed solutions and
00401 *     compute error bounds and backward error estimates for them.
00402 *
00403       CALL CSYRFS( UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B, LDB, X,
00404      $             LDX, FERR, BERR, WORK, RWORK, INFO )
00405 *
00406 *     Set INFO = N+1 if the matrix is singular to working precision.
00407 *
00408       IF( RCOND.LT.SLAMCH( 'Epsilon' ) )
00409      $   INFO = N + 1
00410 *
00411       WORK( 1 ) = LWKOPT
00412 *
00413       RETURN
00414 *
00415 *     End of CSYSVX
00416 *
00417       END
 All Files Functions