c2--567--1---------2---------3---------4---------5---------6---------712
      Subroutine factor(a, px, n, pivot)
c
c     -- Factor/LU factorization/N. Nath/February '92
c
c     -- Function:  This subroutine performs an LU factorization of a 
c        general matrix a of dimension nxn.  If pivot is true, partial
c        pivoting (row interchange) is done:  p*a = l*u; if false, 
c        straightforward factorization is done:  a = l*u.  l and u are 
c        respectively a unit lower triangular matrix, with all diagonal
c        elements = 1, and an upper triangular matrix.  Upon 
c        factorization, a(i,j) is overwritten by l(i,j) for i > j 
c        (l(i,i) = 1 is not written) and by u(i,j) otherwise.  
C        p = e(n-1)...e(1) is a permutation, encoded in the integer
C        vector px(1:n-1); e(j) interchanges rows j and px(j).
c
c        The LU factorization provides an efficient and stable approach
c        in solving linear equations.  In opposition to Gaussian 
c        elimination, this approach does not initially require the 
c        source vector; thus, several linear equations involving the
c        same matrix can be solved simultaneously.  Pivoting increases
c        stability (less round-off errors), but it is not always 
c        needed (e.g., when a is diagonally dominant).  Partial 
c        pivoting (where the kth pivot is the largest entry in the 
c        subcolumn a(k:n,k)), as opposed to complete pivoting (where the
c        largest entry in the submatrix a(k:n,k:n) is permuted in 
c        in the (k,k) position), is quite satisfactory and almost as
c        efficient as no pivoting. 
c
c        To solve a*x = b (or, equivalently to perform ainv*b, where 
c        ainv = inverse of a):  (1) set c = p*b (a permutation of b's
c        components), (2) solve l*y = c (a lower triangular system), 
c        and (3) solve u*x = y (an upper triangular system).  These are 
c        done in a separate routine.  (Note:  if ainv is ever needed,
c        it is obtainable with b = unit matrix.)
c
c     -- Caution:  This routine destroys a (and replaces it by l(i,j) 
c        for i > j and by u(i,j) otherwise)!  Similarly, the adjoining 
c        routine destroys b (and replaces it by the solution)!   
c
            implicit none
c
c     -- Input & Output Variables:
c
            real        a(n,n)!matrix to be factored (@input);
                              !  l(i,j) for i>j, u(i,j) otherwise
                              !  (@output)   
            integer	px(n),!row permutation code (output)
     1                  n     !linear dimension of a (input)
            logical     pivot !true => partial pivoting; false
                              !  => no pivoting (input)
c
c     -- Other Variables:
c
            real        temp, norm, absv, eps
            parameter   (eps = 1e-18)
            integer     i, j, k, m
c
c     -- Ref.:  Golub, G. H., and Van Loan, C. F., Matrix Computations,
c        Johns Hopkins University Press, 1989
c
      if (pivot) go to 200
c
c           Straight factorization:  a = l*u 
c
      do j = 1,n
c
c              Solve l(1:j-1,1:j-1)u(1:j-1,j) = a(1:j-1,j) 
c                       !result:  u(1:j-1,j)
         do k = 1,j-1
            do i = k+1,j-1
               a(i,j) = a(i,j) - a(i,k)*a(k,j)
            end do
         end do
c
c              Compute v(j:n) = a(j:n,j) - l(j:n,1:j-1)u(1:j-1,j)
c				!partial result:  u(j,j) = v(j)
         do k = 1,j-1
            do i = j,n
               a(i,j) = a(i,j) - a(i,k)*a(k,j)
            end do
         end do
c
c              Compute l(j+1:n,j) = v(j+1:n)/v(j)
c                       !result:  l(j+1:n,j)
         if (a(j,j) .eq. 0.0) a(j,j) = eps
                        !bypasses division by zero for singular a
         do i = j+1,n
            a(i,j) = a(i,j)/a(j,j)
         end do

      end do

      return

 200  continue
