c-------------------------------------------------------------    
C
      SUBROUTINE DECOMP (NX,NN,A,SCALES,IPS,D) 
C
C     THIS ROUTINE DOES AN LU DECOMPOSITION ON A SQUARE MATRIX.
C     IT USES ROW SCALING AND ROW PIVOTING TO MAKE THE PROCESS
C     AS ACCURATE AS POSSIBLE.
C     THE PARAMETERS OF THE CALLING SEQUENCE ARE
C     N       THE ORDER OF THE MATRIX TO BE DECOMPOSED
C     NX      THE NUMBER OF LINES OF THE MATRIX A
C	      IN THE CALLING PROGRAM
C     A       THE MATRIX TO BE DECOMPOSED
C	      THE DECOMPOSED MATRIX. THE MAIN DIAGONAL AND EVERYTHING
C	      ABOVE IT MAKE UP THE UPPER TRIANGULAR MATRIX. THE LOWER
C	      TRIANGULAR MATRIX IS MADE UP OF EVERYTHING BELOW THE MAIN
C	      DIAGONAL. THE MAIN DIAGONAL FOR THE LOWER TRIANGULAR
C	      MATRIX IS ASSUMED TO BE ALL ONES.
C     D       THE DETERMINANT OF THE COEFFICIENT MATRIX
C     KERR    ERROR PRINT CONTROL AND RETURN ERROR FLAG
C             IF KERR SET NEGATIVE BY USER NO MESSAGE PRINTED
C             IF SET TO ZERO, ERROR MESSAGE WILL BE PRINTED.
C             ON RETURN, FOLLOWING CODE -
C             =1-MATRIX WITH ZERO ROW
C             =2-SINGULAR MATRIX
C             =3- SOLUTION OBTAINED
C     IPS    VECTOR CONTAINING PERMUTATED ROW NUMBERS OF A 
C
      DIMENSION A(NX,NX), SCALES(NX), IPS(NX)
      N=NN
      KERR=0
C
C     FIND THE SOLUTION IN THE CASE N=1
C
      IF(N.GT.1) GO TO 20 
      D=A(1,1)
      IPS(1)=1
      IF(D.EQ.0.) GO TO 10
      KERR=3 
      RETURN
   10 IF(KERR.EQ.0) CALL SING(1)
      KERR=1 
      RETURN
C
C	INITIALIZE
C
   20 ODD=1.
C
C	FIND THE SCALING FACTOR FOR EACH ROW AND STORE IN SCALES
C
      DO 40 I=1,N 
      IPS(I)=I 
      RMAX=0.
      DO 30 J=1,N
         RMAX=AMAX1(RMAX,ABS(A(I,J))) 
   30 CONTINUE
      IF(RMAX.EQ.0.) GO TO 110 
      SCALES(I)=1./RMAX
   40 CONTINUE
C
C     MAIN LOOP. GAUSSIAN ELIMINATION IS CARRIED OUT. LOOP IS GONE
C     THROUGH ONCE FOR EACH ELEMENT ON THE MAIN DIAGONAL . K IS THE
C     DIAGONAL ELEMENT NUMBER. LAST DIAGONAL ELEMENT IS AN EXCEPTION. 
C
      NM1=N-1
      DO 90 K=1,NM1 
         BIG=0.
         DO 50 I=K,N 
            IP=IPS(I) 
            SIZE=ABS(A(IP,K))*SCALES(IP)
            IF(BIG.GT.SIZE) GO TO 50
            BIG=SIZE
            IPIV=I
   50    CONTINUE
         IF(BIG.EQ.0.) GO TO 120
         IF(IPIV.EQ.K) GO TO 60
         J=IPS(K)
         IPS(K)=IPS(IPIV)
         IPS(IPIV)=J
         ODD=-ODD
   60    CONTINUE 
         KP=IPS(K) 
         PIVOT=A(KP,K)
         KP1=K+1
         DO 80 I=KP1,N
            IP=IPS(I)
            EM=A(IP,K)/PIVOT
            A(IP,K)=EM
            IF(EM.EQ.0.) GO TO 80
            DO 70 J=KP1,N
               A(IP,J)=A(IP,J)-EM*A(KP,J) 
   70       CONTINUE
   80    CONTINUE
   90 CONTINUE
      KP=IPS(N)
      IF(A(KP,N).EQ.0.) GO TO 120
C
C	CALCULATE THE DETERMINANT
C
      D=1.
      DO 100 I=1,N
         IP=IPS(I)
         D=D*A(IP,I)
  100 CONTINUE 
      D=D*ODD 
      KERR=3 
      RETURN
C
C	ERROR EXIST
C
  110 IF(KERR.EQ.0) CALL SING(1) 
      D=0.
      KERR=1
      RETURN
  120 IF(KERR.EQ.0) CALL SING(2) 
      D=0.
      KERR=2 
      RETURN 
      END
C
      SUBROUTINE SING(IWHY)
C
      GO TO (10,30),IWHY
   10 WRITE(UNIT=*, FMT=*) ' *** MATRIX WITH ZERO ROW IN DECOMPOSE'
      RETURN
   30 WRITE(UNIT=*, FMT=*) ' *** SINGULAR MATRIX IN DECOMPOSE.'
      WRITE(UNIT=*, FMT=*) '     ZERO DIVIDE IN SOLVE.'
      RETURN
      END
C
      SUBROUTINE SOLVE (NX,NN,A,B,X,IPS) 
C
C     THIS ROUTINE DOES FORWARD ELIMINATION AND BACK SUBSTITUTION
C     TO SOLVE THE EQUATION LU*X=B
C     THE PARAMETERS OF THE CALLING SEQUENCE ARE
C     N       THE ORDER OF LU
C     NX      THE NUMBER OF LINES OF THE MATRIX LU
C	      IN THE CALLING PROGRAM
C     A       THE DECOMPOSED MATRIX, L AND U ARE PACKED
C	      INTO ONE MATRIX BY DECOMP
C     B       THE RIGHT HAND SIDE OF THE EQUATION
C     IPS     VECTOR CONTAINING PERMUTATED ROW NUMBERS OF A
C	      THIS VECTOR IS BY DECOMP
C
      DIMENSION A(NX,NX), B(NX), X(NX), IPS(NX)
      N=NN
C
C     FIND THE SOLUTION IN THE CASE N=1
C
      IF(N.GT.1) GO TO 10
      X(1)=B(1)/A(1,1)
      RETURN
C
C     FORWARD ELIMINATION USING L AND STORING THE RESULT IN X
C
   10 NP1=N+1
      IP=IPS(1)
      X(1)=B(IP) 
      DO 30 I=2,N 
         IP=IPS(I)
         IM1=I-1
         SUM=0.
         DO 20 J=1,IM1
            SUM=SUM+A(IP,J)*X(J) 
   20    CONTINUE
         X(I)=B(IP)-SUM
   30 CONTINUE
C
C	BACK SUBSTITUTION USING U AND X
C
      IP=IPS(N) 
      X(N)=X(N)/A(IP,N)
      DO 50 IBACK=2,N
         I=NP1-IBACK
         IP=IPS(I)
         IP1=I+1 
         SUM=0.
         DO 40 J=IP1,N
               SUM=SUM+A(IP,J)*X(J)
   40    CONTINUE 
         X(I)=(X(I)-SUM)/A(IP,I)
   50 CONTINUE
      RETURN
      END
