C
C $Id: lsolve.F,v 1.2 1998/07/16 16:40:04 jjv5 Exp arjan $
C
C------------------------------------------------------------------------
      SUBROUTINE LSOLVE(N,A,B,W,THRESH,X,IERROR)
C
C     USES GAUSSIAN ELIMINATION WITH ROW PIVOTING TO SOLVE THE ORDER N
C     LINEAR SYSTEM:
C
C                       A*X = B.
C
C     W IS A WORK VECTOR OF LENGTH N, AND THRESH IS A THRESHOLD FOR
C     ZERO PIVOTAL ELEMENTS OF A.
C
C     ERROR CODES:  IERROR = 0 - SOLUTION FOUND SUCCESSFULLY;
C                   IERROR = 1 - ZERO PIVOT ENCOUNTERED, SINGULAR MATRIX.
C
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(N,N),B(N),W(N),X(N)
      IERROR = 0
      IF(THRESH.LE.0.0D0) THRESH = 1.0D-12
      IF(N.EQ.1)THEN
        IF(ABS(A(1,1)).LT.THRESH)THEN
          IERROR = 1
        ELSE
          X(1) = B(1)/A(1,1)
        ENDIF
        RETURN
      ENDIF
C
C     COMPUTE NORM OF A AS THE AVERAGE ABSOLUTE VALUE.
C
      AVE = 0.0D0
      DO 20 I=1,N
        DO 10 J=1,N
          AVE = AVE + ABS(A(I,J))
 10     CONTINUE
 20   CONTINUE
      AVE = AVE/(N*N)
      AMIN = AVE*THRESH
C
C     BEGIN GAUSSIAN ELIMINATION.
C
      DO 100 K=1,N-1
C
C       CHECK ENTRIES K THROUGH N OF THE KTH COLUMN OF A TO FIND THE
C       LARGEST VALUE.
C
        AIKMAX = 0.0D0
        IROW = K
        DO 30 I=K,N
          AIK = ABS(A(I,K))
          IF(AIK.GT.AIKMAX)THEN
            AIKMAX = AIK
            IROW = I
          ENDIF
 30     CONTINUE
C
C       AIKMAX IS THE ABSOLUTE VALUE OF THE PIVOTAL ELEMENT.  IF
C       AIKMAX IS SMALLER THAN THE THRESHOLD AMIN THEN THE LEADING
C       SUBMATRIX IS NEARLY SINGULAR.  IN THIS CASE, RETURN TO CALLING
C       ROUTINE WITH AN ERROR.
C
        IF(AIKMAX.LT.AMIN)THEN
          IERROR = 1
          RETURN
        ENDIF
C
C       IF IROW IS NOT EQUAL TO K THEN SWAP ROWS K AND IROW OF THE
C       MATRIX A AND SWAP ENTRIES K AND IROW OF THE VECTOR B.
C
        IF(IROW.NE.K)THEN
          DO 60 J=1,N
            W(J) = A(IROW,J)
            A(IROW,J) = A(K,J)
            A(K,J) = W(J)
 60       CONTINUE
          BSWAP = B(IROW)
          B(IROW) = B(K)
          B(K) = BSWAP
        ELSE
          DO 70 J=K+1,N
            W(J) = A(K,J)
 70       CONTINUE
        ENDIF
C
C-RDC C       FOR J GREATER THAN OR EQUAL TO I, OVERWRITE A(I,J) WITH
C-RDC C       U(I,J); FOR I GREATER THAN J OVERWRITE A(I,J) WITH L(I,J).
C
         DO 90 I=K+1,N
           T = A(I,K)/A(K,K)
           A(I,K) = T
           DO 80 J=K+1,N
             A(I,J) = A(I,J) - T*W(J)
 80       CONTINUE
 90     CONTINUE
 100  CONTINUE
       IF(ABS(A(N,N)).LT.AMIN)THEN
         IERROR = 1
         RETURN
       ENDIF
C
C     WE NOW HAVE STORED IN A THE L-U DECOMPOSITION OF P*A, WHERE
C     P IS A PERMUTATION MATRIX OF THE N BY N IDENTITY MATRIX.  IN
C     THE VECTOR B WE HAVE P*B.  WE NOW SOLVE L*U*X = P*B FOR THE
C-RDC C     VECTOR X.  FIRST OVERWRITE B WITH THE SOLUTION TO L*Y = B
C     VIA FORWARD ELIMINATION.
C
      DO 160 I=2,N
        DO 150 J=1,I-1
          B(I) = B(I) - A(I,J)*B(J)
 150    CONTINUE
 160  CONTINUE
C
C     NOW SOLVE U*X = B FOR X VIA BACK SUBSTITUTION.
C
      X(N) = B(N)/A(N,N)
      DO 200 K=2,N
        I = N+1-K
        X(I) = B(I)
        DO 190 J=I+1,N
          X(I) = X(I) - A(I,J)*X(J)
 190    CONTINUE
        X(I) = X(I)/A(I,I)
 200  CONTINUE
      RETURN
      END