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