1 /*****************************************************************************
2 /
3 / SPACE (SPArse Cholesky Elimination) Library: gelim.c
4 /
5 / author J"urgen Schulze, University of Paderborn
6 / created 01jan10
7 /
8 / This file contains functions dealing with the elimination graph object
9 /
10 ******************************************************************************
11
12 Data type: struct gelim
13 graph_t *G; pointer to graph object
14 int maxedges; max number of edges that can be stored
15 int *len; length of v's adjacency list
16 int *elen; number of elements adjacent to v
17 int *parent; parent in front tree / representative of v
18 int *degree; boundary size / (approximate) degree
19 int *score; holds the score of uneliminated vertex v
20 Comments:
21 o Structure used to hold the elimination graphs of a bottom-up ordering
22 o G->totvwght: total weight of all uneliminated vertices
23 o G->xadj[v] = -1 => there is no adjacency list for variable/element v
24 => variable v has degree 0 (in this case G->vwght[v] > 0)
25 => variable v istinguishable/removed by mass elimination
26 or element v has been absorbed (in this case G->vwght[v] = 0)
27 o G->vwght[v]: weight of the princial variable v; if v becomes an element,
28 weight[v] remains unchanged for the rest of the elim. process
29 = 0 => variable v is nonprincipal/removed by mass elimination
30 o len[v], elen[v]: the adjacency list of vertex/element v contains len[v]
31 entries; the first elen[v] entries are elements
32 (if v is an element, then elen[v] = 0 will hold)
33 o parent[v]: for an (absorbed) element, parent[v] points to the parent of
34 element v in the front tree; for an indistinguishable vertex,
35 parent[v] points to its representative vertex (which may have
36 also found to be indistinguishable to another one)
37 o degree[v]: for an uneliminated vertex, the (approximate) degree in Gelim;
38 for an element, the weight of its boundary (i.e. degree[v]
39 gives the exakt degree of v at the time of its elimination)
40 o score[v]: vertices are eliminated according to their score value >= 0;
41 additionally, the score vector is used to represent the status
42 of a node in the actual stage:
43 -1, iff variable v will be eliminated in an upcomming stage
44 -2, iff variable v is nonprincipal/removed by mass elim.
45 -3, iff variable v has been eliminated and now forms an element
46 -4, iff element v has been absorbed
47 Methods in lib/gelim.c
48 - Gelim = newElimGraph(int nvtx, int nedges);
49 - void freeElimGraph(gelim_t *Gelim);
50 - void printElimGraph(gelim_t *Gelim);
51 - Gelim = setupElimGraph(graph_t *G);
52 o allocates memory for the elimination graph by calling newElimGraph and
53 initializes the vectors, i.e. len[u] = xadj[u+1]-xadj[u]; elen[u] = 0;
54 parent[u] = -1; degree[u] = exact (external) degree of vertex u;
55 score[u] = -1; xadj[u] = -1, if len[u] = 0
56 - int crunchElimGraph(gelim_t *Gelim);
57 o tries to compress the adjacency vector
58 on success the function return TRUE, otherwise FALSE
59 - void buildElement(gelim_t *Gelim, int me);
60 o turns variable me into an element; if me is an leaf, the element is
61 constructed in-place, otherwise its adjacency list is appended to G
62 o all relevant vectors are updated, i.e.
63 vwght[me] = 0, degree[me] = |Lme|, score[me] = -3
64 for all neighboring elements: parent[e] = me, score[e] = -4
65 - void updateAdjncy(gelim_t *Gelim, int *reachset, int nreach, int *tmp,
66 int *pflag);
67 o updates the adjacency structure of all vertices in reachset
68 IMPORTANT REQUIREMENTS:
69 (1) all values stored in tmp[u] are smaller than *pflag
70 - void findIndNodes(gelim_t *Gelim, int *reachset, int nreach, int *bin,
71 int *next, int *tmp, int *pflag);
72 o searches reachset for indistinguishable vertices
73 IMPORTANT REQUIREMENTS:
74 (1) the adjacency lists of all vertices in reachset have been updated
75 by a call to updateAdjncy
76 (2) bin[i] = -1 for all 0 <= i < G->nvtx
77 (3) all values stored in tmp[u] are smaller than *pflag
78 o on return bin[i] = -1 holds again
79 - void updateDegree(gelim_t *Gelim, int *reachset, int nreach, int *bin);
80 o computes new approximate degrees for all vertices in reachset
81 IMPORTANT REQUIREMENTS:
82 (1) the adjacency lists of all vertices in reachset have been updated
83 by a call to updateAdjncy
84 (2) the boundary size of each newly formed element has been computed
85 (3) bin[i] = -1 for all 0 <= i < G->nvtx
86 o on return bin[i] = -1 holds again
87 - void updateScore(gelim_t *Gelim, int *reachset, int nreach, int scoretype,
88 int *bin);
89 o updates the score of all vertices in reachset
90 IMPORTANT REQUIREMENTS:
91 (1) the approximate degrees are correctly computed (by updateDegree)
92 (2) bin[i] = -1 for all 0 <= i < G->nvtx
93 o on return bin[i] = -1 holds again
94 - T = extractElimTree(gelim_t *Gelim);
95 o uses the status of the nodes (stored in the score vector) and the
96 parent vector to set up the elimination tree T; vectors T->ncolfactor
97 and T->ncolupdate are initialized using vectors G->vwght and degree
98
99 ******************************************************************************/
100
101 #include <space.h>
102 /* #define DEBUG */
103
104
105 /*****************************************************************************
106 ******************************************************************************/
107 gelim_t*
newElimGraph(PORD_INT nvtx,PORD_INT nedges)108 newElimGraph(PORD_INT nvtx, PORD_INT nedges)
109 { gelim_t *Gelim;
110
111 mymalloc(Gelim, 1, gelim_t);
112 Gelim->G = newGraph(nvtx, nedges);
113 Gelim->maxedges = nedges;
114
115 mymalloc(Gelim->len, nvtx, PORD_INT);
116 mymalloc(Gelim->elen, nvtx, PORD_INT);
117 mymalloc(Gelim->parent, nvtx, PORD_INT);
118 mymalloc(Gelim->degree, nvtx, PORD_INT);
119 mymalloc(Gelim->score, nvtx, PORD_INT);
120
121 return(Gelim);
122 }
123
124
125 /*****************************************************************************
126 ******************************************************************************/
127 void
freeElimGraph(gelim_t * Gelim)128 freeElimGraph(gelim_t *Gelim)
129 {
130 freeGraph(Gelim->G);
131 free(Gelim->len);
132 free(Gelim->elen);
133 free(Gelim->parent);
134 free(Gelim->degree);
135 free(Gelim->score);
136 free(Gelim);
137 }
138
139
140 /*****************************************************************************
141 ******************************************************************************/
142 void
printElimGraph(gelim_t * Gelim)143 printElimGraph(gelim_t *Gelim)
144 { graph_t *G;
145 PORD_INT count, u, v, i, istart;
146
147 G = Gelim->G;
148 for (u = 0; u < G->nvtx; u++)
149 { istart = G->xadj[u];
150
151 /* ---------------------------------------------------------------
152 case 1: u is a principal variable
153 => vwght[u]: weight of all mapped indistinguishable variables
154 => degree[u]: approximate degree
155 ---------------------------------------------------------------- */
156 if ((Gelim->score[u] == -1) || (Gelim->score[u] >= 0))
157 { printf("--- adjacency list of variable %d (weight %d, degree %d, "
158 "score %d):\n", u, G->vwght[u], Gelim->degree[u],
159 Gelim->score[u]);
160 printf("elements:\n");
161 count = 0;
162 for (i = istart; i < istart + Gelim->elen[u]; i++)
163 { printf("%5d", G->adjncy[i]);
164 if ((++count % 16) == 0)
165 printf("\n");
166 }
167 if ((count % 16) != 0)
168 printf("\n");
169 printf("variables:\n");
170 count = 0;
171 for (i = istart + Gelim->elen[u]; i < istart + Gelim->len[u]; i++)
172 { printf("%5d", G->adjncy[i]);
173 if ((++count % 16) == 0)
174 printf("\n");
175 }
176 if ((count % 16) != 0)
177 printf("\n");
178 }
179
180 /* ---------------------------------------------------------------
181 case 2: u is nonprincipal/removed by mass elimination
182 ---------------------------------------------------------------- */
183 else if (Gelim->score[u] == -2)
184 printf("--- variable %d is nonprincipal/removed by mass elim. "
185 "(parent %d)\n", u, Gelim->parent[u]);
186
187 /* -----------------------------------------------
188 case 3: u is an element:
189 => degree[u]: weight of boundary
190 ----------------------------------------------- */
191 else if (Gelim->score[u] == -3)
192 { printf("--- boundary of element %d (degree %d, score %d):"
193 "\n", u, Gelim->degree[u], Gelim->score[u]);
194 count = 0;
195 for (i = istart; i < istart + Gelim->len[u]; i++)
196 { v = G->adjncy[i];
197 if (G->vwght[v] > 0)
198 { printf("%5d", G->adjncy[i]);
199 if ((++count % 16) == 0)
200 printf("\n");
201 }
202 }
203 if ((count % 16) != 0)
204 printf("\n");
205 }
206
207 /* --------------------------------
208 case 4: u is an absorbed element
209 -------------------------------- */
210 else if (Gelim->score[u] == -4)
211 printf("--- element %d has been absorbed (parent %d)\n", u,
212 Gelim->parent[u]);
213
214 /* ----------------------------------------
215 none of the above cases is true => error
216 ---------------------------------------- */
217 else
218 { fprintf(stderr, "\nError in function printElimGraph\n"
219 " node %d has invalid score %d\n", u, Gelim->score[u]);
220 quit();
221 }
222 }
223 }
224
225
226 /*****************************************************************************
227 ******************************************************************************/
228 gelim_t*
setupElimGraph(graph_t * G)229 setupElimGraph(graph_t *G)
230 { gelim_t *Gelim;
231 PORD_INT *xadj, *adjncy, *vwght, *xadjGelim, *adjncyGelim, *vwghtGelim;
232 PORD_INT *len, *elen, *parent, *degree, *score;
233 PORD_INT nvtx, nedges, deg, u, i, istart, istop;
234
235 nvtx = G->nvtx;
236 nedges = G->nedges;
237 xadj = G->xadj;
238 adjncy = G->adjncy;
239 vwght = G->vwght;
240
241 Gelim = newElimGraph(nvtx, nedges+nvtx);
242 xadjGelim = Gelim->G->xadj;
243 adjncyGelim = Gelim->G->adjncy;
244 vwghtGelim = Gelim->G->vwght;
245 len = Gelim->len;
246 elen = Gelim->elen;
247 parent = Gelim->parent;
248 degree = Gelim->degree;
249 score = Gelim->score;
250
251 /* --------------
252 copy the graph
253 -------------- */
254 Gelim->G->type = G->type;
255 Gelim->G->totvwght = G->totvwght;
256 for (u = 0; u < nvtx; u++)
257 { xadjGelim[u] = xadj[u];
258 vwghtGelim[u] = vwght[u];
259 }
260 xadjGelim[nvtx] = xadj[nvtx];
261 for (i = 0; i < nedges; i++)
262 adjncyGelim[i] = adjncy[i];
263 Gelim->G->nedges = nedges;
264
265 /* ----------------------
266 initialize all vectors
267 ---------------------- */
268 for (u = 0; u < nvtx; u++)
269 { istart = xadj[u];
270 istop = xadj[u+1];
271 len[u] = istop - istart;
272 elen[u] = 0;
273 parent[u] = -1;
274 deg = 0;
275
276 switch(Gelim->G->type) /* compute the external degree of u */
277 { case UNWEIGHTED:
278 deg = len[u];
279 break;
280 case WEIGHTED:
281 for (i = istart; i < istop; i++)
282 deg += vwght[adjncy[i]];
283 break;
284 default:
285 fprintf(stderr, "\nError in function setupElimGraph\n"
286 " unrecognized graph type %d\n", Gelim->G->type);
287 }
288 degree[u] = deg;
289
290 if (len[u] == 0) /* len(u) = 0 => adjncy list of u not in use */
291 xadjGelim[u] = -1; /* mark with -1, otherwise crunchElimGraph fails */
292 score[u] = -1;
293 }
294
295 return(Gelim);
296 }
297
298
299 /*****************************************************************************
300 ******************************************************************************/
301 PORD_INT
crunchElimGraph(gelim_t * Gelim)302 crunchElimGraph(gelim_t *Gelim)
303 { PORD_INT *xadj, *adjncy, *len;
304 PORD_INT nvtx, nedges, u, i, isrc, idest;
305
306 nvtx = Gelim->G->nvtx;
307 nedges = Gelim->G->nedges;
308 xadj = Gelim->G->xadj;
309 adjncy = Gelim->G->adjncy;
310 len = Gelim->len;
311
312 /* ---------------------------------------------
313 mark begining of u's adjacency list by -(u+1)
314 --------------------------------------------- */
315 for (u = 0; u < nvtx; u++)
316 { i = xadj[u]; /* is adjacency list of u still in use? */
317 if (i != -1) /* verify that list is non-empty */
318 { if (len[u] == 0)
319 { fprintf(stderr, "\nError in function crunchElimGraph\n"
320 " adjacency list of node %d is empty\n", u);
321 quit();
322 }
323 xadj[u] = adjncy[i]; /* if so, move first item to xadj[u] */
324 adjncy[i] = -(u+1); /* u's adjacency list is headed by -(u+1) */
325 if (len[u] == 0)
326 printf("error: u %d, len %d\n", u, len[u]);
327 }
328 }
329
330 /* --------------------------
331 crunch all adjacency lists
332 -------------------------- */
333 idest = isrc = 0;
334 while (isrc < Gelim->G->nedges)
335 { u = adjncy[isrc++];
336 if (u < 0) /* a new adjacency list starts here */
337 { u = -u - 1; /* it's the adjacency list of u */
338 adjncy[idest] = xadj[u]; /* first item was stored in xadj[u] */
339 xadj[u] = idest++;
340 for (i = 1; i < len[u]; i++)
341 adjncy[idest++] = adjncy[isrc++];
342 }
343 }
344 Gelim->G->nedges = idest;
345
346 /* ------------------
347 was it successful?
348 ------------------ */
349 if (idest < nedges) return(TRUE);
350 else return (FALSE);
351 }
352
353
354 /*****************************************************************************
355 ******************************************************************************/
356 void
buildElement(gelim_t * Gelim,PORD_INT me)357 buildElement(gelim_t *Gelim, PORD_INT me)
358 { graph_t *G;
359 PORD_INT *xadj, *adjncy, *vwght, *len, *elen, *parent, *degree, *score;
360 PORD_INT degme, elenme, vlenme, mesrcptr, medeststart, medeststart2;
361 PORD_INT medestptr, ln, p, i, j, v, e;
362
363 G = Gelim->G;
364 xadj = G->xadj;
365 adjncy = G->adjncy;
366 vwght = G->vwght;
367 len = Gelim->len;
368 elen = Gelim->elen;
369 parent = Gelim->parent;
370 degree = Gelim->degree;
371 score = Gelim->score;
372
373 /* ---------------------------------
374 construct boundary of element Lme
375 --------------------------------- */
376 degme = 0;
377 G->totvwght -= vwght[me]; /* me eliminated => reduce weight of Gelim */
378 vwght[me] = -vwght[me];
379 score[me] = -3; /* variable me becomes an element */
380
381 elenme = elen[me];
382 vlenme = len[me] - elenme;
383 mesrcptr = xadj[me];
384
385 /* -----------------------------------------------------------
386 if me is a leaf => its boundary can be constructed in-place
387 ----------------------------------------------------------- */
388 if (elenme == 0)
389 { medeststart = xadj[me]; /* Lme overwrites old variable */
390 medestptr = medeststart; /* boundary of Lme starts here */
391 for (i = 0; i < vlenme; i++)
392 { v = adjncy[mesrcptr++];
393 if (vwght[v] > 0) /* v not yet placed in boundary */
394 { degme += vwght[v]; /* increase size of Lme */
395 vwght[v] = -vwght[v]; /* flag v as being in Lme */
396 adjncy[medestptr++] = v;
397 }
398 }
399 }
400
401 /* -------------------------------------------------------------------
402 me is not a leaf => its boundary must be constructed in empty space
403 ------------------------------------------------------------------- */
404 else
405 { medeststart = G->nedges; /* Lme appended to graph */
406 medestptr = medeststart; /* boundary of Lme starts here */
407 for (i = 0; i <= elenme; i++)
408 { if (i < elenme) /* working on elements */
409 { len[me]--;
410 e = adjncy[mesrcptr++]; /* merge boundary of element e with Lme */
411 p = xadj[e]; /* adjacency list of e starts here */
412 ln = len[e];
413 }
414 else
415 { e = me; /* merge uncovered variables with Lme */
416 p = mesrcptr; /* variables start here */
417 ln = vlenme;
418 }
419 for (j = 0; j < ln; j++)
420 { len[e]--; /* pick next variable, decrease length */
421 v = adjncy[p++];
422 if (vwght[v] > 0)
423 { degme += vwght[v]; /* increase size of Lme */
424 vwght[v] = -vwght[v]; /* flag v as being in Lme */
425
426 /* ------------------------------------------
427 add v to Lme, compress adjncy if necessary
428 ------------------------------------------ */
429 if (medestptr == Gelim->maxedges)
430 { if (len[me] == 0) xadj[me] = -1;
431 else xadj[me] = mesrcptr;
432 if (len[e] == 0) xadj[e] = -1;
433 else xadj[e] = p;
434
435 /* crunch adjacency list -- !!!we need more memory!!! */
436 if (!crunchElimGraph(Gelim))
437 { fprintf(stderr, "\nError in function buildElement\n"
438 " unable to construct element (not enough memory)\n");
439 quit();
440 }
441
442 /* crunch partially constructed element me */
443 medeststart2 = G->nedges;
444 for (p = medeststart; p < medestptr; p++)
445 adjncy[G->nedges++] = adjncy[p];
446 medeststart = medeststart2;
447 medestptr = G->nedges;
448
449 mesrcptr = xadj[me];
450 p = xadj[e];
451 }
452 adjncy[medestptr++] = v;
453 }
454 }
455
456 /* ----------------------
457 mark absorbed elements
458 ---------------------- */
459 if (e != me)
460 { xadj[e] = -1;
461 parent[e] = me;
462 score[e] = -4;
463 }
464 }
465
466 G->nedges = medestptr; /* new element Lme ends here */
467 }
468
469 /* -----------------------------------
470 element me successfully constructed
471 ----------------------------------- */
472 degree[me] = degme;
473 xadj[me] = medeststart;
474 vwght[me] = -vwght[me];
475 elen[me] = 0;
476 len[me] = medestptr - medeststart;
477 if (len[me] == 0)
478 xadj[me] = -1;
479
480 /* ---------------------------
481 unmark all variables in Lme
482 --------------------------- */
483 mesrcptr = xadj[me];
484 vlenme = len[me];
485 for (i = 0; i < vlenme; i++)
486 { v = adjncy[mesrcptr++];
487 vwght[v] = -vwght[v];
488 }
489 }
490
491
492 /*****************************************************************************
493 ******************************************************************************/
494 void
updateAdjncy(gelim_t * Gelim,PORD_INT * reachset,PORD_INT nreach,PORD_INT * tmp,PORD_INT * pflag)495 updateAdjncy(gelim_t *Gelim, PORD_INT *reachset, PORD_INT nreach, PORD_INT *tmp, PORD_INT *pflag)
496 { PORD_INT *xadj, *adjncy, *vwght, *len, *elen, *parent, *score;
497 PORD_INT u, v, e, me, i, j, jj, jdest, jfirstolde, jfirstv, jstart, jstop;
498 PORD_INT covered, marku;
499
500 xadj = Gelim->G->xadj;
501 adjncy = Gelim->G->adjncy;
502 vwght = Gelim->G->vwght;
503 len = Gelim->len;
504 elen = Gelim->elen;
505 parent = Gelim->parent;
506 score = Gelim->score;
507
508 /* -----------------------------------------------------------------
509 build the new element/variable list for each variable in reachset
510 ----------------------------------------------------------------- */
511 for (i = 0; i < nreach; i++)
512 { u = reachset[i];
513 vwght[u] = -vwght[u]; /* mark all variables in reachset */
514 jstart = xadj[u];
515 jstop = xadj[u] + len[u];
516 jdest = jfirstolde = jstart;
517
518 #ifdef DEBUG
519 printf("Updating adjacency list of node %d\n", u);
520 #endif
521
522 /* --------------------------------------------------------
523 scan the list of elements associated with variable u
524 place newly formed elements at the beginning of the list
525 -------------------------------------------------------- */
526 for (j = jstart; j < jstart + elen[u]; j++)
527 { e = adjncy[j];
528
529 #ifdef DEBUG
530 printf(" >> element %d (score %d, parent %d)\n", e,score[e],parent[e]);
531 #endif
532
533 if (score[e] == -4) /* e has been absorbed in this elim. step */
534 { me = parent[e]; /* me is the newly formed element */
535 if (tmp[me] < *pflag)
536 { adjncy[jdest++] = adjncy[jfirstolde]; /* move 1st old e to end */
537 adjncy[jfirstolde++] = me; /* append me at the beg. */
538 tmp[me] = *pflag;
539 }
540 }
541 else /* e has not been absorbed, i.e. it is */
542 if (tmp[e] < *pflag) /* an old element */
543 { adjncy[jdest++] = e;
544 tmp[e] = *pflag;
545 }
546 }
547 jfirstv = jdest; /* list of variables starts here */
548
549 /* -------------------------------------------------------
550 scan the list of variables associated with variable u
551 place newly formed elements at the begining of the list
552 ------------------------------------------------------- */
553 for (j = jstart + elen[u]; j < jstop; j++)
554 { v = adjncy[j];
555
556 #ifdef DEBUG
557 printf(" >> variable %d (score %d)\n", v, score[v]);
558 #endif
559
560 if (score[v] == -3) /* v has been eliminated in this step */
561 { if (tmp[v] < *pflag) /* and, thus, forms a newly created elem. */
562 { adjncy[jdest++] = adjncy[jfirstv]; /* move 1st var. to end */
563 adjncy[jfirstv++] = adjncy[jfirstolde]; /* move 1st old e to end */
564 adjncy[jfirstolde++] = v; /* append v at the beg. */
565 tmp[v] = *pflag;
566 }
567 }
568 else
569 adjncy[jdest++] = v; /* v is still a variable */
570 }
571 elen[u] = jfirstv - jstart;
572 len[u] = jdest - jstart;
573 (*pflag)++; /* clear tmp for next round */
574
575 #ifdef DEBUG
576 printf(" node %d: neighboring elements:\n", u);
577 for (j = jstart; j < jstart + elen[u]; j++)
578 printf("%5d", adjncy[j]);
579 printf("\n node %d: neighboring variables:\n", u);
580 for (j = jstart + elen[u]; j < jstart + len[u]; j++)
581 printf("%5d", adjncy[j]);
582 printf("\n");
583 #endif
584 }
585
586 /* ---------------------------------------------------------
587 remove from each list all covered edges between variables
588 --------------------------------------------------------- */
589 for (i = 0; i < nreach; i++)
590 { u = reachset[i];
591 jstart = xadj[u];
592 jstop = jstart + len[u];
593 marku = FALSE;
594
595 for (jdest = j = jstart + elen[u]; j < jstop; j++)
596 { v = adjncy[j];
597 if (vwght[v] > 0) /* v does not belong to reachset */
598 adjncy[jdest++] = v; /* edge (u,v) not covered */
599 if (vwght[v] < 0) /* both vertices belong to reachset */
600 { covered = FALSE; /* check for a common element */
601 if (!marku)
602 { for (jj = jstart; jj < jstart + elen[u]; jj++) /* mark elem. */
603 tmp[adjncy[jj]] = *pflag; /* of u */
604 marku = TRUE;
605 }
606 for (jj = xadj[v]; jj < xadj[v] + elen[v]; jj++) /* check elem. */
607 if (tmp[adjncy[jj]] == *pflag) /* of v */
608 { covered = TRUE;
609 break;
610 }
611 if (!covered)
612 adjncy[jdest++] = v;
613 }
614 }
615 len[u] = jdest - jstart;
616 (*pflag)++; /* clear tmp for next round */
617
618 #ifdef DEBUG
619 printf(" node %d: neighboring uncovered variables:\n", u);
620 for (j = jstart + elen[u]; j < jstart + len[u]; j++)
621 printf("%5d", adjncy[j]);
622 printf("\n");
623 #endif
624 }
625
626 /* --------------------------------
627 unmark all variables in reachset
628 -------------------------------- */
629 for (i = 0; i < nreach; i++)
630 { u = reachset[i];
631 vwght[u] = -vwght[u];
632 }
633 }
634
635
636 /*****************************************************************************
637 ******************************************************************************/
638 void
findIndNodes(gelim_t * Gelim,PORD_INT * reachset,PORD_INT nreach,PORD_INT * bin,PORD_INT * next,PORD_INT * tmp,PORD_INT * pflag)639 findIndNodes(gelim_t *Gelim, PORD_INT *reachset, PORD_INT nreach, PORD_INT *bin, PORD_INT *next,
640 PORD_INT *tmp, PORD_INT *pflag)
641 { PORD_INT *xadj, *adjncy, *vwght, *len, *elen, *parent, *score;
642 PORD_INT nvtx, chk, keepon, u, v, w, wlast, i, j, jstart, jstop, jstep, jj, jjstop;
643 nvtx = Gelim->G->nvtx;
644 xadj = Gelim->G->xadj;
645 adjncy = Gelim->G->adjncy;
646 vwght = Gelim->G->vwght;
647 len = Gelim->len;
648 elen = Gelim->elen;
649 parent = Gelim->parent;
650 score = Gelim->score;
651
652 #ifdef DEBUG
653 printf("Checking reachset for indistinguishable variables\n");
654 #endif
655
656 /* -----------------------------------------------------------------------
657 compute checksums for all principal variables on reachset and fill bins
658 NOTE: checksums are stored in parent vector
659 ----------------------------------------------------------------------- */
660 for (i = 0; i < nreach; i++)
661 { u = reachset[i];
662 chk = 0;
663 jstart = xadj[u];
664 jstop = jstart + len[u];
665 /* Modified by JYL: 16 march 2005:
666 * This code was failing in case of
667 * overflow.
668 for (j = jstart; j < jstop; j++)
669 chk += adjncy[j];
670 chk = chk % nvtx;
671 */
672 jstep=max(1000000000/nvtx,1);
673 for (j = jstart; j < jstop; j+=jstep)
674 {
675 jjstop = min(jstop, j+jstep);
676 for (jj = j; jj < jjstop; jj++)
677 chk += adjncy[jj];
678 chk = chk % nvtx;
679 }
680
681 parent[u] = chk;
682 /* JYL: temporary:
683 if (parent[u] < - 10)
684 printf("Probleme %d \n",chk);*/
685 next[u] = bin[chk];
686 bin[chk] = u;
687 }
688
689 /* -----------------------
690 supervariable detection
691 ----------------------- */
692 for (i = 0; i < nreach; i++)
693 { u = reachset[i];
694 if (vwght[u] > 0) /* u is a principal variable */
695 { chk = parent[u]; /* search bin[chk] for ind. nodes */
696 v = bin[chk]; /* okay, v is the first node in this bin */
697 bin[chk] = -1; /* no further examinations of this bin */
698 while (v != -1)
699 { jstart = xadj[v];
700 jstop = xadj[v] + len[v];
701 for (j = jstart; j < jstop; j++)
702 tmp[adjncy[j]] = *pflag;
703 w = next[v]; /* v is principal and w is a potential */
704 wlast = v; /* nonprincipal variable */
705 while (w != -1)
706 { keepon = TRUE;
707 if ((len[w] != len[v]) || (elen[w] != elen[v])
708 || ((score[w] < 0) && (score[v] >= 0))
709 || ((score[w] >= 0) && (score[v] < 0)))
710 keepon = FALSE;
711 if (keepon)
712 { for (jj = xadj[w]; jj < xadj[w] + len[w]; jj++)
713 if (tmp[adjncy[jj]] < *pflag)
714 { keepon = FALSE;
715 break;
716 }
717 }
718 if (keepon) /* found it! mark w as nonprincipal */
719 { parent[w] = v; /* representative of w is v */
720 /* Temporary JY
721 if (parent[w] < - 10)
722 printf("Probleme\n");
723 */
724 #ifdef DEBUG
725 printf(" non-principal variable %d (score %d) mapped onto "
726 "%d (score %d)\n", w, score[w], v, score[v]);
727 #endif
728
729 vwght[v] += vwght[w]; /* add weight of w */
730 vwght[w] = 0;
731 xadj[w] = -1; /* w's adjacency list can be over- */
732 score[w] = -2; /* written during next crunch */
733 w = next[w];
734 next[wlast] = w; /* remove w from bin */
735 }
736 else /* failed */
737 { wlast = w;
738 w = next[w];
739 }
740 }
741 v = next[v]; /* no more variables can be absorbed by v */
742 (*pflag)++; /* clear tmp vector for next round */
743 }
744 }
745 }
746
747 /* -------------------------------------------------------
748 re-initialize parent vector for all principal variables
749 ------------------------------------------------------- */
750 for (i = 0; i < nreach; i++)
751 { u = reachset[i];
752 if (vwght[u] > 0)
753 parent[u] = -1;
754 }
755 }
756
757
758 /*****************************************************************************
759 ******************************************************************************/
760 void
updateDegree(gelim_t * Gelim,PORD_INT * reachset,PORD_INT nreach,PORD_INT * bin)761 updateDegree(gelim_t *Gelim, PORD_INT *reachset, PORD_INT nreach, PORD_INT *bin)
762 { PORD_INT *xadj, *adjncy, *vwght, *len, *elen, *degree;
763 PORD_INT totvwght, deg, vwghtv, u, v, w, e, me, r, i, istart, istop;
764 PORD_INT j, jstart, jstop;
765
766 totvwght = Gelim->G->totvwght;
767 xadj = Gelim->G->xadj;
768 adjncy = Gelim->G->adjncy;
769 vwght = Gelim->G->vwght;
770 len = Gelim->len;
771 elen = Gelim->elen;
772 degree = Gelim->degree;
773
774 /* -------------------------------------------------------------------
775 degree update only for those vertices in reachset that are adjacent
776 to an element
777 ------------------------------------------------------------------- */
778 for (r = 0; r < nreach; r++)
779 { u = reachset[r];
780 if (elen[u] > 0)
781 bin[u] = 1;
782 }
783
784 /* -----------------------------------------
785 and now do the approximate degree updates
786 ----------------------------------------- */
787 for (r = 0; r < nreach; r++)
788 { u = reachset[r];
789 if (bin[u] == 1) /* me is the most recently formed element */
790 { me = adjncy[xadj[u]]; /* in the neighborhood of u */
791
792 #ifdef DEBUG
793 printf("Updating degree of all variables in L(%d) (initiated by %d)\n",
794 me, u);
795 #endif
796
797 /* ----------------------------------------------------------------
798 compute in bin[e] the size of Le\Lme for all unabsorbed elements
799 ---------------------------------------------------------------- */
800 istart = xadj[me];
801 istop = istart + len[me]; /* compute in bin[e] the size */
802 for (i = istart; i < istop; i++) /* of Le/Lme for all elements */
803 { v = adjncy[i]; /* e != me that are adjacent */
804 vwghtv = vwght[v]; /* to a principal var. e Lme */
805 if (vwghtv > 0)
806 { jstart = xadj[v];
807 jstop = jstart + elen[v];
808 for (j = jstart; j < jstop; j++)
809 { e = adjncy[j];
810 if (e != me)
811 { if (bin[e] > 0) bin[e] -= vwghtv;
812 else bin[e] = degree[e] - vwghtv;
813 }
814 }
815 }
816 }
817
818 #ifdef DEBUG
819 for (i = istart; i < istop; i++)
820 { v = adjncy[i];
821 if (vwght[v] > 0)
822 for (j = xadj[v]; j < xadj[v] + elen[v]; j++)
823 { e = adjncy[j];
824 if (e != me)
825 printf(" >> element %d: degree %d, outer degree %d\n", e,
826 degree[e], bin[e]);
827 }
828 }
829 #endif
830
831 /* ------------------------------------------------------
832 update approx. degree for all v in Lme with bin[v] = 1
833 ------------------------------------------------------ */
834 for (i = istart; i < istop; i++)
835 { v = adjncy[i]; /* update the upper bound deg. */
836 vwghtv = vwght[v]; /* of all principal variables */
837 deg = 0; /* in Lme that have not been */
838 if (bin[v] == 1) /* updated yet */
839 { jstart = xadj[v];
840 jstop = jstart + len[v];
841
842 /* scan the element list associated with principal v */
843 for (j = jstart; j < jstart + elen[v]; j++)
844 { e = adjncy[j];
845 if (e != me) deg += bin[e];
846 }
847
848 /* scan the supervariables in the list associated with v */
849 for (j = jstart + elen[v]; j < jstop; j++)
850 { w = adjncy[j];
851 deg += vwght[w];
852 }
853
854 /* compute the external degree of v (add size of Lme) */
855 deg = min(degree[v], deg);
856 degree[v] = max(1, min(deg+degree[me]-vwghtv, totvwght-vwghtv));
857 bin[v] = -1;
858
859 #ifdef DEBUG
860 printf(" >> variable %d (totvwght %d, vwght %d): deg %d, "
861 "degme %d, approx degree %d\n", v, totvwght, vwghtv, deg,
862 degree[me], degree[v]);
863 #endif
864 }
865 }
866
867 /* ------------------------------------
868 clear bin[e] of all elements e != me
869 ------------------------------------ */
870 for (i = istart; i < istop; i++)
871 { v = adjncy[i];
872 vwghtv = vwght[v];
873 if (vwghtv > 0)
874 { jstart = xadj[v];
875 jstop = jstart + elen[v];
876 for (j = jstart; j < jstop; j++)
877 { e = adjncy[j];
878 if (e != me) bin[e] = -1;
879 }
880 }
881 }
882 }
883 }
884 }
885
886
887 /*****************************************************************************
888 ******************************************************************************/
889 void
updateScore(gelim_t * Gelim,PORD_INT * reachset,PORD_INT nreach,PORD_INT scoretype,PORD_INT * bin)890 updateScore(gelim_t *Gelim, PORD_INT *reachset, PORD_INT nreach, PORD_INT scoretype, PORD_INT *bin)
891 { PORD_INT *xadj, *adjncy, *vwght, *len, *elen, *degree, *score;
892 PORD_INT vwghtv, deg, degme, u, v, me, r, i, istart, istop;
893 /* Modified by JYL, 16 march 2005.
894 * scr could overflow for quasi dense rows.
895 * Use a double instead for large degrees
896 * aset it near to MAX_INT in case of problem.
897 */
898 double scr_dbl;
899 PORD_INT scr;
900
901 xadj = Gelim->G->xadj;
902 adjncy = Gelim->G->adjncy;
903 vwght = Gelim->G->vwght;
904 len = Gelim->len;
905 elen = Gelim->elen;
906 degree = Gelim->degree;
907 score = Gelim->score;
908
909 /* ------------------------------------------------------------------
910 score update only for those vertices in reachset that are adjacent
911 to an element
912 ------------------------------------------------------------------ */
913 for (r = 0; r < nreach; r++)
914 { u = reachset[r];
915 if (elen[u] > 0)
916 bin[u] = 1;
917 }
918
919 /* ----------------------------
920 and now do the score updates
921 ---------------------------- */
922 scoretype = scoretype % 10;
923 for (r = 0; r < nreach; r++)
924 { u = reachset[r];
925 if (bin[u] == 1) /* me is the most recently formed element */
926 { me = adjncy[xadj[u]]; /* in the neighborhood of u */
927
928 #ifdef DEBUG
929 printf("Updating score of all variables in L(%d) (initiated by %d)\n",
930 me, u);
931 #endif
932
933 istart = xadj[me];
934 istop = xadj[me] + len[me];
935 for (i = istart; i < istop; i++)
936 { v = adjncy[i]; /* update score of all principal */
937 if (bin[v] == 1) /* variables in Lme that have not */
938 { vwghtv = vwght[v]; /* been updated yet */
939 deg = degree[v];
940 degme = degree[me] - vwghtv;
941 if (deg > 40000 || degme > 40000)
942 {
943 switch(scoretype)
944 { case AMD:
945 scr_dbl = (double)deg;
946 break;
947 case AMF:
948 scr_dbl = (double)deg*(double)(deg-1)/2 - (double)degme*(double)(degme-1)/2;
949 break;
950 case AMMF:
951 scr_dbl = ((double)deg*(double)(deg-1)/2 - (double)degme*(double)(degme-1)/2) / (double)vwghtv;
952 break;
953 case AMIND:
954 scr_dbl = max(0, ((double)deg*(double)(deg-1)/2 - (double)degme*(double)(degme-1)/2)
955 - (double)deg*(double)vwghtv);
956 break;
957 default:
958 fprintf(stderr, "\nError in function updateScore\n"
959 " unrecognized selection strategy %d\n", scoretype);
960 quit();
961 }
962 /* Some buckets have offset nvtx / 2.
963 * Using MAX_INT - nvtx should then be safe */
964 score[v] = (PORD_INT) (min(scr_dbl,MAX_INT-Gelim->G->nvtx));
965 }
966 else
967 {
968 switch(scoretype)
969 { case AMD:
970 scr = deg;
971 break;
972 case AMF:
973 scr = deg*(deg-1)/2 - degme*(degme-1)/2;
974 break;
975 case AMMF:
976 scr = (deg*(deg-1)/2 - degme*(degme-1)/2) / vwghtv;
977 break;
978 case AMIND:
979 scr = max(0, (deg*(deg-1)/2 - degme*(degme-1)/2)
980 - deg*vwghtv);
981 break;
982 default:
983 fprintf(stderr, "\nError in function updateScore\n"
984 " unrecognized selection strategy %d\n", scoretype);
985 quit();
986 }
987 score[v] = scr;
988 }
989 bin[v] = -1;
990
991 #ifdef DEBUG
992 printf(" >> variable %d (me %d): weight %d, (ext)degme %d, "
993 "degree %d, score %d\n", u, me, vwghtv, degme, degree[v],
994 score[v]);
995 #endif
996
997 if (score[v] < 0)
998 { fprintf(stderr, "\nError in function updateScore\n"
999 " score[%d] = %d is negative\n", v, score[v]);
1000 quit();
1001 }
1002 }
1003 }
1004 }
1005 }
1006 }
1007
1008
1009 /*****************************************************************************)
1010 ******************************************************************************/
1011 elimtree_t*
extractElimTree(gelim_t * Gelim)1012 extractElimTree(gelim_t *Gelim)
1013 { elimtree_t *T;
1014 PORD_INT *vwght, *par, *degree, *score, *sib, *fch;
1015 PORD_INT *ncolfactor, *ncolupdate, *parent, *vtx2front;
1016 PORD_INT nvtx, nfronts, root, u, v, front;
1017
1018 nvtx = Gelim->G->nvtx;
1019 vwght = Gelim->G->vwght;
1020 par = Gelim->parent;
1021 degree = Gelim->degree;
1022 score = Gelim->score;
1023
1024 /* ------------------------
1025 allocate working storage
1026 ------------------------ */
1027 mymalloc(sib, nvtx, PORD_INT);
1028 mymalloc(fch, nvtx, PORD_INT);
1029 for (u = 0; u < nvtx; u++)
1030 sib[u] = fch[u] = -1;
1031
1032 /* --------------------------------------------------------------
1033 count fronts and create top-down view of the tree given by par
1034 -------------------------------------------------------------- */
1035 nfronts = 0;
1036 root = -1;
1037 for (u = 0; u < nvtx; u++)
1038 switch(score[u])
1039 { case -2: /* variable u is nonprincipal */
1040 break;
1041 case -3: /* variable u has been elim. and now forms an elem. */
1042 sib[u] = root;
1043 root = u;
1044 nfronts++;
1045 break;
1046 case -4: /* element u has been absorbed by par[u] */
1047 v = par[u];
1048 sib[u] = fch[v];
1049 fch[v] = u;
1050 nfronts++;
1051 break;
1052 default:
1053 fprintf(stderr, "\nError in function extractElimTree\n"
1054 " ordering not complete (score[%d] = %d)\n", u, score[u]);
1055 quit();
1056 }
1057
1058 #ifdef DEBUG
1059 for (u = 0; u < nvtx; u++)
1060 printf("node %d: score %d, par %d, fch %d, sib %d\n", u, score[u],
1061 par[u], fch[u], sib[u]);
1062 #endif
1063
1064 /* --------------------------------------
1065 allocate space for the elimtree object
1066 -------------------------------------- */
1067 T = newElimTree(nvtx, nfronts);
1068 ncolfactor = T->ncolfactor;
1069 ncolupdate = T->ncolupdate;
1070 parent = T->parent;
1071 vtx2front = T->vtx2front;
1072
1073 /* -------------------------------------------------------------
1074 fill the vtx2front vector so that representative vertices are
1075 mapped in a post-order traversal
1076 ------------------------------------------------------------- */
1077 nfronts = 0;
1078 u = root;
1079 while (u != -1)
1080 { while (fch[u] != -1)
1081 u = fch[u];
1082 vtx2front[u] = nfronts++;
1083 while ((sib[u] == -1) && (par[u] != -1))
1084 { u = par[u];
1085 vtx2front[u] = nfronts++;
1086 }
1087 u = sib[u];
1088 }
1089
1090 /* ---------------------------------------------------
1091 fill in the vtx2front map for nonprincipal vertices
1092 --------------------------------------------------- */
1093 for (u = 0; u < nvtx; u++)
1094 if (score[u] == -2)
1095 { v = u;
1096 while ((par[v] != -1) && (score[v] == -2))
1097 v = par[v];
1098 vtx2front[u] = vtx2front[v];
1099 }
1100
1101 /* -------------------------------------------------------------
1102 set up the parent vector of T and fill ncolfactor, ncolupdate
1103 ------------------------------------------------------------- */
1104 for (u = 0; u < nvtx; u++)
1105 { front = vtx2front[u];
1106 if (score[u] == -3)
1107 { parent[front] = -1;
1108 ncolfactor[front] = vwght[u];
1109 ncolupdate[front] = degree[u];
1110 }
1111 if (score[u] == -4)
1112 { parent[front] = vtx2front[par[u]];
1113 ncolfactor[front] = vwght[u];
1114 ncolupdate[front] = degree[u];
1115 }
1116 }
1117
1118 /* ----------------------------
1119 set up all other arrays of T
1120 ---------------------------- */
1121 initFchSilbRoot(T);
1122
1123 /* ----------------------
1124 free memory and return
1125 ---------------------- */
1126 free(sib); free(fch);
1127 return(T);
1128 }
1129
1130