lemon-project-template-glpk

comparison deps/glpk/src/glpqmd.c @ 11:4fc6ad2fb8a6

find changesets by author, revision, files, or words in the commit message
Test GLPK in src/main.cc
author Alpar Juttner <alpar@cs.elte.hu>
date Sun, 06 Nov 2011 21:43:29 +0100
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:bf3ec92398a1
1 /* glpqmd.c (quotient minimum degree algorithm) */
2
3 /***********************************************************************
4 * This code is part of GLPK (GNU Linear Programming Kit).
5 *
6 * THIS CODE IS THE RESULT OF TRANSLATION OF THE FORTRAN SUBROUTINES
7 * GENQMD, QMDRCH, QMDQT, QMDUPD, AND QMDMRG FROM THE BOOK:
8 *
9 * ALAN GEORGE, JOSEPH W-H LIU. COMPUTER SOLUTION OF LARGE SPARSE
10 * POSITIVE DEFINITE SYSTEMS. PRENTICE-HALL, 1981.
11 *
12 * THE TRANSLATION HAS BEEN DONE WITH THE PERMISSION OF THE AUTHORS
13 * OF THE ORIGINAL FORTRAN SUBROUTINES: ALAN GEORGE AND JOSEPH LIU,
14 * UNIVERSITY OF WATERLOO, WATERLOO, ONTARIO, CANADA.
15 *
16 * The translation was made by Andrew Makhorin <mao@gnu.org>.
17 *
18 * GLPK is free software: you can redistribute it and/or modify it
19 * under the terms of the GNU General Public License as published by
20 * the Free Software Foundation, either version 3 of the License, or
21 * (at your option) any later version.
22 *
23 * GLPK is distributed in the hope that it will be useful, but WITHOUT
24 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
25 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
26 * License for more details.
27 *
28 * You should have received a copy of the GNU General Public License
29 * along with GLPK. If not, see <http://www.gnu.org/licenses/>.
30 ***********************************************************************/
31
32 #include "glpqmd.h"
33
34 /***********************************************************************
35 * NAME
36 *
37 * genqmd - GENeral Quotient Minimum Degree algorithm
38 *
39 * SYNOPSIS
40 *
41 * #include "glpqmd.h"
42 * void genqmd(int *neqns, int xadj[], int adjncy[], int perm[],
43 * int invp[], int deg[], int marker[], int rchset[], int nbrhd[],
44 * int qsize[], int qlink[], int *nofsub);
45 *
46 * PURPOSE
47 *
48 * This routine implements the minimum degree algorithm. It makes use
49 * of the implicit representation of the elimination graph by quotient
50 * graphs, and the notion of indistinguishable nodes.
51 *
52 * CAUTION
53 *
54 * The adjancy vector adjncy will be destroyed.
55 *
56 * INPUT PARAMETERS
57 *
58 * neqns - number of equations;
59 * (xadj, adjncy) -
60 * the adjancy structure.
61 *
62 * OUTPUT PARAMETERS
63 *
64 * perm - the minimum degree ordering;
65 * invp - the inverse of perm.
66 *
67 * WORKING PARAMETERS
68 *
69 * deg - the degree vector. deg[i] is negative means node i has been
70 * numbered;
71 * marker - a marker vector, where marker[i] is negative means node i
72 * has been merged with another nodeand thus can be ignored;
73 * rchset - vector used for the reachable set;
74 * nbrhd - vector used for neighborhood set;
75 * qsize - vector used to store the size of indistinguishable
76 * supernodes;
77 * qlink - vector used to store indistinguishable nodes, i, qlink[i],
78 * qlink[qlink[i]], ... are the members of the supernode
79 * represented by i.
80 *
81 * PROGRAM SUBROUTINES
82 *
83 * qmdrch, qmdqt, qmdupd.
84 ***********************************************************************/
85
86 void genqmd(int *_neqns, int xadj[], int adjncy[], int perm[],
87 int invp[], int deg[], int marker[], int rchset[], int nbrhd[],
88 int qsize[], int qlink[], int *_nofsub)
89 { int inode, ip, irch, j, mindeg, ndeg, nhdsze, node, np, num,
90 nump1, nxnode, rchsze, search, thresh;
91 # define neqns (*_neqns)
92 # define nofsub (*_nofsub)
93 /* Initialize degree vector and other working variables. */
94 mindeg = neqns;
95 nofsub = 0;
96 for (node = 1; node <= neqns; node++)
97 { perm[node] = node;
98 invp[node] = node;
99 marker[node] = 0;
100 qsize[node] = 1;
101 qlink[node] = 0;
102 ndeg = xadj[node+1] - xadj[node];
103 deg[node] = ndeg;
104 if (ndeg < mindeg) mindeg = ndeg;
105 }
106 num = 0;
107 /* Perform threshold search to get a node of min degree.
108 Variable search point to where search should start. */
109 s200: search = 1;
110 thresh = mindeg;
111 mindeg = neqns;
112 s300: nump1 = num + 1;
113 if (nump1 > search) search = nump1;
114 for (j = search; j <= neqns; j++)
115 { node = perm[j];
116 if (marker[node] >= 0)
117 { ndeg = deg[node];
118 if (ndeg <= thresh) goto s500;
119 if (ndeg < mindeg) mindeg = ndeg;
120 }
121 }
122 goto s200;
123 /* Node has minimum degree. Find its reachable sets by calling
124 qmdrch. */
125 s500: search = j;
126 nofsub += deg[node];
127 marker[node] = 1;
128 qmdrch(&node, xadj, adjncy, deg, marker, &rchsze, rchset, &nhdsze,
129 nbrhd);
130 /* Eliminate all nodes indistinguishable from node. They are given
131 by node, qlink[node], ... . */
132 nxnode = node;
133 s600: num++;
134 np = invp[nxnode];
135 ip = perm[num];
136 perm[np] = ip;
137 invp[ip] = np;
138 perm[num] = nxnode;
139 invp[nxnode] = num;
140 deg[nxnode] = -1;
141 nxnode = qlink[nxnode];
142 if (nxnode > 0) goto s600;
143 if (rchsze > 0)
144 { /* Update the degrees of the nodes in the reachable set and
145 identify indistinguishable nodes. */
146 qmdupd(xadj, adjncy, &rchsze, rchset, deg, qsize, qlink,
147 marker, &rchset[rchsze+1], &nbrhd[nhdsze+1]);
148 /* Reset marker value of nodes in reach set. Update threshold
149 value for cyclic search. Also call qmdqt to form new
150 quotient graph. */
151 marker[node] = 0;
152 for (irch = 1; irch <= rchsze; irch++)
153 { inode = rchset[irch];
154 if (marker[inode] >= 0)
155 { marker[inode] = 0;
156 ndeg = deg[inode];
157 if (ndeg < mindeg) mindeg = ndeg;
158 if (ndeg <= thresh)
159 { mindeg = thresh;
160 thresh = ndeg;
161 search = invp[inode];
162 }
163 }
164 }
165 if (nhdsze > 0)
166 qmdqt(&node, xadj, adjncy, marker, &rchsze, rchset, nbrhd);
167 }
168 if (num < neqns) goto s300;
169 return;
170 # undef neqns
171 # undef nofsub
172 }
173
174 /***********************************************************************
175 * NAME
176 *
177 * qmdrch - Quotient MD ReaCHable set
178 *
179 * SYNOPSIS
180 *
181 * #include "glpqmd.h"
182 * void qmdrch(int *root, int xadj[], int adjncy[], int deg[],
183 * int marker[], int *rchsze, int rchset[], int *nhdsze,
184 * int nbrhd[]);
185 *
186 * PURPOSE
187 *
188 * This subroutine determines the reachable set of a node through a
189 * given subset. The adjancy structure is assumed to be stored in a
190 * quotient graph format.
191 *
192 * INPUT PARAMETERS
193 *
194 * root - the given node not in the subset;
195 * (xadj, adjncy) -
196 * the adjancy structure pair;
197 * deg - the degree vector. deg[i] < 0 means the node belongs to the
198 * given subset.
199 *
200 * OUTPUT PARAMETERS
201 *
202 * (rchsze, rchset) -
203 * the reachable set;
204 * (nhdsze, nbrhd) -
205 * the neighborhood set.
206 *
207 * UPDATED PARAMETERS
208 *
209 * marker - the marker vector for reach and nbrhd sets. > 0 means the
210 * node is in reach set. < 0 means the node has been merged
211 * with others in the quotient or it is in nbrhd set.
212 ***********************************************************************/
213
214 void qmdrch(int *_root, int xadj[], int adjncy[], int deg[],
215 int marker[], int *_rchsze, int rchset[], int *_nhdsze,
216 int nbrhd[])
217 { int i, istop, istrt, j, jstop, jstrt, nabor, node;
218 # define root (*_root)
219 # define rchsze (*_rchsze)
220 # define nhdsze (*_nhdsze)
221 /* Loop through the neighbors of root in the quotient graph. */
222 nhdsze = 0;
223 rchsze = 0;
224 istrt = xadj[root];
225 istop = xadj[root+1] - 1;
226 if (istop < istrt) return;
227 for (i = istrt; i <= istop; i++)
228 { nabor = adjncy[i];
229 if (nabor == 0) return;
230 if (marker[nabor] == 0)
231 { if (deg[nabor] >= 0)
232 { /* Include nabor into the reachable set. */
233 rchsze++;
234 rchset[rchsze] = nabor;
235 marker[nabor] = 1;
236 goto s600;
237 }
238 /* nabor has been eliminated. Find nodes reachable from
239 it. */
240 marker[nabor] = -1;
241 nhdsze++;
242 nbrhd[nhdsze] = nabor;
243 s300: jstrt = xadj[nabor];
244 jstop = xadj[nabor+1] - 1;
245 for (j = jstrt; j <= jstop; j++)
246 { node = adjncy[j];
247 nabor = - node;
248 if (node < 0) goto s300;
249 if (node == 0) goto s600;
250 if (marker[node] == 0)
251 { rchsze++;
252 rchset[rchsze] = node;
253 marker[node] = 1;
254 }
255 }
256 }
257 s600: ;
258 }
259 return;
260 # undef root
261 # undef rchsze
262 # undef nhdsze
263 }
264
265 /***********************************************************************
266 * NAME
267 *
268 * qmdqt - Quotient MD Quotient graph Transformation
269 *
270 * SYNOPSIS
271 *
272 * #include "glpqmd.h"
273 * void qmdqt(int *root, int xadj[], int adjncy[], int marker[],
274 * int *rchsze, int rchset[], int nbrhd[]);
275 *
276 * PURPOSE
277 *
278 * This subroutine performs the quotient graph transformation after a
279 * node has been eliminated.
280 *
281 * INPUT PARAMETERS
282 *
283 * root - the node just eliminated. It becomes the representative of
284 * the new supernode;
285 * (xadj, adjncy) -
286 * the adjancy structure;
287 * (rchsze, rchset) -
288 * the reachable set of root in the old quotient graph;
289 * nbrhd - the neighborhood set which will be merged with root to form
290 * the new supernode;
291 * marker - the marker vector.
292 *
293 * UPDATED PARAMETERS
294 *
295 * adjncy - becomes the adjncy of the quotient graph.
296 ***********************************************************************/
297
298 void qmdqt(int *_root, int xadj[], int adjncy[], int marker[],
299 int *_rchsze, int rchset[], int nbrhd[])
300 { int inhd, irch, j, jstop, jstrt, link, nabor, node;
301 # define root (*_root)
302 # define rchsze (*_rchsze)
303 irch = 0;
304 inhd = 0;
305 node = root;
306 s100: jstrt = xadj[node];
307 jstop = xadj[node+1] - 2;
308 if (jstop >= jstrt)
309 { /* Place reach nodes into the adjacent list of node. */
310 for (j = jstrt; j <= jstop; j++)
311 { irch++;
312 adjncy[j] = rchset[irch];
313 if (irch >= rchsze) goto s400;
314 }
315 }
316 /* Link to other space provided by the nbrhd set. */
317 link = adjncy[jstop+1];
318 node = - link;
319 if (link >= 0)
320 { inhd++;
321 node = nbrhd[inhd];
322 adjncy[jstop+1] = - node;
323 }
324 goto s100;
325 /* All reachable nodes have been saved. End the adjacent list.
326 Add root to the neighborhood list of each node in the reach
327 set. */
328 s400: adjncy[j+1] = 0;
329 for (irch = 1; irch <= rchsze; irch++)
330 { node = rchset[irch];
331 if (marker[node] >= 0)
332 { jstrt = xadj[node];
333 jstop = xadj[node+1] - 1;
334 for (j = jstrt; j <= jstop; j++)
335 { nabor = adjncy[j];
336 if (marker[nabor] < 0)
337 { adjncy[j] = root;
338 goto s600;
339 }
340 }
341 }
342 s600: ;
343 }
344 return;
345 # undef root
346 # undef rchsze
347 }
348
349 /***********************************************************************
350 * NAME
351 *
352 * qmdupd - Quotient MD UPDate
353 *
354 * SYNOPSIS
355 *
356 * #include "glpqmd.h"
357 * void qmdupd(int xadj[], int adjncy[], int *nlist, int list[],
358 * int deg[], int qsize[], int qlink[], int marker[], int rchset[],
359 * int nbrhd[]);
360 *
361 * PURPOSE
362 *
363 * This routine performs degree update for a set of nodes in the minimum
364 * degree algorithm.
365 *
366 * INPUT PARAMETERS
367 *
368 * (xadj, adjncy) -
369 * the adjancy structure;
370 * (nlist, list) -
371 * the list of nodes whose degree has to be updated.
372 *
373 * UPDATED PARAMETERS
374 *
375 * deg - the degree vector;
376 * qsize - size of indistinguishable supernodes;
377 * qlink - linked list for indistinguishable nodes;
378 * marker - used to mark those nodes in reach/nbrhd sets.
379 *
380 * WORKING PARAMETERS
381 *
382 * rchset - the reachable set;
383 * nbrhd - the neighborhood set.
384 *
385 * PROGRAM SUBROUTINES
386 *
387 * qmdmrg.
388 ***********************************************************************/
389
390 void qmdupd(int xadj[], int adjncy[], int *_nlist, int list[],
391 int deg[], int qsize[], int qlink[], int marker[], int rchset[],
392 int nbrhd[])
393 { int deg0, deg1, il, inhd, inode, irch, j, jstop, jstrt, mark,
394 nabor, nhdsze, node, rchsze;
395 # define nlist (*_nlist)
396 /* Find all eliminated supernodes that are adjacent to some nodes
397 in the given list. Put them into (nhdsze, nbrhd). deg0 contains
398 the number of nodes in the list. */
399 if (nlist <= 0) return;
400 deg0 = 0;
401 nhdsze = 0;
402 for (il = 1; il <= nlist; il++)
403 { node = list[il];
404 deg0 += qsize[node];
405 jstrt = xadj[node];
406 jstop = xadj[node+1] - 1;
407 for (j = jstrt; j <= jstop; j++)
408 { nabor = adjncy[j];
409 if (marker[nabor] == 0 && deg[nabor] < 0)
410 { marker[nabor] = -1;
411 nhdsze++;
412 nbrhd[nhdsze] = nabor;
413 }
414 }
415 }
416 /* Merge indistinguishable nodes in the list by calling the
417 subroutine qmdmrg. */
418 if (nhdsze > 0)
419 qmdmrg(xadj, adjncy, deg, qsize, qlink, marker, &deg0, &nhdsze,
420 nbrhd, rchset, &nbrhd[nhdsze+1]);
421 /* Find the new degrees of the nodes that have not been merged. */
422 for (il = 1; il <= nlist; il++)
423 { node = list[il];
424 mark = marker[node];
425 if (mark == 0 || mark == 1)
426 { marker[node] = 2;
427 qmdrch(&node, xadj, adjncy, deg, marker, &rchsze, rchset,
428 &nhdsze, nbrhd);
429 deg1 = deg0;
430 if (rchsze > 0)
431 { for (irch = 1; irch <= rchsze; irch++)
432 { inode = rchset[irch];
433 deg1 += qsize[inode];
434 marker[inode] = 0;
435 }
436 }
437 deg[node] = deg1 - 1;
438 if (nhdsze > 0)
439 { for (inhd = 1; inhd <= nhdsze; inhd++)
440 { inode = nbrhd[inhd];
441 marker[inode] = 0;
442 }
443 }
444 }
445 }
446 return;
447 # undef nlist
448 }
449
450 /***********************************************************************
451 * NAME
452 *
453 * qmdmrg - Quotient MD MeRGe
454 *
455 * SYNOPSIS
456 *
457 * #include "qmdmrg.h"
458 * void qmdmrg(int xadj[], int adjncy[], int deg[], int qsize[],
459 * int qlink[], int marker[], int *deg0, int *nhdsze, int nbrhd[],
460 * int rchset[], int ovrlp[]);
461 *
462 * PURPOSE
463 *
464 * This routine merges indistinguishable nodes in the minimum degree
465 * ordering algorithm. It also computes the new degrees of these new
466 * supernodes.
467 *
468 * INPUT PARAMETERS
469 *
470 * (xadj, adjncy) -
471 * the adjancy structure;
472 * deg0 - the number of nodes in the given set;
473 * (nhdsze, nbrhd) -
474 * the set of eliminated supernodes adjacent to some nodes in
475 * the set.
476 *
477 * UPDATED PARAMETERS
478 *
479 * deg - the degree vector;
480 * qsize - size of indistinguishable nodes;
481 * qlink - linked list for indistinguishable nodes;
482 * marker - the given set is given by those nodes with marker value set
483 * to 1. Those nodes with degree updated will have marker value
484 * set to 2.
485 *
486 * WORKING PARAMETERS
487 *
488 * rchset - the reachable set;
489 * ovrlp - temp vector to store the intersection of two reachable sets.
490 ***********************************************************************/
491
492 void qmdmrg(int xadj[], int adjncy[], int deg[], int qsize[],
493 int qlink[], int marker[], int *_deg0, int *_nhdsze, int nbrhd[],
494 int rchset[], int ovrlp[])
495 { int deg1, head, inhd, iov, irch, j, jstop, jstrt, link, lnode,
496 mark, mrgsze, nabor, node, novrlp, rchsze, root;
497 # define deg0 (*_deg0)
498 # define nhdsze (*_nhdsze)
499 /* Initialization. */
500 if (nhdsze <= 0) return;
501 for (inhd = 1; inhd <= nhdsze; inhd++)
502 { root = nbrhd[inhd];
503 marker[root] = 0;
504 }
505 /* Loop through each eliminated supernode in the set
506 (nhdsze, nbrhd). */
507 for (inhd = 1; inhd <= nhdsze; inhd++)
508 { root = nbrhd[inhd];
509 marker[root] = -1;
510 rchsze = 0;
511 novrlp = 0;
512 deg1 = 0;
513 s200: jstrt = xadj[root];
514 jstop = xadj[root+1] - 1;
515 /* Determine the reachable set and its intersection with the
516 input reachable set. */
517 for (j = jstrt; j <= jstop; j++)
518 { nabor = adjncy[j];
519 root = - nabor;
520 if (nabor < 0) goto s200;
521 if (nabor == 0) break;
522 mark = marker[nabor];
523 if (mark == 0)
524 { rchsze++;
525 rchset[rchsze] = nabor;
526 deg1 += qsize[nabor];
527 marker[nabor] = 1;
528 }
529 else if (mark == 1)
530 { novrlp++;
531 ovrlp[novrlp] = nabor;
532 marker[nabor] = 2;
533 }
534 }
535 /* From the overlapped set, determine the nodes that can be
536 merged together. */
537 head = 0;
538 mrgsze = 0;
539 for (iov = 1; iov <= novrlp; iov++)
540 { node = ovrlp[iov];
541 jstrt = xadj[node];
542 jstop = xadj[node+1] - 1;
543 for (j = jstrt; j <= jstop; j++)
544 { nabor = adjncy[j];
545 if (marker[nabor] == 0)
546 { marker[node] = 1;
547 goto s1100;
548 }
549 }
550 /* Node belongs to the new merged supernode. Update the
551 vectors qlink and qsize. */
552 mrgsze += qsize[node];
553 marker[node] = -1;
554 lnode = node;
555 s900: link = qlink[lnode];
556 if (link > 0)
557 { lnode = link;
558 goto s900;
559 }
560 qlink[lnode] = head;
561 head = node;
562 s1100: ;
563 }
564 if (head > 0)
565 { qsize[head] = mrgsze;
566 deg[head] = deg0 + deg1 - 1;
567 marker[head] = 2;
568 }
569 /* Reset marker values. */
570 root = nbrhd[inhd];
571 marker[root] = 0;
572 if (rchsze > 0)
573 { for (irch = 1; irch <= rchsze; irch++)
574 { node = rchset[irch];
575 marker[node] = 0;
576 }
577 }
578 }
579 return;
580 # undef deg0
581 # undef nhdsze
582 }
583
584 /* eof */