1!/*****************************************************************************/
2! *
3! *  Elmer, A Finite Element Software for Multiphysical Problems
4! *
5! *  Copyright 1st April 1995 - , CSC - IT Center for Science Ltd., Finland
6! *
7! * This library is free software; you can redistribute it and/or
8! * modify it under the terms of the GNU Lesser General Public
9! * License as published by the Free Software Foundation; either
10! * version 2.1 of the License, or (at your option) any later version.
11! *
12! * This library is distributed in the hope that it will be useful,
13! * but WITHOUT ANY WARRANTY; without even the implied warranty of
14! * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
15! * Lesser General Public License for more details.
16! *
17! * You should have received a copy of the GNU Lesser General Public
18! * License along with this library (in file ../LGPL-2.1); if not, write
19! * to the Free Software Foundation, Inc., 51 Franklin Street,
20! * Fifth Floor, Boston, MA  02110-1301  USA
21! *
22! *****************************************************************************/
23!
24!/******************************************************************************
25! *
26! *  Authors: Juha Ruokolainen
27! *  Email:   Juha.Ruokolainen@csc.fi
28! *  Web:     http://www.csc.fi/elmer
29! *  Address: CSC - IT Center for Science Ltd.
30! *           Keilaranta 14
31! *           02101 Espoo, Finland
32! *
33! *  Original Date: 01 Oct 1996
34! *
35! *****************************************************************************/
36
37!> \ingroup ElmerLib
38!> \}
39
40!------------------------------------------------------------------------------
41!>  LU Decomposition & matrix inverse of small full matrices.
42!> Don't use this for anything big,
43!>  use for example LAPACK routines instead.
44!------------------------------------------------------------------------------
45
46MODULE LUDecomposition
47
48  USE Types
49  IMPLICIT NONE
50
51 CONTAINS
52
53  SUBROUTINE InvertMatrix( A,n )
54
55    INTEGER :: n
56    REAL(KIND=dp) :: A(:,:)
57
58    REAL(KIND=dp) :: s
59    INTEGER :: i,j,k
60    INTEGER :: pivot(n)
61
62   ! /*
63   !  *  AP = LU
64   !  */
65   CALL LUDecomp( a,n,pivot )
66
67   DO i=1,n
68     IF ( ABS(A(i,i)) == 0.0d0 ) THEN
69       CALL Error( 'InvertMatrix', 'Matrix is singular.' )
70       RETURN
71     END IF
72     A(i,i) = 1.0d0 / A(i,i)
73   END DO
74
75   ! /*
76   !  *  INV(U)
77   !  */
78   DO i=N-1,1,-1
79     DO j=N,i+1,-1
80       s = -A(i,j)
81       DO k=i+1,j-1
82         s = s - A(i,k)*A(k,j)
83       END DO
84       A(i,j) = s
85     END DO
86   END DO
87
88   ! /*
89   !  * INV(L)
90   !  */
91   DO i=n-1,1,-1
92     DO j=n,i+1,-1
93       s = 0.D00
94       DO k=i+1,j
95         s = s - A(j,k)*A(k,i)
96       END DO
97       A(j,i) = A(i,i)*s
98     END DO
99   END DO
100
101   ! /*
102   !  * A  = INV(AP)
103   !  */
104   DO i=1,n
105     DO j=1,n
106       s = 0.0D0
107       DO k=MAX(i,j),n
108         IF ( k /= i ) THEN
109           s = s + A(i,k)*A(k,j)
110         ELSE
111           s = s + A(k,j)
112         END IF
113       END DO
114       A(i,j) = s
115     END DO
116   END DO
117
118   ! /*
119   !  * A = INV(A) (at last)
120   !  */
121   DO i=n,1,-1
122     IF ( pivot(i) /= i ) THEN
123       DO j = 1,n
124         s = A(i,j)
125         A(i,j) = A(pivot(i),j)
126         A(pivot(i),j) = s
127       END DO
128     END IF
129   END DO
130
131 END SUBROUTINE InvertMatrix
132
133
134 SUBROUTINE LUSolve( n,A,x )
135   REAL(KIND=dp) :: A(n,n)
136   REAL(KIND=dp) :: x(n)
137   INTEGER :: n
138
139   REAL(KIND=dp) :: s
140   INTEGER :: i,j,k
141   INTEGER :: pivot(n)
142
143   ! /*
144   !  *  AP = LU
145   !  */
146   CALL LUDecomp( A,n,pivot )
147
148   DO i=1,n
149     IF ( ABS(A(i,i)) == 0.0d0 ) THEN
150       CALL Error( 'LUSolve', 'Matrix is singular.' )
151       RETURN
152     END IF
153     A(i,i) = 1.0d0 / A(i,i)
154   END DO
155
156   !
157   ! Forward substitute
158   DO i=1,n
159      s = x(i)
160      DO j=1,i-1
161         s = s - A(i,j) * x(j)
162      END DO
163      x(i) = A(i,i) * s
164   END DO
165
166   !
167   ! Backward substitute (solve x from Ux = z)
168   DO i=n,1,-1
169      s = x(i)
170      DO j=i+1,n
171         s = s - A(i,j) * x(j)
172      END DO
173      x(i) = s
174   END DO
175
176   DO i=n,1,-1
177      IF ( pivot(i) /= i ) THEN
178         s = x(i)
179         x(i) = x(pivot(i))
180         x(pivot(i)) = s
181      END IF
182   END DO
183
184  END SUBROUTINE LUSolve
185
186
187!/*
188! * LU- decomposition by gaussian elimination. Row pivoting is used.
189! *
190! * result : AP = L'U ; L' = LD; pivot[i] is the swapped column number
191! * for column i.
192! *
193! * Result is stored in place of original matrix.
194! *
195! */
196  SUBROUTINE LUDecomp( a,n,pivot )
197
198    REAL(KIND=dp), DIMENSION (:,:) :: a
199    INTEGER :: n
200    INTEGER, DIMENSION (:) :: pivot
201
202    INTEGER :: i,j,k,l
203    REAL(KIND=dp) :: swap
204
205    DO i=1,n
206      j = i
207      DO k=i+1,n
208        IF ( ABS(A(i,k)) > ABS(A(i,j)) ) j = k
209      END DO
210
211      IF ( ABS(A(i,j)) == 0.0d0) THEN
212        CALL Error( 'LUDecomp', 'Matrix is singluar.' )
213        RETURN
214      END IF
215
216      pivot(i) = j
217
218      IF ( j /= i ) THEN
219        DO k=1,i
220          swap = A(k,j)
221          A(k,j) = A(k,i)
222          A(k,i) = swap
223        END DO
224      END IF
225
226      DO k=i+1,n
227        A(i,k) = A(i,k) / A(i,i)
228      END DO
229
230      DO k=i+1,n
231        IF ( j /= i ) THEN
232          swap = A(k,i)
233          A(k,i) = A(k,j)
234          A(k,j) = swap
235        END IF
236
237        DO  l=i+1,n
238          A(k,l) = A(k,l) - A(k,i) * A(i,l)
239        END DO
240      END DO
241    END DO
242
243    pivot(n) = n
244    IF ( ABS(A(n,n)) == 0.0d0 ) THEN
245      CALL Error( 'LUDecomp',  'Matrix is (at least almost) singular.' )
246    END IF
247
248  END SUBROUTINE LUDecomp
249
250
251
252  SUBROUTINE ComplexInvertMatrix( A,n )
253
254    COMPLEx(KIND=dp), DIMENSION(:,:) :: A
255    INTEGER :: n
256
257    COMPLEX(KIND=dp) :: s
258    INTEGER :: i,j,k
259    INTEGER :: pivot(n)
260
261   ! /*
262   !  *  AP = LU
263   !  */
264   CALL ComplexLUDecomp( a,n,pivot )
265
266   DO i=1,n
267     IF ( ABS(A(i,i))==0.0d0 ) THEN
268       CALL Error( 'ComplexInvertMatrix', 'Matrix is singular.' )
269       RETURN
270     END IF
271     A(i,i) = 1.0D0/A(i,i)
272   END DO
273
274   ! /*
275   !  *  INV(U)
276   !  */
277   DO i=N-1,1,-1
278     DO j=N,i+1,-1
279       s = -A(i,j)
280       DO k=i+1,j-1
281         s = s - A(i,k)*A(k,j)
282       END DO
283       A(i,j) = s
284     END DO
285   END DO
286
287   ! /*
288   !  * INV(L)
289   !  */
290   DO i=n-1,1,-1
291     DO j=n,i+1,-1
292       s = 0.D00
293       DO k=i+1,j
294         s = s - A(j,k)*A(k,i)
295       END DO
296       A(j,i) = A(i,i)*s
297     END DO
298   END DO
299
300   ! /*
301   !  * A  = INV(AP)
302   !  */
303   DO i=1,n
304     DO j=1,n
305       s = 0.0D0
306       DO k=MAX(i,j),n
307         IF ( k /= i ) THEN
308           s = s + A(i,k)*A(k,j)
309         ELSE
310           s = s + A(k,j)
311         END IF
312       END DO
313       A(i,j) = s
314     END DO
315   END DO
316
317   ! /*
318   !  * A = INV(A) (at last)
319   !  */
320   DO i=n,1,-1
321     IF ( pivot(i) /= i ) THEN
322       DO j = 1,n
323         s = A(i,j)
324         A(i,j) = A(pivot(i),j)
325         A(pivot(i),j) = s
326       END DO
327     END IF
328   END DO
329
330 END SUBROUTINE ComplexInvertMatrix
331
332!/*
333! * LU- decomposition by gaussian elimination. Row pivoting is used.
334! *
335! * result : AP = L'U ; L' = LD; pivot[i] is the swapped column number
336! * for column i.
337! *
338! * Result is stored in place of original matrix.
339! *
340! */
341
342  SUBROUTINE ComplexLUDecomp( a,n,pivot )
343
344    COMPLEX(KIND=dp), DIMENSION (:,:) :: a
345    INTEGER :: n
346    INTEGER, DIMENSION (:) :: pivot
347
348    INTEGER :: i,j,k,l
349    COMPLEX(KIND=dp) :: swap
350
351    DO i=1,n
352      j = i
353      DO k=i+1,n
354        IF ( ABS(A(i,k)) > ABS(A(i,j)) ) j = k
355      END DO
356
357      IF ( ABS(A(i,j))==0.0d0 ) THEN
358        CALL Error( 'ComplexLUDecomp', 'Matrix is singluar.' )
359        RETURN
360      END IF
361
362      pivot(i) = j
363
364      IF ( j /= i ) THEN
365        DO k=1,i
366          swap = A(k,j)
367          A(k,j) = A(k,i)
368          A(k,i) = swap
369        END DO
370      END IF
371
372      DO k=i+1,n
373        A(i,k) = A(i,k) / A(i,i)
374      END DO
375
376      DO k=i+1,n
377        IF ( j /= i ) THEN
378          swap = A(k,i)
379          A(k,i) = A(k,j)
380          A(k,j) = swap
381        END IF
382
383        DO  l=i+1,n
384          A(k,l) = A(k,l) - A(k,i) * A(i,l)
385        END DO
386      END DO
387    END DO
388
389    pivot(n) = n
390    IF ( ABS(A(n,n))==0.0d0 ) THEN
391      CALL Error( 'ComplexLUDecomp', 'Matrix is (at least almost) singular.' )
392    END IF
393
394  END SUBROUTINE ComplexLUDecomp
395
396
397END MODULE LUDecomposition
398
399!> \} ElmerLib
400