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