--- rpl/lapack/lapack/zcgesv.f 2010/08/13 21:04:01 1.7
+++ rpl/lapack/lapack/zcgesv.f 2012/08/22 09:48:28 1.12
@@ -1,12 +1,211 @@
+*> \brief ZCGESV 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 ZCGESV + dependencies
+*>
+*> [TGZ]
+*>
+*> [ZIP]
+*>
+*> [TXT]
+*> \endhtmlonly
+*
+* Definition:
+* ===========
+*
+* SUBROUTINE ZCGESV( N, NRHS, A, LDA, IPIV, B, LDB, X, LDX, WORK,
+* SWORK, RWORK, ITER, INFO )
+*
+* .. Scalar Arguments ..
+* INTEGER INFO, ITER, LDA, LDB, LDX, N, NRHS
+* ..
+* .. Array Arguments ..
+* INTEGER IPIV( * )
+* DOUBLE PRECISION RWORK( * )
+* COMPLEX SWORK( * )
+* COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( N, * ),
+* $ X( LDX, * )
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*>
+*> ZCGESV computes the solution to a complex system of linear equations
+*> A * X = B,
+*> where A is an N-by-N matrix and X and B are N-by-NRHS matrices.
+*>
+*> ZCGESV first attempts to factorize the matrix in COMPLEX and use this
+*> factorization within an iterative refinement procedure to produce a
+*> solution with COMPLEX*16 normwise backward error quality (see below).
+*> If the approach fails the method switches to a COMPLEX*16
+*> factorization and solve.
+*>
+*> The iterative refinement is not going to be a winning strategy if
+*> the ratio COMPLEX performance over COMPLEX*16 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 COMPLEX*16 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 COMPLEX*16 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 COMPLEX*16 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 COMPLEX*16 array, dimension (N*NRHS)
+*> This array is used to hold the residual vectors.
+*> \endverbatim
+*>
+*> \param[out] SWORK
+*> \verbatim
+*> SWORK is COMPLEX 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] RWORK
+*> \verbatim
+*> RWORK is DOUBLE PRECISION array, dimension (N)
+*> \endverbatim
+*>
+*> \param[out] ITER
+*> \verbatim
+*> ITER is INTEGER
+*> < 0: iterative refinement has failed, COMPLEX*16
+*> 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 CGETRF
+*> -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 COMPLEX*16 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 complex16GEsolve
+*
+* =====================================================================
SUBROUTINE ZCGESV( N, NRHS, A, LDA, IPIV, B, LDB, X, LDX, WORK,
- + SWORK, RWORK, ITER, INFO )
+ $ SWORK, RWORK, ITER, INFO )
*
-* -- LAPACK PROTOTYPE driver routine (version 3.2.2) --
+* -- 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..--
-* January 2007
+* November 2011
*
-* ..
* .. Scalar Arguments ..
INTEGER INFO, ITER, LDA, LDB, LDX, N, NRHS
* ..
@@ -15,117 +214,10 @@
DOUBLE PRECISION RWORK( * )
COMPLEX SWORK( * )
COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( N, * ),
- + X( LDX, * )
+ $ X( LDX, * )
* ..
*
-* Purpose
-* =======
-*
-* ZCGESV computes the solution to a complex system of linear equations
-* A * X = B,
-* where A is an N-by-N matrix and X and B are N-by-NRHS matrices.
-*
-* ZCGESV first attempts to factorize the matrix in COMPLEX and use this
-* factorization within an iterative refinement procedure to produce a
-* solution with COMPLEX*16 normwise backward error quality (see below).
-* If the approach fails the method switches to a COMPLEX*16
-* factorization and solve.
-*
-* The iterative refinement is not going to be a winning strategy if
-* the ratio COMPLEX performance over COMPLEX*16 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) COMPLEX*16 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) COMPLEX*16 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) COMPLEX*16 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) COMPLEX*16 array, dimension (N*NRHS)
-* This array is used to hold the residual vectors.
-*
-* SWORK (workspace) COMPLEX 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.
-*
-* RWORK (workspace) DOUBLE PRECISION array, dimension (N)
-*
-* ITER (output) INTEGER
-* < 0: iterative refinement has failed, COMPLEX*16
-* 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 CGETRF
-* -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 COMPLEX*16 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
@@ -139,7 +231,7 @@
*
COMPLEX*16 NEGONE, ONE
PARAMETER ( NEGONE = ( -1.0D+00, 0.0D+00 ),
- + ONE = ( 1.0D+00, 0.0D+00 ) )
+ $ ONE = ( 1.0D+00, 0.0D+00 ) )
*
* .. Local Scalars ..
INTEGER I, IITER, PTSA, PTSX
@@ -148,7 +240,7 @@
*
* .. External Subroutines ..
EXTERNAL CGETRS, CGETRF, CLAG2Z, XERBLA, ZAXPY, ZGEMM,
- + ZLACPY, ZLAG2C
+ $ ZLACPY, ZLAG2C
* ..
* .. External Functions ..
INTEGER IZAMAX
@@ -190,7 +282,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.
@@ -243,7 +335,7 @@
* Solve the system SA*SX = SB.
*
CALL CGETRS( 'No transpose', N, NRHS, SWORK( PTSA ), N, IPIV,
- + SWORK( PTSX ), N, INFO )
+ $ SWORK( PTSX ), N, INFO )
*
* Convert SX back to double precision
*
@@ -254,7 +346,7 @@
CALL ZLACPY( 'All', N, NRHS, B, LDB, WORK, N )
*
CALL ZGEMM( '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.
@@ -263,7 +355,7 @@
XNRM = CABS1( X( IZAMAX( N, X( 1, I ), 1 ), I ) )
RNRM = CABS1( WORK( IZAMAX( 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
@@ -289,7 +381,7 @@
* Solve the system SA*SX = SR.
*
CALL CGETRS( '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.
@@ -305,7 +397,7 @@
CALL ZLACPY( 'All', N, NRHS, B, LDB, WORK, N )
*
CALL ZGEMM( '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.
@@ -314,7 +406,7 @@
XNRM = CABS1( X( IZAMAX( N, X( 1, I ), 1 ), I ) )
RNRM = CABS1( WORK( IZAMAX( 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
@@ -343,11 +435,11 @@
CALL ZGETRF( N, N, A, LDA, IPIV, INFO )
*
IF( INFO.NE.0 )
- + RETURN
+ $ RETURN
*
CALL ZLACPY( 'All', N, NRHS, B, LDB, X, LDX )
CALL ZGETRS( 'No transpose', N, NRHS, A, LDA, IPIV, X, LDX,
- + INFO )
+ $ INFO )
*
RETURN
*