      SUBROUTINE LNSYS(A,D,M,N,B,X,OPTION,IERR)
C
C   SOLVE AX = B WHERE A IS A MATRIX OF ORDER N. M IS THE NUMBER OF
C   ROWS IN THE DIMENSION STATEMENT FOR A IN THE CALLING PROGRAM.
C
C   OPTION=1   COMPUTE AN LU DECOMPOSITION OF A AND STORE IT IN D.
C              STORE THE PIVOT INDICES IN PIVOT AND SOLVE AX = B.
C   OPTION=2   COMPUTE AN LU DECOMPOSITION OF A AND STORE IT IN D.
C              STORE THE PIVOT INDICES IN PIVOT AND SOLVE AX = B.
C              THEN COMPUTE THE RESIDUAL AND ONE CORRECTION. THE
C              CORRECTION IS STORED IN R, THE NEW VALUE X1 IN X,
C              THE RELATIVE ERROR
C                         NORM(X0-X1)/NORM(X1)
C              IN THE VARIABLE ERROR, AND THE RELATIVE RESIDUAL
C                         NORM(RESIDUAL)/NORM(B)
C              IN THE VARIABLE RELRSD. THESE VALUES CAN BE OBTAINED
C              USING THE COMMON/XXLIN/ GIVEN BELOW.
C   OPTION=3   SAME AS OPTION=1, EXCEPT THAT THE LU DECOMPOSITION
C              HAS ALREADY BEEN STORED IN D AND THE PIVOT INDICES
C              IN PIVOT.
C   OPTION=4   SAME AS OPTION=2, EXCEPT THAT THE LU DECOMPOSITION
C              HAS ALREADY BEEN STORED IN D AND THE PIVOT INDICES
C              IN PIVOT.
C
C   THE LU DECOMPOSITION IS OBTAINED USING SCALED PARTIAL PIVOTING.
C   FOR OPTIONS 1 AND 2, IERR IS A VARIABLE THAT REPORTS THE STATUS
C   OF THE RESULTS. IERR = 0 IF THE LU DECOMPOSITION IS OBTAINED.
C   OTHERWISE, IERR = -K WHEN THE K-TH ROW OF A CONTAINS ONLY ZEROS
C   OR IERR = K WHEN THE K-TH PIVOT ELEMENT IS 0.
C
C   IT IS ASSUMED THAT N .LE. 128. THIS ASSUMPTION MAY BE MODIFIED
C   BY CHANGING THE DIMENSION STATEMENTS FOR THE ARRAYS PIVOT, R,
C   AND SCALE. ALSO MODIFY THE DIMENSION STATEMENT FOR PIVOT IN THE
C   SUBROUTINE IEGS.
C
      REAL A(M,N),D(M,N),B(N),X(N)
      INTEGER OPTION,PIVOT(128)
      REAL NORMX,NORME,NORMB,NORMR,R(128),SCALE(128)
      COMMON /XXLIN/ ERROR,RELRSD,PIVOT
C
      NM1 = N - 1
      ISWIT = 1
      IF (OPTION .GT. 2) GO TO 100
C
      DO 11 I = 1,N
      SCALE(I) = 0.0
        DO 10 J = 1,N
        D(I,J) = A(I,J)
   10   SCALE(I) = SCALE(I) + ABS(D(I,J))
      IF (SCALE(I) .NE. 0.0) GO TO 11
      IERR = -I
      RETURN
   11 CONTINUE
C
C             OBTAIN THE LU DECOMPOSITION OF A
C
      IERR = 0
      DO 43 K = 1,NM1
      C = ABS(D(K,K))/SCALE(K)
      L = K
      KP1 = K + 1
        DO 20 I = KP1,N
        T = ABS(D(I,K))/SCALE(I)
        IF (T .LE. C) GO TO 20
        C = T
        L = I
   20   CONTINUE
C
      IF (C .NE. 0.0) GO TO 30
      IERR = K
      RETURN
C
C     INTERCHANGE ROWS K AND L
C
   30 PIVOT(K) = L
      IF (K .EQ. L) GO TO 40
        DO 31 J = K,N
        T = D(K,J)
        D(K,J) = D(L,J)
   31   D(L,J) = T
      T = SCALE(K)
      SCALE(K) = SCALE(L)
      SCALE(L) = T
C
C     ELIMINATE THE K-TH UNKNOWN BELOW THE DIAGONAL
C
   40   DO 42 I = KP1,N
        D(I,K) = D(I,K)/D(K,K)
        T = D(I,K)
          DO 41 J = KP1,N
   41     D(I,J) = D(I,J) - T*D(K,J)
   42   CONTINUE
   43 CONTINUE
C
      IF (D(N,N) .NE. 0.0) GO TO 100
      IERR = N
      RETURN
C
C               STORE B IN R AND SET X = 0
C
  100 DO 110 I = 1,N
      R(I) = B(I)
  110 X(I) = 0.0
      GO TO 200
C
C             COMPUTE THE RESIDUAL R = B - AX
C
  120 DO 131 I = 1,N
      SUM = 0.0
        DO 130 J = 1,N
  130   SUM = SUM + A(I,J)*X(J)
  131 R(I) = B(I) - SUM
C
      NORMB = 0.0
      NORMR = 0.0
      DO 140 I = 1,N
      NORMB = AMAX1(NORMB,ABS(B(I)))
  140 NORMR = AMAX1(NORMR,ABS(R(I)))
      RELRSD = 0.0
      IF (NORMB .NE. 0.0) RELRSD = NORMR/NORMB
      ISWIT = 2
C
C             SOLVE LZ = R AND STORE Z IN R
C
  200 DO 212 K = 1,NM1
      L = PIVOT(K)
      IF (K .EQ. L) GO TO 210
      T = R(K)
      R(K) = R(L)
      R(L) = T
  210 KP1 = K + 1
        DO 211 I = KP1,N
  211   R(I) = R(I) - D(I,K)*R(K)
  212 CONTINUE
C
C      SOLVE UE = R, STORE E IN R, AND SET X = X + E
C
      R(N) = R(N)/D(N,N)
      X(N) = X(N) + R(N)
      DO 221 NMI = 1,NM1
      I = N - NMI
      IP1 = I + 1
      SUM = 0.0
        DO 220 J = IP1,N
  220   SUM = SUM + D(I,J)*R(J)
      R(I) = (R(I) - SUM)/D(I,I)
  221 X(I) = X(I) + R(I)
C
      GO TO (300,230,300,230),OPTION
  230 IF (ISWIT .EQ. 1) GO TO 120
C
C             CALCULATE THE CORRECTION ERROR
C
      NORMX = 0.0
      NORME = 0.0
      DO 250 I = 1,N
      NORMX = AMAX1(NORMX,ABS(X(I)))
  250 NORME = AMAX1(NORME,ABS(R(I)))
      ERROR = 0.0
      IF (NORMX .NE. 0.0) ERROR = NORME/NORMX
  300 RETURN
      END
