3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
28
29
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
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
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