      SUBROUTINE EQSYSI(N,NDIM,A,B,NEWMAT)
C
C  SOLVES THE EQUATION  SYSIEM  A*X=B.
C  WHEN NEWMAT=TRUE, THE SYSIEM 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 EQSYSI 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 
C  EQSYST WITH 1 IMPROVEMENT
C
      INCLUDE 'PREC'
      PARAMETER (MDIM=2000)
      DIMENSION A(NDIM,NDIM),B(NDIM), AA(MDIM,MDIM),BB(MDIM)
     *,RESID(MDIM)
      INTEGER LASTN(MDIM)
      LOGICAL NEWMAT
      SAVE LASTN
C
      write(*,*) 'eqsysi'
      DO I=1,NDIM
         DO J=1,NDIM
            AA(I,J)=A(I,J)
         ENDDO
         BB(I)=B(I)
      ENDDO


      IF(NEWMAT) THEN
C
C     FIND THE LAST NON-ZERO ELEMENT IN EACH COLUMN
C
        IF(N.GT.MDIM) CALL STOP('EQSYSI: 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
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C IMPROVE
C  

      DO NIMP=1,3
         write(*,*) 'NIMP = ',NIMP
         DO J=1,NDIM
            TMP=0.0
            DO K=1,NDIM
               TMP=AA(J,K)*B(K) +TMP
            ENDDO
            RESID(J)=TMP-BB(J)
         ENDDO
C     
C     BACKSUBSTITUTE
C     
         DO 901 K=N,1,-1
            BK=RESID(K)
            DO 801 L=K+1,N
               BK=-AA(K,L)*RESID(L)+BK
 801        CONTINUE
            RESID(K)=BK/A(K,K)
            write(*,*) K,B(K),RESID(K)
            B(K)=B(K)-RESID(K)
 901     CONTINUE
      enddo
      RETURN
      END
C     
C**********************************************************************
C
