|
5164
|
1 C----------------------------------------------------------------------- |
|
|
2 C AMDBAR: approximate minimum degree, without aggressive absorption |
|
|
3 C----------------------------------------------------------------------- |
|
|
4 |
|
|
5 SUBROUTINE AMDBAR |
|
|
6 $ (N, PE, IW, LEN, IWLEN, PFREE, NV, NEXT, |
|
|
7 $ LAST, HEAD, ELEN, DEGREE, NCMPA, W) |
|
|
8 |
|
|
9 INTEGER N, IWLEN, PFREE, NCMPA, IW (IWLEN), PE (N), |
|
|
10 $ DEGREE (N), NV (N), NEXT (N), LAST (N), HEAD (N), |
|
|
11 $ ELEN (N), W (N), LEN (N) |
|
|
12 |
|
|
13 C Given a representation of the nonzero pattern of a symmetric matrix, |
|
|
14 C A, (excluding the diagonal) perform an approximate minimum |
|
|
15 C (UMFPACK/MA38-style) degree ordering to compute a pivot order |
|
|
16 C such that the introduction of nonzeros (fill-in) in the Cholesky |
|
|
17 C factors A = LL^T are kept low. At each step, the pivot |
|
|
18 C selected is the one with the minimum UMFPACK/MA38-style |
|
|
19 C upper-bound on the external degree. |
|
|
20 C |
|
|
21 C This routine does not do aggresive absorption (as done by AMD). |
|
|
22 |
|
|
23 C ********************************************************************** |
|
|
24 C ***** CAUTION: ARGUMENTS ARE NOT CHECKED FOR ERRORS ON INPUT. ****** |
|
|
25 C ********************************************************************** |
|
|
26 |
|
|
27 C References: |
|
|
28 C |
|
|
29 C [1] Timothy A. Davis and Iain Duff, "An unsymmetric-pattern |
|
|
30 C multifrontal method for sparse LU factorization", SIAM J. |
|
|
31 C Matrix Analysis and Applications, vol. 18, no. 1, pp. |
|
|
32 C 140-158. Discusses UMFPACK / MA38, which first introduced |
|
|
33 C the approximate minimum degree used by this routine. |
|
|
34 C |
|
|
35 C [2] Patrick Amestoy, Timothy A. Davis, and Iain S. Duff, "An |
|
|
36 C approximate degree ordering algorithm," SIAM J. Matrix |
|
|
37 C Analysis and Applications, vol. 17, no. 4, pp. 886-905, |
|
|
38 C 1996. Discusses AMD, AMDBAR, and MC47B. |
|
|
39 C |
|
|
40 C [3] Alan George and Joseph Liu, "The evolution of the minimum |
|
|
41 C degree ordering algorithm," SIAM Review, vol. 31, no. 1, |
|
|
42 C pp. 1-19, 1989. We list below the features mentioned in |
|
|
43 C that paper that this code includes: |
|
|
44 C |
|
|
45 C mass elimination: |
|
|
46 C Yes. MA27 relied on supervariable detection for mass |
|
|
47 C elimination. |
|
|
48 C indistinguishable nodes: |
|
|
49 C Yes (we call these "supervariables"). This was also in |
|
|
50 C the MA27 code - although we modified the method of |
|
|
51 C detecting them (the previous hash was the true degree, |
|
|
52 C which we no longer keep track of). A supervariable is |
|
|
53 C a set of rows with identical nonzero pattern. All |
|
|
54 C variables in a supervariable are eliminated together. |
|
|
55 C Each supervariable has as its numerical name that of |
|
|
56 C one of its variables (its principal variable). |
|
|
57 C quotient graph representation: |
|
|
58 C Yes. We use the term "element" for the cliques formed |
|
|
59 C during elimination. This was also in the MA27 code. |
|
|
60 C The algorithm can operate in place, but it will work |
|
|
61 C more efficiently if given some "elbow room." |
|
|
62 C element absorption: |
|
|
63 C Yes. This was also in the MA27 code. |
|
|
64 C external degree: |
|
|
65 C Yes. The MA27 code was based on the true degree. |
|
|
66 C incomplete degree update and multiple elimination: |
|
|
67 C No. This was not in MA27, either. Our method of |
|
|
68 C degree update within MC47B/BD is element-based, not |
|
|
69 C variable-based. It is thus not well-suited for use |
|
|
70 C with incomplete degree update or multiple elimination. |
|
|
71 |
|
|
72 C----------------------------------------------------------------------- |
|
|
73 C Authors, and Copyright (C) 1995 by: |
|
|
74 C Timothy A. Davis, Patrick Amestoy, Iain S. Duff, & John K. Reid. |
|
|
75 C |
|
|
76 C Acknowledgements: |
|
|
77 C This work (and the UMFPACK package) was supported by the |
|
|
78 C National Science Foundation (ASC-9111263 and DMS-9223088). |
|
|
79 C The UMFPACK/MA38 approximate degree update algorithm, the |
|
|
80 C unsymmetric analog which forms the basis of MC47B/BD, was |
|
|
81 C developed while Tim Davis was supported by CERFACS (Toulouse, |
|
|
82 C France) in a post-doctoral position. |
|
|
83 C |
|
|
84 C Date: September, 1995 |
|
|
85 C----------------------------------------------------------------------- |
|
|
86 |
|
|
87 C----------------------------------------------------------------------- |
|
|
88 C INPUT ARGUMENTS (unaltered): |
|
|
89 C----------------------------------------------------------------------- |
|
|
90 |
|
|
91 C n: The matrix order. |
|
|
92 C |
|
|
93 C Restriction: 1 .le. n .lt. (iovflo/2)-2, where iovflo is |
|
|
94 C the largest positive integer that your computer can represent. |
|
|
95 |
|
|
96 C iwlen: The length of iw (1..iwlen). On input, the matrix is |
|
|
97 C stored in iw (1..pfree-1). However, iw (1..iwlen) should be |
|
|
98 C slightly larger than what is required to hold the matrix, at |
|
|
99 C least iwlen .ge. pfree + n is recommended. Otherwise, |
|
|
100 C excessive compressions will take place. |
|
|
101 C *** We do not recommend running this algorithm with *** |
|
|
102 C *** iwlen .lt. pfree + n. *** |
|
|
103 C *** Better performance will be obtained if *** |
|
|
104 C *** iwlen .ge. pfree + n *** |
|
|
105 C *** or better yet *** |
|
|
106 C *** iwlen .gt. 1.2 * pfree *** |
|
|
107 C *** (where pfree is its value on input). *** |
|
|
108 C The algorithm will not run at all if iwlen .lt. pfree-1. |
|
|
109 C |
|
|
110 C Restriction: iwlen .ge. pfree-1 |
|
|
111 |
|
|
112 C----------------------------------------------------------------------- |
|
|
113 C INPUT/OUPUT ARGUMENTS: |
|
|
114 C----------------------------------------------------------------------- |
|
|
115 |
|
|
116 C pe: On input, pe (i) is the index in iw of the start of row i, or |
|
|
117 C zero if row i has no off-diagonal non-zeros. |
|
|
118 C |
|
|
119 C During execution, it is used for both supervariables and |
|
|
120 C elements: |
|
|
121 C |
|
|
122 C * Principal supervariable i: index into iw of the |
|
|
123 C description of supervariable i. A supervariable |
|
|
124 C represents one or more rows of the matrix |
|
|
125 C with identical nonzero pattern. |
|
|
126 C * Non-principal supervariable i: if i has been absorbed |
|
|
127 C into another supervariable j, then pe (i) = -j. |
|
|
128 C That is, j has the same pattern as i. |
|
|
129 C Note that j might later be absorbed into another |
|
|
130 C supervariable j2, in which case pe (i) is still -j, |
|
|
131 C and pe (j) = -j2. |
|
|
132 C * Unabsorbed element e: the index into iw of the description |
|
|
133 C of element e, if e has not yet been absorbed by a |
|
|
134 C subsequent element. Element e is created when |
|
|
135 C the supervariable of the same name is selected as |
|
|
136 C the pivot. |
|
|
137 C * Absorbed element e: if element e is absorbed into element |
|
|
138 C e2, then pe (e) = -e2. This occurs when the pattern of |
|
|
139 C e (that is, Le) is found to be a subset of the pattern |
|
|
140 C of e2 (that is, Le2). If element e is "null" (it has |
|
|
141 C no nonzeros outside its pivot block), then pe (e) = 0. |
|
|
142 C |
|
|
143 C On output, pe holds the assembly tree/forest, which implicitly |
|
|
144 C represents a pivot order with identical fill-in as the actual |
|
|
145 C order (via a depth-first search of the tree). |
|
|
146 C |
|
|
147 C On output: |
|
|
148 C If nv (i) .gt. 0, then i represents a node in the assembly tree, |
|
|
149 C and the parent of i is -pe (i), or zero if i is a root. |
|
|
150 C If nv (i) = 0, then (i,-pe (i)) represents an edge in a |
|
|
151 C subtree, the root of which is a node in the assembly tree. |
|
|
152 |
|
|
153 C pfree: On input the tail end of the array, iw (pfree..iwlen), |
|
|
154 C is empty, and the matrix is stored in iw (1..pfree-1). |
|
|
155 C During execution, additional data is placed in iw, and pfree |
|
|
156 C is modified so that iw (pfree..iwlen) is always the unused part |
|
|
157 C of iw. On output, pfree is set equal to the size of iw that |
|
|
158 C would have been needed for no compressions to occur. If |
|
|
159 C ncmpa is zero, then pfree (on output) is less than or equal to |
|
|
160 C iwlen, and the space iw (pfree+1 ... iwlen) was not used. |
|
|
161 C Otherwise, pfree (on output) is greater than iwlen, and all the |
|
|
162 C memory in iw was used. |
|
|
163 |
|
|
164 C----------------------------------------------------------------------- |
|
|
165 C INPUT/MODIFIED (undefined on output): |
|
|
166 C----------------------------------------------------------------------- |
|
|
167 |
|
|
168 C len: On input, len (i) holds the number of entries in row i of the |
|
|
169 C matrix, excluding the diagonal. The contents of len (1..n) |
|
|
170 C are undefined on output. |
|
|
171 |
|
|
172 C iw: On input, iw (1..pfree-1) holds the description of each row i |
|
|
173 C in the matrix. The matrix must be symmetric, and both upper |
|
|
174 C and lower triangular parts must be present. The diagonal must |
|
|
175 C not be present. Row i is held as follows: |
|
|
176 C |
|
|
177 C len (i): the length of the row i data structure |
|
|
178 C iw (pe (i) ... pe (i) + len (i) - 1): |
|
|
179 C the list of column indices for nonzeros |
|
|
180 C in row i (simple supervariables), excluding |
|
|
181 C the diagonal. All supervariables start with |
|
|
182 C one row/column each (supervariable i is just |
|
|
183 C row i). |
|
|
184 C if len (i) is zero on input, then pe (i) is ignored |
|
|
185 C on input. |
|
|
186 C |
|
|
187 C Note that the rows need not be in any particular order, |
|
|
188 C and there may be empty space between the rows. |
|
|
189 C |
|
|
190 C During execution, the supervariable i experiences fill-in. |
|
|
191 C This is represented by placing in i a list of the elements |
|
|
192 C that cause fill-in in supervariable i: |
|
|
193 C |
|
|
194 C len (i): the length of supervariable i |
|
|
195 C iw (pe (i) ... pe (i) + elen (i) - 1): |
|
|
196 C the list of elements that contain i. This list |
|
|
197 C is kept short by removing absorbed elements. |
|
|
198 C iw (pe (i) + elen (i) ... pe (i) + len (i) - 1): |
|
|
199 C the list of supervariables in i. This list |
|
|
200 C is kept short by removing nonprincipal |
|
|
201 C variables, and any entry j that is also |
|
|
202 C contained in at least one of the elements |
|
|
203 C (j in Le) in the list for i (e in row i). |
|
|
204 C |
|
|
205 C When supervariable i is selected as pivot, we create an |
|
|
206 C element e of the same name (e=i): |
|
|
207 C |
|
|
208 C len (e): the length of element e |
|
|
209 C iw (pe (e) ... pe (e) + len (e) - 1): |
|
|
210 C the list of supervariables in element e. |
|
|
211 C |
|
|
212 C An element represents the fill-in that occurs when supervariable |
|
|
213 C i is selected as pivot (which represents the selection of row i |
|
|
214 C and all non-principal variables whose principal variable is i). |
|
|
215 C We use the term Le to denote the set of all supervariables |
|
|
216 C in element e. Absorbed supervariables and elements are pruned |
|
|
217 C from these lists when computationally convenient. |
|
|
218 C |
|
|
219 C CAUTION: THE INPUT MATRIX IS OVERWRITTEN DURING COMPUTATION. |
|
|
220 C The contents of iw are undefined on output. |
|
|
221 |
|
|
222 C----------------------------------------------------------------------- |
|
|
223 C OUTPUT (need not be set on input): |
|
|
224 C----------------------------------------------------------------------- |
|
|
225 |
|
|
226 C nv: During execution, abs (nv (i)) is equal to the number of rows |
|
|
227 C that are represented by the principal supervariable i. If i is |
|
|
228 C a nonprincipal variable, then nv (i) = 0. Initially, |
|
|
229 C nv (i) = 1 for all i. nv (i) .lt. 0 signifies that i is a |
|
|
230 C principal variable in the pattern Lme of the current pivot |
|
|
231 C element me. On output, nv (e) holds the true degree of element |
|
|
232 C e at the time it was created (including the diagonal part). |
|
|
233 |
|
|
234 C ncmpa: The number of times iw was compressed. If this is |
|
|
235 C excessive, then the execution took longer than what could have |
|
|
236 C been. To reduce ncmpa, try increasing iwlen to be 10% or 20% |
|
|
237 C larger than the value of pfree on input (or at least |
|
|
238 C iwlen .ge. pfree + n). The fastest performance will be |
|
|
239 C obtained when ncmpa is returned as zero. If iwlen is set to |
|
|
240 C the value returned by pfree on *output*, then no compressions |
|
|
241 C will occur. |
|
|
242 |
|
|
243 C elen: See the description of iw above. At the start of execution, |
|
|
244 C elen (i) is set to zero. During execution, elen (i) is the |
|
|
245 C number of elements in the list for supervariable i. When e |
|
|
246 C becomes an element, elen (e) = -nel is set, where nel is the |
|
|
247 C current step of factorization. elen (i) = 0 is done when i |
|
|
248 C becomes nonprincipal. |
|
|
249 C |
|
|
250 C For variables, elen (i) .ge. 0 holds until just before the |
|
|
251 C permutation vectors are computed. For elements, |
|
|
252 C elen (e) .lt. 0 holds. |
|
|
253 C |
|
|
254 C On output elen (1..n) holds the inverse permutation (the same |
|
|
255 C as the 'INVP' argument in Sparspak). That is, if k = elen (i), |
|
|
256 C then row i is the kth pivot row. Row i of A appears as the |
|
|
257 C (elen(i))-th row in the permuted matrix, PAP^T. |
|
|
258 |
|
|
259 C last: In a degree list, last (i) is the supervariable preceding i, |
|
|
260 C or zero if i is the head of the list. In a hash bucket, |
|
|
261 C last (i) is the hash key for i. last (head (hash)) is also |
|
|
262 C used as the head of a hash bucket if head (hash) contains a |
|
|
263 C degree list (see head, below). |
|
|
264 C |
|
|
265 C On output, last (1..n) holds the permutation (the same as the |
|
|
266 C 'PERM' argument in Sparspak). That is, if i = last (k), then |
|
|
267 C row i is the kth pivot row. Row last (k) of A is the k-th row |
|
|
268 C in the permuted matrix, PAP^T. |
|
|
269 |
|
|
270 C----------------------------------------------------------------------- |
|
|
271 C LOCAL (not input or output - used only during execution): |
|
|
272 C----------------------------------------------------------------------- |
|
|
273 |
|
|
274 C degree: If i is a supervariable, then degree (i) holds the |
|
|
275 C current approximation of the external degree of row i (an upper |
|
|
276 C bound). The external degree is the number of nonzeros in row i, |
|
|
277 C minus abs (nv (i)) (the diagonal part). The bound is equal to |
|
|
278 C the external degree if elen (i) is less than or equal to two. |
|
|
279 C |
|
|
280 C We also use the term "external degree" for elements e to refer |
|
|
281 C to |Le \ Lme|. If e is an element, then degree (e) holds |Le|, |
|
|
282 C which is the degree of the off-diagonal part of the element e |
|
|
283 C (not including the diagonal part). |
|
|
284 |
|
|
285 C head: head is used for degree lists. head (deg) is the first |
|
|
286 C supervariable in a degree list (all supervariables i in a |
|
|
287 C degree list deg have the same approximate degree, namely, |
|
|
288 C deg = degree (i)). If the list deg is empty then |
|
|
289 C head (deg) = 0. |
|
|
290 C |
|
|
291 C During supervariable detection head (hash) also serves as a |
|
|
292 C pointer to a hash bucket. |
|
|
293 C If head (hash) .gt. 0, there is a degree list of degree hash. |
|
|
294 C The hash bucket head pointer is last (head (hash)). |
|
|
295 C If head (hash) = 0, then the degree list and hash bucket are |
|
|
296 C both empty. |
|
|
297 C If head (hash) .lt. 0, then the degree list is empty, and |
|
|
298 C -head (hash) is the head of the hash bucket. |
|
|
299 C After supervariable detection is complete, all hash buckets |
|
|
300 C are empty, and the (last (head (hash)) = 0) condition is |
|
|
301 C restored for the non-empty degree lists. |
|
|
302 |
|
|
303 C next: next (i) is the supervariable following i in a link list, or |
|
|
304 C zero if i is the last in the list. Used for two kinds of |
|
|
305 C lists: degree lists and hash buckets (a supervariable can be |
|
|
306 C in only one kind of list at a time). |
|
|
307 |
|
|
308 C w: The flag array w determines the status of elements and |
|
|
309 C variables, and the external degree of elements. |
|
|
310 C |
|
|
311 C for elements: |
|
|
312 C if w (e) = 0, then the element e is absorbed |
|
|
313 C if w (e) .ge. wflg, then w (e) - wflg is the size of |
|
|
314 C the set |Le \ Lme|, in terms of nonzeros (the |
|
|
315 C sum of abs (nv (i)) for each principal variable i that |
|
|
316 C is both in the pattern of element e and NOT in the |
|
|
317 C pattern of the current pivot element, me). |
|
|
318 C if wflg .gt. w (e) .gt. 0, then e is not absorbed and has |
|
|
319 C not yet been seen in the scan of the element lists in |
|
|
320 C the computation of |Le\Lme| in loop 150 below. |
|
|
321 C |
|
|
322 C for variables: |
|
|
323 C during supervariable detection, if w (j) .ne. wflg then j is |
|
|
324 C not in the pattern of variable i |
|
|
325 C |
|
|
326 C The w array is initialized by setting w (i) = 1 for all i, |
|
|
327 C and by setting wflg = 2. It is reinitialized if wflg becomes |
|
|
328 C too large (to ensure that wflg+n does not cause integer |
|
|
329 C overflow). |
|
|
330 |
|
|
331 C----------------------------------------------------------------------- |
|
|
332 C LOCAL INTEGERS: |
|
|
333 C----------------------------------------------------------------------- |
|
|
334 |
|
|
335 INTEGER DEG, DEGME, DMAX, E, ELENME, ELN, HASH, HMOD, I, |
|
|
336 $ ILAST, INEXT, J, JLAST, JNEXT, K, KNT1, KNT2, KNT3, |
|
|
337 $ LENJ, LN, MAXMEM, ME, MEM, MINDEG, NEL, NEWMEM, |
|
|
338 $ NLEFT, NVI, NVJ, NVPIV, SLENME, WE, WFLG, WNVI, X |
|
|
339 |
|
|
340 C deg: the degree of a variable or element |
|
|
341 C degme: size, |Lme|, of the current element, me (= degree (me)) |
|
|
342 C dext: external degree, |Le \ Lme|, of some element e |
|
|
343 C dmax: largest |Le| seen so far |
|
|
344 C e: an element |
|
|
345 C elenme: the length, elen (me), of element list of pivotal var. |
|
|
346 C eln: the length, elen (...), of an element list |
|
|
347 C hash: the computed value of the hash function |
|
|
348 C hmod: the hash function is computed modulo hmod = max (1,n-1) |
|
|
349 C i: a supervariable |
|
|
350 C ilast: the entry in a link list preceding i |
|
|
351 C inext: the entry in a link list following i |
|
|
352 C j: a supervariable |
|
|
353 C jlast: the entry in a link list preceding j |
|
|
354 C jnext: the entry in a link list, or path, following j |
|
|
355 C k: the pivot order of an element or variable |
|
|
356 C knt1: loop counter used during element construction |
|
|
357 C knt2: loop counter used during element construction |
|
|
358 C knt3: loop counter used during compression |
|
|
359 C lenj: len (j) |
|
|
360 C ln: length of a supervariable list |
|
|
361 C maxmem: amount of memory needed for no compressions |
|
|
362 C me: current supervariable being eliminated, and the |
|
|
363 C current element created by eliminating that |
|
|
364 C supervariable |
|
|
365 C mem: memory in use assuming no compressions have occurred |
|
|
366 C mindeg: current minimum degree |
|
|
367 C nel: number of pivots selected so far |
|
|
368 C newmem: amount of new memory needed for current pivot element |
|
|
369 C nleft: n - nel, the number of nonpivotal rows/columns remaining |
|
|
370 C nvi: the number of variables in a supervariable i (= nv (i)) |
|
|
371 C nvj: the number of variables in a supervariable j (= nv (j)) |
|
|
372 C nvpiv: number of pivots in current element |
|
|
373 C slenme: number of variables in variable list of pivotal variable |
|
|
374 C we: w (e) |
|
|
375 C wflg: used for flagging the w array. See description of iw. |
|
|
376 C wnvi: wflg - nv (i) |
|
|
377 C x: either a supervariable or an element |
|
|
378 |
|
|
379 C----------------------------------------------------------------------- |
|
|
380 C LOCAL POINTERS: |
|
|
381 C----------------------------------------------------------------------- |
|
|
382 |
|
|
383 INTEGER P, P1, P2, P3, PDST, PEND, PJ, PME, PME1, PME2, PN, PSRC |
|
|
384 |
|
|
385 C Any parameter (pe (...) or pfree) or local variable |
|
|
386 C starting with "p" (for Pointer) is an index into iw, |
|
|
387 C and all indices into iw use variables starting with |
|
|
388 C "p." The only exception to this rule is the iwlen |
|
|
389 C input argument. |
|
|
390 |
|
|
391 C p: pointer into lots of things |
|
|
392 C p1: pe (i) for some variable i (start of element list) |
|
|
393 C p2: pe (i) + elen (i) - 1 for some var. i (end of el. list) |
|
|
394 C p3: index of first supervariable in clean list |
|
|
395 C pdst: destination pointer, for compression |
|
|
396 C pend: end of memory to compress |
|
|
397 C pj: pointer into an element or variable |
|
|
398 C pme: pointer into the current element (pme1...pme2) |
|
|
399 C pme1: the current element, me, is stored in iw (pme1...pme2) |
|
|
400 C pme2: the end of the current element |
|
|
401 C pn: pointer into a "clean" variable, also used to compress |
|
|
402 C psrc: source pointer, for compression |
|
|
403 |
|
|
404 C----------------------------------------------------------------------- |
|
|
405 C FUNCTIONS CALLED: |
|
|
406 C----------------------------------------------------------------------- |
|
|
407 |
|
|
408 INTRINSIC MAX, MIN, MOD |
|
|
409 |
|
|
410 C======================================================================= |
|
|
411 C INITIALIZATIONS |
|
|
412 C======================================================================= |
|
|
413 |
|
|
414 WFLG = 2 |
|
|
415 MINDEG = 1 |
|
|
416 NCMPA = 0 |
|
|
417 NEL = 0 |
|
|
418 HMOD = MAX (1, N-1) |
|
|
419 DMAX = 0 |
|
|
420 MEM = PFREE - 1 |
|
|
421 MAXMEM = MEM |
|
|
422 ME = 0 |
|
|
423 |
|
|
424 DO 10 I = 1, N |
|
|
425 LAST (I) = 0 |
|
|
426 HEAD (I) = 0 |
|
|
427 NV (I) = 1 |
|
|
428 W (I) = 1 |
|
|
429 ELEN (I) = 0 |
|
|
430 DEGREE (I) = LEN (I) |
|
|
431 10 CONTINUE |
|
|
432 |
|
|
433 C ---------------------------------------------------------------- |
|
|
434 C initialize degree lists and eliminate rows with no off-diag. nz. |
|
|
435 C ---------------------------------------------------------------- |
|
|
436 |
|
|
437 DO 20 I = 1, N |
|
|
438 |
|
|
439 DEG = DEGREE (I) |
|
|
440 |
|
|
441 IF (DEG .GT. 0) THEN |
|
|
442 |
|
|
443 C ---------------------------------------------------------- |
|
|
444 C place i in the degree list corresponding to its degree |
|
|
445 C ---------------------------------------------------------- |
|
|
446 |
|
|
447 INEXT = HEAD (DEG) |
|
|
448 IF (INEXT .NE. 0) LAST (INEXT) = I |
|
|
449 NEXT (I) = INEXT |
|
|
450 HEAD (DEG) = I |
|
|
451 |
|
|
452 ELSE |
|
|
453 |
|
|
454 C ---------------------------------------------------------- |
|
|
455 C we have a variable that can be eliminated at once because |
|
|
456 C there is no off-diagonal non-zero in its row. |
|
|
457 C ---------------------------------------------------------- |
|
|
458 |
|
|
459 NEL = NEL + 1 |
|
|
460 ELEN (I) = -NEL |
|
|
461 PE (I) = 0 |
|
|
462 W (I) = 0 |
|
|
463 |
|
|
464 ENDIF |
|
|
465 |
|
|
466 20 CONTINUE |
|
|
467 |
|
|
468 C======================================================================= |
|
|
469 C WHILE (selecting pivots) DO |
|
|
470 C======================================================================= |
|
|
471 |
|
|
472 30 CONTINUE |
|
|
473 IF (NEL .LT. N) THEN |
|
|
474 |
|
|
475 C======================================================================= |
|
|
476 C GET PIVOT OF MINIMUM DEGREE |
|
|
477 C======================================================================= |
|
|
478 |
|
|
479 C ------------------------------------------------------------- |
|
|
480 C find next supervariable for elimination |
|
|
481 C ------------------------------------------------------------- |
|
|
482 |
|
|
483 DO 40 DEG = MINDEG, N |
|
|
484 ME = HEAD (DEG) |
|
|
485 IF (ME .GT. 0) GOTO 50 |
|
|
486 40 CONTINUE |
|
|
487 50 CONTINUE |
|
|
488 MINDEG = DEG |
|
|
489 |
|
|
490 C ------------------------------------------------------------- |
|
|
491 C remove chosen variable from link list |
|
|
492 C ------------------------------------------------------------- |
|
|
493 |
|
|
494 INEXT = NEXT (ME) |
|
|
495 IF (INEXT .NE. 0) LAST (INEXT) = 0 |
|
|
496 HEAD (DEG) = INEXT |
|
|
497 |
|
|
498 C ------------------------------------------------------------- |
|
|
499 C me represents the elimination of pivots nel+1 to nel+nv(me). |
|
|
500 C place me itself as the first in this set. It will be moved |
|
|
501 C to the nel+nv(me) position when the permutation vectors are |
|
|
502 C computed. |
|
|
503 C ------------------------------------------------------------- |
|
|
504 |
|
|
505 ELENME = ELEN (ME) |
|
|
506 ELEN (ME) = - (NEL + 1) |
|
|
507 NVPIV = NV (ME) |
|
|
508 NEL = NEL + NVPIV |
|
|
509 |
|
|
510 C======================================================================= |
|
|
511 C CONSTRUCT NEW ELEMENT |
|
|
512 C======================================================================= |
|
|
513 |
|
|
514 C ------------------------------------------------------------- |
|
|
515 C At this point, me is the pivotal supervariable. It will be |
|
|
516 C converted into the current element. Scan list of the |
|
|
517 C pivotal supervariable, me, setting tree pointers and |
|
|
518 C constructing new list of supervariables for the new element, |
|
|
519 C me. p is a pointer to the current position in the old list. |
|
|
520 C ------------------------------------------------------------- |
|
|
521 |
|
|
522 C flag the variable "me" as being in Lme by negating nv (me) |
|
|
523 NV (ME) = -NVPIV |
|
|
524 DEGME = 0 |
|
|
525 |
|
|
526 IF (ELENME .EQ. 0) THEN |
|
|
527 |
|
|
528 C ---------------------------------------------------------- |
|
|
529 C construct the new element in place |
|
|
530 C ---------------------------------------------------------- |
|
|
531 |
|
|
532 PME1 = PE (ME) |
|
|
533 PME2 = PME1 - 1 |
|
|
534 |
|
|
535 DO 60 P = PME1, PME1 + LEN (ME) - 1 |
|
|
536 I = IW (P) |
|
|
537 NVI = NV (I) |
|
|
538 IF (NVI .GT. 0) THEN |
|
|
539 |
|
|
540 C ---------------------------------------------------- |
|
|
541 C i is a principal variable not yet placed in Lme. |
|
|
542 C store i in new list |
|
|
543 C ---------------------------------------------------- |
|
|
544 |
|
|
545 DEGME = DEGME + NVI |
|
|
546 C flag i as being in Lme by negating nv (i) |
|
|
547 NV (I) = -NVI |
|
|
548 PME2 = PME2 + 1 |
|
|
549 IW (PME2) = I |
|
|
550 |
|
|
551 C ---------------------------------------------------- |
|
|
552 C remove variable i from degree list. |
|
|
553 C ---------------------------------------------------- |
|
|
554 |
|
|
555 ILAST = LAST (I) |
|
|
556 INEXT = NEXT (I) |
|
|
557 IF (INEXT .NE. 0) LAST (INEXT) = ILAST |
|
|
558 IF (ILAST .NE. 0) THEN |
|
|
559 NEXT (ILAST) = INEXT |
|
|
560 ELSE |
|
|
561 C i is at the head of the degree list |
|
|
562 HEAD (DEGREE (I)) = INEXT |
|
|
563 ENDIF |
|
|
564 |
|
|
565 ENDIF |
|
|
566 60 CONTINUE |
|
|
567 C this element takes no new memory in iw: |
|
|
568 NEWMEM = 0 |
|
|
569 |
|
|
570 ELSE |
|
|
571 |
|
|
572 C ---------------------------------------------------------- |
|
|
573 C construct the new element in empty space, iw (pfree ...) |
|
|
574 C ---------------------------------------------------------- |
|
|
575 |
|
|
576 P = PE (ME) |
|
|
577 PME1 = PFREE |
|
|
578 SLENME = LEN (ME) - ELENME |
|
|
579 |
|
|
580 DO 120 KNT1 = 1, ELENME + 1 |
|
|
581 |
|
|
582 IF (KNT1 .GT. ELENME) THEN |
|
|
583 C search the supervariables in me. |
|
|
584 E = ME |
|
|
585 PJ = P |
|
|
586 LN = SLENME |
|
|
587 ELSE |
|
|
588 C search the elements in me. |
|
|
589 E = IW (P) |
|
|
590 P = P + 1 |
|
|
591 PJ = PE (E) |
|
|
592 LN = LEN (E) |
|
|
593 ENDIF |
|
|
594 |
|
|
595 C ------------------------------------------------------- |
|
|
596 C search for different supervariables and add them to the |
|
|
597 C new list, compressing when necessary. this loop is |
|
|
598 C executed once for each element in the list and once for |
|
|
599 C all the supervariables in the list. |
|
|
600 C ------------------------------------------------------- |
|
|
601 |
|
|
602 DO 110 KNT2 = 1, LN |
|
|
603 I = IW (PJ) |
|
|
604 PJ = PJ + 1 |
|
|
605 NVI = NV (I) |
|
|
606 IF (NVI .GT. 0) THEN |
|
|
607 |
|
|
608 C ------------------------------------------------- |
|
|
609 C compress iw, if necessary |
|
|
610 C ------------------------------------------------- |
|
|
611 |
|
|
612 IF (PFREE .GT. IWLEN) THEN |
|
|
613 C prepare for compressing iw by adjusting |
|
|
614 C pointers and lengths so that the lists being |
|
|
615 C searched in the inner and outer loops contain |
|
|
616 C only the remaining entries. |
|
|
617 |
|
|
618 PE (ME) = P |
|
|
619 LEN (ME) = LEN (ME) - KNT1 |
|
|
620 IF (LEN (ME) .EQ. 0) THEN |
|
|
621 C nothing left of supervariable me |
|
|
622 PE (ME) = 0 |
|
|
623 ENDIF |
|
|
624 PE (E) = PJ |
|
|
625 LEN (E) = LN - KNT2 |
|
|
626 IF (LEN (E) .EQ. 0) THEN |
|
|
627 C nothing left of element e |
|
|
628 PE (E) = 0 |
|
|
629 ENDIF |
|
|
630 |
|
|
631 NCMPA = NCMPA + 1 |
|
|
632 C store first item in pe |
|
|
633 C set first entry to -item |
|
|
634 DO 70 J = 1, N |
|
|
635 PN = PE (J) |
|
|
636 IF (PN .GT. 0) THEN |
|
|
637 PE (J) = IW (PN) |
|
|
638 IW (PN) = -J |
|
|
639 ENDIF |
|
|
640 70 CONTINUE |
|
|
641 |
|
|
642 C psrc/pdst point to source/destination |
|
|
643 PDST = 1 |
|
|
644 PSRC = 1 |
|
|
645 PEND = PME1 - 1 |
|
|
646 |
|
|
647 C while loop: |
|
|
648 80 CONTINUE |
|
|
649 IF (PSRC .LE. PEND) THEN |
|
|
650 C search for next negative entry |
|
|
651 J = -IW (PSRC) |
|
|
652 PSRC = PSRC + 1 |
|
|
653 IF (J .GT. 0) THEN |
|
|
654 IW (PDST) = PE (J) |
|
|
655 PE (J) = PDST |
|
|
656 PDST = PDST + 1 |
|
|
657 C copy from source to destination |
|
|
658 LENJ = LEN (J) |
|
|
659 DO 90 KNT3 = 0, LENJ - 2 |
|
|
660 IW (PDST + KNT3) = IW (PSRC + KNT3) |
|
|
661 90 CONTINUE |
|
|
662 PDST = PDST + LENJ - 1 |
|
|
663 PSRC = PSRC + LENJ - 1 |
|
|
664 ENDIF |
|
|
665 GOTO 80 |
|
|
666 ENDIF |
|
|
667 |
|
|
668 C move the new partially-constructed element |
|
|
669 P1 = PDST |
|
|
670 DO 100 PSRC = PME1, PFREE - 1 |
|
|
671 IW (PDST) = IW (PSRC) |
|
|
672 PDST = PDST + 1 |
|
|
673 100 CONTINUE |
|
|
674 PME1 = P1 |
|
|
675 PFREE = PDST |
|
|
676 PJ = PE (E) |
|
|
677 P = PE (ME) |
|
|
678 ENDIF |
|
|
679 |
|
|
680 C ------------------------------------------------- |
|
|
681 C i is a principal variable not yet placed in Lme |
|
|
682 C store i in new list |
|
|
683 C ------------------------------------------------- |
|
|
684 |
|
|
685 DEGME = DEGME + NVI |
|
|
686 C flag i as being in Lme by negating nv (i) |
|
|
687 NV (I) = -NVI |
|
|
688 IW (PFREE) = I |
|
|
689 PFREE = PFREE + 1 |
|
|
690 |
|
|
691 C ------------------------------------------------- |
|
|
692 C remove variable i from degree link list |
|
|
693 C ------------------------------------------------- |
|
|
694 |
|
|
695 ILAST = LAST (I) |
|
|
696 INEXT = NEXT (I) |
|
|
697 IF (INEXT .NE. 0) LAST (INEXT) = ILAST |
|
|
698 IF (ILAST .NE. 0) THEN |
|
|
699 NEXT (ILAST) = INEXT |
|
|
700 ELSE |
|
|
701 C i is at the head of the degree list |
|
|
702 HEAD (DEGREE (I)) = INEXT |
|
|
703 ENDIF |
|
|
704 |
|
|
705 ENDIF |
|
|
706 110 CONTINUE |
|
|
707 |
|
|
708 IF (E .NE. ME) THEN |
|
|
709 C set tree pointer and flag to indicate element e is |
|
|
710 C absorbed into new element me (the parent of e is me) |
|
|
711 PE (E) = -ME |
|
|
712 W (E) = 0 |
|
|
713 ENDIF |
|
|
714 120 CONTINUE |
|
|
715 |
|
|
716 PME2 = PFREE - 1 |
|
|
717 C this element takes newmem new memory in iw (possibly zero) |
|
|
718 NEWMEM = PFREE - PME1 |
|
|
719 MEM = MEM + NEWMEM |
|
|
720 MAXMEM = MAX (MAXMEM, MEM) |
|
|
721 ENDIF |
|
|
722 |
|
|
723 C ------------------------------------------------------------- |
|
|
724 C me has now been converted into an element in iw (pme1..pme2) |
|
|
725 C ------------------------------------------------------------- |
|
|
726 |
|
|
727 C degme holds the external degree of new element |
|
|
728 DEGREE (ME) = DEGME |
|
|
729 PE (ME) = PME1 |
|
|
730 LEN (ME) = PME2 - PME1 + 1 |
|
|
731 |
|
|
732 C ------------------------------------------------------------- |
|
|
733 C make sure that wflg is not too large. With the current |
|
|
734 C value of wflg, wflg+n must not cause integer overflow |
|
|
735 C ------------------------------------------------------------- |
|
|
736 |
|
|
737 IF (WFLG + N .LE. WFLG) THEN |
|
|
738 DO 130 X = 1, N |
|
|
739 IF (W (X) .NE. 0) W (X) = 1 |
|
|
740 130 CONTINUE |
|
|
741 WFLG = 2 |
|
|
742 ENDIF |
|
|
743 |
|
|
744 C======================================================================= |
|
|
745 C COMPUTE (w (e) - wflg) = |Le\Lme| FOR ALL ELEMENTS |
|
|
746 C======================================================================= |
|
|
747 |
|
|
748 C ------------------------------------------------------------- |
|
|
749 C Scan 1: compute the external degrees of previous elements |
|
|
750 C with respect to the current element. That is: |
|
|
751 C (w (e) - wflg) = |Le \ Lme| |
|
|
752 C for each element e that appears in any supervariable in Lme. |
|
|
753 C The notation Le refers to the pattern (list of |
|
|
754 C supervariables) of a previous element e, where e is not yet |
|
|
755 C absorbed, stored in iw (pe (e) + 1 ... pe (e) + iw (pe (e))). |
|
|
756 C The notation Lme refers to the pattern of the current element |
|
|
757 C (stored in iw (pme1..pme2)). If (w (e) - wflg) becomes |
|
|
758 C zero, then the element e will be absorbed in scan 2. |
|
|
759 C ------------------------------------------------------------- |
|
|
760 |
|
|
761 DO 150 PME = PME1, PME2 |
|
|
762 I = IW (PME) |
|
|
763 ELN = ELEN (I) |
|
|
764 IF (ELN .GT. 0) THEN |
|
|
765 C note that nv (i) has been negated to denote i in Lme: |
|
|
766 NVI = -NV (I) |
|
|
767 WNVI = WFLG - NVI |
|
|
768 DO 140 P = PE (I), PE (I) + ELN - 1 |
|
|
769 E = IW (P) |
|
|
770 WE = W (E) |
|
|
771 IF (WE .GE. WFLG) THEN |
|
|
772 C unabsorbed element e has been seen in this loop |
|
|
773 WE = WE - NVI |
|
|
774 ELSE IF (WE .NE. 0) THEN |
|
|
775 C e is an unabsorbed element |
|
|
776 C this is the first we have seen e in all of Scan 1 |
|
|
777 WE = DEGREE (E) + WNVI |
|
|
778 ENDIF |
|
|
779 W (E) = WE |
|
|
780 140 CONTINUE |
|
|
781 ENDIF |
|
|
782 150 CONTINUE |
|
|
783 |
|
|
784 C======================================================================= |
|
|
785 C DEGREE UPDATE AND ELEMENT ABSORPTION |
|
|
786 C======================================================================= |
|
|
787 |
|
|
788 C ------------------------------------------------------------- |
|
|
789 C Scan 2: for each i in Lme, sum up the degree of Lme (which |
|
|
790 C is degme), plus the sum of the external degrees of each Le |
|
|
791 C for the elements e appearing within i, plus the |
|
|
792 C supervariables in i. Place i in hash list. |
|
|
793 C ------------------------------------------------------------- |
|
|
794 |
|
|
795 DO 180 PME = PME1, PME2 |
|
|
796 I = IW (PME) |
|
|
797 P1 = PE (I) |
|
|
798 P2 = P1 + ELEN (I) - 1 |
|
|
799 PN = P1 |
|
|
800 HASH = 0 |
|
|
801 DEG = 0 |
|
|
802 |
|
|
803 C ---------------------------------------------------------- |
|
|
804 C scan the element list associated with supervariable i |
|
|
805 C ---------------------------------------------------------- |
|
|
806 |
|
|
807 C UMFPACK/MA38-style approximate degree: |
|
|
808 DO 160 P = P1, P2 |
|
|
809 E = IW (P) |
|
|
810 WE = W (E) |
|
|
811 IF (WE .NE. 0) THEN |
|
|
812 C e is an unabsorbed element |
|
|
813 DEG = DEG + WE - WFLG |
|
|
814 IW (PN) = E |
|
|
815 PN = PN + 1 |
|
|
816 HASH = HASH + E |
|
|
817 ENDIF |
|
|
818 160 CONTINUE |
|
|
819 |
|
|
820 C count the number of elements in i (including me): |
|
|
821 ELEN (I) = PN - P1 + 1 |
|
|
822 |
|
|
823 C ---------------------------------------------------------- |
|
|
824 C scan the supervariables in the list associated with i |
|
|
825 C ---------------------------------------------------------- |
|
|
826 |
|
|
827 P3 = PN |
|
|
828 DO 170 P = P2 + 1, P1 + LEN (I) - 1 |
|
|
829 J = IW (P) |
|
|
830 NVJ = NV (J) |
|
|
831 IF (NVJ .GT. 0) THEN |
|
|
832 C j is unabsorbed, and not in Lme. |
|
|
833 C add to degree and add to new list |
|
|
834 DEG = DEG + NVJ |
|
|
835 IW (PN) = J |
|
|
836 PN = PN + 1 |
|
|
837 HASH = HASH + J |
|
|
838 ENDIF |
|
|
839 170 CONTINUE |
|
|
840 |
|
|
841 C ---------------------------------------------------------- |
|
|
842 C update the degree and check for mass elimination |
|
|
843 C ---------------------------------------------------------- |
|
|
844 |
|
|
845 IF (ELEN (I) .EQ. 1 .AND. P3 .EQ. PN) THEN |
|
|
846 |
|
|
847 C ------------------------------------------------------- |
|
|
848 C mass elimination |
|
|
849 C ------------------------------------------------------- |
|
|
850 |
|
|
851 C There is nothing left of this node except for an |
|
|
852 C edge to the current pivot element. elen (i) is 1, |
|
|
853 C and there are no variables adjacent to node i. |
|
|
854 C Absorb i into the current pivot element, me. |
|
|
855 |
|
|
856 PE (I) = -ME |
|
|
857 NVI = -NV (I) |
|
|
858 DEGME = DEGME - NVI |
|
|
859 NVPIV = NVPIV + NVI |
|
|
860 NEL = NEL + NVI |
|
|
861 NV (I) = 0 |
|
|
862 ELEN (I) = 0 |
|
|
863 |
|
|
864 ELSE |
|
|
865 |
|
|
866 C ------------------------------------------------------- |
|
|
867 C update the upper-bound degree of i |
|
|
868 C ------------------------------------------------------- |
|
|
869 |
|
|
870 C the following degree does not yet include the size |
|
|
871 C of the current element, which is added later: |
|
|
872 DEGREE (I) = MIN (DEGREE (I), DEG) |
|
|
873 |
|
|
874 C ------------------------------------------------------- |
|
|
875 C add me to the list for i |
|
|
876 C ------------------------------------------------------- |
|
|
877 |
|
|
878 C move first supervariable to end of list |
|
|
879 IW (PN) = IW (P3) |
|
|
880 C move first element to end of element part of list |
|
|
881 IW (P3) = IW (P1) |
|
|
882 C add new element to front of list. |
|
|
883 IW (P1) = ME |
|
|
884 C store the new length of the list in len (i) |
|
|
885 LEN (I) = PN - P1 + 1 |
|
|
886 |
|
|
887 C ------------------------------------------------------- |
|
|
888 C place in hash bucket. Save hash key of i in last (i). |
|
|
889 C ------------------------------------------------------- |
|
|
890 |
|
|
891 HASH = MOD (HASH, HMOD) + 1 |
|
|
892 J = HEAD (HASH) |
|
|
893 IF (J .LE. 0) THEN |
|
|
894 C the degree list is empty, hash head is -j |
|
|
895 NEXT (I) = -J |
|
|
896 HEAD (HASH) = -I |
|
|
897 ELSE |
|
|
898 C degree list is not empty |
|
|
899 C use last (head (hash)) as hash head |
|
|
900 NEXT (I) = LAST (J) |
|
|
901 LAST (J) = I |
|
|
902 ENDIF |
|
|
903 LAST (I) = HASH |
|
|
904 ENDIF |
|
|
905 180 CONTINUE |
|
|
906 |
|
|
907 DEGREE (ME) = DEGME |
|
|
908 |
|
|
909 C ------------------------------------------------------------- |
|
|
910 C Clear the counter array, w (...), by incrementing wflg. |
|
|
911 C ------------------------------------------------------------- |
|
|
912 |
|
|
913 DMAX = MAX (DMAX, DEGME) |
|
|
914 WFLG = WFLG + DMAX |
|
|
915 |
|
|
916 C make sure that wflg+n does not cause integer overflow |
|
|
917 IF (WFLG + N .LE. WFLG) THEN |
|
|
918 DO 190 X = 1, N |
|
|
919 IF (W (X) .NE. 0) W (X) = 1 |
|
|
920 190 CONTINUE |
|
|
921 WFLG = 2 |
|
|
922 ENDIF |
|
|
923 C at this point, w (1..n) .lt. wflg holds |
|
|
924 |
|
|
925 C======================================================================= |
|
|
926 C SUPERVARIABLE DETECTION |
|
|
927 C======================================================================= |
|
|
928 |
|
|
929 DO 250 PME = PME1, PME2 |
|
|
930 I = IW (PME) |
|
|
931 IF (NV (I) .LT. 0) THEN |
|
|
932 C i is a principal variable in Lme |
|
|
933 |
|
|
934 C ------------------------------------------------------- |
|
|
935 C examine all hash buckets with 2 or more variables. We |
|
|
936 C do this by examing all unique hash keys for super- |
|
|
937 C variables in the pattern Lme of the current element, me |
|
|
938 C ------------------------------------------------------- |
|
|
939 |
|
|
940 HASH = LAST (I) |
|
|
941 C let i = head of hash bucket, and empty the hash bucket |
|
|
942 J = HEAD (HASH) |
|
|
943 IF (J .EQ. 0) GOTO 250 |
|
|
944 IF (J .LT. 0) THEN |
|
|
945 C degree list is empty |
|
|
946 I = -J |
|
|
947 HEAD (HASH) = 0 |
|
|
948 ELSE |
|
|
949 C degree list is not empty, restore last () of head |
|
|
950 I = LAST (J) |
|
|
951 LAST (J) = 0 |
|
|
952 ENDIF |
|
|
953 IF (I .EQ. 0) GOTO 250 |
|
|
954 |
|
|
955 C while loop: |
|
|
956 200 CONTINUE |
|
|
957 IF (NEXT (I) .NE. 0) THEN |
|
|
958 |
|
|
959 C ---------------------------------------------------- |
|
|
960 C this bucket has one or more variables following i. |
|
|
961 C scan all of them to see if i can absorb any entries |
|
|
962 C that follow i in hash bucket. Scatter i into w. |
|
|
963 C ---------------------------------------------------- |
|
|
964 |
|
|
965 LN = LEN (I) |
|
|
966 ELN = ELEN (I) |
|
|
967 C do not flag the first element in the list (me) |
|
|
968 DO 210 P = PE (I) + 1, PE (I) + LN - 1 |
|
|
969 W (IW (P)) = WFLG |
|
|
970 210 CONTINUE |
|
|
971 |
|
|
972 C ---------------------------------------------------- |
|
|
973 C scan every other entry j following i in bucket |
|
|
974 C ---------------------------------------------------- |
|
|
975 |
|
|
976 JLAST = I |
|
|
977 J = NEXT (I) |
|
|
978 |
|
|
979 C while loop: |
|
|
980 220 CONTINUE |
|
|
981 IF (J .NE. 0) THEN |
|
|
982 |
|
|
983 C ------------------------------------------------- |
|
|
984 C check if j and i have identical nonzero pattern |
|
|
985 C ------------------------------------------------- |
|
|
986 |
|
|
987 IF (LEN (J) .NE. LN) THEN |
|
|
988 C i and j do not have same size data structure |
|
|
989 GOTO 240 |
|
|
990 ENDIF |
|
|
991 IF (ELEN (J) .NE. ELN) THEN |
|
|
992 C i and j do not have same number of adjacent el |
|
|
993 GOTO 240 |
|
|
994 ENDIF |
|
|
995 C do not flag the first element in the list (me) |
|
|
996 DO 230 P = PE (J) + 1, PE (J) + LN - 1 |
|
|
997 IF (W (IW (P)) .NE. WFLG) THEN |
|
|
998 C an entry (iw(p)) is in j but not in i |
|
|
999 GOTO 240 |
|
|
1000 ENDIF |
|
|
1001 230 CONTINUE |
|
|
1002 |
|
|
1003 C ------------------------------------------------- |
|
|
1004 C found it! j can be absorbed into i |
|
|
1005 C ------------------------------------------------- |
|
|
1006 |
|
|
1007 PE (J) = -I |
|
|
1008 C both nv (i) and nv (j) are negated since they |
|
|
1009 C are in Lme, and the absolute values of each |
|
|
1010 C are the number of variables in i and j: |
|
|
1011 NV (I) = NV (I) + NV (J) |
|
|
1012 NV (J) = 0 |
|
|
1013 ELEN (J) = 0 |
|
|
1014 C delete j from hash bucket |
|
|
1015 J = NEXT (J) |
|
|
1016 NEXT (JLAST) = J |
|
|
1017 GOTO 220 |
|
|
1018 |
|
|
1019 C ------------------------------------------------- |
|
|
1020 240 CONTINUE |
|
|
1021 C j cannot be absorbed into i |
|
|
1022 C ------------------------------------------------- |
|
|
1023 |
|
|
1024 JLAST = J |
|
|
1025 J = NEXT (J) |
|
|
1026 GOTO 220 |
|
|
1027 ENDIF |
|
|
1028 |
|
|
1029 C ---------------------------------------------------- |
|
|
1030 C no more variables can be absorbed into i |
|
|
1031 C go to next i in bucket and clear flag array |
|
|
1032 C ---------------------------------------------------- |
|
|
1033 |
|
|
1034 WFLG = WFLG + 1 |
|
|
1035 I = NEXT (I) |
|
|
1036 IF (I .NE. 0) GOTO 200 |
|
|
1037 ENDIF |
|
|
1038 ENDIF |
|
|
1039 250 CONTINUE |
|
|
1040 |
|
|
1041 C======================================================================= |
|
|
1042 C RESTORE DEGREE LISTS AND REMOVE NONPRINCIPAL SUPERVAR. FROM ELEMENT |
|
|
1043 C======================================================================= |
|
|
1044 |
|
|
1045 P = PME1 |
|
|
1046 NLEFT = N - NEL |
|
|
1047 DO 260 PME = PME1, PME2 |
|
|
1048 I = IW (PME) |
|
|
1049 NVI = -NV (I) |
|
|
1050 IF (NVI .GT. 0) THEN |
|
|
1051 C i is a principal variable in Lme |
|
|
1052 C restore nv (i) to signify that i is principal |
|
|
1053 NV (I) = NVI |
|
|
1054 |
|
|
1055 C ------------------------------------------------------- |
|
|
1056 C compute the external degree (add size of current elem) |
|
|
1057 C ------------------------------------------------------- |
|
|
1058 |
|
|
1059 DEG = MAX (1, MIN (DEGREE (I) + DEGME-NVI, NLEFT-NVI)) |
|
|
1060 |
|
|
1061 C ------------------------------------------------------- |
|
|
1062 C place the supervariable at the head of the degree list |
|
|
1063 C ------------------------------------------------------- |
|
|
1064 |
|
|
1065 INEXT = HEAD (DEG) |
|
|
1066 IF (INEXT .NE. 0) LAST (INEXT) = I |
|
|
1067 NEXT (I) = INEXT |
|
|
1068 LAST (I) = 0 |
|
|
1069 HEAD (DEG) = I |
|
|
1070 |
|
|
1071 C ------------------------------------------------------- |
|
|
1072 C save the new degree, and find the minimum degree |
|
|
1073 C ------------------------------------------------------- |
|
|
1074 |
|
|
1075 MINDEG = MIN (MINDEG, DEG) |
|
|
1076 DEGREE (I) = DEG |
|
|
1077 |
|
|
1078 C ------------------------------------------------------- |
|
|
1079 C place the supervariable in the element pattern |
|
|
1080 C ------------------------------------------------------- |
|
|
1081 |
|
|
1082 IW (P) = I |
|
|
1083 P = P + 1 |
|
|
1084 ENDIF |
|
|
1085 260 CONTINUE |
|
|
1086 |
|
|
1087 C======================================================================= |
|
|
1088 C FINALIZE THE NEW ELEMENT |
|
|
1089 C======================================================================= |
|
|
1090 |
|
|
1091 NV (ME) = NVPIV + DEGME |
|
|
1092 C nv (me) is now the degree of pivot (including diagonal part) |
|
|
1093 C save the length of the list for the new element me |
|
|
1094 LEN (ME) = P - PME1 |
|
|
1095 IF (LEN (ME) .EQ. 0) THEN |
|
|
1096 C there is nothing left of the current pivot element |
|
|
1097 PE (ME) = 0 |
|
|
1098 W (ME) = 0 |
|
|
1099 ENDIF |
|
|
1100 IF (NEWMEM .NE. 0) THEN |
|
|
1101 C element was not constructed in place: deallocate part |
|
|
1102 C of it (final size is less than or equal to newmem, |
|
|
1103 C since newly nonprincipal variables have been removed). |
|
|
1104 PFREE = P |
|
|
1105 MEM = MEM - NEWMEM + LEN (ME) |
|
|
1106 ENDIF |
|
|
1107 |
|
|
1108 C======================================================================= |
|
|
1109 C END WHILE (selecting pivots) |
|
|
1110 GOTO 30 |
|
|
1111 ENDIF |
|
|
1112 C======================================================================= |
|
|
1113 |
|
|
1114 C======================================================================= |
|
|
1115 C COMPUTE THE PERMUTATION VECTORS |
|
|
1116 C======================================================================= |
|
|
1117 |
|
|
1118 C ---------------------------------------------------------------- |
|
|
1119 C The time taken by the following code is O(n). At this |
|
|
1120 C point, elen (e) = -k has been done for all elements e, |
|
|
1121 C and elen (i) = 0 has been done for all nonprincipal |
|
|
1122 C variables i. At this point, there are no principal |
|
|
1123 C supervariables left, and all elements are absorbed. |
|
|
1124 C ---------------------------------------------------------------- |
|
|
1125 |
|
|
1126 C ---------------------------------------------------------------- |
|
|
1127 C compute the ordering of unordered nonprincipal variables |
|
|
1128 C ---------------------------------------------------------------- |
|
|
1129 |
|
|
1130 DO 290 I = 1, N |
|
|
1131 IF (ELEN (I) .EQ. 0) THEN |
|
|
1132 |
|
|
1133 C ---------------------------------------------------------- |
|
|
1134 C i is an un-ordered row. Traverse the tree from i until |
|
|
1135 C reaching an element, e. The element, e, was the |
|
|
1136 C principal supervariable of i and all nodes in the path |
|
|
1137 C from i to when e was selected as pivot. |
|
|
1138 C ---------------------------------------------------------- |
|
|
1139 |
|
|
1140 J = -PE (I) |
|
|
1141 C while (j is a variable) do: |
|
|
1142 270 CONTINUE |
|
|
1143 IF (ELEN (J) .GE. 0) THEN |
|
|
1144 J = -PE (J) |
|
|
1145 GOTO 270 |
|
|
1146 ENDIF |
|
|
1147 E = J |
|
|
1148 |
|
|
1149 C ---------------------------------------------------------- |
|
|
1150 C get the current pivot ordering of e |
|
|
1151 C ---------------------------------------------------------- |
|
|
1152 |
|
|
1153 K = -ELEN (E) |
|
|
1154 |
|
|
1155 C ---------------------------------------------------------- |
|
|
1156 C traverse the path again from i to e, and compress the |
|
|
1157 C path (all nodes point to e). Path compression allows |
|
|
1158 C this code to compute in O(n) time. Order the unordered |
|
|
1159 C nodes in the path, and place the element e at the end. |
|
|
1160 C ---------------------------------------------------------- |
|
|
1161 |
|
|
1162 J = I |
|
|
1163 C while (j is a variable) do: |
|
|
1164 280 CONTINUE |
|
|
1165 IF (ELEN (J) .GE. 0) THEN |
|
|
1166 JNEXT = -PE (J) |
|
|
1167 PE (J) = -E |
|
|
1168 IF (ELEN (J) .EQ. 0) THEN |
|
|
1169 C j is an unordered row |
|
|
1170 ELEN (J) = K |
|
|
1171 K = K + 1 |
|
|
1172 ENDIF |
|
|
1173 J = JNEXT |
|
|
1174 GOTO 280 |
|
|
1175 ENDIF |
|
|
1176 C leave elen (e) negative, so we know it is an element |
|
|
1177 ELEN (E) = -K |
|
|
1178 ENDIF |
|
|
1179 290 CONTINUE |
|
|
1180 |
|
|
1181 C ---------------------------------------------------------------- |
|
|
1182 C reset the inverse permutation (elen (1..n)) to be positive, |
|
|
1183 C and compute the permutation (last (1..n)). |
|
|
1184 C ---------------------------------------------------------------- |
|
|
1185 |
|
|
1186 DO 300 I = 1, N |
|
|
1187 K = ABS (ELEN (I)) |
|
|
1188 LAST (K) = I |
|
|
1189 ELEN (I) = K |
|
|
1190 300 CONTINUE |
|
|
1191 |
|
|
1192 C======================================================================= |
|
|
1193 C RETURN THE MEMORY USAGE IN IW |
|
|
1194 C======================================================================= |
|
|
1195 |
|
|
1196 C If maxmem is less than or equal to iwlen, then no compressions |
|
|
1197 C occurred, and iw (maxmem+1 ... iwlen) was unused. Otherwise |
|
|
1198 C compressions did occur, and iwlen would have had to have been |
|
|
1199 C greater than or equal to maxmem for no compressions to occur. |
|
|
1200 C Return the value of maxmem in the pfree argument. |
|
|
1201 |
|
|
1202 PFREE = MAXMEM |
|
|
1203 |
|
|
1204 RETURN |
|
|
1205 END |
|
|
1206 |
