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