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