--- rpl/lapack/lapack/dsgesv.f 2010/04/21 13:45:24 1.2
+++ rpl/lapack/lapack/dsgesv.f 2020/05/21 21:46:01 1.21
@@ -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 = 0 and ITER >= 0, see description below), then A is
+*> unchanged, if double precision factorization has been used
+*> (INFO = 0 and ITER < 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 = 0 and ITER >= 0) or the double precision
+*> factorization (if INFO = 0 and ITER < 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 successfully 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 June 2016
+*
+*> \ingroup doubleGEsolve
+*
+* =====================================================================
SUBROUTINE DSGESV( N, NRHS, A, LDA, IPIV, B, LDB, X, LDX, WORK,
- + SWORK, ITER, INFO )
+ $ SWORK, ITER, INFO )
*
-* -- LAPACK PROTOTYPE driver routine (version 3.2) --
+* -- LAPACK driver routine (version 3.8.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-* February 2007
+* June 2016
*
-* ..
* .. Scalar Arguments ..
INTEGER INFO, ITER, LDA, LDB, LDX, N, NRHS
* ..
@@ -14,115 +207,10 @@
INTEGER IPIV( * )
REAL SWORK( * )
DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( N, * ),
- + X( LDX, * )
+ $ 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 or input/ouptut) 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 ..
LOGICAL DOITREF
@@ -142,8 +230,8 @@
DOUBLE PRECISION ANRM, CTE, EPS, RNRM, XNRM
*
* .. External Subroutines ..
- EXTERNAL DAXPY, DGEMM, DLACPY, DLAG2S, SLAG2D, SGETRF,
- + SGETRS, XERBLA
+ EXTERNAL DAXPY, DGEMM, DLACPY, DLAG2S, DGETRF, DGETRS,
+ $ SGETRF, SGETRS, SLAG2D, XERBLA
* ..
* .. External Functions ..
INTEGER IDAMAX
@@ -179,7 +267,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.
@@ -232,7 +320,7 @@
* Solve the system SA*SX = SB.
*
CALL SGETRS( 'No transpose', N, NRHS, SWORK( PTSA ), N, IPIV,
- + SWORK( PTSX ), N, INFO )
+ $ SWORK( PTSX ), N, INFO )
*
* Convert SX back to double precision
*
@@ -243,7 +331,7 @@
CALL DLACPY( 'All', N, NRHS, B, LDB, WORK, N )
*
CALL DGEMM( 'No Transpose', 'No Transpose', N, NRHS, N, NEGONE, A,
- + LDA, X, LDX, ONE, WORK, N )
+ $ LDA, X, LDX, ONE, WORK, N )
*
* Check whether the NRHS normwise backward errors satisfy the
* stopping criterion. If yes, set ITER=0 and return.
@@ -252,7 +340,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
@@ -278,7 +366,7 @@
* Solve the system SA*SX = SR.
*
CALL SGETRS( 'No transpose', N, NRHS, SWORK( PTSA ), N, IPIV,
- + SWORK( PTSX ), N, INFO )
+ $ SWORK( PTSX ), N, INFO )
*
* Convert SX back to double precision and update the current
* iterate.
@@ -294,7 +382,7 @@
CALL DLACPY( 'All', N, NRHS, B, LDB, WORK, N )
*
CALL DGEMM( 'No Transpose', 'No Transpose', N, NRHS, N, NEGONE,
- + A, LDA, X, LDX, ONE, WORK, N )
+ $ A, LDA, X, LDX, ONE, WORK, N )
*
* Check whether the NRHS normwise backward errors satisfy the
* stopping criterion. If yes, set ITER=IITER>0 and return.
@@ -303,7 +391,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
@@ -318,7 +406,7 @@
30 CONTINUE
*
* If we are at this place of the code, this is because we have
-* performed ITER=ITERMAX iterations and never satisified the
+* performed ITER=ITERMAX iterations and never satisfied the
* stopping criterion, set up the ITER flag accordingly and follow up
* on double precision routine.
*
@@ -332,11 +420,11 @@
CALL DGETRF( N, N, A, LDA, IPIV, INFO )
*
IF( INFO.NE.0 )
- + RETURN
+ $ RETURN
*
CALL DLACPY( 'All', N, NRHS, B, LDB, X, LDX )
CALL DGETRS( 'No transpose', N, NRHS, A, LDA, IPIV, X, LDX,
- + INFO )
+ $ INFO )
*
RETURN
*