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