1
subroutine nroc (n, ic, ia, ja, a, jar, ar, p, flag)
4
c ----------------------------------------------------------------
6
c yale sparse matrix package - nonsymmetric codes
7
c solving the system of equations mx = b
10
c the coefficient matrix can be processed by an ordering routine
11
c (e.g., to reduce fillin or ensure numerical stability) before using
12
c the remaining subroutines. if no reordering is done, then set
13
c r(i) = c(i) = ic(i) = i for i=1,...,n. if an ordering subroutine
14
c is used, then nroc should be used to reorder the coefficient matrix
15
c the calling sequence is --
16
c ( (matrix ordering))
17
c (nroc (matrix reordering))
18
c nsfc (symbolic factorization to determine where fillin will
19
c occur during numeric factorization)
20
c nnfc (numeric factorization into product ldu of unit lower
21
c triangular matrix l, diagonal matrix d, and unit
22
c upper triangular matrix u, and solution of linear
24
c nnsc (solution of linear system for additional right-hand
25
c side using ldu factorization from nnfc)
26
c (if only one system of equations is to be solved, then the
27
c subroutine trk should be used.)
29
c ii. storage of sparse matrices
30
c the nonzero entries of the coefficient matrix m are stored
31
c row-by-row in the array a. to identify the individual nonzero
32
c entries in each row, we need to know in which column each entry
33
c lies. the column indices which correspond to the nonzero entries
34
c of m are stored in the array ja. i.e., if a(k) = m(i,j), then
35
c ja(k) = j. in addition, we need to know where each row starts and
36
c how long it is. the index positions in ja and a where the rows of
37
c m begin are stored in the array ia. i.e., if m(i,j) is the first
38
c (leftmost) entry in the i-th row and a(k) = m(i,j), then
39
c ia(i) = k. moreover, the index in ja and a of the first location
40
c following the last element in the last row is stored in ia(n+1).
41
c thus, the number of entries in the i-th row is given by
42
c ia(i+1) - ia(i), the nonzero entries of the i-th row are stored
44
c a(ia(i)), a(ia(i)+1), ..., a(ia(i+1)-1),
45
c and the corresponding column indices are stored consecutively in
46
c ja(ia(i)), ja(ia(i)+1), ..., ja(ia(i+1)-1).
47
c for example, the 5 by 5 matrix
50
c m = ( 0. 4. 5. 6. 0.)
55
c ---+--------------------------
57
c ja - 1 3 2 2 3 4 4 4 5
58
c a - 1. 2. 3. 4. 5. 6. 7. 8. 9. .
60
c the strict upper (lower) triangular portion of the matrix
61
c u (l) is stored in a similar fashion using the arrays iu, ju, u
62
c (il, jl, l) except that an additional array iju (ijl) is used to
63
c compress storage of ju (jl) by allowing some sequences of column
64
c (row) indices to used for more than one row (column) (n.b., l is
65
c stored by columns). iju(k) (ijl(k)) points to the starting
66
c location in ju (jl) of entries for the kth row (column).
67
c compression in ju (jl) occurs in two ways. first, if a row
68
c (column) i was merged into the current row (column) k, and the
69
c number of elements merged in from (the tail portion of) row
70
c (column) i is the same as the final length of row (column) k, then
71
c the kth row (column) and the tail of row (column) i are identical
72
c and iju(k) (ijl(k)) points to the start of the tail. second, if
73
c some tail portion of the (k-1)st row (column) is identical to the
74
c head of the kth row (column), then iju(k) (ijl(k)) points to the
75
c start of that tail portion. for example, the nonzero structure of
76
c the strict upper triangular part of the matrix
82
c would be represented as
88
c the diagonal entries of l and u are assumed to be equal to one and
89
c are not stored. the array d contains the reciprocals of the
90
c diagonal entries of the matrix d.
92
c iii. additional storage savings
93
c in nsfc, r and ic can be the same array in the calling
94
c sequence if no reordering of the coefficient matrix has been done.
95
c in nnfc, r, c, and ic can all be the same array if no
96
c reordering has been done. if only the rows have been reordered,
97
c then c and ic can be the same array. if the row and column
98
c orderings are the same, then r and c can be the same array. z and
99
c row can be the same array.
100
c in nnsc or nntc, r and c can be the same array if no
101
c reordering has been done or if the row and column orderings are the
102
c same. z and b can be the same array. however, then b will be
106
c following is a list of parameters to the programs. names are
107
c uniform among the various subroutines. class abbreviations are --
108
c n - integer variable
110
c v - supplies a value to a subroutine
111
c r - returns a result from a subroutine
112
c i - used internally by a subroutine
117
c fva - a - nonzero entries of the coefficient matrix m, stored
119
c - size = number of nonzero entries in m.
120
c fva - b - right-hand side b.
122
c nva - c - ordering of the columns of m.
124
c fvra - d - reciprocals of the diagonal entries of the matrix d.
126
c nr - flag - error flag. values and their meanings are --
127
c - 0 no errors detected
128
c - n+k null row in a -- row = k
129
c - 2n+k duplicate entry in a -- row = k
130
c - 3n+k insufficient storage for jl -- row = k
131
c - 4n+1 insufficient storage for l
132
c - 5n+k null pivot -- row = k
133
c - 6n+k insufficient storage for ju -- row = k
134
c - 7n+1 insufficient storage for u
135
c - 8n+k zero pivot -- row = k
136
c nva - ia - pointers to delimit the rows of a.
138
c nvra - ijl - pointers to the first element in each column in jl,
139
c - used to compress storage in jl.
141
c nvra - iju - pointers to the first element in each row in ju, used
142
c - to compress storage in ju.
144
c nvra - il - pointers to delimit the columns of l.
146
c nvra - iu - pointers to delimit the rows of u.
148
c nva - ja - column numbers corresponding to the elements of a.
149
c - size = size of a.
150
c nvra - jl - row numbers corresponding to the elements of l.
152
c nv - jlmax - declared dimension of jl. jlmax must be larger than
153
c - the number of nonzeros in the strict lower triangle
154
c - of m plus fillin minus compression.
155
c nvra - ju - column numbers corresponding to the elements of u.
157
c nv - jumax - declared dimension of ju. jumax must be larger than
158
c - the number of nonzeros in the strict upper triangle
159
c - of m plus fillin minus compression.
160
c fvra - l - nonzero entries in the strict lower triangular portion
161
c - of the matrix l, stored by columns.
163
c nv - lmax - declared dimension of l. lmax must be larger than
164
c - the number of nonzeros in the strict lower triangle
165
c - of m plus fillin (il(n+1)-1 after nsfc).
166
c nv - n - number of variables/equations.
167
c nva - r - ordering of the rows of m.
169
c fvra - u - nonzero entries in the strict upper triangular portion
170
c - of the matrix u, stored by rows.
172
c nv - umax - declared dimension of u. umax must be larger than
173
c - the number of nonzeros in the strict upper triangle
174
c - of m plus fillin (iu(n+1)-1 after nsfc).
175
c fra - z - solution x.
178
c ----------------------------------------------------------------
181
c*** reorders rows of a, leaving row order unchanged
184
c input parameters.. n, ic, ia, ja, a
185
c output parameters.. ja, a, flag
187
c parameters used internally..
188
c nia - p - at the kth step, p is a linked list of the reordered
189
c - column indices of the kth row of a. p(n+1) points
190
c - to the first entry in the list.
192
c nia - jar - at the kth step,jar contains the elements of the
193
c - reordered column indices of a.
195
c fia - ar - at the kth step, ar contains the elements of the
196
c - reordered row of a.
199
integer ic(1), ia(1), ja(1), jar(1), p(1), flag
200
double precision a(1), ar(1)
202
c ****** for each nonempty row *******************************
206
if(jmin .gt. jmax) go to 5
208
c ****** insert each element in the list *********************
212
1 if(p(i) .ge. newj) go to 2
215
2 if(p(i) .eq. newj) go to 102
221
c ****** replace old row in ja and a *************************
231
c ** error.. duplicate entry in a
added
removed
Lib/integrate/odepack/nroc.f
