      SUBROUTINE LUIMP(A, KA, N, Q, KQ, IPVT, B, X, R, IND)
C ----------------------------------------------------------------------
C     PURPOSE
C       GIVEN AN APPROXIMATE SOLUTION X OF A LINEAR SYSTEM AX = B
C       OBTAINED USING SGECO OR SGEFA. LUIMP ATTEMPTS TO COMPUTE
C       AN IMPROVED SOLUTION CORRECT TO MACHINE PRECISION.
C
C     PARAMETERS
C
C       A    AN ARRAY OF DIMENSION (KA,N) CONTAINING THE MATRIX
C            A OF ORDER N.
C       Q    AN ARRAY OF DIMENSION (KQ,N) CONTAINING THE LU
C            DECOMPOSITION OF A PRODUCED BY SGECO OR SGEFA.
C       IPVT AN ARRAY OF DIMENSION N CONTAINING THE PERMUTATION
C            INFORMATION GIVEN BY SGECO OR SGEFA.
C       B    THE RIGHT HAND SIDE OF THE EQUATION AX = B.
C       X    ON INPUT X IS THE APPROXIMATE SOLUTION OF AX = B TO
C            BE IMPROVED. ON OUTPUT X IS THE SOLUTION OBTAINED.
C       R    AN ARRAY FOR INTERNAL USE BY THE ROUTINE.
C       IND  VARIABLE THAT REPORTS THE STATUS OF THE RESULTS.
C            IND = 0 IF IMPROVEMENT OF X IS SUCCESSFUL WITH A
C            GAIN IN ACCURACY OF AT LEAST 50 PER CENT EACH
C            ITERATION. OTHERWISE IND = 1.
C
C     METHOD
C       LUIMP EXECUTES THE ITERATION CYCLE
C            (1)  AR = B - AX
C            (2)  X = X + R
C       WITH AN INITIAL GIVEN X. THE RESIDUAL VECTOR B - AX IS
C       COMPUTED TO HIGH ACCURACY USING DOUBLE PRECISION. SGESL
C       IS THEN USED TO SOLVE (1).
C
C ----------------------------------------------------------------------
      DIMENSION A(KA,N), Q(KQ,N), IPVT(N), B(N), X(N), R(N)
      DOUBLE PRECISION DSUM
C
C     ********** EPS IS A MACHINE DEPENDENT PARAMETER. ASSIGN EPS
C                THE VALUE U WHERE U IS THE SMALLEST FLOATING POINT
C                NUMBER SUCH THAT 1.0 + U .GT. 1.0.
C
      EPS = SPMPAR(1)
C
      IND = 0
      XNRM = 0.0
      DO 10 I = 1,N
   10 XNRM = XNRM + X(I)*X(I)
      IF (XNRM .EQ. 0.0) RETURN
      EPS2 = EPS*EPS
      RATIO = 1.0
C
C                  COMPUTE THE RESIDUAL VECTOR
C
   20 DO 22 I = 1,N
      DSUM = B(I)
        DO 21 J = 1,N
   21   DSUM = DSUM - DBLE(A(I,J))*DBLE(X(J))
   22 R(I) = DSUM
C
C                  FIND THE CORRECTION VECTOR
C
      CALL SGESL(Q, KQ, N, IPVT, R, 0)
      RNRM = 0.0
      DO 30 I = 1,N
   30 RNRM = RNRM + R(I)*R(I)
      IF (RNRM .LE. EPS2*XNRM) RETURN
C
C                FORM A NEW APPROXIMATE SOLUTION
C
      DO 40 I = 1,N
   40 X(I) = X(I) + R(I)
      XNRM = 0.0
      DO 41 I = 1,N
   41 XNRM = XNRM + X(I)*X(I)
C
      IF (XNRM .EQ. 0.0) RETURN
      RAT = RATIO
      RATIO = RNRM/XNRM
      IF (RATIO .LE. 0.25*RAT) GO TO 20
C
      IF (RATIO .GT. AMIN1(RAT,4.0*EPS2)) IND = 1
      RETURN
      END
