SUBROUTINE DSGESV(N,NRHS,A,LDA,IPIV,X,LDX,WORK,SWORK,ITER,INFO)
*
* -- LAPACK driver routine (version @(version)) --* @(who)
* @(date)
*
* .. Scalar Arguments ..
**INTEGER** INFO,ITER,LDA,LDX,N,NRHS
**DOUBLE PRECISION** BWD
* ..

* .. Array Arguments ..
**INTEGER** IPIV(*****)
** REAL** SWORK(*****)
**DOUBLE PRECISION** A(LDA,*****),WORK(*****),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 tries 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 switch to a
* DOUBLE PRECISION factorization and solve.
*
* The iterative refinement process is stopped if
* ITER < ITERMAX = 30

* or
* BWD = RNRM /
* ( XNRM * ANRM * EPS * MIN(4.0E00,SQRT(FLOAT(N))/6E00) )
* < 1
* where
* o ITER is the number of iteration in the iterative refinement
* process
* o RNRM is the 2-norm of the residual

* o XNRM is the 2-norm of the solution
* o ANRM is the Frobenius-norm of the matrix A
* o EPS is the relative machine precision returned by DLAMCH
*
* Arguments
* =========
*
* N (input) INTEGER
* The number of linear equations, i.e., the order of the

* matrix A. M >= 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 M-by-N coefficient matrix A.
* On exit, 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,M).
*
* 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).
*
* X (input/output) DOUBLE PRECISION array, dimension (LDX,NRHS)
* On entry, the N-by-NRHS matrix of right hand side matrix B.
* On exit, 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 (2*N*NRHS)
*
* SWORK (workspace) SINGLE PRECISION array, dimension (N*(N+NRHS))
*
* ITER (output) INTEGER
* < 0: iterative refinement has failed
* -1 : SPDRAT is too small, it is not worth working in

* SINGLE PRECISION
* -2 : overflow of an entry when moving from double to
* SINGLE PRECISION
* -3 : failure of SGETRF
* -4 : failure of SGETRS
* -31: stop the iterative refinement after the 30th
* iterations
* > 0: iterative refinement has been used successfully.

* Returns the number of iteration
*
*
* 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 ..
**DOUBLE PRECISION** MONE,ONE
**PARAMETER** (MONE**=-**1.0D0,ONE**=**1.0D0)

*
* .. Local Scalars ..
**INTEGER** I,IITER,ITERMAX,DPTB,DPTR,OK,SPTSA,SPTSX
**DOUBLE PRECISION** ANRM,BWDCTE,EPS,RNRM,SPDRAT,XNRM
*
* .. External Subroutines ..
**EXTERNAL** DAXPY,DGEMM,DLACPY,DSLACONVERTGED2S,DSLACONVERTGES2D,
+ SGETRF,SGETRS,XERBLA

* ..
* .. External Functions ..
**INTEGER** ILAENV
**DOUBLE PRECISION** DLAMCH,DLANGE,DNRM2
**EXTERNAL** DLAMCH,DLANGE,DNRM2,ILAENV
* ..
* .. Intrinsic Functions ..

**INTRINSIC** MAX,SQRT
* ..
* .. Executable Statements ..
*
* Test the input parameters.
*

ITERMAX **=** 30
SPDRAT **=** 100.00E+00
*
BWD **=** 10.00E00

ITER **=** 0
INFO **=** 0
*
**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** (LDX**.LT.**MAX(1,N)) **THEN**

INFO **=** **-**7
**END IF**
**IF** (INFO**.NE.**0) **THEN**

CALL XERBLA('DSGESV ',**-**INFO)
**RETURN**
**END IF**
*
* Quick return if N .EQ. 0

*
**IF** (N**.EQ.**0) **RETURN**
*
* The iterative refinement is not going to be a winning strategy if
* SPDRAT, the speed ratio low_prec_perf/high_prec_perf is too small.

* (will need to take the number of right-hand sides into account as
* well)
*
**IF** (SPDRAT**.LT.**1.5) **THEN**
ITER **=** **-**1

