ROMS
Loading...
Searching...
No Matches
ludcmp.F File Reference
#include "cppdefs.h"
Include dependency graph for ludcmp.F:

Go to the source code of this file.

Functions/Subroutines

subroutine ludcmp (a, n, np, indx, d)
 

Function/Subroutine Documentation

◆ ludcmp()

subroutine ludcmp ( real(r8), dimension(np,np), intent(inout) a,
integer, intent(in) n,
integer, intent(in) np,
integer, dimension(n), intent(out) indx,
real(r8), intent(out) d )

Definition at line 2 of file ludcmp.F.

3!
4!git $Id$
5!================================================== Hernan G. Arango ===
6! Copyright (c) 2002-2025 The ROMS Group !
7! Licensed under a MIT/X style license !
8! See License_ROMS.md !
9!=======================================================================
10! !
11! Given an N x N matrix A, with physical dimension NP, this routine !
12! replaces it by the LU decomposition of a rowwise permutation of !
13! itself. A and N are input. A is output, arranged according to !
14! Crout algorithm; INDX is an output vector which records the row !
15! permutation effected by the partial pivoting; D is output as +1 !
16! or -1 depending on whether the number of row interchanges was !
17! even or odd, respectively. IER is output as 1 or 0 depending on !
18! whether the matrix A is singular or not, respectively. It is used !
19! in combination with LUBKSB to solve linear equations or invert a !
20! matrix. !
21! !
22! Reference: Press, W.H, et al., 1989: Numerical Recipes, The Art !
23! of Scientific Computing, pp 31-37. !
24! !
25!=======================================================================
26!
27 USE mod_kinds
28!
29! Imported variable declarations.
30!
31 integer, intent(in) :: n, np
32
33 integer, intent(out) :: indx(n)
34
35 real(r8), intent(inout) :: a(np,np)
36 real(r8), intent(out) :: d
37!
38! Local variable declarations.
39!
40 integer i, ier, imax, j, k
41
42 real(r8), parameter :: tiny = 1.0e-20_r8
43
44 real(r8), dimension(n) :: vv
45
46 real(r8) :: aamax, dum, MySum
47!
48!-----------------------------------------------------------------------
49! Replace matrix A by its LU decomposition.
50!-----------------------------------------------------------------------
51!
52 ier=0
53 d=1.0_r8
54 DO i=1,n
55 aamax=0.0_r8
56 DO j=1,n
57 IF (abs(a(i,j)).gt.aamax) aamax=abs(a(i,j))
58 END DO
59 IF (aamax.eq.0.0_r8) THEN
60 ier=1
61 ELSE
62 vv(i)=1.0_r8/aamax
63 END IF
64 END DO
65 IF (ier.eq.1) RETURN
66 DO j=1,n
67 IF (j.gt.1) THEN
68 DO i=1,j-1
69 mysum=a(i,j)
70 IF (i.gt.1) THEN
71 DO k=1,i-1
72 mysum=mysum-a(i,k)*a(k,j)
73 END DO
74 a(i,j)=mysum
75 END IF
76 END DO
77 END IF
78 aamax=0.0_r8
79 DO i=j,n
80 mysum=a(i,j)
81 IF (j.gt.1) THEN
82 DO k=1,j-1
83 mysum=mysum-a(i,k)*a(k,j)
84 END DO
85 a(i,j)=mysum
86 END IF
87 dum=vv(i)*abs(mysum)
88 IF (dum.ge.aamax) THEN
89 imax=i
90 aamax=dum
91 END IF
92 END DO
93 IF (j.ne.imax) THEN
94 DO k=1,n
95 dum=a(imax,k)
96 a(imax,k)=a(j,k)
97 a(j,k)=dum
98 END DO
99 d=-d
100 vv(imax)=vv(j)
101 END IF
102 indx(j)=imax
103 IF (j.ne.n) THEN
104 IF (a(j,j).eq.0.0_r8) a(j,j)=tiny
105 dum=1./a(j,j)
106 DO i=j+1,n
107 a(i,j)=a(i,j)*dum
108 END DO
109 END IF
110 END DO
111 IF (a(n,n).eq.0.0_r8) a(n,n)=tiny
112
113 RETURN

Referenced by set_avg_mod::set_detide_tile().

Here is the caller graph for this function: