MPI-AMRVAC  3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
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 !*******************************************************
14 MODULE mod_lu
15  implicit none
16 
17 CONTAINS
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 
144 END MODULE mod_lu
Definition: mod_lu.t:14
subroutine ludcmp(A, N, INDX, D, CODE)
Definition: mod_lu.t:30
subroutine lubksb(A, N, INDX, B)
Definition: mod_lu.t:110