--- rpl/lapack/lapack/dsposv.f 2010/08/07 13:22:25 1.2
+++ rpl/lapack/lapack/dsposv.f 2012/12/14 14:22:40 1.10
@@ -1,11 +1,209 @@
+*> \brief DSPOSV computes the solution to system of linear equations A * X = B for PO matrices
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+*> \htmlonly
+*> Download DSPOSV + dependencies
+*>
+*> [TGZ]
+*>
+*> [ZIP]
+*>
+*> [TXT]
+*> \endhtmlonly
+*
+* Definition:
+* ===========
+*
+* SUBROUTINE DSPOSV( UPLO, N, NRHS, A, LDA, B, LDB, X, LDX, WORK,
+* SWORK, ITER, INFO )
+*
+* .. Scalar Arguments ..
+* CHARACTER UPLO
+* INTEGER INFO, ITER, LDA, LDB, LDX, N, NRHS
+* ..
+* .. Array Arguments ..
+* REAL SWORK( * )
+* DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( N, * ),
+* $ X( LDX, * )
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*>
+*> DSPOSV computes the solution to a real system of linear equations
+*> A * X = B,
+*> where A is an N-by-N symmetric positive definite matrix and X and B
+*> are N-by-NRHS matrices.
+*>
+*> DSPOSV first attempts to factorize the matrix in SINGLE PRECISION
+*> and use this factorization within an iterative refinement procedure
+*> to produce a solution with DOUBLE PRECISION normwise backward error
+*> quality (see below). If the approach fails the method switches to a
+*> DOUBLE PRECISION factorization and solve.
+*>
+*> The iterative refinement is not going to be a winning strategy if
+*> the ratio SINGLE PRECISION performance over DOUBLE PRECISION
+*> performance is too small. A reasonable strategy should take the
+*> number of right-hand sides and the size of the matrix into account.
+*> This might be done with a call to ILAENV in the future. Up to now, we
+*> always try iterative refinement.
+*>
+*> The iterative refinement process is stopped if
+*> ITER > ITERMAX
+*> or for all the RHS we have:
+*> RNRM < SQRT(N)*XNRM*ANRM*EPS*BWDMAX
+*> where
+*> o ITER is the number of the current iteration in the iterative
+*> refinement process
+*> o RNRM is the infinity-norm of the residual
+*> o XNRM is the infinity-norm of the solution
+*> o ANRM is the infinity-operator-norm of the matrix A
+*> o EPS is the machine epsilon returned by DLAMCH('Epsilon')
+*> The value ITERMAX and BWDMAX are fixed to 30 and 1.0D+00
+*> respectively.
+*> \endverbatim
+*
+* Arguments:
+* ==========
+*
+*> \param[in] UPLO
+*> \verbatim
+*> UPLO is CHARACTER*1
+*> = 'U': Upper triangle of A is stored;
+*> = 'L': Lower triangle of A is stored.
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*> N is INTEGER
+*> The number of linear equations, i.e., the order of the
+*> matrix A. N >= 0.
+*> \endverbatim
+*>
+*> \param[in] NRHS
+*> \verbatim
+*> NRHS is INTEGER
+*> The number of right hand sides, i.e., the number of columns
+*> of the matrix B. NRHS >= 0.
+*> \endverbatim
+*>
+*> \param[in,out] A
+*> \verbatim
+*> A is DOUBLE PRECISION array,
+*> dimension (LDA,N)
+*> On entry, the symmetric matrix A. If UPLO = 'U', the leading
+*> N-by-N upper triangular part of A contains the upper
+*> triangular part of the matrix A, and the strictly lower
+*> triangular part of A is not referenced. If UPLO = 'L', the
+*> leading N-by-N lower triangular part of A contains the lower
+*> triangular part of the matrix A, and the strictly upper
+*> triangular part of A is not referenced.
+*> On exit, if iterative refinement has been successfully used
+*> (INFO.EQ.0 and ITER.GE.0, see description below), then A is
+*> unchanged, if double precision factorization has been used
+*> (INFO.EQ.0 and ITER.LT.0, see description below), then the
+*> array A contains the factor U or L from the Cholesky
+*> factorization A = U**T*U or A = L*L**T.
+*> \endverbatim
+*>
+*> \param[in] LDA
+*> \verbatim
+*> LDA is INTEGER
+*> The leading dimension of the array A. LDA >= max(1,N).
+*> \endverbatim
+*>
+*> \param[in] B
+*> \verbatim
+*> B is DOUBLE PRECISION array, dimension (LDB,NRHS)
+*> The N-by-NRHS right hand side matrix B.
+*> \endverbatim
+*>
+*> \param[in] LDB
+*> \verbatim
+*> LDB is INTEGER
+*> The leading dimension of the array B. LDB >= max(1,N).
+*> \endverbatim
+*>
+*> \param[out] X
+*> \verbatim
+*> X is DOUBLE PRECISION array, dimension (LDX,NRHS)
+*> If INFO = 0, the N-by-NRHS solution matrix X.
+*> \endverbatim
+*>
+*> \param[in] LDX
+*> \verbatim
+*> LDX is INTEGER
+*> The leading dimension of the array X. LDX >= max(1,N).
+*> \endverbatim
+*>
+*> \param[out] WORK
+*> \verbatim
+*> WORK is DOUBLE PRECISION array, dimension (N,NRHS)
+*> This array is used to hold the residual vectors.
+*> \endverbatim
+*>
+*> \param[out] SWORK
+*> \verbatim
+*> SWORK is REAL array, dimension (N*(N+NRHS))
+*> This array is used to use the single precision matrix and the
+*> right-hand sides or solutions in single precision.
+*> \endverbatim
+*>
+*> \param[out] ITER
+*> \verbatim
+*> ITER is INTEGER
+*> < 0: iterative refinement has failed, double precision
+*> factorization has been performed
+*> -1 : the routine fell back to full precision for
+*> implementation- or machine-specific reasons
+*> -2 : narrowing the precision induced an overflow,
+*> the routine fell back to full precision
+*> -3 : failure of SPOTRF
+*> -31: stop the iterative refinement after the 30th
+*> iterations
+*> > 0: iterative refinement has been sucessfully used.
+*> Returns the number of iterations
+*> \endverbatim
+*>
+*> \param[out] INFO
+*> \verbatim
+*> INFO is INTEGER
+*> = 0: successful exit
+*> < 0: if INFO = -i, the i-th argument had an illegal value
+*> > 0: if INFO = i, the leading minor of order i of (DOUBLE
+*> PRECISION) A is not positive definite, so the
+*> factorization could not be completed, and the solution
+*> has not been computed.
+*> \endverbatim
+*
+* Authors:
+* ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \date November 2011
+*
+*> \ingroup doublePOsolve
+*
+* =====================================================================
SUBROUTINE DSPOSV( UPLO, N, NRHS, A, LDA, B, LDB, X, LDX, WORK,
- + SWORK, ITER, INFO )
+ $ SWORK, ITER, INFO )
*
-* -- LAPACK PROTOTYPE driver routine (version 3.1.2) --
-* Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..
-* May 2007
+* -- LAPACK driver routine (version 3.4.0) --
+* -- LAPACK is a software package provided by Univ. of Tennessee, --
+* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+* November 2011
*
-* ..
* .. Scalar Arguments ..
CHARACTER UPLO
INTEGER INFO, ITER, LDA, LDB, LDX, N, NRHS
@@ -13,120 +211,10 @@
* .. Array Arguments ..
REAL SWORK( * )
DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( N, * ),
- + X( LDX, * )
+ $ X( LDX, * )
* ..
*
-* Purpose
-* =======
-*
-* DSPOSV computes the solution to a real system of linear equations
-* A * X = B,
-* where A is an N-by-N symmetric positive definite matrix and X and B
-* are N-by-NRHS matrices.
-*
-* DSPOSV first attempts to factorize the matrix in SINGLE PRECISION
-* and use this factorization within an iterative refinement procedure
-* to produce a solution with DOUBLE PRECISION normwise backward error
-* quality (see below). If the approach fails the method switches to a
-* DOUBLE PRECISION factorization and solve.
-*
-* The iterative refinement is not going to be a winning strategy if
-* the ratio SINGLE PRECISION performance over DOUBLE PRECISION
-* performance is too small. A reasonable strategy should take the
-* number of right-hand sides and the size of the matrix into account.
-* This might be done with a call to ILAENV in the future. Up to now, we
-* always try iterative refinement.
-*
-* The iterative refinement process is stopped if
-* ITER > ITERMAX
-* or for all the RHS we have:
-* RNRM < SQRT(N)*XNRM*ANRM*EPS*BWDMAX
-* where
-* o ITER is the number of the current iteration in the iterative
-* refinement process
-* o RNRM is the infinity-norm of the residual
-* o XNRM is the infinity-norm of the solution
-* o ANRM is the infinity-operator-norm of the matrix A
-* o EPS is the machine epsilon returned by DLAMCH('Epsilon')
-* The value ITERMAX and BWDMAX are fixed to 30 and 1.0D+00
-* respectively.
-*
-* Arguments
-* =========
-*
-* UPLO (input) CHARACTER
-* = 'U': Upper triangle of A is stored;
-* = 'L': Lower triangle of A is stored.
-*
-* N (input) INTEGER
-* The number of linear equations, i.e., the order of the
-* matrix A. N >= 0.
-*
-* NRHS (input) INTEGER
-* The number of right hand sides, i.e., the number of columns
-* of the matrix B. NRHS >= 0.
-*
-* A (input/output) DOUBLE PRECISION array,
-* dimension (LDA,N)
-* On entry, the symmetric matrix A. If UPLO = 'U', the leading
-* N-by-N upper triangular part of A contains the upper
-* triangular part of the matrix A, and the strictly lower
-* triangular part of A is not referenced. If UPLO = 'L', the
-* leading N-by-N lower triangular part of A contains the lower
-* triangular part of the matrix A, and the strictly upper
-* triangular part of A is not referenced.
-* On exit, if iterative refinement has been successfully used
-* (INFO.EQ.0 and ITER.GE.0, see description below), then A is
-* unchanged, if double precision factorization has been used
-* (INFO.EQ.0 and ITER.LT.0, see description below), then the
-* array A contains the factor U or L from the Cholesky
-* factorization A = U**T*U or A = L*L**T.
-*
-*
-* LDA (input) INTEGER
-* The leading dimension of the array A. LDA >= max(1,N).
-*
-* B (input) DOUBLE PRECISION array, dimension (LDB,NRHS)
-* The N-by-NRHS right hand side matrix B.
-*
-* LDB (input) INTEGER
-* The leading dimension of the array B. LDB >= max(1,N).
-*
-* X (output) DOUBLE PRECISION array, dimension (LDX,NRHS)
-* If INFO = 0, the N-by-NRHS solution matrix X.
-*
-* LDX (input) INTEGER
-* The leading dimension of the array X. LDX >= max(1,N).
-*
-* WORK (workspace) DOUBLE PRECISION array, dimension (N,NRHS)
-* This array is used to hold the residual vectors.
-*
-* SWORK (workspace) REAL array, dimension (N*(N+NRHS))
-* This array is used to use the single precision matrix and the
-* right-hand sides or solutions in single precision.
-*
-* ITER (output) INTEGER
-* < 0: iterative refinement has failed, double precision
-* factorization has been performed
-* -1 : the routine fell back to full precision for
-* implementation- or machine-specific reasons
-* -2 : narrowing the precision induced an overflow,
-* the routine fell back to full precision
-* -3 : failure of SPOTRF
-* -31: stop the iterative refinement after the 30th
-* iterations
-* > 0: iterative refinement has been sucessfully used.
-* Returns the number of iterations
-*
-* INFO (output) INTEGER
-* = 0: successful exit
-* < 0: if INFO = -i, the i-th argument had an illegal value
-* > 0: if INFO = i, the leading minor of order i of (DOUBLE
-* PRECISION) A is not positive definite, so the
-* factorization could not be completed, and the solution
-* has not been computed.
-*
-* =========
+* =====================================================================
*
* .. Parameters ..
LOGICAL DOITREF
@@ -147,7 +235,7 @@
*
* .. External Subroutines ..
EXTERNAL DAXPY, DSYMM, DLACPY, DLAT2S, DLAG2S, SLAG2D,
- + SPOTRF, SPOTRS, XERBLA
+ $ SPOTRF, SPOTRS, XERBLA
* ..
* .. External Functions ..
INTEGER IDAMAX
@@ -186,7 +274,7 @@
* Quick return if (N.EQ.0).
*
IF( N.EQ.0 )
- + RETURN
+ $ RETURN
*
* Skip single precision iterative refinement if a priori slower
* than double precision factorization.
@@ -239,7 +327,7 @@
* Solve the system SA*SX = SB.
*
CALL SPOTRS( UPLO, N, NRHS, SWORK( PTSA ), N, SWORK( PTSX ), N,
- + INFO )
+ $ INFO )
*
* Convert SX back to double precision
*
@@ -250,7 +338,7 @@
CALL DLACPY( 'All', N, NRHS, B, LDB, WORK, N )
*
CALL DSYMM( 'Left', UPLO, N, NRHS, NEGONE, A, LDA, X, LDX, ONE,
- + WORK, N )
+ $ WORK, N )
*
* Check whether the NRHS normwise backward errors satisfy the
* stopping criterion. If yes, set ITER=0 and return.
@@ -259,7 +347,7 @@
XNRM = ABS( X( IDAMAX( N, X( 1, I ), 1 ), I ) )
RNRM = ABS( WORK( IDAMAX( N, WORK( 1, I ), 1 ), I ) )
IF( RNRM.GT.XNRM*CTE )
- + GO TO 10
+ $ GO TO 10
END DO
*
* If we are here, the NRHS normwise backward errors satisfy the
@@ -285,7 +373,7 @@
* Solve the system SA*SX = SR.
*
CALL SPOTRS( UPLO, N, NRHS, SWORK( PTSA ), N, SWORK( PTSX ), N,
- + INFO )
+ $ INFO )
*
* Convert SX back to double precision and update the current
* iterate.
@@ -301,7 +389,7 @@
CALL DLACPY( 'All', N, NRHS, B, LDB, WORK, N )
*
CALL DSYMM( 'L', UPLO, N, NRHS, NEGONE, A, LDA, X, LDX, ONE,
- + WORK, N )
+ $ WORK, N )
*
* Check whether the NRHS normwise backward errors satisfy the
* stopping criterion. If yes, set ITER=IITER>0 and return.
@@ -310,7 +398,7 @@
XNRM = ABS( X( IDAMAX( N, X( 1, I ), 1 ), I ) )
RNRM = ABS( WORK( IDAMAX( N, WORK( 1, I ), 1 ), I ) )
IF( RNRM.GT.XNRM*CTE )
- + GO TO 20
+ $ GO TO 20
END DO
*
* If we are here, the NRHS normwise backward errors satisfy the
@@ -339,7 +427,7 @@
CALL DPOTRF( UPLO, N, A, LDA, INFO )
*
IF( INFO.NE.0 )
- + RETURN
+ $ RETURN
*
CALL DLACPY( 'All', N, NRHS, B, LDB, X, LDX )
CALL DPOTRS( UPLO, N, NRHS, A, LDA, X, LDX, INFO )