LAPACK 3.3.1 Linear Algebra PACKage

# dsposv.f

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