diff options
Diffstat (limited to 'SRC/zcgesv.f')
| -rw-r--r-- | SRC/zcgesv.f | 280 |
1 files changed, 145 insertions, 135 deletions
diff --git a/SRC/zcgesv.f b/SRC/zcgesv.f index 5a6d5c28..a6e6aa8e 100644 --- a/SRC/zcgesv.f +++ b/SRC/zcgesv.f @@ -1,49 +1,46 @@ SUBROUTINE ZCGESV( N, NRHS, A, LDA, IPIV, B, LDB, X, LDX, WORK, - + SWORK, ITER, INFO) + + SWORK, RWORK, ITER, INFO ) * -* -- LAPACK PROTOTYPE driver routine (version 3.1.1) -- +* -- LAPACK PROTOTYPE driver routine (version 3.2) -- * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. * January 2007 * * .. -* .. WARNING: PROTOTYPE .. -* This is an LAPACK PROTOTYPE routine which means that the -* interface of this routine is likely to be changed in the future -* based on community feedback. -* -* .. * .. Scalar Arguments .. - INTEGER INFO,ITER,LDA,LDB,LDX,N,NRHS + INTEGER INFO, ITER, LDA, LDB, LDX, N, NRHS * .. * .. Array Arguments .. - INTEGER IPIV(*) - COMPLEX SWORK(*) - COMPLEX*16 A(LDA,*),B(LDB,*),WORK(N,*),X(LDX,*) + INTEGER IPIV( * ) + DOUBLE PRECISION RWORK( * ) + COMPLEX SWORK( * ) + COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( N, * ), + + X( LDX, * ) * .. * * Purpose * ======= * -* ZCGESV computes the solution to a real system of linear equations +* 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 SINGLE COMPLEX PRECISION -* and use this factorization within an iterative refinement procedure to -* produce a solution with DOUBLE COMPLEX PRECISION normwise backward error -* quality (see below). If the approach fails the method switches to a -* DOUBLE COMPLEX PRECISION factorization and solve. +* 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 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 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 +* RNRM < SQRT(N)*XNRM*ANRM*EPS*BWDMAX * where * o ITER is the number of the current iteration in the iterative * refinement process @@ -51,7 +48,8 @@ * 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. +* The value ITERMAX and BWDMAX are fixed to 30 and 1.0D+00 +* respectively. * * Arguments * ========= @@ -80,12 +78,12 @@ * 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 +* 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 matrix of right hand side matrix B. +* The N-by-NRHS right hand side matrix B. * * LDB (input) INTEGER * The leading dimension of the array B. LDB >= max(1,N). @@ -100,17 +98,19 @@ * 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 +* 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, double precision +* < 0: iterative refinement has failed, COMPLEX*16 * factorization has been performed -* -1 : taking into account machine parameters, N, NRHS, it -* is a priori not worth working in SINGLE PRECISION -* -2 : overflow of an entry when moving from double to -* SINGLE PRECISION -* -3 : failure of SGETRF +* -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. @@ -119,34 +119,43 @@ * 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 +* > 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 .. - COMPLEX*16 NEGONE,ONE - PARAMETER (NEGONE=(-1.0D+00,0.0D+00),ONE=(1.0D+00,0.0D+00)) + LOGICAL DOITREF + PARAMETER ( DOITREF = .TRUE. ) +* + INTEGER ITERMAX + PARAMETER ( ITERMAX = 30 ) +* + DOUBLE PRECISION BWDMAX + PARAMETER ( BWDMAX = 1.0E+00 ) +* + COMPLEX*16 NEGONE, ONE + PARAMETER ( NEGONE = ( -1.0D+00, 0.0D+00 ), + + ONE = ( 1.0D+00, 0.0D+00 ) ) * * .. Local Scalars .. - LOGICAL DOITREF - INTEGER I,IITER,ITERMAX,OK,PTSA,PTSX - DOUBLE PRECISION ANRM,BWDMAX,CTE,EPS,RNRM,XNRM - COMPLEX*16 ZDUM + INTEGER I, IITER, PTSA, PTSX + DOUBLE PRECISION ANRM, CTE, EPS, RNRM, XNRM + COMPLEX*16 ZDUM * * .. External Subroutines .. - EXTERNAL CGETRS,CGETRF,CLAG2Z,XERBLA,ZAXPY, - $ ZGEMM,ZLACPY,ZLAG2C + EXTERNAL CGETRS, CGETRF, CLAG2Z, XERBLA, ZAXPY, ZGEMM, + + ZLACPY, ZLAG2C * .. * .. External Functions .. - INTEGER IZAMAX - DOUBLE PRECISION DLAMCH,ZLANGE - EXTERNAL IZAMAX,DLAMCH,ZLANGE + INTEGER IZAMAX + DOUBLE PRECISION DLAMCH, ZLANGE + EXTERNAL IZAMAX, DLAMCH, ZLANGE * .. * .. Intrinsic Functions .. - INTRINSIC ABS,DBLE,MAX,SQRT + INTRINSIC ABS, DBLE, MAX, SQRT * .. * .. Statement Functions .. DOUBLE PRECISION CABS1 @@ -156,51 +165,47 @@ * .. * .. Executable Statements .. * - ITERMAX = 30 - BWDMAX = 1.0E+00 - DOITREF = .TRUE. -* - OK = 0 INFO = 0 ITER = 0 * * Test the input parameters. * - IF (N.LT.0) THEN - INFO = -1 - ELSE IF (NRHS.LT.0) THEN - INFO = -2 - ELSE IF (LDA.LT.MAX(1,N)) THEN - INFO = -4 - ELSE IF (LDB.LT.MAX(1,N)) THEN - INFO = -7 - ELSE IF (LDX.LT.MAX(1,N)) THEN - INFO = -9 + IF( N.LT.0 ) THEN + INFO = -1 + ELSE IF( NRHS.LT.0 ) THEN + INFO = -2 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -4 + ELSE IF( LDB.LT.MAX( 1, N ) ) THEN + INFO = -7 + ELSE IF( LDX.LT.MAX( 1, N ) ) THEN + INFO = -9 END IF - IF (INFO.NE.0) THEN - CALL XERBLA('ZCGESV',-INFO) - RETURN + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'ZCGESV', -INFO ) + RETURN END IF * * Quick return if (N.EQ.0). * - IF (N.EQ.0) RETURN + IF( N.EQ.0 ) + + RETURN * * Skip single precision iterative refinement if a priori slower * than double precision factorization. * - IF (.NOT.DOITREF) THEN - ITER = -1 - GO TO 40 + IF( .NOT.DOITREF ) THEN + ITER = -1 + GO TO 40 END IF * * Compute some constants. * - ANRM = ZLANGE('I',N,N,A,LDA,WORK) - EPS = DLAMCH('Epsilon') - CTE = ANRM*EPS*SQRT(DBLE(N))*BWDMAX + ANRM = ZLANGE( 'I', N, N, A, LDA, RWORK ) + EPS = DLAMCH( 'Epsilon' ) + CTE = ANRM*EPS*SQRT( DBLE( N ) )*BWDMAX * -* Set the pointers PTSA, PTSX for referencing SA and SX in SWORK. +* Set the indices PTSA, PTSX for referencing SA and SX in SWORK. * PTSA = 1 PTSX = PTSA + N*N @@ -208,58 +213,59 @@ * Convert B from double precision to single precision and store the * result in SX. * - CALL ZLAG2C(N,NRHS,B,LDB,SWORK(PTSX),N,INFO) + CALL ZLAG2C( N, NRHS, B, LDB, SWORK( PTSX ), N, INFO ) * - IF (INFO.NE.0) THEN - ITER = -2 - GO TO 40 + IF( INFO.NE.0 ) THEN + ITER = -2 + GO TO 40 END IF * * Convert A from double precision to single precision and store the * result in SA. * - CALL ZLAG2C(N,N,A,LDA,SWORK(PTSA),N,INFO) + CALL ZLAG2C( N, N, A, LDA, SWORK( PTSA ), N, INFO ) * - IF (INFO.NE.0) THEN - ITER = -2 - GO TO 40 + IF( INFO.NE.0 ) THEN + ITER = -2 + GO TO 40 END IF * * Compute the LU factorization of SA. * - CALL CGETRF(N,N,SWORK(PTSA),N,IPIV,INFO) + CALL CGETRF( N, N, SWORK( PTSA ), N, IPIV, INFO ) * - IF (INFO.NE.0) THEN - ITER = -3 - GO TO 40 + IF( INFO.NE.0 ) THEN + ITER = -3 + GO TO 40 END IF * * Solve the system SA*SX = SB. * - CALL CGETRS('No transpose',N,NRHS,SWORK(PTSA),N,IPIV, - + SWORK(PTSX),N,INFO) + CALL CGETRS( 'No transpose', N, NRHS, SWORK( PTSA ), N, IPIV, + + SWORK( PTSX ), N, INFO ) * * Convert SX back to double precision * - CALL CLAG2Z(N,NRHS,SWORK(PTSX),N,X,LDX,INFO) + CALL CLAG2Z( N, NRHS, SWORK( PTSX ), N, X, LDX, INFO ) * * Compute R = B - AX (R is WORK). * - CALL ZLACPY('All',N,NRHS,B,LDB,WORK,N) + 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) + CALL ZGEMM( 'No Transpose', 'No Transpose', N, NRHS, N, NEGONE, A, + + LDA, X, LDX, ONE, WORK, N ) * -* Check whether the NRHS normwised backward errors satisfy the +* Check whether the NRHS normwise backward errors satisfy the * stopping criterion. If yes, set ITER=0 and return. * - DO I = 1,NRHS - 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) GOTO 10 + DO I = 1, NRHS + 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 END DO * -* If we are here, the NRHS normwised backward errors satisfy the +* If we are here, the NRHS normwise backward errors satisfy the * stopping criterion. We are good to exit. * ITER = 0 @@ -267,54 +273,57 @@ * 10 CONTINUE * - DO 30 IITER = 1,ITERMAX + DO 30 IITER = 1, ITERMAX * -* Convert R (in WORK) from double precision to single precision -* and store the result in SX. +* Convert R (in WORK) from double precision to single precision +* and store the result in SX. * - CALL ZLAG2C(N,NRHS,WORK,N,SWORK(PTSX),N,INFO) + CALL ZLAG2C( N, NRHS, WORK, N, SWORK( PTSX ), N, INFO ) * - IF (INFO.NE.0) THEN - ITER = -2 - GO TO 40 - END IF + IF( INFO.NE.0 ) THEN + ITER = -2 + GO TO 40 + END IF * -* Solve the system SA*SX = SR. +* Solve the system SA*SX = SR. * - CALL CGETRS('No transpose',N,NRHS,SWORK(PTSA),N,IPIV, - + SWORK(PTSX),N,INFO) + CALL CGETRS( 'No transpose', N, NRHS, SWORK( PTSA ), N, IPIV, + + SWORK( PTSX ), N, INFO ) * -* Convert SX back to double precision and update the current -* iterate. +* Convert SX back to double precision and update the current +* iterate. * - CALL CLAG2Z(N,NRHS,SWORK(PTSX),N,WORK,N,INFO) + CALL CLAG2Z( N, NRHS, SWORK( PTSX ), N, WORK, N, INFO ) * - CALL ZAXPY(N*NRHS,ONE,WORK,1,X,1) + DO I = 1, NRHS + CALL ZAXPY( N, ONE, WORK( 1, I ), 1, X( 1, I ), 1 ) + END DO * -* Compute R = B - AX (R is WORK). +* Compute R = B - AX (R is WORK). * - CALL ZLACPY('All',N,NRHS,B,LDB,WORK,N) + 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) + CALL ZGEMM( 'No Transpose', 'No Transpose', N, NRHS, N, NEGONE, + + A, LDA, X, LDX, ONE, WORK, N ) * -* Check whether the NRHS normwised backward errors satisfy the -* stopping criterion. If yes, set ITER=IITER>0 and return. +* Check whether the NRHS normwise backward errors satisfy the +* stopping criterion. If yes, set ITER=IITER>0 and return. * - DO I = 1,NRHS - 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) GOTO 20 - END DO + DO I = 1, NRHS + 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 + END DO * -* If we are here, the NRHS normwised backward errors satisfy the -* stopping criterion, we are good to exit. +* If we are here, the NRHS normwise backward errors satisfy the +* stopping criterion, we are good to exit. * - ITER = IITER + ITER = IITER * - RETURN + RETURN * - 20 CONTINUE + 20 CONTINUE * 30 CONTINUE * @@ -330,13 +339,14 @@ * Single-precision iterative refinement failed to converge to a * satisfactory solution, so we resort to double precision. * - CALL ZGETRF(N,N,A,LDA,IPIV,INFO) + CALL ZGETRF( N, N, A, LDA, IPIV, INFO ) * - CALL ZLACPY('All',N,NRHS,B,LDB,X,LDX) + IF( INFO.NE.0 ) + + RETURN * - IF (INFO.EQ.0) THEN - CALL ZGETRS('No transpose',N,NRHS,A,LDA,IPIV,X,LDX,INFO) - END IF + CALL ZLACPY( 'All', N, NRHS, B, LDB, X, LDX ) + CALL ZGETRS( 'No transpose', N, NRHS, A, LDA, IPIV, X, LDX, + + INFO ) * RETURN * |