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