C ----------------------------------------------------------------------
C           Subroutine LUFAC.FOR by www.numerical-methods.com           |
C ----------------------------------------------------------------------
C
C Subroutine LUfac returns the components of the matrices P, L and U in
C  the LU factorisation of a real matrix A
C
C                            PA  = LU
C
C where A is an n by n real matrix and b is a given vector
C
C SUBROUTINE LUFAC(A, N, PERM, LFAIL)
C real        a: on input the nxn matrix A, on output L and U
C integer     n: the dimension of the matrix/vector
C integer  perm: an n-vector, the column index of the permutation matrix P
C logical lfail: returns 'true' if the method fails, otherwise 'false'
C
C As a result of the subroutine, the input array a is over-written. The
C lower triangle of the returned matrix with a set of 1s down the
C diagonal is the matrix L. The upper triangle and the diagonal of the
C returned matrix is the matrx U. perm records the row interchanges that
C result from pivotting; perm(i)=j means that the P(i,j)=1.
C
C Source of the code: http://www.numerical-methods.com/fortran/LUFAC.FOR
C Source of the user-guide: http://www.numerical-methods.com/fortran/
C  lufac.htm
C
C Licence: This is 'open source'; the software may be used and applied
C  within other systems as long as its provenance is appropriately
C  acknowledged. See the GNU Licence http://www.gnu.org/licenses/lgpl.txt
C  for more information or contact webmaster@numerical-methods.com
 
      SUBROUTINE LUFAC(MAXN, N, A, PERM, LFAIL, ROWNORMS)

C Input parameters
C ----------------
C The limit on the dimension of the matrices A and B 
      INTEGER    MAXN
C The dimension of the matrices
      INTEGER    N
C The matrix A
      REAL*8     A(MAXN,MAXN)
C The vector PERM
      INTEGER    PERM(MAXN)
C Failure flag
      LOGICAL    LFAIL
C      
      REAL*8     ROWNORMS(MAXN)
      
C Local variables
      REAL*8     TINY
      INTEGER    I,IMAX, J, K
      REAL*8     UIICOLMAX, DUM, SUM, RELDIVISOR
      REAL*8     ANORM

C Initialisation
      TINY = 0.0000000001
      LFAIL=.FALSE.

C  Initialise the permutation matrix, as the identity matrix

      DO 10 I = 1,N
        PERM(I)=I
10    CONTINUE

      
C Set the row norms and matrix norm ANORM
      ANORM=0.0D0
      DO 30 I=1,N
        ROWNORMS(I)=0.0D0
        DO 40 J=1,N
          IF (ABS(A(I,J)).GT.ROWNORMS(I)) THEN
            ROWNORMS(I)=ABS(A(I,J))
          END IF
40      CONTINUE            
        IF (ROWNORMS(I).LT.TINY) THEN
          WRITE(*,*) 'Matrix is Singular or very Ill-conditioned'
          WRITE(*,*) 'LU factorisation is abandoned'
          LFAIL=.TRUE.
       END IF
        IF (ROWNORMS(I).GT.ANORM) THEN 
          ANORM=ROWNORMS(I)
        END IF
30    CONTINUE
 
C The main loop is a counter effectively down the diagonal
      DO 70 I=1,N

c  The pivotting method
c   Pass down from the diagonal to check which row would have the
C    largest divisor.
      UIICOLMAX=0.0D0
        DO 80 J=I,N
          SUM=0.0D0
          DO 90 K=1,I-1
            SUM=SUM+A(J,K)*A(K,I)
90        CONTINUE
          RELDIVISOR=(ABS(A(J,I)-SUM))/ROWNORMS(PERM(I))
          IF (RELDIVISOR.GT.UIICOLMAX) THEN
            UIICOLMAX=RELDIVISOR
            IMAX=J
          END IF
80      CONTINUE
  
C  If the largest divisor is zero or very small then the matrix is
C   singular or ill-conditioned and the factorisation is abandoned
        IF (UIICOLMAX.LT.TINY*ANORM) THEN
          WRITE(*,*) 'Matrix is Singular or very Ill-conditioned'
          WRITE(*,*) 'LU factorisation is abandoned'
          LFAIL=.TRUE.
        END IF
  
C  The rows are swapped and the exchange is recorded in perm
        IF (I.NE.IMAX) THEN 
          DO 100 K=1,N
            DUM=A(IMAX,K)
            A(IMAX,K)=A(I,K)
            A(I,K)=DUM
 100      CONTINUE  
          PERMI = PERM(I)
          PERM(I) = PERM(IMAX)
          PERM(IMAX) = PERMI
        END IF
        
C  The first inner loop counts j=i..n, so that a(i,j) addresses the diagonal and super-
C   diagonal elements of a on row i. That is ith row of the U matrix is formed.
        DO 200 J=I,N
          SUM=0.0D0
          DO 210 K=1,I-1
            SUM=SUM+A(I,K)*A(K,J)
210       CONTINUE
          A(I,J)=A(I,J)-SUM
200     CONTINUE       

C  The second  inner loop counts j=i..n, so that a(j,i) addresses the sub-
C   diagonal elements of A on column i. That is the ith column of the L matrix is formed.
        DO 300 J=I+1,N
          SUM=0.0D0
          DO 310 K=1,I-1
            SUM=SUM+A(J,K)*A(K,I)
310       CONTINUE
          A(J,I)=(A(J,I)-SUM)/A(I,I)
300     CONTINUE            

70    CONTINUE

      END