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
 *