subroutine nroc (n, ic, ia, ja, a, jar, ar, p, flag)
c
c ----------------------------------------------------------------
c
c yale sparse matrix package - nonsymmetric codes
c solving the system of equations mx = b
c
c i. calling sequences
c the coefficient matrix can be processed by an ordering routine
c (e.g., to reduce fillin or ensure numerical stability) before using
c the remaining subroutines. if no reordering is done, then set
c r(i) = c(i) = ic(i) = i for i=1,...,n. if an ordering subroutine
c is used, then nroc should be used to reorder the coefficient matrix
c the calling sequence is --
c ( (matrix ordering))
c (nroc (matrix reordering))
c nsfc (symbolic factorization to determine where fillin will
c occur during numeric factorization)
c nnfc (numeric factorization into product ldu of unit lower
c triangular matrix l, diagonal matrix d, and unit
c upper triangular matrix u, and solution of linear
c system)
c nnsc (solution of linear system for additional right-hand
c side using ldu factorization from nnfc)
c (if only one system of equations is to be solved, then the
c subroutine trk should be used.)
c
c ii. storage of sparse matrices
c the nonzero entries of the coefficient matrix m are stored
c row-by-row in the array a. to identify the individual nonzero
c entries in each row, we need to know in which column each entry
c lies. the column indices which correspond to the nonzero entries
c of m are stored in the array ja. i.e., if a(k) = m(i,j), then
c ja(k) = j. in addition, we need to know where each row starts and
c how long it is. the index positions in ja and a where the rows of
c m begin are stored in the array ia. i.e., if m(i,j) is the first
c (leftmost) entry in the i-th row and a(k) = m(i,j), then
c ia(i) = k. moreover, the index in ja and a of the first location
c following the last element in the last row is stored in ia(n+1).
c thus, the number of entries in the i-th row is given by
c ia(i+1) - ia(i), the nonzero entries of the i-th row are stored
c consecutively in
c a(ia(i)), a(ia(i)+1), ..., a(ia(i+1)-1),
c and the corresponding column indices are stored consecutively in
c ja(ia(i)), ja(ia(i)+1), ..., ja(ia(i+1)-1).
c for example, the 5 by 5 matrix
c ( 1. 0. 2. 0. 0.)
c ( 0. 3. 0. 0. 0.)
c m = ( 0. 4. 5. 6. 0.)
c ( 0. 0. 0. 7. 0.)
c ( 0. 0. 0. 8. 9.)
c would be stored as
c - 1 2 3 4 5 6 7 8 9
c ---+--------------------------
c ia - 1 3 4 7 8 10
c ja - 1 3 2 2 3 4 4 4 5
c a - 1. 2. 3. 4. 5. 6. 7. 8. 9. .
c
c the strict upper (lower) triangular portion of the matrix
c u (l) is stored in a similar fashion using the arrays iu, ju, u
c (il, jl, l) except that an additional array iju (ijl) is used to
c compress storage of ju (jl) by allowing some sequences of column
c (row) indices to used for more than one row (column) (n.b., l is
c stored by columns). iju(k) (ijl(k)) points to the starting
c location in ju (jl) of entries for the kth row (column).
c compression in ju (jl) occurs in two ways. first, if a row
c (column) i was merged into the current row (column) k, and the
c number of elements merged in from (the tail portion of) row
c (column) i is the same as the final length of row (column) k, then
c the kth row (column) and the tail of row (column) i are identical
c and iju(k) (ijl(k)) points to the start of the tail. second, if
c some tail portion of the (k-1)st row (column) is identical to the
c head of the kth row (column), then iju(k) (ijl(k)) points to the
c start of that tail portion. for example, the nonzero structure of
c the strict upper triangular part of the matrix
c d 0 x x x
c 0 d 0 x x
c 0 0 d x 0
c 0 0 0 d x
c 0 0 0 0 d
c would be represented as
c - 1 2 3 4 5 6
c ----+------------
c iu - 1 4 6 7 8 8
c ju - 3 4 5 4
c iju - 1 2 4 3 .
c the diagonal entries of l and u are assumed to be equal to one and
c are not stored. the array d contains the reciprocals of the
c diagonal entries of the matrix d.
c
c iii. additional storage savings
c in nsfc, r and ic can be the same array in the calling
c sequence if no reordering of the coefficient matrix has been done.
c in nnfc, r, c, and ic can all be the same array if no
c reordering has been done. if only the rows have been reordered,
c then c and ic can be the same array. if the row and column
c orderings are the same, then r and c can be the same array. z and
c row can be the same array.
c in nnsc or nntc, r and c can be the same array if no
c reordering has been done or if the row and column orderings are the
c same. z and b can be the same array. however, then b will be
c destroyed.
c
c iv. parameters
c following is a list of parameters to the programs. names are
c uniform among the various subroutines. class abbreviations are --
c n - integer variable
c f - real variable
c v - supplies a value to a subroutine
c r - returns a result from a subroutine
c i - used internally by a subroutine
c a - array
c
c class - parameter
c ------+----------
c fva - a - nonzero entries of the coefficient matrix m, stored
c - by rows.
c - size = number of nonzero entries in m.
c fva - b - right-hand side b.
c - size = n.
c nva - c - ordering of the columns of m.
c - size = n.
c fvra - d - reciprocals of the diagonal entries of the matrix d.
c - size = n.
c nr - flag - error flag. values and their meanings are --
c - 0 no errors detected
c - n+k null row in a -- row = k
c - 2n+k duplicate entry in a -- row = k
c - 3n+k insufficient storage for jl -- row = k
c - 4n+1 insufficient storage for l
c - 5n+k null pivot -- row = k
c - 6n+k insufficient storage for ju -- row = k
c - 7n+1 insufficient storage for u
c - 8n+k zero pivot -- row = k
c nva - ia - pointers to delimit the rows of a.
c - size = n+1.
c nvra - ijl - pointers to the first element in each column in jl,
c - used to compress storage in jl.
c - size = n.
c nvra - iju - pointers to the first element in each row in ju, used
c - to compress storage in ju.
c - size = n.
c nvra - il - pointers to delimit the columns of l.
c - size = n+1.
c nvra - iu - pointers to delimit the rows of u.
c - size = n+1.
c nva - ja - column numbers corresponding to the elements of a.
c - size = size of a.
c nvra - jl - row numbers corresponding to the elements of l.
c - size = jlmax.
c nv - jlmax - declared dimension of jl. jlmax must be larger than
c - the number of nonzeros in the strict lower triangle
c - of m plus fillin minus compression.
c nvra - ju - column numbers corresponding to the elements of u.
c - size = jumax.
c nv - jumax - declared dimension of ju. jumax must be larger than
c - the number of nonzeros in the strict upper triangle
c - of m plus fillin minus compression.
c fvra - l - nonzero entries in the strict lower triangular portion
c - of the matrix l, stored by columns.
c - size = lmax.
c nv - lmax - declared dimension of l. lmax must be larger than
c - the number of nonzeros in the strict lower triangle
c - of m plus fillin (il(n+1)-1 after nsfc).
c nv - n - number of variables/equations.
c nva - r - ordering of the rows of m.
c - size = n.
c fvra - u - nonzero entries in the strict upper triangular portion
c - of the matrix u, stored by rows.
c - size = umax.
c nv - umax - declared dimension of u. umax must be larger than
c - the number of nonzeros in the strict upper triangle
c - of m plus fillin (iu(n+1)-1 after nsfc).
c fra - z - solution x.
c - size = n.
c
c ----------------------------------------------------------------
c
c*** subroutine nroc
c*** reorders rows of a, leaving row order unchanged
c
c
c input parameters.. n, ic, ia, ja, a
c output parameters.. ja, a, flag
c
c parameters used internally..
c nia - p - at the kth step, p is a linked list of the reordered
c - column indices of the kth row of a. p(n+1) points
c - to the first entry in the list.
c - size = n+1.
c nia - jar - at the kth step,jar contains the elements of the
c - reordered column indices of a.
c - size = n.
c fia - ar - at the kth step, ar contains the elements of the
c - reordered row of a.
c - size = n.
c
integer ic(*), ia(*), ja(*), jar(*), p(*), flag
c real a(*), ar(*)
double precision a(*), ar(*)
c
c ****** for each nonempty row *******************************
do 5 k=1,n
jmin = ia(k)
jmax = ia(k+1) - 1
if(jmin .gt. jmax) go to 5
p(n+1) = n + 1
c ****** insert each element in the list *********************
do 3 j=jmin,jmax
newj = ic(ja(j))
i = n + 1
1 if(p(i) .ge. newj) go to 2
i = p(i)
go to 1
2 if(p(i) .eq. newj) go to 102
p(newj) = p(i)
p(i) = newj
jar(newj) = ja(j)
ar(newj) = a(j)
3 continue
c ****** replace old row in ja and a *************************
i = n + 1
do 4 j=jmin,jmax
i = p(i)
ja(j) = jar(i)
4 a(j) = ar(i)
5 continue
flag = 0
return
c
c ** error.. duplicate entry in a
102 flag = n + k
return
end