**GO TO** 160
**END IF**
*
ANRM **=** DLANGE('Frobenius',N,N,A,LDA,WORK)
EPS **=** DLAMCH('P')
BWDCTE **=** ANRM*****EPS*****MIN(4.0E00,SQRT(FLOAT(N))**/**6E00)

*
SPTSA **=** 1
SPTSX **=** SPTSA **+** N*****N

*
DPTB **=** 1
DPTR **=** DPTB **+** N*****NRHS

*
CALL DLACPY('All',N,NRHS,X,LDX,WORK(DPTB),N)
*
CALL DSLACONVERTGED2S(N,N,A,LDA,SWORK(SPTSA),N,INFO)
*
**IF** (INFO**.NE.**0) **THEN**

ITER **=** **-**2
**GO TO** 160
**END IF**
*

* Compute the LU factorization of SA.
*
CALL SGETRF(N,N,SWORK(SPTSA),N,IPIV,INFO)
*
**IF** (INFO**.NE.**0) **THEN**

ITER **=** **-**3
**GO TO** 160
**END IF**
*

* Solve the system A*X = B, overwriting B with X.
*
CALL DSLACONVERTGED2S(N,NRHS,WORK(DPTB),N,SWORK(SPTSX),N,INFO)
*
**IF** (INFO**.NE.**0) **THEN**

ITER **=** **-**2
**GO TO** 160
**END IF**
*

CALL SGETRS('No transpose',N,NRHS,SWORK(SPTSA),N,IPIV,
+ SWORK(SPTSX),N,INFO)
*
**IF** (INFO**.NE.**0) **THEN**

ITER **=** **-**4
**GO TO** 160
**END IF**
*

CALL DSLACONVERTGES2D(N,NRHS,SWORK(SPTSX),N,X,LDX,INFO)
*
CALL DLACPY('All',N,NRHS,WORK(DPTB),N,WORK(DPTR),N)
*
CALL DGEMM('No Transpose','No Transpose',N,NRHS,N,MONE,A,LDA,X,
+ LDX,ONE,WORK(DPTR),N)

*
OK **=** 1
**DO** I **=** 1,NRHS
XNRM **=** DNRM2(N,X(1,I),1)
RNRM **=** DNRM2(N,WORK(DPTR**+** (I**-**1)*****N),1)
BWD **=** RNRM**/**XNRM**/**BWDCTE
**IF** (BWD**.GT.**1) OK **=** 0

**END DO**
*
ITER **=** 1
*
**DO** 90 IITER **=** 1,ITERMAX

*
**IF** (OK**.EQ.**1) **GO TO** 160
*
CALL DSLACONVERTGED2S(N,NRHS,WORK(DPTR),N,SWORK(SPTSX),N,INFO)

*
**IF** (INFO**.NE.**0) **THEN**
ITER **=** **-**2

**GO TO** 160
**END IF**
*
CALL SGETRS('No transpose',N,NRHS,SWORK(SPTSA),N,IPIV,
+ SWORK(SPTSX),N,INFO)

*
**IF** (INFO**.NE.**0) **THEN**
ITER **=** **-**4

**GO TO** 160
**END IF**
*
CALL DSLACONVERTGES2D(N,NRHS,SWORK(SPTSX),N,WORK(DPTR),N,INFO)
*
CALL DAXPY(N*****NRHS,ONE,WORK(DPTR),1,X,1)

*
CALL DLACPY('All',N,NRHS,WORK(DPTB),N,WORK(DPTR),N)
*
CALL DGEMM('No Transpose','No Transpose',N,NRHS,N,MONE,A,LDA,
+ X,LDX,ONE,WORK(DPTR),N)

*
*
OK **=** 1
**DO** I **=** 1,NRHS
XNRM **=** DNRM2(N,X(1,I),1)
RNRM **=** DNRM2(N,WORK(DPTR**+** (I**-**1)*****N),1)
BWD **=** RNRM**/**XNRM**/**BWDCTE
**IF** (BWD**.GT.**1) OK **=** 0

**END DO**
*
ITER **=** ITER **+** 1
*
90 **CONTINUE**

160 **CONTINUE**
*
**IF** (BWD**.LT.**1) **RETURN**
*

BWD **=** 0.0E0
**IF** (ITER**.GE.**0) **THEN**
ITER **=** **-**ITERMAX **-** 1

**END IF**
*
CALL DGETRF(N,N,A,LDA,IPIV,INFO)
*
CALL DLACPY('All',N,NRHS,WORK(DPTB),N,X,LDX)

*
**IF** (INFO**.EQ.**0) **THEN**
*
* Solve the system A*X = B, overwriting B with X.
*

CALL DGETRS('No transpose',N,NRHS,A,LDA,IPIV,X,LDX,INFO)
**END IF**
**RETURN**
*
* End of DGESV
*

END