1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
|
subroutine nnfc
* (n, r,c,ic, ia,ja,a, z, b,
* lmax,il,jl,ijl,l, d, umax,iu,ju,iju,u,
* row, tmp, irl,jrl, flag)
clll. optimize
c*** subroutine nnfc
c*** numerical ldu-factorization of sparse nonsymmetric matrix and
c solution of system of linear equations (compressed pointer
c storage)
c
c
c input variables.. n, r, c, ic, ia, ja, a, b,
c il, jl, ijl, lmax, iu, ju, iju, umax
c output variables.. z, l, d, u, flag
c
c parameters used internally..
c nia - irl, - vectors used to find the rows of l. at the kth step
c nia - jrl of the factorization, jrl(k) points to the head
c - of a linked list in jrl of column indices j
c - such j .lt. k and l(k,j) is nonzero. zero
c - indicates the end of the list. irl(j) (j.lt.k)
c - points to the smallest i such that i .ge. k and
c - l(i,j) is nonzero.
c - size of each = n.
c fia - row - holds intermediate values in calculation of u and l.
c - size = n.
c fia - tmp - holds new right-hand side b* for solution of the
c - equation ux = b*.
c - size = n.
c
c internal variables..
c jmin, jmax - indices of the first and last positions in a row to
c be examined.
c sum - used in calculating tmp.
c
integer rk,umax
integer r(1), c(1), ic(1), ia(1), ja(1), il(1), jl(1), ijl(1)
integer iu(1), ju(1), iju(1), irl(1), jrl(1), flag
double precision a(1), l(1), d(1), u(1), z(1), b(1), row(1)
double precision tmp(1), lki, sum, dk
c
c ****** initialize pointers and test storage ***********************
if(il(n+1)-1 .gt. lmax) go to 104
if(iu(n+1)-1 .gt. umax) go to 107
do 1 k=1,n
irl(k) = il(k)
jrl(k) = 0
1 continue
c
c ****** for each row ***********************************************
do 19 k=1,n
c ****** reverse jrl and zero row where kth row of l will fill in ***
row(k) = 0
i1 = 0
if (jrl(k) .eq. 0) go to 3
i = jrl(k)
2 i2 = jrl(i)
jrl(i) = i1
i1 = i
row(i) = 0
i = i2
if (i .ne. 0) go to 2
c ****** set row to zero where u will fill in ***********************
3 jmin = iju(k)
jmax = jmin + iu(k+1) - iu(k) - 1
if (jmin .gt. jmax) go to 5
do 4 j=jmin,jmax
4 row(ju(j)) = 0
c ****** place kth row of a in row **********************************
5 rk = r(k)
jmin = ia(rk)
jmax = ia(rk+1) - 1
do 6 j=jmin,jmax
row(ic(ja(j))) = a(j)
6 continue
c ****** initialize sum, and link through jrl ***********************
sum = b(rk)
i = i1
if (i .eq. 0) go to 10
c ****** assign the kth row of l and adjust row, sum ****************
7 lki = -row(i)
c ****** if l is not required, then comment out the following line **
l(irl(i)) = -lki
sum = sum + lki * tmp(i)
jmin = iu(i)
jmax = iu(i+1) - 1
if (jmin .gt. jmax) go to 9
mu = iju(i) - jmin
do 8 j=jmin,jmax
8 row(ju(mu+j)) = row(ju(mu+j)) + lki * u(j)
9 i = jrl(i)
if (i .ne. 0) go to 7
c
c ****** assign kth row of u and diagonal d, set tmp(k) *************
10 if (row(k) .eq. 0.0d0) go to 108
dk = 1.0d0 / row(k)
d(k) = dk
tmp(k) = sum * dk
if (k .eq. n) go to 19
jmin = iu(k)
jmax = iu(k+1) - 1
if (jmin .gt. jmax) go to 12
mu = iju(k) - jmin
do 11 j=jmin,jmax
11 u(j) = row(ju(mu+j)) * dk
12 continue
c
c ****** update irl and jrl, keeping jrl in decreasing order ********
i = i1
if (i .eq. 0) go to 18
14 irl(i) = irl(i) + 1
i1 = jrl(i)
if (irl(i) .ge. il(i+1)) go to 17
ijlb = irl(i) - il(i) + ijl(i)
j = jl(ijlb)
15 if (i .gt. jrl(j)) go to 16
j = jrl(j)
go to 15
16 jrl(i) = jrl(j)
jrl(j) = i
17 i = i1
if (i .ne. 0) go to 14
18 if (irl(k) .ge. il(k+1)) go to 19
j = jl(ijl(k))
jrl(k) = jrl(j)
jrl(j) = k
19 continue
c
c ****** solve ux = tmp by back substitution **********************
k = n
do 22 i=1,n
sum = tmp(k)
jmin = iu(k)
jmax = iu(k+1) - 1
if (jmin .gt. jmax) go to 21
mu = iju(k) - jmin
do 20 j=jmin,jmax
20 sum = sum - u(j) * tmp(ju(mu+j))
21 tmp(k) = sum
z(c(k)) = sum
22 k = k-1
flag = 0
return
c
c ** error.. insufficient storage for l
104 flag = 4*n + 1
return
c ** error.. insufficient storage for u
107 flag = 7*n + 1
return
c ** error.. zero pivot
108 flag = 8*n + k
return
end
|
