C
C PURPOSE: SOLVES THE EQUATION  SYSTEM  A*X=B.
C
C INPUTS:
C     N      SIZE OF SYSTEM TO BE SOLVED (N,N)
C     NDIM   DECLARED SIZE OF A (NDIM,NDIM) AND B (NDIM)
C     A      MATRIX 
C     B      INPUT VECTOR
C     NEWMAT IF .TRUE. PERFORM LU DECOMPOSITION
C OUTPUTS:
C     A      MODIFIED (LU DECOMPOSED) MATRIX, CAN BE USED IN LATER CALLS
C     B      SOLUTION VECTOR
C COMMON:
C
C COMMENTS: OCTOBER 6, 1999, P. JUDGE
C 
      SUBROUTINE EQSYST(N,NDIM,A,B,NEWMAT)
C
C  
C  WHEN NEWMAT=TRUE, THE SYSTEM IS REARRANGED INTO U*X=L*B, WHERE U
C  IS UPPER AND L IS LOWER TRIANGULAR. THESE ARE THEN REUSED IN LATER
C  CALLS WITH NEWMAT=FALSE AND NEW RIGHT HAND SIDES B. THE SOLUTION
C  VECTOR IS RETURNED IN B. NO PIVOTING, I.E. THE MATRIX A IS ASSUMED
C  TO HAVE NONZERO DIAGONAL ELEMENTS.
C
C  CODED BY: A. NORDLUND (FEB-1979)
C
C  THIS IS A MODIFIED VERSION OF EQSYST WHICH TESTS FOR ZERO ELEMENTS
C  BELOW THE DIAGONAL AND ALSO STOPS AT THE LAST NON-ZERO ELEMENT ABOVE
C  THE DIAGONAL. CONSIDERABLE SAVINGS ARE OBTAINED FOR LOOSE MATRICES.
C
C  THIS IS A COLUMN ORIENTED VERSION (M. CARLSSON JAN-1986)
C  TEMPORARY SCALARS BL, ALL, ALM AND BK ARE USED TO SHOW THE COMPILER
C  THAT THERE IS NO VECTOR DEPENDENCY IN THE INNERMOST DO-LOOP
C
      INCLUDE 'PREC'
      PARAMETER (MDIM=2000)
      DIMENSION A(NDIM,NDIM),B(NDIM)
      INTEGER LASTN(MDIM)
      LOGICAL NEWMAT
      SAVE LASTN
C
      IF(NEWMAT) THEN
C
C     FIND THE LAST NON-ZERO ELEMENT IN EACH COLUMN
C
        IF(N.GT.MDIM) CALL STOP('EQSYST: N.GT.MDIM')
        DO 120 L=1,N
          DO 100 K=N,L+1,-1
            IF(A(K,L).NE.0.0) GOTO 110
  100     CONTINUE
          K=L
  110     LASTN(L)=K
  120   CONTINUE
C
C  COLUMN LOOP: ELIMINATE ELEMENTS BELOW THE DIAGONAL IN COLUMN L.
C
        DO 500 L=1,N-1
C
C  STORE -A(K,L)/A(L,L) IN ELEMENT A(K,L)
C  MULTIPLY RIGHT HAND SIDE WITH -A(K,L)/A(L,L)
C
          ALL=A(L,L)
          BL=B(L)
          DO 200 K=L+1,LASTN(L)
            A(K,L)=-A(K,L)/ALL
            B(K)=A(K,L)*BL+B(K)
  200     CONTINUE
C
C  ADD FRACTION -A(K,L)/A(L,L) OF ROW L TO ROW K.
C  IN EACH COLUMN GO THROUGH ALL ROWS
C
          DO 400 M=L+1,N
            IF(A(L,M).NE.0.0) THEN
              ALM=A(L,M)
              LASTN(M)=MAX(LASTN(L),LASTN(M))
              DO 300 K=L+1,LASTN(L)
                A(K,M)=A(K,L)*ALM+A(K,M)
  300         CONTINUE
            END IF
  400     CONTINUE
  500   CONTINUE
      ELSE
        DO 700 L=1,N-1
          BL=B(L)
          DO 600 K=L+1,LASTN(L)
            B(K)=A(K,L)*BL+B(K)
  600     CONTINUE
  700   CONTINUE
      END IF
C
C BACKSUBSTITUTE
C
      DO 900 K=N,1,-1
        BK=B(K)
        DO 800 L=K+1,N
          BK=-A(K,L)*B(L)+BK
  800   CONTINUE
        B(K)=BK/A(K,K)
  900 CONTINUE
C
      RETURN
      END
C
C**********************************************************************
C