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