LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
sposvx.f
Go to the documentation of this file.
00001 *> \brief <b> SPOSVX computes the solution to system of linear equations A * X = B for PO 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 SPOSVX + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sposvx.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sposvx.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sposvx.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE SPOSVX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, EQUED,
00022 *                          S, B, LDB, X, LDX, RCOND, FERR, BERR, WORK,
00023 *                          IWORK, INFO )
00024 * 
00025 *       .. Scalar Arguments ..
00026 *       CHARACTER          EQUED, FACT, UPLO
00027 *       INTEGER            INFO, LDA, LDAF, LDB, LDX, N, NRHS
00028 *       REAL               RCOND
00029 *       ..
00030 *       .. Array Arguments ..
00031 *       INTEGER            IWORK( * )
00032 *       REAL               A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
00033 *      $                   BERR( * ), FERR( * ), S( * ), WORK( * ),
00034 *      $                   X( LDX, * )
00035 *       ..
00036 *  
00037 *
00038 *> \par Purpose:
00039 *  =============
00040 *>
00041 *> \verbatim
00042 *>
00043 *> SPOSVX uses the Cholesky factorization A = U**T*U or A = L*L**T to
00044 *> compute the solution to a real system of linear equations
00045 *>    A * X = B,
00046 *> where A is an N-by-N symmetric positive definite matrix and X and B
00047 *> are N-by-NRHS matrices.
00048 *>
00049 *> Error bounds on the solution and a condition estimate are also
00050 *> provided.
00051 *> \endverbatim
00052 *
00053 *> \par Description:
00054 *  =================
00055 *>
00056 *> \verbatim
00057 *>
00058 *> The following steps are performed:
00059 *>
00060 *> 1. If FACT = 'E', real scaling factors are computed to equilibrate
00061 *>    the system:
00062 *>       diag(S) * A * diag(S) * inv(diag(S)) * X = diag(S) * B
00063 *>    Whether or not the system will be equilibrated depends on the
00064 *>    scaling of the matrix A, but if equilibration is used, A is
00065 *>    overwritten by diag(S)*A*diag(S) and B by diag(S)*B.
00066 *>
00067 *> 2. If FACT = 'N' or 'E', the Cholesky decomposition is used to
00068 *>    factor the matrix A (after equilibration if FACT = 'E') as
00069 *>       A = U**T* U,  if UPLO = 'U', or
00070 *>       A = L * L**T,  if UPLO = 'L',
00071 *>    where U is an upper triangular matrix and L is a lower triangular
00072 *>    matrix.
00073 *>
00074 *> 3. If the leading i-by-i principal minor is not positive definite,
00075 *>    then the routine returns with INFO = i. Otherwise, the factored
00076 *>    form of A is used to estimate the condition number of the matrix
00077 *>    A.  If the reciprocal of the condition number is less than machine
00078 *>    precision, INFO = N+1 is returned as a warning, but the routine
00079 *>    still goes on to solve for X and compute error bounds as
00080 *>    described below.
00081 *>
00082 *> 4. The system of equations is solved for X using the factored form
00083 *>    of A.
00084 *>
00085 *> 5. Iterative refinement is applied to improve the computed solution
00086 *>    matrix and calculate error bounds and backward error estimates
00087 *>    for it.
00088 *>
00089 *> 6. If equilibration was used, the matrix X is premultiplied by
00090 *>    diag(S) so that it solves the original system before
00091 *>    equilibration.
00092 *> \endverbatim
00093 *
00094 *  Arguments:
00095 *  ==========
00096 *
00097 *> \param[in] FACT
00098 *> \verbatim
00099 *>          FACT is CHARACTER*1
00100 *>          Specifies whether or not the factored form of the matrix A is
00101 *>          supplied on entry, and if not, whether the matrix A should be
00102 *>          equilibrated before it is factored.
00103 *>          = 'F':  On entry, AF contains the factored form of A.
00104 *>                  If EQUED = 'Y', the matrix A has been equilibrated
00105 *>                  with scaling factors given by S.  A and AF will not
00106 *>                  be modified.
00107 *>          = 'N':  The matrix A will be copied to AF and factored.
00108 *>          = 'E':  The matrix A will be equilibrated if necessary, then
00109 *>                  copied to AF and factored.
00110 *> \endverbatim
00111 *>
00112 *> \param[in] UPLO
00113 *> \verbatim
00114 *>          UPLO is CHARACTER*1
00115 *>          = 'U':  Upper triangle of A is stored;
00116 *>          = 'L':  Lower triangle of A is stored.
00117 *> \endverbatim
00118 *>
00119 *> \param[in] N
00120 *> \verbatim
00121 *>          N is INTEGER
00122 *>          The number of linear equations, i.e., the order of the
00123 *>          matrix A.  N >= 0.
00124 *> \endverbatim
00125 *>
00126 *> \param[in] NRHS
00127 *> \verbatim
00128 *>          NRHS is INTEGER
00129 *>          The number of right hand sides, i.e., the number of columns
00130 *>          of the matrices B and X.  NRHS >= 0.
00131 *> \endverbatim
00132 *>
00133 *> \param[in,out] A
00134 *> \verbatim
00135 *>          A is REAL array, dimension (LDA,N)
00136 *>          On entry, the symmetric matrix A, except if FACT = 'F' and
00137 *>          EQUED = 'Y', then A must contain the equilibrated matrix
00138 *>          diag(S)*A*diag(S).  If UPLO = 'U', the leading
00139 *>          N-by-N upper triangular part of A contains the upper
00140 *>          triangular part of the matrix A, and the strictly lower
00141 *>          triangular part of A is not referenced.  If UPLO = 'L', the
00142 *>          leading N-by-N lower triangular part of A contains the lower
00143 *>          triangular part of the matrix A, and the strictly upper
00144 *>          triangular part of A is not referenced.  A is not modified if
00145 *>          FACT = 'F' or 'N', or if FACT = 'E' and EQUED = 'N' on exit.
00146 *>
00147 *>          On exit, if FACT = 'E' and EQUED = 'Y', A is overwritten by
00148 *>          diag(S)*A*diag(S).
00149 *> \endverbatim
00150 *>
00151 *> \param[in] LDA
00152 *> \verbatim
00153 *>          LDA is INTEGER
00154 *>          The leading dimension of the array A.  LDA >= max(1,N).
00155 *> \endverbatim
00156 *>
00157 *> \param[in,out] AF
00158 *> \verbatim
00159 *>          AF is REAL array, dimension (LDAF,N)
00160 *>          If FACT = 'F', then AF is an input argument and on entry
00161 *>          contains the triangular factor U or L from the Cholesky
00162 *>          factorization A = U**T*U or A = L*L**T, in the same storage
00163 *>          format as A.  If EQUED .ne. 'N', then AF is the factored form
00164 *>          of the equilibrated matrix diag(S)*A*diag(S).
00165 *>
00166 *>          If FACT = 'N', then AF is an output argument and on exit
00167 *>          returns the triangular factor U or L from the Cholesky
00168 *>          factorization A = U**T*U or A = L*L**T of the original
00169 *>          matrix A.
00170 *>
00171 *>          If FACT = 'E', then AF is an output argument and on exit
00172 *>          returns the triangular factor U or L from the Cholesky
00173 *>          factorization A = U**T*U or A = L*L**T of the equilibrated
00174 *>          matrix A (see the description of A for the form of the
00175 *>          equilibrated matrix).
00176 *> \endverbatim
00177 *>
00178 *> \param[in] LDAF
00179 *> \verbatim
00180 *>          LDAF is INTEGER
00181 *>          The leading dimension of the array AF.  LDAF >= max(1,N).
00182 *> \endverbatim
00183 *>
00184 *> \param[in,out] EQUED
00185 *> \verbatim
00186 *>          EQUED is CHARACTER*1
00187 *>          Specifies the form of equilibration that was done.
00188 *>          = 'N':  No equilibration (always true if FACT = 'N').
00189 *>          = 'Y':  Equilibration was done, i.e., A has been replaced by
00190 *>                  diag(S) * A * diag(S).
00191 *>          EQUED is an input argument if FACT = 'F'; otherwise, it is an
00192 *>          output argument.
00193 *> \endverbatim
00194 *>
00195 *> \param[in,out] S
00196 *> \verbatim
00197 *>          S is REAL array, dimension (N)
00198 *>          The scale factors for A; not accessed if EQUED = 'N'.  S is
00199 *>          an input argument if FACT = 'F'; otherwise, S is an output
00200 *>          argument.  If FACT = 'F' and EQUED = 'Y', each element of S
00201 *>          must be positive.
00202 *> \endverbatim
00203 *>
00204 *> \param[in,out] B
00205 *> \verbatim
00206 *>          B is REAL array, dimension (LDB,NRHS)
00207 *>          On entry, the N-by-NRHS right hand side matrix B.
00208 *>          On exit, if EQUED = 'N', B is not modified; if EQUED = 'Y',
00209 *>          B is overwritten by diag(S) * B.
00210 *> \endverbatim
00211 *>
00212 *> \param[in] LDB
00213 *> \verbatim
00214 *>          LDB is INTEGER
00215 *>          The leading dimension of the array B.  LDB >= max(1,N).
00216 *> \endverbatim
00217 *>
00218 *> \param[out] X
00219 *> \verbatim
00220 *>          X is REAL array, dimension (LDX,NRHS)
00221 *>          If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X to
00222 *>          the original system of equations.  Note that if EQUED = 'Y',
00223 *>          A and B are modified on exit, and the solution to the
00224 *>          equilibrated system is inv(diag(S))*X.
00225 *> \endverbatim
00226 *>
00227 *> \param[in] LDX
00228 *> \verbatim
00229 *>          LDX is INTEGER
00230 *>          The leading dimension of the array X.  LDX >= max(1,N).
00231 *> \endverbatim
00232 *>
00233 *> \param[out] RCOND
00234 *> \verbatim
00235 *>          RCOND is REAL
00236 *>          The estimate of the reciprocal condition number of the matrix
00237 *>          A after equilibration (if done).  If RCOND is less than the
00238 *>          machine precision (in particular, if RCOND = 0), the matrix
00239 *>          is singular to working precision.  This condition is
00240 *>          indicated by a return code of INFO > 0.
00241 *> \endverbatim
00242 *>
00243 *> \param[out] FERR
00244 *> \verbatim
00245 *>          FERR is REAL array, dimension (NRHS)
00246 *>          The estimated forward error bound for each solution vector
00247 *>          X(j) (the j-th column of the solution matrix X).
00248 *>          If XTRUE is the true solution corresponding to X(j), FERR(j)
00249 *>          is an estimated upper bound for the magnitude of the largest
00250 *>          element in (X(j) - XTRUE) divided by the magnitude of the
00251 *>          largest element in X(j).  The estimate is as reliable as
00252 *>          the estimate for RCOND, and is almost always a slight
00253 *>          overestimate of the true error.
00254 *> \endverbatim
00255 *>
00256 *> \param[out] BERR
00257 *> \verbatim
00258 *>          BERR is REAL array, dimension (NRHS)
00259 *>          The componentwise relative backward error of each solution
00260 *>          vector X(j) (i.e., the smallest relative change in
00261 *>          any element of A or B that makes X(j) an exact solution).
00262 *> \endverbatim
00263 *>
00264 *> \param[out] WORK
00265 *> \verbatim
00266 *>          WORK is REAL array, dimension (3*N)
00267 *> \endverbatim
00268 *>
00269 *> \param[out] IWORK
00270 *> \verbatim
00271 *>          IWORK is INTEGER array, dimension (N)
00272 *> \endverbatim
00273 *>
00274 *> \param[out] INFO
00275 *> \verbatim
00276 *>          INFO is INTEGER
00277 *>          = 0: successful exit
00278 *>          < 0: if INFO = -i, the i-th argument had an illegal value
00279 *>          > 0: if INFO = i, and i is
00280 *>                <= N:  the leading minor of order i of A is
00281 *>                       not positive definite, so the factorization
00282 *>                       could not be completed, and the solution has not
00283 *>                       been computed. RCOND = 0 is returned.
00284 *>                = N+1: U is nonsingular, but RCOND is less than machine
00285 *>                       precision, meaning that the matrix is singular
00286 *>                       to working precision.  Nevertheless, the
00287 *>                       solution and error bounds are computed because
00288 *>                       there are a number of situations where the
00289 *>                       computed solution can be more accurate than the
00290 *>                       value of RCOND would suggest.
00291 *> \endverbatim
00292 *
00293 *  Authors:
00294 *  ========
00295 *
00296 *> \author Univ. of Tennessee 
00297 *> \author Univ. of California Berkeley 
00298 *> \author Univ. of Colorado Denver 
00299 *> \author NAG Ltd. 
00300 *
00301 *> \date April 2012
00302 *
00303 *> \ingroup realPOsolve
00304 *
00305 *  =====================================================================
00306       SUBROUTINE SPOSVX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, EQUED,
00307      $                   S, B, LDB, X, LDX, RCOND, FERR, BERR, WORK,
00308      $                   IWORK, INFO )
00309 *
00310 *  -- LAPACK driver routine (version 3.4.1) --
00311 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00312 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00313 *     April 2012
00314 *
00315 *     .. Scalar Arguments ..
00316       CHARACTER          EQUED, FACT, UPLO
00317       INTEGER            INFO, LDA, LDAF, LDB, LDX, N, NRHS
00318       REAL               RCOND
00319 *     ..
00320 *     .. Array Arguments ..
00321       INTEGER            IWORK( * )
00322       REAL               A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
00323      $                   BERR( * ), FERR( * ), S( * ), WORK( * ),
00324      $                   X( LDX, * )
00325 *     ..
00326 *
00327 *  =====================================================================
00328 *
00329 *     .. Parameters ..
00330       REAL               ZERO, ONE
00331       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00332 *     ..
00333 *     .. Local Scalars ..
00334       LOGICAL            EQUIL, NOFACT, RCEQU
00335       INTEGER            I, INFEQU, J
00336       REAL               AMAX, ANORM, BIGNUM, SCOND, SMAX, SMIN, SMLNUM
00337 *     ..
00338 *     .. External Functions ..
00339       LOGICAL            LSAME
00340       REAL               SLAMCH, SLANSY
00341       EXTERNAL           LSAME, SLAMCH, SLANSY
00342 *     ..
00343 *     .. External Subroutines ..
00344       EXTERNAL           SLACPY, SLAQSY, SPOCON, SPOEQU, SPORFS, SPOTRF,
00345      $                   SPOTRS, XERBLA
00346 *     ..
00347 *     .. Intrinsic Functions ..
00348       INTRINSIC          MAX, MIN
00349 *     ..
00350 *     .. Executable Statements ..
00351 *
00352       INFO = 0
00353       NOFACT = LSAME( FACT, 'N' )
00354       EQUIL = LSAME( FACT, 'E' )
00355       IF( NOFACT .OR. EQUIL ) THEN
00356          EQUED = 'N'
00357          RCEQU = .FALSE.
00358       ELSE
00359          RCEQU = LSAME( EQUED, 'Y' )
00360          SMLNUM = SLAMCH( 'Safe minimum' )
00361          BIGNUM = ONE / SMLNUM
00362       END IF
00363 *
00364 *     Test the input parameters.
00365 *
00366       IF( .NOT.NOFACT .AND. .NOT.EQUIL .AND. .NOT.LSAME( FACT, 'F' ) )
00367      $     THEN
00368          INFO = -1
00369       ELSE IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) )
00370      $          THEN
00371          INFO = -2
00372       ELSE IF( N.LT.0 ) THEN
00373          INFO = -3
00374       ELSE IF( NRHS.LT.0 ) THEN
00375          INFO = -4
00376       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00377          INFO = -6
00378       ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN
00379          INFO = -8
00380       ELSE IF( LSAME( FACT, 'F' ) .AND. .NOT.
00381      $         ( RCEQU .OR. LSAME( EQUED, 'N' ) ) ) THEN
00382          INFO = -9
00383       ELSE
00384          IF( RCEQU ) THEN
00385             SMIN = BIGNUM
00386             SMAX = ZERO
00387             DO 10 J = 1, N
00388                SMIN = MIN( SMIN, S( J ) )
00389                SMAX = MAX( SMAX, S( J ) )
00390    10       CONTINUE
00391             IF( SMIN.LE.ZERO ) THEN
00392                INFO = -10
00393             ELSE IF( N.GT.0 ) THEN
00394                SCOND = MAX( SMIN, SMLNUM ) / MIN( SMAX, BIGNUM )
00395             ELSE
00396                SCOND = ONE
00397             END IF
00398          END IF
00399          IF( INFO.EQ.0 ) THEN
00400             IF( LDB.LT.MAX( 1, N ) ) THEN
00401                INFO = -12
00402             ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
00403                INFO = -14
00404             END IF
00405          END IF
00406       END IF
00407 *
00408       IF( INFO.NE.0 ) THEN
00409          CALL XERBLA( 'SPOSVX', -INFO )
00410          RETURN
00411       END IF
00412 *
00413       IF( EQUIL ) THEN
00414 *
00415 *        Compute row and column scalings to equilibrate the matrix A.
00416 *
00417          CALL SPOEQU( N, A, LDA, S, SCOND, AMAX, INFEQU )
00418          IF( INFEQU.EQ.0 ) THEN
00419 *
00420 *           Equilibrate the matrix.
00421 *
00422             CALL SLAQSY( UPLO, N, A, LDA, S, SCOND, AMAX, EQUED )
00423             RCEQU = LSAME( EQUED, 'Y' )
00424          END IF
00425       END IF
00426 *
00427 *     Scale the right hand side.
00428 *
00429       IF( RCEQU ) THEN
00430          DO 30 J = 1, NRHS
00431             DO 20 I = 1, N
00432                B( I, J ) = S( I )*B( I, J )
00433    20       CONTINUE
00434    30    CONTINUE
00435       END IF
00436 *
00437       IF( NOFACT .OR. EQUIL ) THEN
00438 *
00439 *        Compute the Cholesky factorization A = U**T *U or A = L*L**T.
00440 *
00441          CALL SLACPY( UPLO, N, N, A, LDA, AF, LDAF )
00442          CALL SPOTRF( UPLO, N, AF, LDAF, INFO )
00443 *
00444 *        Return if INFO is non-zero.
00445 *
00446          IF( INFO.GT.0 )THEN
00447             RCOND = ZERO
00448             RETURN
00449          END IF
00450       END IF
00451 *
00452 *     Compute the norm of the matrix A.
00453 *
00454       ANORM = SLANSY( '1', UPLO, N, A, LDA, WORK )
00455 *
00456 *     Compute the reciprocal of the condition number of A.
00457 *
00458       CALL SPOCON( UPLO, N, AF, LDAF, ANORM, RCOND, WORK, IWORK, INFO )
00459 *
00460 *     Compute the solution matrix X.
00461 *
00462       CALL SLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
00463       CALL SPOTRS( UPLO, N, NRHS, AF, LDAF, X, LDX, INFO )
00464 *
00465 *     Use iterative refinement to improve the computed solution and
00466 *     compute error bounds and backward error estimates for it.
00467 *
00468       CALL SPORFS( UPLO, N, NRHS, A, LDA, AF, LDAF, B, LDB, X, LDX,
00469      $             FERR, BERR, WORK, IWORK, INFO )
00470 *
00471 *     Transform the solution matrix X to a solution of the original
00472 *     system.
00473 *
00474       IF( RCEQU ) THEN
00475          DO 50 J = 1, NRHS
00476             DO 40 I = 1, N
00477                X( I, J ) = S( I )*X( I, J )
00478    40       CONTINUE
00479    50    CONTINUE
00480          DO 60 J = 1, NRHS
00481             FERR( J ) = FERR( J ) / SCOND
00482    60    CONTINUE
00483       END IF
00484 *
00485 *     Set INFO = N+1 if the matrix is singular to working precision.
00486 *
00487       IF( RCOND.LT.SLAMCH( 'Epsilon' ) )
00488      $   INFO = N + 1
00489 *
00490       RETURN
00491 *
00492 *     End of SPOSVX
00493 *
00494       END
 All Files Functions