var/da/da_tools/da_ludcmp.inc

   1 subroutine da_ludcmp(n, np, indx, a, d)
   2
   3    !-----------------------------------------------------------------------
   4    ! Purpose: Adapted Numerical Recipes routine to solve the set of n linear
   5    ! equations
   6    ! A.X=B. Routine takes in to account possibility that B will begin with many
   7    ! zero elements, so it is efficient for matrix inversion.
   8    !-----------------------------------------------------------------------
   9
  10    implicit none
  11
  12    integer, intent(in)    :: n           ! Logical size of array.
  13    integer, intent(in)    :: np          ! Physical size of array.
  14    integer, intent(out)   :: indx(1:n)   ! Permutation vector returned by LUDCMP
  15    real,    intent(inout) :: a(1:np,1:np)! LU decomposition of matrix A in A.x=B
  16    real,    intent(out)   :: d           ! On input = B, on output = x.
  17
  18    real, parameter      :: tiny = 1.0e-20
  19    real                 :: aamax , dum , sum
  20    integer              :: i , imax , j , k
  21    real                 :: vv(1:np)
  22
  23    if (trace_use) call da_trace_entry("da_ludcmp")
  24
  25    d = 1.0
  26    do i = 1 , n
  27       aamax = 0.0
  28
  29       do j = 1 , n
  30          if (abs(a(i,j)) > aamax) aamax = abs(a(i,j))
  31       end do
  32       if (aamax == 0.0) then
  33          call da_error(__FILE__,__LINE__,(/"Singular matrix"/))
  34       end if
  35       vv(i) = 1.0 / aamax
  36    end do
  37
  38    do j = 1 , n
  39       if (j > 1) then
  40          do i = 1 , j - 1
  41             sum = a(i,j)
  42             if (i > 1) then
  43                do k = 1 , i - 1
  44                   sum = sum - a(i,k) * a(k,j)
  45                end do
  46                a(i,j) = sum
  47             end if
  48          end do
  49       end if
  50
  51       aamax = 0.0
  52       do i = j , n
  53          sum = a(i,j)
  54          if (j > 1) then
  55             do k = 1 , j - 1
  56                sum = sum - a(i,k) * a(k,j)
  57             end do
  58             a(i,j) = sum
  59          end if
  60          dum = vv(i) * abs(sum)
  61          if (dum >= aamax) then
  62             imax = i
  63             aamax = dum
  64          end if
  65       end do
  66
  67       if (j /= imax) then
  68          do k = 1 , n
  69             dum = a(imax,k)
  70             a(imax,k) = a(j,k)
  71             a(j,k) = dum
  72          end do
  73          d = -d
  74          vv(imax) = vv(j)
  75       end if
  76
  77       indx(j) = imax
  78       if (j /= n) then
  79          if (a(j,j) == 0.0) a(j,j) = tiny
  80          dum = 1.0 / a(j,j)
  81          do i = j + 1 , n
  82             a(i,j) = a(i,j) * dum
  83          end do
  84       end if
  85    end do
  86
  87    if (a(n,n) == 0.0) a(n,n) = tiny
  88
  89    if (trace_use) call da_trace_exit("da_ludcmp")
  90
  91 end subroutine da_ludcmp
  92
  93