c
c           Factorization with partial pivoting:  p*a = l*u
c
      do j = 1,n
c
c              Compute a( ,j) = e(j-1)..e(1)a( ,j); i.e., swap
c                a(1:j-1,j) <-> a(px(1:j-1),j)
c
         do i = 1,j-1
            temp = a(i,j)
            a(i,j) = a(px(i),j)
                        !px(i<j) is defined later
            a(px(i),j) = temp
         end do
c
c              Solve l(1:j-1,1:j-1)u(1:j-1,j) = a(1:j-1,j)
c                       !result:  u(1:j-1,j)
         do k = 1,j-1
            do i = k+1,j-1
               a(i,j) = a(i,j) - a(i,k)*a(k,j)
            end do
         end do
c
c              Compute v(j:n) = a(j:n,j) - l(j:n,1:j-1)u(1:j-1,j)
c
         do k = 1,j-1
            do i = j,n
               a(i,j) = a(i,j) - a(i,k)*a(k,j)
            end do
         end do
c
c              Determine e(j) and apply to v & l( ,1:j-1)
c
c              .. Find m such that abs(v(m)) = max(abs(v(j:n)))
c
         if (j .eq. n) return

         m = j
         norm = abs(a(j,j))
         do i = j+1,n
            absv = abs(a(i,j))
            if (absv .gt. norm) then
               m = i
               norm = absv
            end if
         end do
c
c              .. Set px(j) = m
c
         px(j) = m
c
c              .. Swap v(j) <-> v(m), l(j,1:j-1) <-> l(m,1:j-1)
c                       !result:  u(j,j) = v(j); l(j,1:j-1) 
         do k = 1,j
            temp = a(j,k)
            a(j,k) = a(m,k) 
            a(m,k) = temp
         end do
c
c              Compute l(j+1:n,j) = v(j+1:n)/v(j)
c                       !result:  l(j+1:n,j)
         if (a(j,j) .eq. 0.0) a(j,j) = eps
                        !bypasses division by zero for singular a
         do i = j+1,n
            a(i,j) = a(i,j)/a(j,j)
         end do

      end do
 
      return
      end

c     ------
      Subroutine inverse(a, b, px, n, nc, pivot)
c
c     -- Inverse/Solve a*x = b/N. Nath/February '92
c
c     -- Function:  This subroutine solves a*x = b (or performs ainv*b,
c        where ainv is inverse of a), as outlined in the adjoining LU 
c        factorization routine.  The solution replaces b.  a is an nxn
c        matrix; b is nxnc; px and other details are as in the adjoining
c        routine.
c
c     -- Caution:  This routine destroys b (and replaces it by the 
c        solution)!  
c
c     -- Input & Output Variables:
c
            real  	a(n,n),!lu-factored a:  l(i,j) for i>j,
                              !  u(i,j) otherwise (input)
     1                  b(n,nc)!source (or operand) (@input);
				      !  solution (@output)
            integer     px(n),!row permutation code (input)
     1                  n,    !linear dimension of a (input)
     2                  nc    !# columns in b (input)
            logical     pivot !true => partial pivoting; false
                              !  => no pivoting (input)
c
c     -- Other Variables:
c
            real        temp
            integer     i, j, k
c
       do j = 1,nc
c
c              Compute c = p*b, where p = e(n-1)...e(1)
c
         if (pivot) then
            do i = 1,n-1
               temp = b(i,j)
               b(i,j) = b(px(i),j)
               b(px(i),j) = temp
            end do
         end if
c
c              Solve l*y = c
c
         do k = 1,n-1
            do i = k+1,n
               b(i,j) = b(i,j) - a(i,k)*b(k,j)
            end do
         end do
c
c              Solve u*x = y
c
         do k = n,2,-1
            b(k,j) = b(k,j)/a(k,k)
            do i = 1,k-1
               b(i,j) = b(i,j) - a(i,k)*b(k,j)
            end do
         end do
         b(1,j) = b(1,j)/a(1,1)

      end do

      return
      end