LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dsposv.f
Go to the documentation of this file.
00001 *> \brief <b> DSPOSV 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 DSPOSV + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsposv.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsposv.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsposv.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE DSPOSV( UPLO, N, NRHS, A, LDA, B, LDB, X, LDX, WORK,
00022 *                          SWORK, ITER, INFO )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       CHARACTER          UPLO
00026 *       INTEGER            INFO, ITER, LDA, LDB, LDX, N, NRHS
00027 *       ..
00028 *       .. Array Arguments ..
00029 *       REAL               SWORK( * )
00030 *       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), WORK( N, * ),
00031 *      $                   X( LDX, * )
00032 *       ..
00033 *  
00034 *
00035 *> \par Purpose:
00036 *  =============
00037 *>
00038 *> \verbatim
00039 *>
00040 *> DSPOSV computes the solution to a real system of linear equations
00041 *>    A * X = B,
00042 *> where A is an N-by-N symmetric positive definite matrix and X and B
00043 *> are N-by-NRHS matrices.
00044 *>
00045 *> DSPOSV first attempts to factorize the matrix in SINGLE PRECISION
00046 *> and use this factorization within an iterative refinement procedure
00047 *> to produce a solution with DOUBLE PRECISION normwise backward error
00048 *> quality (see below). If the approach fails the method switches to a
00049 *> DOUBLE PRECISION factorization and solve.
00050 *>
00051 *> The iterative refinement is not going to be a winning strategy if
00052 *> the ratio SINGLE PRECISION performance over DOUBLE PRECISION
00053 *> performance is too small. A reasonable strategy should take the
00054 *> number of right-hand sides and the size of the matrix into account.
00055 *> This might be done with a call to ILAENV in the future. Up to now, we
00056 *> always try iterative refinement.
00057 *>
00058 *> The iterative refinement process is stopped if
00059 *>     ITER > ITERMAX
00060 *> or for all the RHS we have:
00061 *>     RNRM < SQRT(N)*XNRM*ANRM*EPS*BWDMAX
00062 *> where
00063 *>     o ITER is the number of the current iteration in the iterative
00064 *>       refinement process
00065 *>     o RNRM is the infinity-norm of the residual
00066 *>     o XNRM is the infinity-norm of the solution
00067 *>     o ANRM is the infinity-operator-norm of the matrix A
00068 *>     o EPS is the machine epsilon returned by DLAMCH('Epsilon')
00069 *> The value ITERMAX and BWDMAX are fixed to 30 and 1.0D+00
00070 *> respectively.
00071 *> \endverbatim
00072 *
00073 *  Arguments:
00074 *  ==========
00075 *
00076 *> \param[in] UPLO
00077 *> \verbatim
00078 *>          UPLO is CHARACTER*1
00079 *>          = 'U':  Upper triangle of A is stored;
00080 *>          = 'L':  Lower triangle of A is stored.
00081 *> \endverbatim
00082 *>
00083 *> \param[in] N
00084 *> \verbatim
00085 *>          N is INTEGER
00086 *>          The number of linear equations, i.e., the order of the
00087 *>          matrix A.  N >= 0.
00088 *> \endverbatim
00089 *>
00090 *> \param[in] NRHS
00091 *> \verbatim
00092 *>          NRHS is INTEGER
00093 *>          The number of right hand sides, i.e., the number of columns
00094 *>          of the matrix B.  NRHS >= 0.
00095 *> \endverbatim
00096 *>
00097 *> \param[in,out] A
00098 *> \verbatim
00099 *>          A is DOUBLE PRECISION array,
00100 *>          dimension (LDA,N)
00101 *>          On entry, the symmetric matrix A.  If UPLO = 'U', the leading
00102 *>          N-by-N upper triangular part of A contains the upper
00103 *>          triangular part of the matrix A, and the strictly lower
00104 *>          triangular part of A is not referenced.  If UPLO = 'L', the
00105 *>          leading N-by-N lower triangular part of A contains the lower
00106 *>          triangular part of the matrix A, and the strictly upper
00107 *>          triangular part of A is not referenced.
00108 *>          On exit, if iterative refinement has been successfully used
00109 *>          (INFO.EQ.0 and ITER.GE.0, see description below), then A is
00110 *>          unchanged, if double precision factorization has been used
00111 *>          (INFO.EQ.0 and ITER.LT.0, see description below), then the
00112 *>          array A contains the factor U or L from the Cholesky
00113 *>          factorization A = U**T*U or A = L*L**T.
00114 *> \endverbatim
00115 *>
00116 *> \param[in] LDA
00117 *> \verbatim
00118 *>          LDA is INTEGER
00119 *>          The leading dimension of the array A.  LDA >= max(1,N).
00120 *> \endverbatim
00121 *>
00122 *> \param[in] B
00123 *> \verbatim
00124 *>          B is DOUBLE PRECISION array, dimension (LDB,NRHS)
00125 *>          The N-by-NRHS right hand side matrix B.
00126 *> \endverbatim
00127 *>
00128 *> \param[in] LDB
00129 *> \verbatim
00130 *>          LDB is INTEGER
00131 *>          The leading dimension of the array B.  LDB >= max(1,N).
00132 *> \endverbatim
00133 *>
00134 *> \param[out] X
00135 *> \verbatim
00136 *>          X is DOUBLE PRECISION array, dimension (LDX,NRHS)
00137 *>          If INFO = 0, the N-by-NRHS solution matrix X.
00138 *> \endverbatim
00139 *>
00140 *> \param[in] LDX
00141 *> \verbatim
00142 *>          LDX is INTEGER
00143 *>          The leading dimension of the array X.  LDX >= max(1,N).
00144 *> \endverbatim
00145 *>
00146 *> \param[out] WORK
00147 *> \verbatim
00148 *>          WORK is DOUBLE PRECISION array, dimension (N,NRHS)
00149 *>          This array is used to hold the residual vectors.
00150 *> \endverbatim
00151 *>
00152 *> \param[out] SWORK
00153 *> \verbatim
00154 *>          SWORK is REAL array, dimension (N*(N+NRHS))
00155 *>          This array is used to use the single precision matrix and the
00156 *>          right-hand sides or solutions in single precision.
00157 *> \endverbatim
00158 *>
00159 *> \param[out] ITER
00160 *> \verbatim
00161 *>          ITER is INTEGER
00162 *>          < 0: iterative refinement has failed, double precision
00163 *>               factorization has been performed
00164 *>               -1 : the routine fell back to full precision for
00165 *>                    implementation- or machine-specific reasons
00166 *>               -2 : narrowing the precision induced an overflow,
00167 *>                    the routine fell back to full precision
00168 *>               -3 : failure of SPOTRF
00169 *>               -31: stop the iterative refinement after the 30th
00170 *>                    iterations
00171 *>          > 0: iterative refinement has been sucessfully used.
00172 *>               Returns the number of iterations
00173 *> \endverbatim
00174 *>
00175 *> \param[out] INFO
00176 *> \verbatim
00177 *>          INFO is INTEGER
00178 *>          = 0:  successful exit
00179 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
00180 *>          > 0:  if INFO = i, the leading minor of order i of (DOUBLE
00181 *>                PRECISION) A is not positive definite, so the
00182 *>                factorization could not be completed, and the solution
00183 *>                has not been computed.
00184 *> \endverbatim
00185 *
00186 *  Authors:
00187 *  ========
00188 *
00189 *> \author Univ. of Tennessee 
00190 *> \author Univ. of California Berkeley 
00191 *> \author Univ. of Colorado Denver 
00192 *> \author NAG Ltd. 
00193 *
00194 *> \date November 2011
00195 *
00196 *> \ingroup doublePOsolve
00197 *
00198 *  =====================================================================
00199       SUBROUTINE DSPOSV( UPLO, N, NRHS, A, LDA, B, LDB, X, LDX, WORK,
00200      $                   SWORK, ITER, INFO )
00201 *
00202 *  -- LAPACK driver routine (version 3.4.0) --
00203 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00204 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00205 *     November 2011
00206 *
00207 *     .. Scalar Arguments ..
00208       CHARACTER          UPLO
00209       INTEGER            INFO, ITER, LDA, LDB, LDX, N, NRHS
00210 *     ..
00211 *     .. Array Arguments ..
00212       REAL               SWORK( * )
00213       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), WORK( N, * ),
00214      $                   X( LDX, * )
00215 *     ..
00216 *
00217 *  =====================================================================
00218 *
00219 *     .. Parameters ..
00220       LOGICAL            DOITREF
00221       PARAMETER          ( DOITREF = .TRUE. )
00222 *
00223       INTEGER            ITERMAX
00224       PARAMETER          ( ITERMAX = 30 )
00225 *
00226       DOUBLE PRECISION   BWDMAX
00227       PARAMETER          ( BWDMAX = 1.0E+00 )
00228 *
00229       DOUBLE PRECISION   NEGONE, ONE
00230       PARAMETER          ( NEGONE = -1.0D+0, ONE = 1.0D+0 )
00231 *
00232 *     .. Local Scalars ..
00233       INTEGER            I, IITER, PTSA, PTSX
00234       DOUBLE PRECISION   ANRM, CTE, EPS, RNRM, XNRM
00235 *
00236 *     .. External Subroutines ..
00237       EXTERNAL           DAXPY, DSYMM, DLACPY, DLAT2S, DLAG2S, SLAG2D,
00238      $                   SPOTRF, SPOTRS, XERBLA
00239 *     ..
00240 *     .. External Functions ..
00241       INTEGER            IDAMAX
00242       DOUBLE PRECISION   DLAMCH, DLANSY
00243       LOGICAL            LSAME
00244       EXTERNAL           IDAMAX, DLAMCH, DLANSY, LSAME
00245 *     ..
00246 *     .. Intrinsic Functions ..
00247       INTRINSIC          ABS, DBLE, MAX, SQRT
00248 *     ..
00249 *     .. Executable Statements ..
00250 *
00251       INFO = 0
00252       ITER = 0
00253 *
00254 *     Test the input parameters.
00255 *
00256       IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00257          INFO = -1
00258       ELSE IF( N.LT.0 ) THEN
00259          INFO = -2
00260       ELSE IF( NRHS.LT.0 ) THEN
00261          INFO = -3
00262       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00263          INFO = -5
00264       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00265          INFO = -7
00266       ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
00267          INFO = -9
00268       END IF
00269       IF( INFO.NE.0 ) THEN
00270          CALL XERBLA( 'DSPOSV', -INFO )
00271          RETURN
00272       END IF
00273 *
00274 *     Quick return if (N.EQ.0).
00275 *
00276       IF( N.EQ.0 )
00277      $   RETURN
00278 *
00279 *     Skip single precision iterative refinement if a priori slower
00280 *     than double precision factorization.
00281 *
00282       IF( .NOT.DOITREF ) THEN
00283          ITER = -1
00284          GO TO 40
00285       END IF
00286 *
00287 *     Compute some constants.
00288 *
00289       ANRM = DLANSY( 'I', UPLO, N, A, LDA, WORK )
00290       EPS = DLAMCH( 'Epsilon' )
00291       CTE = ANRM*EPS*SQRT( DBLE( N ) )*BWDMAX
00292 *
00293 *     Set the indices PTSA, PTSX for referencing SA and SX in SWORK.
00294 *
00295       PTSA = 1
00296       PTSX = PTSA + N*N
00297 *
00298 *     Convert B from double precision to single precision and store the
00299 *     result in SX.
00300 *
00301       CALL DLAG2S( N, NRHS, B, LDB, SWORK( PTSX ), N, INFO )
00302 *
00303       IF( INFO.NE.0 ) THEN
00304          ITER = -2
00305          GO TO 40
00306       END IF
00307 *
00308 *     Convert A from double precision to single precision and store the
00309 *     result in SA.
00310 *
00311       CALL DLAT2S( UPLO, N, A, LDA, SWORK( PTSA ), N, INFO )
00312 *
00313       IF( INFO.NE.0 ) THEN
00314          ITER = -2
00315          GO TO 40
00316       END IF
00317 *
00318 *     Compute the Cholesky factorization of SA.
00319 *
00320       CALL SPOTRF( UPLO, N, SWORK( PTSA ), N, INFO )
00321 *
00322       IF( INFO.NE.0 ) THEN
00323          ITER = -3
00324          GO TO 40
00325       END IF
00326 *
00327 *     Solve the system SA*SX = SB.
00328 *
00329       CALL SPOTRS( UPLO, N, NRHS, SWORK( PTSA ), N, SWORK( PTSX ), N,
00330      $             INFO )
00331 *
00332 *     Convert SX back to double precision
00333 *
00334       CALL SLAG2D( N, NRHS, SWORK( PTSX ), N, X, LDX, INFO )
00335 *
00336 *     Compute R = B - AX (R is WORK).
00337 *
00338       CALL DLACPY( 'All', N, NRHS, B, LDB, WORK, N )
00339 *
00340       CALL DSYMM( 'Left', UPLO, N, NRHS, NEGONE, A, LDA, X, LDX, ONE,
00341      $            WORK, N )
00342 *
00343 *     Check whether the NRHS normwise backward errors satisfy the
00344 *     stopping criterion. If yes, set ITER=0 and return.
00345 *
00346       DO I = 1, NRHS
00347          XNRM = ABS( X( IDAMAX( N, X( 1, I ), 1 ), I ) )
00348          RNRM = ABS( WORK( IDAMAX( N, WORK( 1, I ), 1 ), I ) )
00349          IF( RNRM.GT.XNRM*CTE )
00350      $      GO TO 10
00351       END DO
00352 *
00353 *     If we are here, the NRHS normwise backward errors satisfy the
00354 *     stopping criterion. We are good to exit.
00355 *
00356       ITER = 0
00357       RETURN
00358 *
00359    10 CONTINUE
00360 *
00361       DO 30 IITER = 1, ITERMAX
00362 *
00363 *        Convert R (in WORK) from double precision to single precision
00364 *        and store the result in SX.
00365 *
00366          CALL DLAG2S( N, NRHS, WORK, N, SWORK( PTSX ), N, INFO )
00367 *
00368          IF( INFO.NE.0 ) THEN
00369             ITER = -2
00370             GO TO 40
00371          END IF
00372 *
00373 *        Solve the system SA*SX = SR.
00374 *
00375          CALL SPOTRS( UPLO, N, NRHS, SWORK( PTSA ), N, SWORK( PTSX ), N,
00376      $                INFO )
00377 *
00378 *        Convert SX back to double precision and update the current
00379 *        iterate.
00380 *
00381          CALL SLAG2D( N, NRHS, SWORK( PTSX ), N, WORK, N, INFO )
00382 *
00383          DO I = 1, NRHS
00384             CALL DAXPY( N, ONE, WORK( 1, I ), 1, X( 1, I ), 1 )
00385          END DO
00386 *
00387 *        Compute R = B - AX (R is WORK).
00388 *
00389          CALL DLACPY( 'All', N, NRHS, B, LDB, WORK, N )
00390 *
00391          CALL DSYMM( 'L', UPLO, N, NRHS, NEGONE, A, LDA, X, LDX, ONE,
00392      $               WORK, N )
00393 *
00394 *        Check whether the NRHS normwise backward errors satisfy the
00395 *        stopping criterion. If yes, set ITER=IITER>0 and return.
00396 *
00397          DO I = 1, NRHS
00398             XNRM = ABS( X( IDAMAX( N, X( 1, I ), 1 ), I ) )
00399             RNRM = ABS( WORK( IDAMAX( N, WORK( 1, I ), 1 ), I ) )
00400             IF( RNRM.GT.XNRM*CTE )
00401      $         GO TO 20
00402          END DO
00403 *
00404 *        If we are here, the NRHS normwise backward errors satisfy the
00405 *        stopping criterion, we are good to exit.
00406 *
00407          ITER = IITER
00408 *
00409          RETURN
00410 *
00411    20    CONTINUE
00412 *
00413    30 CONTINUE
00414 *
00415 *     If we are at this place of the code, this is because we have
00416 *     performed ITER=ITERMAX iterations and never satisified the
00417 *     stopping criterion, set up the ITER flag accordingly and follow
00418 *     up on double precision routine.
00419 *
00420       ITER = -ITERMAX - 1
00421 *
00422    40 CONTINUE
00423 *
00424 *     Single-precision iterative refinement failed to converge to a
00425 *     satisfactory solution, so we resort to double precision.
00426 *
00427       CALL DPOTRF( UPLO, N, A, LDA, INFO )
00428 *
00429       IF( INFO.NE.0 )
00430      $   RETURN
00431 *
00432       CALL DLACPY( 'All', N, NRHS, B, LDB, X, LDX )
00433       CALL DPOTRS( UPLO, N, NRHS, A, LDA, X, LDX, INFO )
00434 *
00435       RETURN
00436 *
00437 *     End of DSPOSV.
00438 *
00439       END
 All Files Functions