1 /* amdatr.f -- translated by f2c (version of 23 April 1993 18:34:30).
2 You must link the resulting object file with the libraries:
3 -lf2c -lm (in that order)
4 */
5
6 #include "f2c.h"
7
amdatr_(n,pe,iw,len,iwlen,pfree,nv,next,last,head,elen,degree,ncmpa,w,iovflo)8 /* Subroutine */ int amdatr_(n, pe, iw, len, iwlen, pfree, nv, next, last,
9 head, elen, degree, ncmpa, w, iovflo)
10 integer *n, *pe, *iw, *len, *iwlen, *pfree, *nv, *next, *last, *head, *elen, *
11 degree, *ncmpa, *w, *iovflo;
12 {
13 /* System generated locals */
14 integer i__1, i__2, i__3;
15
16 /* Local variables */
17 static integer hash, pend, hmod, lenj, dmax_, wbig, wflg, psrc, pdst,
18 wnvi, e, i, j, k, p, degme, x, nleft, ilast, jlast, inext, jnext,
19 p1, nvpiv, p2, p3, me, ln, we, pj, pn, mindeg, elenme, slenme,
20 maxmem, newmem, deg, eln, mem, nel, pme, nvi, nvj, pme1, pme2,
21 knt1, knt2, knt3;
22
23 /* -----------------------------------------------------------------------
24 */
25 /* The MC47 / AMD suite of minimum degree ordering algorithms. */
26
27 /* This code is one of seven variations of a single algorithm: */
28 /* the primary routine (MC47B/BD, only available in the Harwell */
29 /* Subroutine Library), and 6 variations that differ only in */
30 /* how they compute the degree (available in NETLIB). */
31
32 /* For information on the Harwell Subroutine Library, contact */
33 /* John Harding, Harwell Subroutine Library, B 552, AEA Technology, */
34 /* Harwell, Didcot, Oxon OX11 0RA, telephone (44) 1235 434573, */
35 /* fax (44) 1235 434340, email john.harding@aeat.co.uk, who will */
36 /* provide details of price and conditions of use. */
37 /* -----------------------------------------------------------------------
38 */
39 /* ***********************************************************************
40 */
41 /* NOTICE: "The AMD routines (AMDEXA, AMDBAR, AMDHAF, AMDHAT, AMDTRU, */
42 /* and AMDATR) may be used SOLELY for educational, research, and */
43 /* benchmarking purposes by non-profit organizations and the U.S. */
44 /* government. Commercial and other organizations may make use of the */
45 /* AMD routines SOLELY for benchmarking purposes only. The AMD */
46 /* routines may be modified by or on behalf of the User for such */
47 /* use but at no time shall the AMD routines or any such modified */
48 /* version of them become the property of the User. The AMD routines */
49 /* are provided without warranty of any kind, either expressed or */
50 /* implied. Neither the Authors nor their employers shall be liable */
51 /* for any direct or consequential loss or damage whatsoever arising */
52 /* out of the use or misuse of the AMD routines by the User. The AMD */
53 /* routines must not be sold. You may make copies of the AMD routines, */
54 /* but this NOTICE and the Copyright notice must appear in all copies. */
55 /* Any other use of the AMD routines requires written permission. */
56 /* Your use of the AMD routines is an implicit agreement to these */
57 /* conditions." */
58 /* ***********************************************************************
59 */
60 /* -----------------------------------------------------------------------
61 */
62 /* AMDatr: Approximate Minimum (UMFPACK/MA38-style, true) Degree */
63 /* ordering algorithm, but without aggresive absorption */
64 /* -----------------------------------------------------------------------
65 */
66 /* Variation 6: MC47-style approximate true degree, but with no */
67 /* aggresive absorption. Note that some of the comments in the code */
68 /* below reflect the approximate *external* degree, whereas this */
69 /* code actually uses the approximate *true* degree. */
70
71 /* We recommend using MC47B/BD instead of this routine since MC47B/BD */
72 /* gives better results in about the same time. */
73 /* -----------------------------------------------------------------------
74 */
75 /* Given a representation of the nonzero pattern of a symmetric matrix, */
76 /* A, (excluding the diagonal) perform an approximate minimum */
77 /* (UMFPACK/MA38-style) degree ordering to compute a pivot order */
78 /* such that the introduction of nonzeros (fill-in) in the Cholesky
79 */
80 /* factors A = LL^T are kept low. At each step, the pivot */
81 /* selected is the one with the minimum UMFAPACK/MA38-style */
82 /* upper-bound on the true degree. This routine does not */
83 /* perform aggresive absorption (as done by MC47B/BD). Aggresive */
84 /* absorption in MC47B/BD is used to tighten the bound on the */
85 /* degree. This can result an significant improvement in the */
86 /* quality of the ordering for some matrices. */
87
88 /* The approximate degree algorithm implemented here is the */
89 /* symmetric analog of the degree update algorithm in MA38 and */
90 /* UMFPACK (the Unsymmetric-pattern MultiFrontal PACKage, both by */
91 /* Davis and Duff, available for academic users in NETLIB as */
92 /* linalg/umfpack.shar or via anonymous ftp to */
93 /* ftp.cis.ufl.edu:pub/umfpack). Non-academic users must use */
94 /* MA38 in the Harwell Subroutine Library instead of UMPFACK. */
95 /* **********************************************************************
96 */
97 /* ***** CAUTION: ARGUMENTS ARE NOT CHECKED FOR ERRORS ON INPUT. ******
98 */
99 /* **********************************************************************
100 */
101 /* ** If you want error checking, a more versatile input format, and a **
102 */
103 /* ** simpler user interface, then use MC47A/AD in the Harwell **
104 */
105 /* ** Subroutine Library, which checks for errors, transforms the **
106 */
107 /* ** input, and calls MC47B/BD. **
108 */
109 /* **********************************************************************
110 */
111 /* References: (UF Tech Reports are available via anonymous ftp */
112 /* to ftp.cis.ufl.edu:cis/tech-reports). */
113
114 /* [1] Timothy A. Davis and Iain Duff, "An unsymmetric-pattern */
115 /* multifrontal method for sparse LU factorization", */
116 /* SIAM J. Matrix Analysis and Applications, to appear. */
117 /* also Univ. of Florida Technical Report TR-94-038. */
118 /* Discusses UMFPACK / MA38. */
119
120 /* [2] Patrick Amestoy, Timothy A. Davis, and Iain S. Duff, */
121 /* "An approximate minimum degree ordering algorithm," */
122 /* SIAM J. Matrix Analysis and Applications (to appear), */
123 /* also Univ. of Florida Technical Report TR-94-039. */
124 /* Discusses this routine. */
125
126 /* [3] Alan George and Joseph Liu, "The evolution of the */
127 /* minimum degree ordering algorithm," SIAM Review, vol. */
128 /* 31, no. 1, pp. 1-19, March 1989. We list below the */
129 /* features mentioned in that paper that this code */
130 /* includes: */
131
132 /* mass elimination: */
133 /* Yes. MA27 relied on supervariable detection for mass */
134 /* elimination. */
135 /* indistinguishable nodes: */
136 /* Yes (we call these "supervariables"). This was also in
137 */
138 /* the MA27 code - although we modified the method of */
139 /* detecting them (the previous hash was the true degree, */
140 /* which we no longer keep track of). A supervariable is */
141 /* a set of rows with identical nonzero pattern. All */
142 /* variables in a supervariable are eliminated together. */
143 /* Each supervariable has as its numerical name that of */
144 /* one of its variables (its principal variable). */
145 /* quotient graph representation: */
146 /* Yes. We use the term "element" for the cliques formed */
147 /* during elimination. This was also in the MA27 code. */
148 /* The algorithm can operate in place, but it will work */
149 /* more efficiently if given some "elbow room." */
150 /* element absorption: */
151 /* Yes. This was also in the MA27 code. */
152 /* external degree: */
153 /* Yes. The MA27 code was based on the true degree. */
154 /* incomplete degree update and multiple elimination: */
155 /* No. This was not in MA27, either. Our method of */
156 /* degree update within MC47B/BD is element-based, not */
157 /* variable-based. It is thus not well-suited for use */
158 /* with incomplete degree update or multiple elimination. */
159 /* -----------------------------------------------------------------------
160 */
161 /* Authors, and Copyright (C) 1995 by: */
162 /* Timothy A. Davis, Patrick Amestoy, Iain S. Duff, & John K. Reid.
163 */
164
165 /* Acknowledgements: */
166 /* This work (and the UMFPACK package) was supported by the */
167 /* National Science Foundation (ASC-9111263 and DMS-9223088). */
168 /* The UMFPACK/MA38 approximate degree update algorithm, the */
169 /* unsymmetric analog which forms the basis of MC47B/BD, was */
170 /* developed while Tim Davis was supported by CERFACS (Toulouse, */
171 /* France) in a post-doctoral position. */
172
173 /* Date: September, 1995 */
174 /* -----------------------------------------------------------------------
175 */
176 /* -----------------------------------------------------------------------
177 */
178 /* INPUT ARGUMENTS (unaltered): */
179 /* -----------------------------------------------------------------------
180 */
181 /* n: The matrix order. */
182
183 /* Restriction: 1 .le. n .lt. (iovflo/2)-2 */
184 /* iwlen: The length of iw (1..iwlen). On input, the matrix is */
185 /* stored in iw (1..pfree-1). However, iw (1..iwlen) should be */
186 /* slightly larger than what is required to hold the matrix, at */
187 /* least iwlen .ge. pfree + n is recommended. Otherwise, */
188 /* excessive compressions will take place. */
189 /* *** We do not recommend running this algorithm with *** */
190 /* *** iwlen .lt. pfree + n. *** */
191 /* *** Better performance will be obtained if *** */
192 /* *** iwlen .ge. pfree + n *** */
193 /* *** or better yet *** */
194 /* *** iwlen .gt. 1.2 * pfree *** */
195 /* *** (where pfree is its value on input). *** */
196 /* The algorithm will not run at all if iwlen .lt. pfree-1. */
197
198 /* Restriction: iwlen .ge. pfree-1 */
199 /* iovflo: The largest positive integer that your computer can */
200 /* represent (-iovflo should also be representable). On a 32-bit */
201 /* computer with 2's-complement arithmetic, */
202 /* iovflo = (2^31)-1 = 2,147,483,648. */
203 /* -----------------------------------------------------------------------
204 */
205 /* INPUT/OUPUT ARGUMENTS: */
206 /* -----------------------------------------------------------------------
207 */
208 /* pe: On input, pe (i) is the index in iw of the start of row i, or */
209 /* zero if row i has no off-diagonal non-zeros. */
210
211 /* During execution, it is used for both supervariables and */
212 /* elements: */
213
214 /* * Principal supervariable i: index into iw of the */
215 /* description of supervariable i. A supervariable */
216 /* represents one or more rows of the matrix */
217 /* with identical nonzero pattern. */
218 /* * Non-principal supervariable i: if i has been absorbed */
219 /* into another supervariable j, then pe (i) = -j. */
220 /* That is, j has the same pattern as i. */
221 /* Note that j might later be absorbed into another */
222 /* supervariable j2, in which case pe (i) is still -j, */
223 /* and pe (j) = -j2. */
224 /* * Unabsorbed element e: the index into iw of the description */
225 /* of element e, if e has not yet been absorbed by a */
226 /* subsequent element. Element e is created when */
227 /* the supervariable of the same name is selected as */
228 /* the pivot. */
229 /* * Absorbed element e: if element e is absorbed into element */
230 /* e2, then pe (e) = -e2. This occurs when the pattern of
231 */
232 /* e (that is, Le) is found to be a subset of the pattern */
233 /* of e2 (that is, Le2). If element e is "null" (it has */
234 /* no nonzeros outside its pivot block), then pe (e) = 0. */
235
236 /* On output, pe holds the assembly tree/forest, which implicitly */
237 /* represents a pivot order with identical fill-in as the actual */
238 /* order (via a depth-first search of the tree). */
239
240 /* On output: */
241 /* If nv (i) .gt. 0, then i represents a node in the assembly tree,
242 */
243 /* and the parent of i is -pe (i), or zero if i is a root. */
244 /* If nv (i) = 0, then (i,-pe (i)) represents an edge in a */
245 /* subtree, the root of which is a node in the assembly tree. */
246 /* pfree: On input the tail end of the array, iw (pfree..iwlen), */
247 /* is empty, and the matrix is stored in iw (1..pfree-1). */
248 /* During execution, additional data is placed in iw, and pfree */
249 /* is modified so that iw (pfree..iwlen) is always the unused part
250 */
251 /* of iw. On output, pfree is set equal to the size of iw that */
252 /* would have been needed for no compressions to occur. If */
253 /* ncmpa is zero, then pfree (on output) is less than or equal to */
254 /* iwlen, and the space iw (pfree+1 ... iwlen) was not used. */
255 /* Otherwise, pfree (on output) is greater than iwlen, and all the
256 */
257 /* memory in iw was used. */
258 /* -----------------------------------------------------------------------
259 */
260 /* INPUT/MODIFIED (undefined on output): */
261 /* -----------------------------------------------------------------------
262 */
263 /* len: On input, len (i) holds the number of entries in row i of the */
264 /* matrix, excluding the diagonal. The contents of len (1..n) */
265 /* are undefined on output. */
266 /* iw: On input, iw (1..pfree-1) holds the description of each row i */
267 /* in the matrix. The matrix must be symmetric, and both upper */
268 /* and lower triangular parts must be present. The diagonal must */
269 /* not be present. Row i is held as follows: */
270
271 /* len (i): the length of the row i data structure */
272 /* iw (pe (i) ... pe (i) + len (i) - 1): */
273 /* the list of column indices for nonzeros */
274 /* in row i (simple supervariables), excluding */
275 /* the diagonal. All supervariables start with */
276 /* one row/column each (supervariable i is just */
277 /* row i). */
278 /* if len (i) is zero on input, then pe (i) is ignored */
279 /* on input. */
280
281 /* Note that the rows need not be in any particular order,
282 */
283 /* and there may be empty space between the rows. */
284
285 /* During execution, the supervariable i experiences fill-in. */
286 /* This is represented by placing in i a list of the elements */
287 /* that cause fill-in in supervariable i: */
288
289 /* len (i): the length of supervariable i */
290 /* iw (pe (i) ... pe (i) + elen (i) - 1): */
291 /* the list of elements that contain i. This list
292 */
293 /* is kept short by removing absorbed elements. */
294 /* iw (pe (i) + elen (i) ... pe (i) + len (i) - 1): */
295 /* the list of supervariables in i. This list */
296 /* is kept short by removing nonprincipal */
297 /* variables, and any entry j that is also */
298 /* contained in at least one of the elements */
299 /* (j in Le) in the list for i (e in row i). */
300
301 /* When supervariable i is selected as pivot, we create an */
302 /* element e of the same name (e=i): */
303
304 /* len (e): the length of element e */
305 /* iw (pe (e) ... pe (e) + len (e) - 1): */
306 /* the list of supervariables in element e. */
307
308 /* An element represents the fill-in that occurs when supervariable
309 */
310 /* i is selected as pivot (which represents the selection of row i
311 */
312 /* and all non-principal variables whose principal variable is i).
313 */
314 /* We use the term Le to denote the set of all supervariables */
315 /* in element e. Absorbed supervariables and elements are pruned */
316 /* from these lists when computationally convenient. */
317
318 /* CAUTION: THE INPUT MATRIX IS OVERWRITTEN DURING COMPUTATION. */
319 /* The contents of iw are undefined on output. */
320 /* -----------------------------------------------------------------------
321 */
322 /* OUTPUT (need not be set on input): */
323 /* -----------------------------------------------------------------------
324 */
325 /* nv: During execution, abs (nv (i)) is equal to the number of rows */
326 /* that are represented by the principal supervariable i. If i is
327 */
328 /* a nonprincipal variable, then nv (i) = 0. Initially, */
329 /* nv (i) = 1 for all i. nv (i) .lt. 0 signifies that i is a */
330 /* principal variable in the pattern Lme of the current pivot */
331 /* element me. On output, nv (e) holds the true degree of element
332 */
333 /* e at the time it was created (including the diagonal part). */
334 /* ncmpa: The number of times iw was compressed. If this is */
335 /* excessive, then the execution took longer than what could have */
336 /* been. To reduce ncmpa, try increasing iwlen to be 10% or 20% */
337 /* larger than the value of pfree on input (or at least */
338 /* iwlen .ge. pfree + n). The fastest performance will be */
339 /* obtained when ncmpa is returned as zero. If iwlen is set to */
340 /* the value returned by pfree on *output*, then no compressions */
341 /* will occur. */
342 /* elen: See the description of iw above. At the start of execution, */
343 /* elen (i) is set to zero. During execution, elen (i) is the */
344 /* number of elements in the list for supervariable i. When e */
345 /* becomes an element, elen (e) = -nel is set, where nel is the */
346 /* current step of factorization. elen (i) = 0 is done when i */
347 /* becomes nonprincipal. */
348
349 /* For variables, elen (i) .ge. 0 holds until just before the */
350 /* permutation vectors are computed. For elements, */
351 /* elen (e) .lt. 0 holds. */
352
353 /* On output elen (1..n) holds the inverse permutation (the same */
354 /* as the 'INVP' argument in Sparspak). That is, if k = elen (i),
355 */
356 /* then row i is the kth pivot row. Row i of A appears as the */
357 /* (elen(i))-th row in the permuted matrix, PAP^T. */
358 /* last: In a degree list, last (i) is the supervariable preceding i, */
359 /* or zero if i is the head of the list. In a hash bucket, */
360 /* last (i) is the hash key for i. last (head (hash)) is also */
361 /* used as the head of a hash bucket if head (hash) contains a */
362 /* degree list (see head, below). */
363
364 /* On output, last (1..n) holds the permutation (the same as the */
365 /* 'PERM' argument in Sparspak). That is, if i = last (k), then */
366 /* row i is the kth pivot row. Row last (k) of A is the k-th row */
367 /* in the permuted matrix, PAP^T. */
368 /* -----------------------------------------------------------------------
369 */
370 /* LOCAL (not input or output - used only during execution): */
371 /* -----------------------------------------------------------------------
372 */
373 /* degree: If i is a supervariable, then degree (i) holds the */
374 /* current approximation of the external degree of row i (an upper
375 */
376 /* bound). The external degree is the number of nonzeros in row i,
377 */
378 /* minus abs (nv (i)) (the diagonal part). The bound is equal to */
379 /* the external degree if elen (i) is less than or equal to two. */
380
381 /* We also use the term "external degree" for elements e to refer */
382 /* to |Le \ Lme|. If e is an element, then degree (e) holds |Le|,
383 */
384 /* which is the degree of the off-diagonal part of the element e */
385 /* (not including the diagonal part). */
386 /* head: head is used for degree lists. head (deg) is the first */
387 /* supervariable in a degree list (all supervariables i in a */
388 /* degree list deg have the same approximate degree, namely, */
389 /* deg = degree (i)). If the list deg is empty then */
390 /* head (deg) = 0. */
391
392 /* During supervariable detection head (hash) also serves as a */
393 /* pointer to a hash bucket. */
394 /* If head (hash) .gt. 0, there is a degree list of degree hash. */
395 /* The hash bucket head pointer is last (head (hash)). */
396 /* If head (hash) = 0, then the degree list and hash bucket are */
397 /* both empty. */
398 /* If head (hash) .lt. 0, then the degree list is empty, and */
399 /* -head (hash) is the head of the hash bucket. */
400 /* After supervariable detection is complete, all hash buckets */
401 /* are empty, and the (last (head (hash)) = 0) condition is */
402 /* restored for the non-empty degree lists. */
403 /* next: next (i) is the supervariable following i in a link list, or */
404 /* zero if i is the last in the list. Used for two kinds of */
405 /* lists: degree lists and hash buckets (a supervariable can be */
406 /* in only one kind of list at a time). */
407 /* w: The flag array w determines the status of elements and */
408 /* variables, and the external degree of elements. */
409
410 /* for elements: */
411 /* if w (e) = 0, then the element e is absorbed */
412 /* if w (e) .ge. wflg, then w (e) - wflg is the size of */
413 /* the set |Le \ Lme|, in terms of nonzeros (the */
414 /* sum of abs (nv (i)) for each principal variable i that */
415 /* is both in the pattern of element e and NOT in the */
416 /* pattern of the current pivot element, me). */
417 /* if wflg .gt. w (e) .gt. 0, then e is not absorbed and has */
418 /* not yet been seen in the scan of the element lists in */
419 /* the computation of |Le\Lme| in loop 150 below. */
420
421 /* for variables: */
422 /* during supervariable detection, if w (j) .ne. wflg then j is
423 */
424 /* not in the pattern of variable i */
425
426 /* The w array is initialized by setting w (i) = 1 for all i, */
427 /* and by setting wflg = 2. It is reinitialized if wflg becomes */
428 /* too large (to ensure that wflg+n does not cause integer */
429 /* overflow). */
430 /* -----------------------------------------------------------------------
431 */
432 /* LOCAL INTEGERS: */
433 /* -----------------------------------------------------------------------
434 */
435 /* deg: the degree of a variable or element */
436 /* degme: size, |Lme|, of the current element, me (= degree (me))
437 */
438 /* dext: external degree, |Le \ Lme|, of some element e */
439 /* dmax: largest |Le| seen so far */
440 /* e: an element */
441 /* elenme: the length, elen (me), of element list of pivotal var. */
442 /* eln: the length, elen (...), of an element list */
443 /* hash: the computed value of the hash function */
444 /* hmod: the hash function is computed modulo hmod = max (1,n-1)
445 */
446 /* i: a supervariable */
447 /* ilast: the entry in a link list preceding i */
448 /* inext: the entry in a link list following i */
449 /* j: a supervariable */
450 /* jlast: the entry in a link list preceding j */
451 /* jnext: the entry in a link list, or path, following j */
452 /* k: the pivot order of an element or variable */
453 /* knt1: loop counter used during element construction */
454 /* knt2: loop counter used during element construction */
455 /* knt3: loop counter used during compression */
456 /* lenj: len (j) */
457 /* ln: length of a supervariable list */
458 /* maxmem: amount of memory needed for no compressions */
459 /* me: current supervariable being eliminated, and the */
460 /* current element created by eliminating that */
461 /* supervariable */
462 /* mem: memory in use assuming no compressions have occurred */
463 /* mindeg: current minimum degree */
464 /* nel: number of pivots selected so far */
465 /* newmem: amount of new memory needed for current pivot element */
466 /* nleft: n - nel, the number of nonpivotal rows/columns remaining
467 */
468 /* nvi: the number of variables in a supervariable i (= nv (i))
469 */
470 /* nvj: the number of variables in a supervariable j (= nv (j))
471 */
472 /* nvpiv: number of pivots in current element */
473 /* slenme: number of variables in variable list of pivotal variable
474 */
475 /* wbig: = iovflo - n. wflg is not allowed to be .ge. wbig. */
476 /* we: w (e) */
477 /* wflg: used for flagging the w array. See description of iw. */
478 /* wnvi: wflg - nv (i) */
479 /* x: either a supervariable or an element */
480 /* -----------------------------------------------------------------------
481 */
482 /* LOCAL POINTERS: */
483 /* -----------------------------------------------------------------------
484 */
485 /* Any parameter (pe (...) or pfree) or local variable */
486 /* starting with "p" (for Pointer) is an index into iw, */
487 /* and all indices into iw use variables starting with */
488 /* "p." The only exception to this rule is the iwlen */
489 /* input argument. */
490 /* p: pointer into lots of things */
491 /* p1: pe (i) for some variable i (start of element list) */
492 /* p2: pe (i) + elen (i) - 1 for some var. i (end of el. list)
493 */
494 /* p3: index of first supervariable in clean list */
495 /* pdst: destination pointer, for compression */
496 /* pend: end of memory to compress */
497 /* pj: pointer into an element or variable */
498 /* pme: pointer into the current element (pme1...pme2) */
499 /* pme1: the current element, me, is stored in iw (pme1...pme2) */
500 /* pme2: the end of the current element */
501 /* pn: pointer into a "clean" variable, also used to compress */
502 /* psrc: source pointer, for compression */
503 /* -----------------------------------------------------------------------
504 */
505 /* FUNCTIONS CALLED: */
506 /* -----------------------------------------------------------------------
507 */
508 /* =======================================================================
509 */
510 /* INITIALIZATIONS */
511 /* =======================================================================
512 */
513 /* Parameter adjustments */
514 --w;
515 --degree;
516 --elen;
517 --head;
518 --last;
519 --next;
520 --nv;
521 --len;
522 --iw;
523 --pe;
524
525 /* Function Body */
526 wflg = 2;
527 mindeg = 1;
528 *ncmpa = 0;
529 nel = 0;
530 /* Computing MAX */
531 i__1 = 1, i__2 = *n - 1;
532 hmod = max(i__1,i__2);
533 dmax_ = 0;
534 wbig = *iovflo - *n;
535 mem = *pfree - 1;
536 maxmem = mem;
537 i__1 = *n;
538 for (i = 1; i <= i__1; ++i) {
539 last[i] = 0;
540 head[i] = 0;
541 nv[i] = 1;
542 w[i] = 1;
543 elen[i] = 0;
544 degree[i] = len[i];
545 /* L10: */
546 }
547 /* ----------------------------------------------------------------
548 */
549 /* initialize degree lists and eliminate rows with no off-diag. nz.
550 */
551 /* ----------------------------------------------------------------
552 */
553 i__1 = *n;
554 for (i = 1; i <= i__1; ++i) {
555 deg = degree[i];
556 /* include the diagonal in the true degree */
557 ++deg;
558 degree[i] = deg;
559 if (deg > 1) {
560 /* --------------------------------------------------
561 -------- */
562 /* place i in the degree list corresponding to its deg
563 ree */
564 /* --------------------------------------------------
565 -------- */
566 inext = head[deg];
567 if (inext != 0) {
568 last[inext] = i;
569 }
570 next[i] = inext;
571 head[deg] = i;
572 } else {
573 /* --------------------------------------------------
574 -------- */
575 /* we have a variable that can be eliminated at once b
576 ecause */
577 /* there is no off-diagonal non-zero in its row. */
578 /* --------------------------------------------------
579 -------- */
580 degree[i] = 0;
581 ++nel;
582 elen[i] = -nel;
583 pe[i] = 0;
584 w[i] = 0;
585 }
586 /* L20: */
587 }
588 /* =======================================================================
589 */
590 /* WHILE (selecting pivots) DO */
591 /* =======================================================================
592 */
593 L30:
594 if (nel < *n) {
595 /* ==================================================================
596 ===== */
597 /* GET PIVOT OF MINIMUM DEGREE */
598 /* ==================================================================
599 ===== */
600 /* ---------------------------------------------------------
601 ---- */
602 /* find next supervariable for elimination */
603 /* ---------------------------------------------------------
604 ---- */
605 i__1 = *n;
606 for (deg = mindeg; deg <= i__1; ++deg) {
607 me = head[deg];
608 if (me > 0) {
609 goto L50;
610 }
611 /* L40: */
612 }
613 L50:
614 mindeg = deg;
615 /* ---------------------------------------------------------
616 ---- */
617 /* remove chosen variable from link list */
618 /* ---------------------------------------------------------
619 ---- */
620 inext = next[me];
621 if (inext != 0) {
622 last[inext] = 0;
623 }
624 head[deg] = inext;
625 /* ---------------------------------------------------------
626 ---- */
627 /* me represents the elimination of pivots nel+1 to nel+nv(me
628 ). */
629 /* place me itself as the first in this set. It will be move
630 d */
631 /* to the nel+nv(me) position when the permutation vectors ar
632 e */
633 /* computed. */
634 /* ---------------------------------------------------------
635 ---- */
636 elenme = elen[me];
637 elen[me] = -(nel + 1);
638 nvpiv = nv[me];
639 nel += nvpiv;
640 /* ==================================================================
641 ===== */
642 /* CONSTRUCT NEW ELEMENT */
643 /* ==================================================================
644 ===== */
645 /* ---------------------------------------------------------
646 ---- */
647 /* At this point, me is the pivotal supervariable. It will b
648 e */
649 /* converted into the current element. Scan list of the */
650 /* pivotal supervariable, me, setting tree pointers and */
651 /* constructing new list of supervariables for the new elemen
652 t, */
653 /* me. p is a pointer to the current position in the old lis
654 t. */
655 /* ---------------------------------------------------------
656 ---- */
657 /* flag the variable "me" as being in Lme by negating nv (me)
658 */
659 nv[me] = -nvpiv;
660 degme = 0;
661 if (elenme == 0) {
662 /* --------------------------------------------------
663 -------- */
664 /* construct the new element in place */
665 /* --------------------------------------------------
666 -------- */
667 pme1 = pe[me];
668 pme2 = pme1 - 1;
669 i__1 = pme1 + len[me] - 1;
670 for (p = pme1; p <= i__1; ++p) {
671 i = iw[p];
672 nvi = nv[i];
673 if (nvi > 0) {
674 /* ------------------------------------
675 ---------------- */
676 /* i is a principal variable not yet pla
677 ced in Lme. */
678 /* store i in new list */
679 /* ------------------------------------
680 ---------------- */
681 degme += nvi;
682 /* flag i as being in Lme by negating nv
683 (i) */
684 nv[i] = -nvi;
685 ++pme2;
686 iw[pme2] = i;
687 /* ------------------------------------
688 ---------------- */
689 /* remove variable i from degree list.
690 */
691 /* ------------------------------------
692 ---------------- */
693 ilast = last[i];
694 inext = next[i];
695 if (inext != 0) {
696 last[inext] = ilast;
697 }
698 if (ilast != 0) {
699 next[ilast] = inext;
700 } else {
701 /* i is at the head of the degree
702 list */
703 head[degree[i]] = inext;
704 }
705 }
706 /* L60: */
707 }
708 /* this element takes no new memory in iw: */
709 newmem = 0;
710 } else {
711 /* --------------------------------------------------
712 -------- */
713 /* construct the new element in empty space, iw (pfree
714 ...) */
715 /* --------------------------------------------------
716 -------- */
717 p = pe[me];
718 pme1 = *pfree;
719 slenme = len[me] - elenme;
720 i__1 = elenme + 1;
721 for (knt1 = 1; knt1 <= i__1; ++knt1) {
722 if (knt1 > elenme) {
723 /* search the supervariables in me. */
724 e = me;
725 pj = p;
726 ln = slenme;
727 } else {
728 /* search the elements in me. */
729 e = iw[p];
730 ++p;
731 pj = pe[e];
732 ln = len[e];
733 }
734 /* -------------------------------------------
735 ------------ */
736 /* search for different supervariables and add
737 them to the */
738 /* new list, compressing when necessary. this l
739 oop is */
740 /* executed once for each element in the list a
741 nd once for */
742 /* all the supervariables in the list. */
743 /* -------------------------------------------
744 ------------ */
745 i__2 = ln;
746 for (knt2 = 1; knt2 <= i__2; ++knt2) {
747 i = iw[pj];
748 ++pj;
749 nvi = nv[i];
750 if (nvi > 0) {
751 /* -----------------------------
752 -------------------- */
753 /* compress iw, if necessary */
754 /* -----------------------------
755 -------------------- */
756 if (*pfree > *iwlen) {
757 /* prepare for compressing
758 iw by adjusting */
759 /* pointers and lengths so
760 that the lists being */
761 /* searched in the inner a
762 nd outer loops contain */
763 /* only the remaining entr
764 ies. */
765 pe[me] = p;
766 len[me] -= knt1;
767 if (len[me] == 0) {
768 /* nothing left of
769 supervariable me */
770 pe[me] = 0;
771 }
772 pe[e] = pj;
773 len[e] = ln - knt2;
774 if (len[e] == 0) {
775 /* nothing left of
776 element e */
777 pe[e] = 0;
778 }
779 ++(*ncmpa);
780 /* store first item in pe
781 */
782 /* set first entry to -ite
783 m */
784 i__3 = *n;
785 for (j = 1; j <= i__3; ++j) {
786 pn = pe[j];
787 if (pn > 0) {
788 pe[j] = iw[pn];
789 iw[pn] = -j;
790 }
791 /* L70: */
792 }
793 /* psrc/pdst point to sour
794 ce/destination */
795 pdst = 1;
796 psrc = 1;
797 pend = pme1 - 1;
798 /* while loop: */
799 L80:
800 if (psrc <= pend) {
801 /* search for next
802 negative entry */
803 j = -iw[psrc];
804 ++psrc;
805 if (j > 0) {
806 iw[pdst] = pe[j];
807 pe[j] = pdst;
808 ++pdst;
809 /* copy from
810 source to destination */
811 lenj = len[j];
812 i__3 = lenj - 2;
813 for (knt3 = 0; knt3 <= i__3; ++knt3) {
814 iw[pdst + knt3] = iw[psrc + knt3];
815 /* L90: */
816 }
817 pdst = pdst + lenj - 1;
818 psrc = psrc + lenj - 1;
819 }
820 goto L80;
821 }
822 /* move the new partially-
823 constructed element */
824 p1 = pdst;
825 i__3 = *pfree - 1;
826 for (psrc = pme1; psrc <= i__3; ++psrc) {
827 iw[pdst] = iw[psrc];
828 ++pdst;
829 /* L100: */
830 }
831 pme1 = p1;
832 *pfree = pdst;
833 pj = pe[e];
834 p = pe[me];
835 }
836 /* -----------------------------
837 -------------------- */
838 /* i is a principal variable not
839 yet placed in Lme */
840 /* store i in new list */
841 /* -----------------------------
842 -------------------- */
843 degme += nvi;
844 /* flag i as being in Lme by nega
845 ting nv (i) */
846 nv[i] = -nvi;
847 iw[*pfree] = i;
848 ++(*pfree);
849 /* -----------------------------
850 -------------------- */
851 /* remove variable i from degree
852 link list */
853 /* -----------------------------
854 -------------------- */
855 ilast = last[i];
856 inext = next[i];
857 if (inext != 0) {
858 last[inext] = ilast;
859 }
860 if (ilast != 0) {
861 next[ilast] = inext;
862 } else {
863 /* i is at the head of the
864 degree list */
865 head[degree[i]] = inext;
866 }
867 }
868 /* L110: */
869 }
870 if (e != me) {
871 /* set tree pointer and flag to indicate
872 element e is */
873 /* absorbed into new element me (the par
874 ent of e is me) */
875 pe[e] = -me;
876 w[e] = 0;
877 }
878 /* L120: */
879 }
880 pme2 = *pfree - 1;
881 /* this element takes newmem new memory in iw (possibl
882 y zero) */
883 newmem = *pfree - pme1;
884 mem += newmem;
885 maxmem = max(maxmem,mem);
886 }
887 /* ---------------------------------------------------------
888 ---- */
889 /* me has now been converted into an element in iw (pme1..pme
890 2) */
891 /* ---------------------------------------------------------
892 ---- */
893 /* degme holds the external degree of new element */
894 degree[me] = degme;
895 pe[me] = pme1;
896 len[me] = pme2 - pme1 + 1;
897 /* ---------------------------------------------------------
898 ---- */
899 /* make sure that wflg is not too large. With the current */
900 /* value of wflg, wflg+n must not cause integer overflow */
901 /* ---------------------------------------------------------
902 ---- */
903 if (wflg >= wbig) {
904 i__1 = *n;
905 for (x = 1; x <= i__1; ++x) {
906 if (w[x] != 0) {
907 w[x] = 1;
908 }
909 /* L130: */
910 }
911 wflg = 2;
912 }
913 /* ==================================================================
914 ===== */
915 /* COMPUTE (w (e) - wflg) = |Le\Lme| FOR ALL ELEMENTS */
916 /* ==================================================================
917 ===== */
918 /* ---------------------------------------------------------
919 ---- */
920 /* Scan 1: compute the external degrees of previous elements
921 */
922 /* with respect to the current element. That is: */
923 /* (w (e) - wflg) = |Le \ Lme| */
924 /* for each element e that appears in any supervariable in Lm
925 e. */
926 /* The notation Le refers to the pattern (list of */
927 /* supervariables) of a previous element e, where e is not ye
928 t */
929 /* absorbed, stored in iw (pe (e) + 1 ... pe (e) + iw (pe (e)
930 )). */
931 /* The notation Lme refers to the pattern of the current elem
932 ent */
933 /* (stored in iw (pme1..pme2)). If (w (e) - wflg) becomes
934 */
935 /* zero, then the element e will be absorbed in scan 2. */
936 /* ---------------------------------------------------------
937 ---- */
938 i__1 = pme2;
939 for (pme = pme1; pme <= i__1; ++pme) {
940 i = iw[pme];
941 eln = elen[i];
942 if (eln > 0) {
943 /* note that nv (i) has been negated to denote
944 i in Lme: */
945 nvi = -nv[i];
946 wnvi = wflg - nvi;
947 i__2 = pe[i] + eln - 1;
948 for (p = pe[i]; p <= i__2; ++p) {
949 e = iw[p];
950 we = w[e];
951 if (we >= wflg) {
952 /* unabsorbed element e has been
953 seen in this loop */
954 we -= nvi;
955 } else if (we != 0) {
956 /* e is an unabsorbed element */
957 /* this is the first we have seen
958 e in all of Scan 1 */
959 we = degree[e] + wnvi;
960 }
961 w[e] = we;
962 /* L140: */
963 }
964 }
965 /* L150: */
966 }
967 /* ==================================================================
968 ===== */
969 /* DEGREE UPDATE AND ELEMENT ABSORPTION */
970 /* ==================================================================
971 ===== */
972 /* ---------------------------------------------------------
973 ---- */
974 /* Scan 2: for each i in Lme, sum up the degree of Lme (whic
975 h */
976 /* is degme), plus the sum of the external degrees of each Le
977 */
978 /* for the elements e appearing within i, plus the */
979 /* supervariables in i. Place i in hash list. */
980 /* ---------------------------------------------------------
981 ---- */
982 i__1 = pme2;
983 for (pme = pme1; pme <= i__1; ++pme) {
984 i = iw[pme];
985 p1 = pe[i];
986 p2 = p1 + elen[i] - 1;
987 pn = p1;
988 hash = 0;
989 deg = 0;
990 /* --------------------------------------------------
991 -------- */
992 /* scan the element list associated with supervariable
993 i */
994 /* --------------------------------------------------
995 -------- */
996 /* UMFPACK/MA38-style approximate degree: */
997 i__2 = p2;
998 for (p = p1; p <= i__2; ++p) {
999 e = iw[p];
1000 we = w[e];
1001 if (we != 0) {
1002 /* e is an unabsorbed element */
1003 deg = deg + we - wflg;
1004 iw[pn] = e;
1005 ++pn;
1006 hash += e;
1007 }
1008 /* L160: */
1009 }
1010 /* count the number of elements in i (including me):
1011 */
1012 elen[i] = pn - p1 + 1;
1013 /* --------------------------------------------------
1014 -------- */
1015 /* scan the supervariables in the list associated with
1016 i */
1017 /* --------------------------------------------------
1018 -------- */
1019 p3 = pn;
1020 i__2 = p1 + len[i] - 1;
1021 for (p = p2 + 1; p <= i__2; ++p) {
1022 j = iw[p];
1023 nvj = nv[j];
1024 if (nvj > 0) {
1025 /* j is unabsorbed, and not in Lme. */
1026 /* add to degree and add to new list */
1027 deg += nvj;
1028 iw[pn] = j;
1029 ++pn;
1030 hash += j;
1031 }
1032 /* L170: */
1033 }
1034 /* --------------------------------------------------
1035 -------- */
1036 /* update the degree and check for mass elimination */
1037 /* --------------------------------------------------
1038 -------- */
1039 if (elen[i] == 1 && p3 == pn) {
1040 /* -------------------------------------------
1041 ------------ */
1042 /* mass elimination */
1043 /* -------------------------------------------
1044 ------------ */
1045 /* There is nothing left of this node except fo
1046 r an */
1047 /* edge to the current pivot element. elen (i)
1048 is 1, */
1049 /* and there are no variables adjacent to node
1050 i. */
1051 /* Absorb i into the current pivot element, me.
1052 */
1053 pe[i] = -me;
1054 nvi = -nv[i];
1055 degme -= nvi;
1056 nvpiv += nvi;
1057 nel += nvi;
1058 nv[i] = 0;
1059 elen[i] = 0;
1060 } else {
1061 /* -------------------------------------------
1062 ------------ */
1063 /* update the upper-bound degree of i */
1064 /* -------------------------------------------
1065 ------------ */
1066 /* the following degree does not yet include th
1067 e size */
1068 /* of the current element, which is added later
1069 : */
1070 /* Computing MIN */
1071 i__2 = degree[i];
1072 degree[i] = min(i__2,deg);
1073 /* -------------------------------------------
1074 ------------ */
1075 /* add me to the list for i */
1076 /* -------------------------------------------
1077 ------------ */
1078 /* move first supervariable to end of list */
1079 iw[pn] = iw[p3];
1080 /* move first element to end of element part of
1081 list */
1082 iw[p3] = iw[p1];
1083 /* add new element to front of list. */
1084 iw[p1] = me;
1085 /* store the new length of the list in len (i)
1086 */
1087 len[i] = pn - p1 + 1;
1088 /* -------------------------------------------
1089 ------------ */
1090 /* place in hash bucket. Save hash key of i in
1091 last (i). */
1092 /* -------------------------------------------
1093 ------------ */
1094 hash = hash % hmod + 1;
1095 j = head[hash];
1096 if (j <= 0) {
1097 /* the degree list is empty, hash head i
1098 s -j */
1099 next[i] = -j;
1100 head[hash] = -i;
1101 } else {
1102 /* degree list is not empty */
1103 /* use last (head (hash)) as hash head
1104 */
1105 next[i] = last[j];
1106 last[j] = i;
1107 }
1108 last[i] = hash;
1109 }
1110 /* L180: */
1111 }
1112 degree[me] = degme;
1113 /* ---------------------------------------------------------
1114 ---- */
1115 /* Clear the counter array, w (...), by incrementing wflg. */
1116 /* ---------------------------------------------------------
1117 ---- */
1118 dmax_ = max(dmax_,degme);
1119 wflg += dmax_;
1120 /* make sure that wflg+n does not cause integer overflow */
1121 if (wflg >= wbig) {
1122 i__1 = *n;
1123 for (x = 1; x <= i__1; ++x) {
1124 if (w[x] != 0) {
1125 w[x] = 1;
1126 }
1127 /* L190: */
1128 }
1129 wflg = 2;
1130 }
1131 /* at this point, w (1..n) .lt. wflg holds */
1132 /* ==================================================================
1133 ===== */
1134 /* SUPERVARIABLE DETECTION */
1135 /* ==================================================================
1136 ===== */
1137 i__1 = pme2;
1138 for (pme = pme1; pme <= i__1; ++pme) {
1139 i = iw[pme];
1140 if (nv[i] < 0) {
1141 /* i is a principal variable in Lme */
1142 /* -------------------------------------------
1143 ------------ */
1144 /* examine all hash buckets with 2 or more vari
1145 ables. We */
1146 /* do this by examing all unique hash keys for
1147 super- */
1148 /* variables in the pattern Lme of the current
1149 element, me */
1150 /* -------------------------------------------
1151 ------------ */
1152 hash = last[i];
1153 /* let i = head of hash bucket, and empty the h
1154 ash bucket */
1155 j = head[hash];
1156 if (j == 0) {
1157 goto L250;
1158 }
1159 if (j < 0) {
1160 /* degree list is empty */
1161 i = -j;
1162 head[hash] = 0;
1163 } else {
1164 /* degree list is not empty, restore las
1165 t () of head */
1166 i = last[j];
1167 last[j] = 0;
1168 }
1169 if (i == 0) {
1170 goto L250;
1171 }
1172 /* while loop: */
1173 L200:
1174 if (next[i] != 0) {
1175 /* ------------------------------------
1176 ---------------- */
1177 /* this bucket has one or more variables
1178 following i. */
1179 /* scan all of them to see if i can abso
1180 rb any entries */
1181 /* that follow i in hash bucket. Scatte
1182 r i into w. */
1183 /* ------------------------------------
1184 ---------------- */
1185 ln = len[i];
1186 eln = elen[i];
1187 /* do not flag the first element in the
1188 list (me) */
1189 i__2 = pe[i] + ln - 1;
1190 for (p = pe[i] + 1; p <= i__2; ++p) {
1191 w[iw[p]] = wflg;
1192 /* L210: */
1193 }
1194 /* ------------------------------------
1195 ---------------- */
1196 /* scan every other entry j following i
1197 in bucket */
1198 /* ------------------------------------
1199 ---------------- */
1200 jlast = i;
1201 j = next[i];
1202 /* while loop: */
1203 L220:
1204 if (j != 0) {
1205 /* -----------------------------
1206 -------------------- */
1207 /* check if j and i have identica
1208 l nonzero pattern */
1209 /* -----------------------------
1210 -------------------- */
1211 if (len[j] != ln) {
1212 /* i and j do not have sam
1213 e size data structure */
1214 goto L240;
1215 }
1216 if (elen[j] != eln) {
1217 /* i and j do not have sam
1218 e number of adjacent el */
1219 goto L240;
1220 }
1221 /* do not flag the first element
1222 in the list (me) */
1223 i__2 = pe[j] + ln - 1;
1224 for (p = pe[j] + 1; p <= i__2; ++p) {
1225 if (w[iw[p]] != wflg) {
1226 /* an entry (iw(p))
1227 is in j but not in i */
1228 goto L240;
1229 }
1230 /* L230: */
1231 }
1232 /* -----------------------------
1233 -------------------- */
1234 /* found it! j can be absorbed i
1235 nto i */
1236 /* -----------------------------
1237 -------------------- */
1238 pe[j] = -i;
1239 /* both nv (i) and nv (j) are neg
1240 ated since they */
1241 /* are in Lme, and the absolute v
1242 alues of each */
1243 /* are the number of variables in
1244 i and j: */
1245 nv[i] += nv[j];
1246 nv[j] = 0;
1247 elen[j] = 0;
1248 /* delete j from hash bucket */
1249 j = next[j];
1250 next[jlast] = j;
1251 goto L220;
1252 /* -----------------------------
1253 -------------------- */
1254 L240:
1255 /* j cannot be absorbed into i */
1256 /* -----------------------------
1257 -------------------- */
1258 jlast = j;
1259 j = next[j];
1260 goto L220;
1261 }
1262 /* ------------------------------------
1263 ---------------- */
1264 /* no more variables can be absorbed int
1265 o i */
1266 /* go to next i in bucket and clear flag
1267 array */
1268 /* ------------------------------------
1269 ---------------- */
1270 ++wflg;
1271 i = next[i];
1272 if (i != 0) {
1273 goto L200;
1274 }
1275 }
1276 }
1277 L250:
1278 ;
1279 }
1280 /* ==================================================================
1281 ===== */
1282 /* RESTORE DEGREE LISTS AND REMOVE NONPRINCIPAL SUPERVAR. FROM ELEMEN
1283 T */
1284 /* ==================================================================
1285 ===== */
1286 p = pme1;
1287 nleft = *n - nel;
1288 i__1 = pme2;
1289 for (pme = pme1; pme <= i__1; ++pme) {
1290 i = iw[pme];
1291 nvi = -nv[i];
1292 if (nvi > 0) {
1293 /* i is a principal variable in Lme */
1294 /* restore nv (i) to signify that i is principa
1295 l */
1296 nv[i] = nvi;
1297 /* -------------------------------------------
1298 ------------ */
1299 /* compute the true degree (add size of current
1300 element) */
1301 /* -------------------------------------------
1302 ------------ */
1303 /* Computing MIN */
1304 i__2 = degree[i] + degme;
1305 deg = min(i__2,nleft);
1306 /* -------------------------------------------
1307 ------------ */
1308 /* place the supervariable at the head of the d
1309 egree list */
1310 /* -------------------------------------------
1311 ------------ */
1312 inext = head[deg];
1313 if (inext != 0) {
1314 last[inext] = i;
1315 }
1316 next[i] = inext;
1317 last[i] = 0;
1318 head[deg] = i;
1319 /* -------------------------------------------
1320 ------------ */
1321 /* save the new degree, and find the minimum de
1322 gree */
1323 /* -------------------------------------------
1324 ------------ */
1325 mindeg = min(mindeg,deg);
1326 degree[i] = deg;
1327 /* -------------------------------------------
1328 ------------ */
1329 /* place the supervariable in the element patte
1330 rn */
1331 /* -------------------------------------------
1332 ------------ */
1333 iw[p] = i;
1334 ++p;
1335 }
1336 /* L260: */
1337 }
1338 /* ==================================================================
1339 ===== */
1340 /* FINALIZE THE NEW ELEMENT */
1341 /* ==================================================================
1342 ===== */
1343 nv[me] = nvpiv + degme;
1344 /* nv (me) is now the degree of pivot (including diagonal par
1345 t) */
1346 /* save the length of the list for the new element me */
1347 len[me] = p - pme1;
1348 if (len[me] == 0) {
1349 /* there is nothing left of the current pivot element
1350 */
1351 pe[me] = 0;
1352 w[me] = 0;
1353 }
1354 if (newmem != 0) {
1355 /* element was not constructed in place: deallocate pa
1356 rt */
1357 /* of it (final size is less than or equal to newmem,
1358 */
1359 /* since newly nonprincipal variables have been remove
1360 d). */
1361 *pfree = p;
1362 mem = mem - newmem + len[me];
1363 }
1364 /* ==================================================================
1365 ===== */
1366 /* END WHILE (selecting pivots) */
1367 goto L30;
1368 }
1369 /* =======================================================================
1370 */
1371 /* =======================================================================
1372 */
1373 /* COMPUTE THE PERMUTATION VECTORS */
1374 /* =======================================================================
1375 */
1376 /* ----------------------------------------------------------------
1377 */
1378 /* The time taken by the following code is O(n). At this */
1379 /* point, elen (e) = -k has been done for all elements e, */
1380 /* and elen (i) = 0 has been done for all nonprincipal */
1381 /* variables i. At this point, there are no principal */
1382 /* supervariables left, and all elements are absorbed. */
1383 /* ----------------------------------------------------------------
1384 */
1385 /* ----------------------------------------------------------------
1386 */
1387 /* compute the ordering of unordered nonprincipal variables */
1388 /* ----------------------------------------------------------------
1389 */
1390 i__1 = *n;
1391 for (i = 1; i <= i__1; ++i) {
1392 if (elen[i] == 0) {
1393 /* --------------------------------------------------
1394 -------- */
1395 /* i is an un-ordered row. Traverse the tree from i u
1396 ntil */
1397 /* reaching an element, e. The element, e, was the */
1398 /* principal supervariable of i and all nodes in the p
1399 ath */
1400 /* from i to when e was selected as pivot. */
1401 /* --------------------------------------------------
1402 -------- */
1403 j = -pe[i];
1404 /* while (j is a variable) do: */
1405 L270:
1406 if (elen[j] >= 0) {
1407 j = -pe[j];
1408 goto L270;
1409 }
1410 e = j;
1411 /* --------------------------------------------------
1412 -------- */
1413 /* get the current pivot ordering of e */
1414 /* --------------------------------------------------
1415 -------- */
1416 k = -elen[e];
1417 /* --------------------------------------------------
1418 -------- */
1419 /* traverse the path again from i to e, and compress t
1420 he */
1421 /* path (all nodes point to e). Path compression allo
1422 ws */
1423 /* this code to compute in O(n) time. Order the unord
1424 ered */
1425 /* nodes in the path, and place the element e at the e
1426 nd. */
1427 /* --------------------------------------------------
1428 -------- */
1429 j = i;
1430 /* while (j is a variable) do: */
1431 L280:
1432 if (elen[j] >= 0) {
1433 jnext = -pe[j];
1434 pe[j] = -e;
1435 if (elen[j] == 0) {
1436 /* j is an unordered row */
1437 elen[j] = k;
1438 ++k;
1439 }
1440 j = jnext;
1441 goto L280;
1442 }
1443 /* leave elen (e) negative, so we know it is an elemen
1444 t */
1445 elen[e] = -k;
1446 }
1447 /* L290: */
1448 }
1449 /* ----------------------------------------------------------------
1450 */
1451 /* reset the inverse permutation (elen (1..n)) to be positive, */
1452 /* and compute the permutation (last (1..n)). */
1453 /* ----------------------------------------------------------------
1454 */
1455 i__1 = *n;
1456 for (i = 1; i <= i__1; ++i) {
1457 k = (i__2 = elen[i], abs(i__2));
1458 last[k] = i;
1459 elen[i] = k;
1460 /* L300: */
1461 }
1462 /* =======================================================================
1463 */
1464 /* RETURN THE MEMORY USAGE IN IW */
1465 /* =======================================================================
1466 */
1467 /* If maxmem is less than or equal to iwlen, then no compressions */
1468 /* occurred, and iw (maxmem+1 ... iwlen) was unused. Otherwise */
1469 /* compressions did occur, and iwlen would have had to have been */
1470 /* greater than or equal to maxmem for no compressions to occur. */
1471 /* Return the value of maxmem in the pfree argument. */
1472 *pfree = maxmem;
1473 return 0;
1474 } /* amdatr_ */
1475
1476