MPI-AMRVAC 3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
Loading...
Searching...
No Matches
mod_lu.t
Go to the documentation of this file.
1!*******************************************************
2!* LU decomposition routines used by test_lu.f90 *
3!* *
4!* F90 version by J-P Moreau, Paris *
5!* (www.jpmoreau.fr) *
6!* --------------------------------------------------- *
7!* Reference: *
8!* *
9!* "Numerical Recipes By W.H. Press, B. P. Flannery, *
10!* S.A. Teukolsky and W.T. Vetterling, Cambridge *
11!* University Press, 1986" [BIBLI 08]. *
12!* *
13!*******************************************************
14MODULE mod_lu
15 implicit none
16
17CONTAINS
18
19! ***************************************************************
20! * Given an N x N matrix A, this routine replaces it by the LU *
21! * decomposition of a rowwise permutation of itself. A and N *
22! * are input. INDX is an output vector which records the row *
23! * permutation effected by the partial pivoting; D is output *
24! * as -1 or 1, depending on whether the number of row inter- *
25! * changes was even or odd, respectively. This routine is used *
26! * in combination with LUBKSB to solve linear equations or to *
27! * invert a matrix. Return code is 1, if matrix is singular. *
28! ***************************************************************
29 Subroutine ludcmp(A,N,INDX,D,CODE)
30 integer, parameter :: NMAX=10
31 double precision, parameter :: TINY=1.5d-16
32 double precision :: AMAX,DUM, SUM, A(N,N),VV(NMAX)
33 INTEGER :: CODE, D, INDX(N)
34 integer :: N
35 ! .. local ..
36 integer :: I, J, K, IMAX
37
38 d=1; code=0
39
40 DO i=1,n
41 amax=0.d0
42 DO j=1,n
43 IF (dabs(a(i,j)).GT.amax) amax=dabs(a(i,j))
44 END DO ! j loop
45 IF(amax.LT.tiny) THEN
46 code = 1
47 RETURN
48 END IF
49 vv(i) = 1.d0 / amax
50 END DO ! i loop
51
52 DO j=1,n
53 DO i=1,j-1
54 sum = a(i,j)
55 DO k=1,i-1
56 sum = sum - a(i,k)*a(k,j)
57 END DO ! k loop
58 a(i,j) = sum
59 END DO ! i loop
60 amax = 0.d0
61 DO i=j,n
62 sum = a(i,j)
63 DO k=1,j-1
64 sum = sum - a(i,k)*a(k,j)
65 END DO ! k loop
66 a(i,j) = sum
67 dum = vv(i)*dabs(sum)
68 IF(dum.GE.amax) THEN
69 imax = i
70 amax = dum
71 END IF
72 END DO ! i loop
73
74 IF(j.NE.imax) THEN
75 DO k=1,n
76 dum = a(imax,k)
77 a(imax,k) = a(j,k)
78 a(j,k) = dum
79 END DO ! k loop
80 d = -d
81 vv(imax) = vv(j)
82 END IF
83
84 indx(j) = imax
85 IF(dabs(a(j,j)) < tiny) a(j,j) = tiny
86
87 IF(j.NE.n) THEN
88 dum = 1.d0 / a(j,j)
89 DO i=j+1,n
90 a(i,j) = a(i,j)*dum
91 END DO ! i loop
92 END IF
93 END DO ! j loop
94
95 RETURN
96 END subroutine ludcmp
97
98
99! ******************************************************************
100! * Solves the set of N linear equations A . X = B. Here A is *
101! * input, not as the matrix A but rather as its LU decomposition, *
102! * determined by the routine LUDCMP. INDX is input as the permuta-*
103! * tion vector returned by LUDCMP. B is input as the right-hand *
104! * side vector B, and returns with the solution vector X. A, N and*
105! * INDX are not modified by this routine and can be used for suc- *
106! * cessive calls with different right-hand sides. This routine is *
107! * also efficient for plain matrix inversion. *
108! ******************************************************************
109 Subroutine lubksb(A,N,INDX,B)
110 double precision :: SUM, A(N,N),B(N)
111 INTEGER :: INDX(N), N
112 ! .. local ..
113 integer :: II, LL, I, J
114
115 ii = 0
116
117 DO i=1,n
118 ll = indx(i)
119 sum = b(ll)
120 b(ll) = b(i)
121 IF(ii.NE.0) THEN
122 DO j=ii,i-1
123 sum = sum - a(i,j)*b(j)
124 END DO ! j loop
125 ELSE IF(sum.NE.0.d0) THEN
126 ii = i
127 END IF
128 b(i) = sum
129 END DO ! i loop
130
131 DO i=n,1,-1
132 sum = b(i)
133 IF(i < n) THEN
134 DO j=i+1,n
135 sum = sum - a(i,j)*b(j)
136 END DO ! j loop
137 END IF
138 b(i) = sum / a(i,i)
139 END DO ! i loop
140
141 RETURN
142 END subroutine lubksb
143
144END MODULE mod_lu
subroutine lubksb(a, n, indx, b)
Definition mod_lu.t:110
subroutine ludcmp(a, n, indx, d, code)
Definition mod_lu.t:30