1 files changed, 136 insertions, 130 deletions

diff --git a/SRC/dsgesv.f b/SRC/dsgesv.f
index 5be14625..1bf5a8df 100644
--- a/SRC/dsgesv.f
+++ b/SRC/dsgesv.f
@@ -1,24 +1,19 @@
       SUBROUTINE DSGESV( N, NRHS, A, LDA, IPIV, B, LDB, X, LDX, WORK,
-     +                   SWORK, ITER, INFO)
+     +                   SWORK, 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..
 *     February 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(*)
-      REAL SWORK(*)
-      DOUBLE PRECISION A(LDA,*),B(LDB,*),WORK(N,*),X(LDX,*)
+      INTEGER            IPIV( * )
+      REAL               SWORK( * )
+      DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), WORK( N, * ),
+     +                   X( LDX, * )
 *     ..
 *
 *  Purpose
@@ -28,22 +23,23 @@
 *     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
+*  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 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 
+*      RNRM < SQRT(N)*XNRM*ANRM*EPS*BWDMAX
 *  where
 *      o ITER is the number of the current iteration in the iterative
 *        refinement process
@@ -51,7 +47,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 +77,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) DOUBLE PRECISION 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,16 +97,16 @@
 *          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 
+*          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 : 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
+*               -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
@@ -127,73 +124,77 @@
 *  =========
 *
 *     .. Parameters ..
-      DOUBLE PRECISION NEGONE,ONE
-      PARAMETER (NEGONE=-1.0D+0,ONE=1.0D+0)
+      LOGICAL            DOITREF
+      PARAMETER          ( DOITREF = .TRUE. )
+*
+      INTEGER            ITERMAX
+      PARAMETER          ( ITERMAX = 30 )
+*
+      DOUBLE PRECISION   BWDMAX
+      PARAMETER          ( BWDMAX = 1.0E+00 )
+*
+      DOUBLE PRECISION   NEGONE, ONE
+      PARAMETER          ( NEGONE = -1.0D+0, ONE = 1.0D+0 )
 *
 *     .. Local Scalars ..
-      LOGICAL DOITREF
-      INTEGER I,IITER,ITERMAX,OK,PTSA,PTSX
-      DOUBLE PRECISION ANRM,BWDMAX,CTE,EPS,RNRM,XNRM
+      INTEGER            I, IITER, PTSA, PTSX
+      DOUBLE PRECISION   ANRM, CTE, EPS, RNRM, XNRM
 *
 *     .. External Subroutines ..
-      EXTERNAL DAXPY,DGEMM,DLACPY,DLAG2S,SLAG2D,
-     +         SGETRF,SGETRS,XERBLA
+      EXTERNAL           DAXPY, DGEMM, DLACPY, DLAG2S, SLAG2D, SGETRF,
+     +                   SGETRS, XERBLA
 *     ..
 *     .. External Functions ..
-      INTEGER IDAMAX
-      DOUBLE PRECISION DLAMCH,DLANGE
-      EXTERNAL IDAMAX,DLAMCH,DLANGE
+      INTEGER            IDAMAX
+      DOUBLE PRECISION   DLAMCH, DLANGE
+      EXTERNAL           IDAMAX, DLAMCH, DLANGE
 *     ..
 *     .. Intrinsic Functions ..
-      INTRINSIC ABS,DBLE,MAX,SQRT
+      INTRINSIC          ABS, DBLE, MAX, SQRT
 *     ..
 *     .. 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('DSGESV',-INFO)
