--- rpl/lapack/lapack/dsgesv.f 2011/11/22 10:35:58 1.10 +++ rpl/lapack/lapack/dsgesv.f 2012/07/31 11:06:37 1.11 @@ -1,12 +1,205 @@ +*> \brief DSGESV computes the solution to system of linear equations A * X = B for GE matrices (mixed precision with iterative refinement) +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download DSGESV + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE DSGESV( N, NRHS, A, LDA, IPIV, B, LDB, X, LDX, WORK, +* SWORK, ITER, INFO ) +* +* .. Scalar Arguments .. +* INTEGER INFO, ITER, LDA, LDB, LDX, N, NRHS +* .. +* .. Array Arguments .. +* INTEGER IPIV( * ) +* REAL SWORK( * ) +* DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( N, * ), +* $ X( LDX, * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DSGESV computes the solution to a real system of linear equations +*> A * X = B, +*> where A is an N-by-N matrix and X and B are N-by-NRHS matrices. +*> +*> DSGESV 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] 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 N-by-N coefficient matrix A. +*> 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 factors L and U from the factorization +*> A = P*L*U; the unit diagonal elements of L are not stored. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,N). +*> \endverbatim +*> +*> \param[out] IPIV +*> \verbatim +*> IPIV is INTEGER array, dimension (N) +*> The pivot indices that define the permutation matrix P; +*> row i of the matrix was interchanged with row IPIV(i). +*> Corresponds either to the single precision factorization +*> (if INFO.EQ.0 and ITER.GE.0) or the double precision +*> factorization (if INFO.EQ.0 and ITER.LT.0). +*> \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 SGETRF +*> -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, U(i,i) computed in DOUBLE PRECISION is +*> exactly zero. The factorization has been completed, +*> but the factor U is exactly singular, so the solution +*> could not be computed. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup doubleGEsolve +* +* ===================================================================== SUBROUTINE DSGESV( N, NRHS, A, LDA, IPIV, B, LDB, X, LDX, WORK, $ SWORK, ITER, INFO ) * -* -- LAPACK PROTOTYPE driver routine (version 3.3.1) -- +* -- 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..-- -* -- April 2011 -- +* November 2011 * -* .. * .. Scalar Arguments .. INTEGER INFO, ITER, LDA, LDB, LDX, N, NRHS * .. @@ -17,111 +210,6 @@ $ X( LDX, * ) * .. * -* Purpose -* ======= -* -* DSGESV computes the solution to a real system of linear equations -* A * X = B, -* where A is an N-by-N matrix and X and B are N-by-NRHS matrices. -* -* DSGESV 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 -* ========= -* -* 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 N-by-N coefficient matrix A. -* 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 factors L and U from the factorization -* A = P*L*U; the unit diagonal elements of L are not stored. -* -* LDA (input) INTEGER -* The leading dimension of the array A. LDA >= max(1,N). -* -* IPIV (output) INTEGER array, dimension (N) -* The pivot indices that define the permutation matrix P; -* row i of the matrix was interchanged with row IPIV(i). -* Corresponds either to the single precision factorization -* (if INFO.EQ.0 and ITER.GE.0) or the double precision -* factorization (if INFO.EQ.0 and ITER.LT.0). -* -* 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 SGETRF -* -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, U(i,i) computed in DOUBLE PRECISION is -* exactly zero. The factorization has been completed, -* but the factor U is exactly singular, so the solution -* could not be computed. -* * ===================================================================== * * .. Parameters ..