var/da/da_tools/da_lubksb.inc

   1 subroutine da_lubksb(n, np, indx, a, b)
   2
   3    !-----------------------------------------------------------------------
   4    ! Purpose: Adapted Numerical Recipes routine to solve the set of n linear
   5    ! equations A.X=B.
   6    ! Routine takes in to account possibility that B will begin with many zero elements,
   7    ! 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(in)    :: indx(1:n)      ! Permutation vector returned by LUDCMP.
  15    real,    intent(in)    :: a(1:np,1:np)   ! LU decomposition of matrix A in A.x=B.
  16    real,    intent(inout) :: b(1:n)         ! On input = B, on output = x.
  17
  18    integer :: i , ii , j , ll
  19    real    :: sum
  20
  21    if (trace_use) call da_trace_entry("da_lubksb")
  22
  23    ii = 0
  24    do i = 1 , n
  25       ll = indx(i)
  26       sum = b(ll)
  27       b(ll) = b(i)
  28
  29       if (ii /= 0) then
  30          do j = ii , i - 1
  31             sum = sum - a(i,j) * b(j)
  32          end do
  33       else if (sum /= 0.0) then
  34          ii = i
  35       end if
  36       b(i) = sum
  37    end do
  38
  39    do i = n , 1 , -1
  40       sum = b(i)
  41       if (i < n) then
  42          do j = i + 1 , n
  43             sum = sum - a(i,j) * b(j)
  44          end do
  45       end if
  46       b(i) = sum / a(i,i)
  47    end do
  48
  49    if (trace_use) call da_trace_exit("da_lubksb")
  50
  51 end subroutine da_lubksb
  52
  53