mod_lu.t

Go to the documentation of this file.

!*******************************************************
!*    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" [BIBLI 08].                *
!*                                                     * 
!*******************************************************
MODULE mod_lu
  implicit none
 
CONTAINS
 
!  ***************************************************************
!  * 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 inter-   *
!  * changes 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.   *
!  ***************************************************************
  Subroutine ludcmp(A,N,INDX,D,CODE)
 integer, parameter           :: NMAX=10
 double precision, parameter  :: TINY=1.5d-16
 double precision             :: AMAX,DUM, SUM, A(N,N),VV(NMAX)
 INTEGER                      :: CODE, D, INDX(N)
 integer                      :: N
 ! .. local ..
 integer                      :: I, J, K, IMAX
 
 d=1; code=0
 
 DO i=1,n
   amax=0.d0
   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
 
 RETURN
 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)
 double precision        :: SUM, A(N,N),B(N)
 INTEGER                 :: INDX(N), N
 ! .. local ..
 integer                 :: II, LL, I, J
 
 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
 
 RETURN
 END subroutine lubksb
 
END MODULE mod_lu