matrix_inversion.f90

!WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
!    LU decomposition routines used by test_lu.f90
!
!                 F90 version by J-P Moreau, Paris (www.jpmoreau.fr)
!
! Reference:
! Numerical Recipes By W.H. Press, B. P. Flannery, S.A. Teukolsky
! and W.T. Vetterling, Cambridge University Press, 1986
!
!WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
MODULE lu

implicit none

CONTAINS

!WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
! Given an N x N matrix A, this routine replaces it by the LU decomposition of
! a rowwise permutation of itself. A and N are input. INDX is an output vector
! which records the row permutation effected by the partial pivoting; D is
! output as -1 or 1, depending on whether the number of row interchanges was
! even or odd, respectively. This routine is used in combination with LUBKSB to
! solve linear equations or to invert a matrix. Return code is 1, if matrix is
! singular.
!WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW

  Subroutine ludcmp(A,N,INDX,D,CODE)

    parameter(nmax=100,tiny=1.5d-16)
    INTEGER                                  :: code, d, indx(n)
    REAL                                     :: amax,dum, sum, a(n,n),vv(nmax)

    d = 1
    code = 0

    DO i = 1, n
       amax = 0.0
       DO j = 1, n
          IF (dabs(a(i,j)).GT.amax) amax=dabs(a(i,j))
       END DO ! j loop
       IF(amax.LT.tiny) THEN
          code = 1
          RETURN
       END IF
       vv(i) = 1.d0 / amax
    END DO ! i loop

    DO j=1,n
       DO i=1,j-1
          sum = a(i,j)
          DO k=1,i-1
             sum = sum - a(i,k)*a(k,j) 
          END DO ! k loop
          a(i,j) = sum
       END DO ! i loop
       amax = 0.d0
       DO i=j,n
          sum = a(i,j)
          DO k=1,j-1
             sum = sum - a(i,k)*a(k,j) 
          END DO ! k loop
          a(i,j) = sum
          dum = vv(i)*dabs(sum)
          IF(dum.GE.amax) THEN
             imax = i
             amax = dum
          END IF
       END DO ! i loop  

       IF(j.NE.imax) THEN
          DO k=1,n
             dum = a(imax,k)
             a(imax,k) = a(j,k)
             a(j,k) = dum
          END DO ! k loop
          d = -d
          vv(imax) = vv(j)
       END IF

       indx(j) = imax
       IF(dabs(a(j,j)) < tiny) a(j,j) = tiny

       IF(j.NE.n) THEN
          dum = 1.d0 / a(j,j)
          DO i=j+1,n
             a(i,j) = a(i,j)*dum
          END DO ! i loop
       END IF 
    END DO ! j loop

  END subroutine ludcmp


!  ******************************************************************
!  * Solves the set of N linear equations A . X = B.  Here A is     *
!  * input, not as the matrix A but rather as its LU decomposition, *
!  * determined by the routine LUDCMP. INDX is input as the permuta-*
!  * tion vector returned by LUDCMP. B is input as the right-hand   *
!  * side vector B, and returns with the solution vector X. A, N and*
!  * INDX are not modified by this routine and can be used for suc- *
!  * cessive calls with different right-hand sides. This routine is *
!  * also efficient for plain matrix inversion.                     *
!  ******************************************************************
  Subroutine lubksb(A,N,INDX,B)
    REAL*8  sum, a(n,n),b(n)
    INTEGER indx(n)

    ii = 0

    DO i=1,n
       ll = indx(i)
       sum = b(ll)
       b(ll) = b(i)
       IF(ii.NE.0) THEN
          DO j=ii,i-1
             sum = sum - a(i,j)*b(j)
          END DO ! j loop
       ELSE IF(sum.NE.0.d0) THEN
          ii = i
       END IF
       b(i) = sum
    END DO ! i loop

    DO i=n,1,-1
       sum = b(i)
       IF(i < n) THEN
          DO j=i+1,n
             sum = sum - a(i,j)*b(j)
          END DO ! j loop
       END IF
       b(i) = sum / a(i,i)
    END DO ! i loop

  END subroutine lubksb

END MODULE lu