subroutine preord(n, isym, ja, a, list, index, ierr)
c
c the routines fillin, factor, and solve in this package are
c prototype routines based on sparse elimination algorithms
c and data structures described in:
c
c General Sparse Elimination Requires No Permanent Integer Storage
c by Randolph E. Bank and R. Kent Smith, vol. 8, pp.574-584, SIAM
c Journal of Scientific and Statistical Computing, 1987.
c
c a more complete discussion of these procedures can be found there.
c to keep the code as simple as possible, we have assumed that the
c identity ordering is used. if a reordering algorithm such as
c minimum degree is to be used, one must either physically
c reorder the matrix, or systematically modify references
c to the arrays ja, a, x, and b to reflect the permutation matrix.
c this is a shortcoming of the code, which was written for mainly
c for illustrative purposes, and not of the underlying algorithms.
c
c fillin and factor assume a particular order for the row indices
c in ja. this routine is included to physically reorder ja and
c a in the required way. it is not an essential part of the package,
c and should be called only if necessary.
c
integer ja(*), list(*), index(*)
real a(*), ltemp, utemp
c
c put row indices in ja in decreasing order
c
c input: n, isym, ja, a
c output: ja, a, ierr
c workspace: list(n), index(n)
c
c remarks:
c
c 1. isym = 0 means nonsymmetric storage is used for a
c isym = 1 means symmetric storage is used for a
c
c 2. ierr = 0 is a normal return
c ierr = i means a row index greater than i was
c encountered in column i
c
c 3. if a is symmetric, the reordering of the lower
c triangle in the do 40 loop is redundant
c
ierr = 0
c
c compute offset for a
c
if (isym .eq. 1) then
lshift = 0
else
lshift = ja(n+1) - ja(1)
end if
c
c the main loop
c
do 50 i = 1, n
if (ja(i)+1 .ge. ja(i+1)) go to 50
c
c sort seed indices in increasing order
c
list(i) = i
do 30 indexj = ja(i), ja(i+1)-1
j = ja(indexj)
if (j .ge. i) go to 100
index(j) = indexj
last = i
k = list(i)
10 if (j .lt. k) go to 20
last = k
k = list(k)
go to 10
20 list(last) = j
list(j) = k
30 continue
c
c in place reordering of ja and a
c
k = list(i)
do 40 indexj = ja(i+1)-1, ja(i), -1
j = ja(indexj)
if (j .eq. k) go to 40
c
indexk = index(k)
ja(indexj) = k
ja(indexk) = j
index(j) = indexk
c
ltemp = a(indexj+lshift)
utemp = a(indexj)
a(indexj+lshift) = a(indexk+lshift)
a(indexj) = a(indexk)
a(indexk+lshift) = ltemp
a(indexk) = utemp
c
40 k = list(k)
50 continue
c
return
c
c error return
c
100 ierr = i
return
end
subroutine fillin (n, ja, jl, list)
integer ja(*), jl(*), list(*), uptr
c
c compute jl
c
c input: n, ja
c output: jl
c workspace: list(n)
c
c remarks:
c
c 1. for convenience, jl(n) is the size of the strict
c upper triangular part of the factored matrix
c
c 2. the row indices for each column in ja are assumed
c to be in decreasing order
c
c initialize
c
do 10 i = 1, n
jl(i) = 0
10 list(i) = 0
c
c the main loop
c
jl(n+1) = n + 2
do 50 i = 1, n
if (ja(i) .ge. ja(i+1)) go to 50
c
c loop over seed indices in decreasing order
c
list(i) = i
length = 1
do 30 iseed = ja(i), ja(i+1)-1
k = ja(iseed)
uptr = jl(iseed-1)
c
c add a new entry to list
c
20 list(k) = list(i)
list(i) = k
uptr = uptr + 1
length = length + 1
if (jl(k) .eq. 0) jl(k) = i
k = jl(k)
if (list(k) .eq. 0) go to 20
c
30 jl(iseed) = uptr
c
c clean up loop for list array
c
k = i
do 40 j = 1, length
ksave = k
k = list(k)
40 list(ksave) = 0
c
50 continue
c
c compute size of upper triangle
c
lenjl = n + 1 + ja(n+1) - ja(1)
jl(n) = jl(lenjl) - jl(n+1)
c
return
end
subroutine factor(n,isym,ja,a,jl,u,index,expres,ierr)
integer ja(*), jl(*), index(*), ashift, ushift ,
1 expres(*), expk
real a(*), u(*), dsum, usum, lsum
c
c compute u
c
c input: n, isym, ja, a, jl
c output: u, ierr
c workspace: index(n), expres(n)
c
c remarks:
c
c 1. isym = 0 means nonsymmetric storage for a and u
c isym = 1 means symmetric storage for a and u
c
c 2. ierr = 0 is a normal return
c ierr = i means the i-th diagonal element of the
c factored matrix was zero
c
c 3. for symmetric matrices, the computation of lsum in the
c do 60 loop, and the computation of u(indexk+ushift) in
c the do 80 loop are redundant
c
c 4. the express vector modification was added in September,
c 1986.
c
ierr = 0
c
c compute offsets for a and u
c
if (isym .eq. 1) then
ashift = 0
ushift = 0
lenu = n + 1 + jl(n)
else
ashift = ja(n+1) - ja(1)
ushift = jl(n)
lenu = n + 1 + 2 * jl(n)
end if
c
c initialize
c
do 10 i = 1, lenu
10 u(i) = 0.0e0
expk = n
do 20 i = n, 1, -1
index(i) = 0
if (jl(i) .eq. i+1) then
expres(i) = expk
else
expres(i) = 0
expk = i
end if
20 continue
expres(n) = 0
c
c the main loop
c
do 100 i = 1, n
dsum = 0.0
if (ja(i) .ge. ja(i+1)) go to 90
c
c loop over seed indices for column i in increasing order
c
do 60 iseed = ja(i+1)-1, ja(i), -1
k = ja(iseed)
c
c move off diagonal entries from a to u
c
u(jl(iseed-1)) = a(iseed)
u(jl(iseed-1)+ushift) = a(iseed+ashift)
c
do 50 indexk = jl(iseed-1), jl(iseed)-1
index(k) = indexk
expk = expres(k)
if (expk .gt. 0) then
if (expres(expk) .eq. 0) expres(expk) = -k
end if
if (ja(k) .ge. ja(k+1)) go to 50
c
c loop over seed indices for column k in decreasing order
c
lsum = u(indexk+ushift)
usum = u(indexk)
c
do 40 kseed = ja(k), ja(k+1)-1
j = ja(kseed)
jstart = jl(kseed-1)
c
c find start index for inner product loop
c
25 if (jstart .ge. jl(kseed)) go to 40
if (index(j) .eq. 0) then
if (expres(j) .gt. 0) then
jold = j
j = expres(j)
if (expres(j) .lt. 0) j = -expres(j)
jstart = jstart + j - jold
else
j = jl(j)
jstart = jstart + 1
end if
go to 25
end if
c
c form inner products
c
do 30 indexj = jstart, jl(kseed)-1
lsum = lsum - u(indexj) * u(index(j)+ushift)
usum = usum - u(indexj+ushift) * u(index(j))
30 j = jl(j)
40 continue
c
c update u(indexk) and u(indexk+ushift)
c
u(indexk+ushift) = lsum
u(indexk) = usum
50 k = jl(k)
60 continue
c
c clean up loop for index array, compute diagonal
c
do 80 iseed = ja(i), ja(i+1)-1
k = ja(iseed)
do 70 indexk = jl(iseed-1), jl(iseed)-1
c
usave = u(indexk+ushift)
u(indexk) = u(indexk) * u(k)
u(indexk+ushift) = usave * u(k)
dsum = dsum + u(indexk) * usave
c
index(k) = 0
expk = expres(k)
if (expk .gt. 0) expres(expk) = 0
c
70 k = jl(k)
80 continue
c
c compute the diagonal
c
90 u(i) = a(i) - dsum
if (u(i) .eq. 0.0e0) go to 200
u(i) = 1.0e0 / u(i)
100 continue
c
return
c
c error return
c
200 ierr = i
return
end
subroutine solve (n, isym, ja, jl, u, x, b)
integer ja(*), jl(*), ushift, lshift
real u(*), x(*), b(*)
c
c compute the solution x of a * x = b
c
c input: n, isym, ja, jl, u, b
c output: x
c
c remark: isym = 1 means symmetric storage is used for u
c isym = 0 means nonsymmetric storage is used for u
c isym = -1 means nonsymmetric storage is used, and the
c problem is to be solved with the transpose of a
c
c compute offsets for u
c
if (isym .eq. 0) then
lshift = jl(n)
else
lshift = 0
end if
if (isym .eq. -1) then
ushift = jl(n)
else
ushift = 0
end if
c
c lower triangular system
c
do 30 i = 1, n
sum = 0.0e0
if (ja(i) .ge. ja(i+1)) go to 30
do 20 iseed = ja(i), ja(i+1)-1
k = ja(iseed)
do 10 j = jl(iseed-1), jl(iseed)-1
sum = sum + u(j+lshift) * x(k)
10 k = jl(k)
20 continue
30 x(i) = b(i) - sum
c
c diagonal system
c
do 40 i = 1, n
40 x(i) = u(i) * x(i)
c
c upper triangular system
c
do 70 i = n, 1, -1
if (ja(i) .ge. ja(i+1)) go to 70
do 60 iseed = ja(i), ja(i+1)-1
k = ja(iseed)
do 50 j = jl(iseed-1), jl(iseed)-1
x(k) = x(k) - u(j+ushift) * x(i)
50 k = jl(k)
60 continue
70 continue
c
return
end