da_ludcmp.inc
References to this file elsewhere.
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 d = 1.0
24 do i = 1 , n
25 aamax = 0.0
26
27 do j = 1 , n
28 if (abs(a(i,j)) > aamax) aamax = abs(a(i,j))
29 end do
30 if (aamax == 0.0) then
31 call da_error(__FILE__,__LINE__,(/"Singular matrix"/))
32 end if
33 vv(i) = 1.0 / aamax
34 end do
35
36 do j = 1 , n
37 if (j > 1) then
38 do i = 1 , j - 1
39 sum = a(i,j)
40 if (i > 1) then
41 do k = 1 , i - 1
42 sum = sum - a(i,k) * a(k,j)
43 end do
44 a(i,j) = sum
45 end if
46 end do
47 end if
48
49 aamax = 0.0
50 do i = j , n
51 sum = a(i,j)
52 if (j > 1) then
53 do k = 1 , j - 1
54 sum = sum - a(i,k) * a(k,j)
55 end do
56 a(i,j) = sum
57 end if
58 dum = vv(i) * abs(sum)
59 if (dum >= aamax) then
60 imax = i
61 aamax = dum
62 end if
63 end do
64
65 if (j /= imax) then
66 do k = 1 , n
67 dum = a(imax,k)
68 a(imax,k) = a(j,k)
69 a(j,k) = dum
70 end do
71 d = -d
72 vv(imax) = vv(j)
73 end if
74
75 indx(j) = imax
76 if (j /= n) then
77 if (a(j,j) == 0.0) a(j,j) = tiny
78 dum = 1.0 / a(j,j)
79 do i = j + 1 , n
80 a(i,j) = a(i,j) * dum
81 end do
82 end if
83 end do
84
85 if (a(n,n) == 0.0) a(n,n) = tiny
86
87 end subroutine da_ludcmp
88
89