Annotation of rpl/lapack/lapack/dsposv.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE DSPOSV( UPLO, N, NRHS, A, LDA, B, LDB, X, LDX, WORK,
! 2: + SWORK, ITER, INFO )
! 3: *
! 4: * -- LAPACK PROTOTYPE driver routine (version 3.1.2) --
! 5: * Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..
! 6: * May 2007
! 7: *
! 8: * ..
! 9: * .. Scalar Arguments ..
! 10: CHARACTER UPLO
! 11: INTEGER INFO, ITER, LDA, LDB, LDX, N, NRHS
! 12: * ..
! 13: * .. Array Arguments ..
! 14: REAL SWORK( * )
! 15: DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( N, * ),
! 16: + X( LDX, * )
! 17: * ..
! 18: *
! 19: * Purpose
! 20: * =======
! 21: *
! 22: * DSPOSV computes the solution to a real system of linear equations
! 23: * A * X = B,
! 24: * where A is an N-by-N symmetric positive definite matrix and X and B
! 25: * are N-by-NRHS matrices.
! 26: *
! 27: * DSPOSV first attempts to factorize the matrix in SINGLE PRECISION
! 28: * and use this factorization within an iterative refinement procedure
! 29: * to produce a solution with DOUBLE PRECISION normwise backward error
! 30: * quality (see below). If the approach fails the method switches to a
! 31: * DOUBLE PRECISION factorization and solve.
! 32: *
! 33: * The iterative refinement is not going to be a winning strategy if
! 34: * the ratio SINGLE PRECISION performance over DOUBLE PRECISION
! 35: * performance is too small. A reasonable strategy should take the
! 36: * number of right-hand sides and the size of the matrix into account.
! 37: * This might be done with a call to ILAENV in the future. Up to now, we
! 38: * always try iterative refinement.
! 39: *
! 40: * The iterative refinement process is stopped if
! 41: * ITER > ITERMAX
! 42: * or for all the RHS we have:
! 43: * RNRM < SQRT(N)*XNRM*ANRM*EPS*BWDMAX
! 44: * where
! 45: * o ITER is the number of the current iteration in the iterative
! 46: * refinement process
! 47: * o RNRM is the infinity-norm of the residual
! 48: * o XNRM is the infinity-norm of the solution
! 49: * o ANRM is the infinity-operator-norm of the matrix A
! 50: * o EPS is the machine epsilon returned by DLAMCH('Epsilon')
! 51: * The value ITERMAX and BWDMAX are fixed to 30 and 1.0D+00
! 52: * respectively.
! 53: *
! 54: * Arguments
! 55: * =========
! 56: *
! 57: * UPLO (input) CHARACTER
! 58: * = 'U': Upper triangle of A is stored;
! 59: * = 'L': Lower triangle of A is stored.
! 60: *
! 61: * N (input) INTEGER
! 62: * The number of linear equations, i.e., the order of the
! 63: * matrix A. N >= 0.
! 64: *
! 65: * NRHS (input) INTEGER
! 66: * The number of right hand sides, i.e., the number of columns
! 67: * of the matrix B. NRHS >= 0.
! 68: *
! 69: * A (input/output) DOUBLE PRECISION array,
! 70: * dimension (LDA,N)
! 71: * On entry, the symmetric matrix A. If UPLO = 'U', the leading
! 72: * N-by-N upper triangular part of A contains the upper
! 73: * triangular part of the matrix A, and the strictly lower
! 74: * triangular part of A is not referenced. If UPLO = 'L', the
! 75: * leading N-by-N lower triangular part of A contains the lower
! 76: * triangular part of the matrix A, and the strictly upper
! 77: * triangular part of A is not referenced.
! 78: * On exit, if iterative refinement has been successfully used
! 79: * (INFO.EQ.0 and ITER.GE.0, see description below), then A is
! 80: * unchanged, if double precision factorization has been used
! 81: * (INFO.EQ.0 and ITER.LT.0, see description below), then the
! 82: * array A contains the factor U or L from the Cholesky
! 83: * factorization A = U**T*U or A = L*L**T.
! 84: *
! 85: *
! 86: * LDA (input) INTEGER
! 87: * The leading dimension of the array A. LDA >= max(1,N).
! 88: *
! 89: * B (input) DOUBLE PRECISION array, dimension (LDB,NRHS)
! 90: * The N-by-NRHS right hand side matrix B.
! 91: *
! 92: * LDB (input) INTEGER
! 93: * The leading dimension of the array B. LDB >= max(1,N).
! 94: *
! 95: * X (output) DOUBLE PRECISION array, dimension (LDX,NRHS)
! 96: * If INFO = 0, the N-by-NRHS solution matrix X.
! 97: *
! 98: * LDX (input) INTEGER
! 99: * The leading dimension of the array X. LDX >= max(1,N).
! 100: *
! 101: * WORK (workspace) DOUBLE PRECISION array, dimension (N,NRHS)
! 102: * This array is used to hold the residual vectors.
! 103: *
! 104: * SWORK (workspace) REAL array, dimension (N*(N+NRHS))
! 105: * This array is used to use the single precision matrix and the
! 106: * right-hand sides or solutions in single precision.
! 107: *
! 108: * ITER (output) INTEGER
! 109: * < 0: iterative refinement has failed, double precision
! 110: * factorization has been performed
! 111: * -1 : the routine fell back to full precision for
! 112: * implementation- or machine-specific reasons
! 113: * -2 : narrowing the precision induced an overflow,
! 114: * the routine fell back to full precision
! 115: * -3 : failure of SPOTRF
! 116: * -31: stop the iterative refinement after the 30th
! 117: * iterations
! 118: * > 0: iterative refinement has been sucessfully used.
! 119: * Returns the number of iterations
! 120: *
! 121: * INFO (output) INTEGER
! 122: * = 0: successful exit
! 123: * < 0: if INFO = -i, the i-th argument had an illegal value
! 124: * > 0: if INFO = i, the leading minor of order i of (DOUBLE
! 125: * PRECISION) A is not positive definite, so the
! 126: * factorization could not be completed, and the solution
! 127: * has not been computed.
! 128: *
! 129: * =========
! 130: *
! 131: * .. Parameters ..
! 132: LOGICAL DOITREF
! 133: PARAMETER ( DOITREF = .TRUE. )
! 134: *
! 135: INTEGER ITERMAX
! 136: PARAMETER ( ITERMAX = 30 )
! 137: *
! 138: DOUBLE PRECISION BWDMAX
! 139: PARAMETER ( BWDMAX = 1.0E+00 )
! 140: *
! 141: DOUBLE PRECISION NEGONE, ONE
! 142: PARAMETER ( NEGONE = -1.0D+0, ONE = 1.0D+0 )
! 143: *
! 144: * .. Local Scalars ..
! 145: INTEGER I, IITER, PTSA, PTSX
! 146: DOUBLE PRECISION ANRM, CTE, EPS, RNRM, XNRM
! 147: *
! 148: * .. External Subroutines ..
! 149: EXTERNAL DAXPY, DSYMM, DLACPY, DLAT2S, DLAG2S, SLAG2D,
! 150: + SPOTRF, SPOTRS, XERBLA
! 151: * ..
! 152: * .. External Functions ..
! 153: INTEGER IDAMAX
! 154: DOUBLE PRECISION DLAMCH, DLANSY
! 155: LOGICAL LSAME
! 156: EXTERNAL IDAMAX, DLAMCH, DLANSY, LSAME
! 157: * ..
! 158: * .. Intrinsic Functions ..
! 159: INTRINSIC ABS, DBLE, MAX, SQRT
! 160: * ..
! 161: * .. Executable Statements ..
! 162: *
! 163: INFO = 0
! 164: ITER = 0
! 165: *
! 166: * Test the input parameters.
! 167: *
! 168: IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
! 169: INFO = -1
! 170: ELSE IF( N.LT.0 ) THEN
! 171: INFO = -2
! 172: ELSE IF( NRHS.LT.0 ) THEN
! 173: INFO = -3
! 174: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
! 175: INFO = -5
! 176: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
! 177: INFO = -7
! 178: ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
! 179: INFO = -9
! 180: END IF
! 181: IF( INFO.NE.0 ) THEN
! 182: CALL XERBLA( 'DSPOSV', -INFO )
! 183: RETURN
! 184: END IF
! 185: *
! 186: * Quick return if (N.EQ.0).
! 187: *
! 188: IF( N.EQ.0 )
! 189: + RETURN
! 190: *
! 191: * Skip single precision iterative refinement if a priori slower
! 192: * than double precision factorization.
! 193: *
! 194: IF( .NOT.DOITREF ) THEN
! 195: ITER = -1
! 196: GO TO 40
! 197: END IF
! 198: *
! 199: * Compute some constants.
! 200: *
! 201: ANRM = DLANSY( 'I', UPLO, N, A, LDA, WORK )
! 202: EPS = DLAMCH( 'Epsilon' )
! 203: CTE = ANRM*EPS*SQRT( DBLE( N ) )*BWDMAX
! 204: *
! 205: * Set the indices PTSA, PTSX for referencing SA and SX in SWORK.
! 206: *
! 207: PTSA = 1
! 208: PTSX = PTSA + N*N
! 209: *
! 210: * Convert B from double precision to single precision and store the
! 211: * result in SX.
! 212: *
! 213: CALL DLAG2S( N, NRHS, B, LDB, SWORK( PTSX ), N, INFO )
! 214: *
! 215: IF( INFO.NE.0 ) THEN
! 216: ITER = -2
! 217: GO TO 40
! 218: END IF
! 219: *
! 220: * Convert A from double precision to single precision and store the
! 221: * result in SA.
! 222: *
! 223: CALL DLAT2S( UPLO, N, A, LDA, SWORK( PTSA ), N, INFO )
! 224: *
! 225: IF( INFO.NE.0 ) THEN
! 226: ITER = -2
! 227: GO TO 40
! 228: END IF
! 229: *
! 230: * Compute the Cholesky factorization of SA.
! 231: *
! 232: CALL SPOTRF( UPLO, N, SWORK( PTSA ), N, INFO )
! 233: *
! 234: IF( INFO.NE.0 ) THEN
! 235: ITER = -3
! 236: GO TO 40
! 237: END IF
! 238: *
! 239: * Solve the system SA*SX = SB.
! 240: *
! 241: CALL SPOTRS( UPLO, N, NRHS, SWORK( PTSA ), N, SWORK( PTSX ), N,
! 242: + INFO )
! 243: *
! 244: * Convert SX back to double precision
! 245: *
! 246: CALL SLAG2D( N, NRHS, SWORK( PTSX ), N, X, LDX, INFO )
! 247: *
! 248: * Compute R = B - AX (R is WORK).
! 249: *
! 250: CALL DLACPY( 'All', N, NRHS, B, LDB, WORK, N )
! 251: *
! 252: CALL DSYMM( 'Left', UPLO, N, NRHS, NEGONE, A, LDA, X, LDX, ONE,
! 253: + WORK, N )
! 254: *
! 255: * Check whether the NRHS normwise backward errors satisfy the
! 256: * stopping criterion. If yes, set ITER=0 and return.
! 257: *
! 258: DO I = 1, NRHS
! 259: XNRM = ABS( X( IDAMAX( N, X( 1, I ), 1 ), I ) )
! 260: RNRM = ABS( WORK( IDAMAX( N, WORK( 1, I ), 1 ), I ) )
! 261: IF( RNRM.GT.XNRM*CTE )
! 262: + GO TO 10
! 263: END DO
! 264: *
! 265: * If we are here, the NRHS normwise backward errors satisfy the
! 266: * stopping criterion. We are good to exit.
! 267: *
! 268: ITER = 0
! 269: RETURN
! 270: *
! 271: 10 CONTINUE
! 272: *
! 273: DO 30 IITER = 1, ITERMAX
! 274: *
! 275: * Convert R (in WORK) from double precision to single precision
! 276: * and store the result in SX.
! 277: *
! 278: CALL DLAG2S( N, NRHS, WORK, N, SWORK( PTSX ), N, INFO )
! 279: *
! 280: IF( INFO.NE.0 ) THEN
! 281: ITER = -2
! 282: GO TO 40
! 283: END IF
! 284: *
! 285: * Solve the system SA*SX = SR.
! 286: *
! 287: CALL SPOTRS( UPLO, N, NRHS, SWORK( PTSA ), N, SWORK( PTSX ), N,
! 288: + INFO )
! 289: *
! 290: * Convert SX back to double precision and update the current
! 291: * iterate.
! 292: *
! 293: CALL SLAG2D( N, NRHS, SWORK( PTSX ), N, WORK, N, INFO )
! 294: *
! 295: DO I = 1, NRHS
! 296: CALL DAXPY( N, ONE, WORK( 1, I ), 1, X( 1, I ), 1 )
! 297: END DO
! 298: *
! 299: * Compute R = B - AX (R is WORK).
! 300: *
! 301: CALL DLACPY( 'All', N, NRHS, B, LDB, WORK, N )
! 302: *
! 303: CALL DSYMM( 'L', UPLO, N, NRHS, NEGONE, A, LDA, X, LDX, ONE,
! 304: + WORK, N )
! 305: *
! 306: * Check whether the NRHS normwise backward errors satisfy the
! 307: * stopping criterion. If yes, set ITER=IITER>0 and return.
! 308: *
! 309: DO I = 1, NRHS
! 310: XNRM = ABS( X( IDAMAX( N, X( 1, I ), 1 ), I ) )
! 311: RNRM = ABS( WORK( IDAMAX( N, WORK( 1, I ), 1 ), I ) )
! 312: IF( RNRM.GT.XNRM*CTE )
! 313: + GO TO 20
! 314: END DO
! 315: *
! 316: * If we are here, the NRHS normwise backward errors satisfy the
! 317: * stopping criterion, we are good to exit.
! 318: *
! 319: ITER = IITER
! 320: *
! 321: RETURN
! 322: *
! 323: 20 CONTINUE
! 324: *
! 325: 30 CONTINUE
! 326: *
! 327: * If we are at this place of the code, this is because we have
! 328: * performed ITER=ITERMAX iterations and never satisified the
! 329: * stopping criterion, set up the ITER flag accordingly and follow
! 330: * up on double precision routine.
! 331: *
! 332: ITER = -ITERMAX - 1
! 333: *
! 334: 40 CONTINUE
! 335: *
! 336: * Single-precision iterative refinement failed to converge to a
! 337: * satisfactory solution, so we resort to double precision.
! 338: *
! 339: CALL DPOTRF( UPLO, N, A, LDA, INFO )
! 340: *
! 341: IF( INFO.NE.0 )
! 342: + RETURN
! 343: *
! 344: CALL DLACPY( 'All', N, NRHS, B, LDB, X, LDX )
! 345: CALL DPOTRS( UPLO, N, NRHS, A, LDA, X, LDX, INFO )
! 346: *
! 347: RETURN
! 348: *
! 349: * End of DSPOSV.
! 350: *
! 351: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dsposv.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack