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)
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(*), c(*), ic(*), ia(*), ja(*), il(*), jl(*), ijl(*)
integer iu(*), ju(*), iju(*), irl(*), jrl(*), flag
c real a(*), l(*), d(*), u(*), z(*), b(*), row(*)
c real tmp(*), lki, sum, dk
double precision a(*), l(*), d(*), u(*), z(*), b(*), row(*)
double precision tmp(*), 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