subroutine nroc (n, ic, ia, ja, a, jar, ar, p, flag) clll. optimize 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(1), ia(1), ja(1), jar(1), p(1), flag double precision a(1), ar(1) 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