-          RETURN
+      IF( INFO.NE.0 ) THEN
+         CALL XERBLA( 'DSGESV', -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 = DLANGE('I',N,N,A,LDA,WORK)
-      EPS = DLAMCH('Epsilon')
-      CTE = ANRM*EPS*SQRT(DBLE(N))*BWDMAX
+      ANRM = DLANGE( 'I', N, N, A, LDA, WORK )
+      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
@@ -201,120 +202,124 @@
 *     Convert B from double precision to single precision and store the
 *     result in SX.
 *
-      CALL DLAG2S(N,NRHS,B,LDB,SWORK(PTSX),N,INFO)
+      CALL DLAG2S( 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 DLAG2S(N,N,A,LDA,SWORK(PTSA),N,INFO)
+      CALL DLAG2S( 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 SGETRF(N,N,SWORK(PTSA),N,IPIV,INFO)
+      CALL SGETRF( 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 SGETRS('No transpose',N,NRHS,SWORK(PTSA),N,IPIV,
-     +            SWORK(PTSX),N,INFO)
+      CALL SGETRS( 'No transpose', N, NRHS, SWORK( PTSA ), N, IPIV,
+     +             SWORK( PTSX ), N, INFO )
 *
 *     Convert SX back to double precision
 *
-      CALL SLAG2D(N,NRHS,SWORK(PTSX),N,X,LDX,INFO)
+      CALL SLAG2D( N, NRHS, SWORK( PTSX ), N, X, LDX, INFO )
 *
 *     Compute R = B - AX (R is WORK).
 *
-      CALL DLACPY('All',N,NRHS,B,LDB,WORK,N)
+      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)
+      CALL DGEMM( '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 = ABS(X(IDAMAX(N,X(1,I),1),I))
-          RNRM = ABS(WORK(IDAMAX(N,WORK(1,I),1),I))
-          IF (RNRM.GT.XNRM*CTE) GOTO 10
+      DO I = 1, NRHS
+         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
       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
       RETURN
 *
- 10   CONTINUE
+   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 DLAG2S(N,NRHS,WORK,N,SWORK(PTSX),N,INFO)
+         CALL DLAG2S( 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 SGETRS('No transpose',N,NRHS,SWORK(PTSA),N,IPIV,
-     +                SWORK(PTSX),N,INFO)
+         CALL SGETRS( '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 SLAG2D(N,NRHS,SWORK(PTSX),N,WORK,N,INFO)
+         CALL SLAG2D( N, NRHS, SWORK( PTSX ), N, WORK, N, INFO )
 *
-          CALL DAXPY(N*NRHS,ONE,WORK,1,X,1)
+         DO I = 1, NRHS
+            CALL DAXPY( 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 DLACPY('All',N,NRHS,B,LDB,WORK,N)
+         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)
+         CALL DGEMM( '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 = ABS(X(IDAMAX(N,X(1,I),1),I))
-              RNRM = ABS(WORK(IDAMAX(N,WORK(1,I),1),I))
-              IF (RNRM.GT.XNRM*CTE) GOTO 20
-          END DO
+         DO I = 1, NRHS
+            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
+         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
 *
 *     If we are at this place of the code, this is because we have
-*     performed ITER=ITERMAX iterations and never satisified the stopping
-*     criterion, set up the ITER flag accordingly and follow up on double
-*     precision routine.
+*     performed ITER=ITERMAX iterations and never satisified the
+*     stopping criterion, set up the ITER flag accordingly and follow up
+*     on double precision routine.
 *
       ITER = -ITERMAX - 1
 *
@@ -323,13 +328,14 @@
 *     Single-precision iterative refinement failed to converge to a
 *     satisfactory solution, so we resort to double precision.
 *
-      CALL DGETRF(N,N,A,LDA,IPIV,INFO)
+      CALL DGETRF( N, N, A, LDA, IPIV, INFO )
 *
-      CALL DLACPY('All',N,NRHS,B,LDB,X,LDX)
+      IF( INFO.NE.0 )
+     +   RETURN
 *
-      IF (INFO.EQ.0) THEN
-          CALL DGETRS('No transpose',N,NRHS,A,LDA,IPIV,X,LDX,INFO)
-      END IF
+      CALL DLACPY( 'All', N, NRHS, B, LDB, X, LDX )
+      CALL DGETRS( 'No transpose', N, NRHS, A, LDA, IPIV, X, LDX,
+     +             INFO )
 *
       RETURN
 *