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    ii = 0
22    do i = 1 , n
23       ll = indx(i)
24       sum = b(ll)
25       b(ll) = b(i)
26 
27       if (ii /= 0) then
28          do j = ii , i - 1
29             sum = sum - a(i,j) * b(j)
30          end do
31       else if (sum /= 0.0) then
32          ii = i
33       end if
34       b(i) = sum
35    end do
36 
37    do i = n , 1 , -1
38       sum = b(i)
39       if (i < n) then
40          do j = i + 1 , n
41             sum = sum - a(i,j) * b(j)
42          end do
43       end if
44       b(i) = sum / a(i,i)
45    end do
46 
47 end subroutine da_lubksb
48 
49