      SUBROUTINE ORIMP(M, N, A, KA, QR, MS, S, IP, B, X, R, EXIT)
C ----------------------------------------------------------------------
C      PURPOSE
C        GIVEN AN APPROXIMATE LEAST SQUARES SOLUTION X OF A LINEAR
C        SYSTEM AX = B OBTAINED USING ORSOL. ORIMP ATTEMPTS TO COMPUTE
C        AN IMPROVED SOLUTION CORRECT TO MACHINE PRECISION.
C
C      CONTROL
C
C        DIMENSION A(KA,N), QR(MS,N), B(M), X(M), R(M), S(N), IP(N)
C        LOGICAL EXIT
C        .
C        .
C        .
C        CALL ORIMP(M, N, A, KA, QR, MS, S, IP, B, X, R, EXIT)
C
C        WHERE
C        M    IS AN INTEGER INPUT VARIABLE, THE NUMBER OF ROWS OF A.
C        N    IS AN INTEGER INPUT VARIABLE, THE NUMBER OF COLUMNS OF A
C             (1 .LT. N .LE. M).
C        A    IS A REAL INPUT ARRAY, THE GIVEN M BY N MATRIX.
C        QR   IS A REAL INPUT ARRAY, THE ORTHOGONAL AND TRIANGULAR
C             FACTORS OF A PRODUCED BY ORTHO.
C        B    IS A REAL INPUT ARRAY, THE RIGHT HAND SIDE OF AX = B.
C        X    AS A REAL INPUT ARRAY IS THE APPROXIMATE LEAST SQUARES
C             SOLUTION TOGETHER WITH THE RESIDUAL INFORMATION PRODUCED
C             BY ORSOL.
C        X    AS A REAL OUTPUT ARRAY IS THE IMPROVED LEAST SQUARES
C             SOLUTION WITH A RESIDUAL OF MINIMUM LENGTH.
C        R    IS A REAL OUTPUT ARRAY, THE CORRECTION VECTOR ADDED TO X.
C        KA   IS AN INTEGER INPUT VARIABLE, THE LEADING DIMENSION OF
C             OF A IN THE CALLING PROGRAM.
C        MS   IS AN INTEGER INPUT VARIABLE, THE LEADING DIMENSION OF
C             OF QR IN THE CALLING PROGRAM.
C        S    IS A REAL INPUT ARRAY, A RELEVANT PART OF Q PRODUCED BY
C             ORTHO.
C        IP   IS AN INTEGER INPUT ARRAY, THE PERMUTATION INFORMATION
C             PRODUCED BY ORTHO.
C        EXIT IS SET TO THE VALUE .TRUE. IF IMPROVEMENT OF X IS SUCCESS-
C             FUL WITH A GAIN IN ACCURACY OF AT LEAST 50 PER CENT EACH
C             ITERATION AND .FALSE. OTHERWISE.
C
C      METHOD
C        ORIMP EXECUTES THE ITERATION CYCLE
C              (1)  AR = B - AX
C              (2)  X = X + R
C        WITH A GIVEN INITIAL X.  THE RESIDUAL VECTOR B - AX IS COMPUTED
C        TO HIGH ACCURACY BY DOUBLE PRECISION.  ORSOL IS THEN USED TO
C        SOLVE (1).
C
C ----------------------------------------------------------------------
C
      DIMENSION A(KA,N), QR(MS,N), B(M), X(M), R(M), S(N), IP(N)
      LOGICAL EXIT
      DOUBLE PRECISION DSUM
      DATA ZERO/0.0/, ONE/1.0/, FOUR/4.0/, FOURTH/0.25/
C
C     ********** EPS IS A MACHINE DEPENDENT PARAMETER. ASSIGN EPS
C                THE VALUE U WHERE U IS THE SMALLEST POSITIVE FLOATING
C                POINT NUMBER SUCH THAT 1.0 + U .GT. 1.0.
C
      EPS = SPMPAR(1)
C
      EXIT = .TRUE.
      NP1 = N + 1
      EPS2 = EPS*EPS
C
      XNRM2 = ZERO
      DO 10 I = 1,N
   10 XNRM2 = XNRM2 + X(I)*X(I)
      IF (XNRM2 .EQ. ZERO) RETURN
      RATIO = ONE
C
C                                       FIND THE RESIDUAL VECTOR.
C
   20 DO 22 K = 1,M
      DSUM = B(K)
      DO 21 J = 1,N
   21 DSUM = DSUM - DBLE(A(K,J))*DBLE(X(J))
   22 R(K) = DSUM
C
C                                       FIND THE CORRECTION VECTOR.
C
      CALL ORSOL(M, N, QR, MS, S, IP, R)
      RNRM2 = ZERO
      DO 30 K = 1,N
   30 RNRM2 = RNRM2 + R(K)*R(K)
      IF (RNRM2 .LE. EPS2*XNRM2) RETURN
C
C                                       FORM NEW APPROXIMATE SOLUTION.
C
      DO 40 K = 1,N
   40 X(K) = X(K) + R(K)
      XNRM2 = ZERO
      DO 41 K = 1,N
   41 XNRM2 = XNRM2 + X(K)*X(K)
      IF (M .EQ. N) GO TO 50
      DO 42 K = NP1,M
   42 X(K) = R(K)
C
   50 IF (XNRM2 .EQ. ZERO) RETURN
      RAT = RATIO
      RATIO = RNRM2/XNRM2
      IF (RATIO .LE. FOURTH*RAT) GO TO 20
C
      IF (RATIO .LE. AMIN1(RAT,FOUR*EPS2)) RETURN
      EXIT = .FALSE.
      RETURN
      END
