1 /*****************************************************************************
2 /
3 / SPACE (SPArse Cholesky Elimination) Library: gbipart.c
4 /
5 / author J"urgen Schulze, University of Paderborn
6 / created 00dec26
7 /
8 / This file contains functions dealing with bipartite graphs
9 /
10 ******************************************************************************
11
12 Data type: struct gbipart
13 graph_t *G; pointer to graph object with E c X x Y
14 int nX; the vertices 0,...,nX-1 belong to X
15 int nY; the vertices nX,...,nX+nY-1 belong to Y
16 Comments:
17 o Structure used to smooth a separator computed for a subgraph Gbisect.
18 The separator is paired with the border vertices in black/white partition,
19 thus, resulting in a bipartite graph.
20 Methods in lib/gbipart.c:
21 - Gbipart = newBipartiteGraph(int nX, int nY, int nedges);
22 - void freeBipartiteGraph(gbipart_t *Gbipart);
23 - void printGbipart(gbipart_t *Gbipart);
24 - Gbipart = setupBipartiteGraph(graph_t *G, int *bipartvertex, int nX, int nY,
25 int *vtxmap)
26 o Gbipart is induced by the vertices in bipartvertex. The first
27 nX vertices are the vertices 0...nX-1 and the last nY vertices
28 are the vertices nX...nX+nY-1 of Gbipart. Vector vtxmap maps the
29 vertices in bipartvertex to the vertices of the bipartite graph.
30 - void maximumMatching(gbipart_t *Gbipart, int *matching);
31 - void maximumFlow(gbipart_t *Gbipart, int *flow, int *rc)
32 o flow[i] stores the flow over the edge in adjncy[i] of Gbipart. It is
33 positive, if the edge is from X to Y, otherwise flow is negative.
34 o rc[u] stores the residual capacity of edge (source,u), u e X,
35 respectively (u,sink), u e Y. All edges between X and Y have
36 infinite capacity, therefore, no rc value must be computed for them.
37 - void DMviaMatching(gbipart_t *Gbipart, int *matching, int *dmflag,
38 int *dmwght);
39 o on return. vector dmflag is filled with the following values:
40 / SI, iff x e X is reachable via exposed node e X
41 dmflag[x] = < SX, iff x e X is reachable via exposed node e Y
42 \ SR, iff x e X - (SI u SX)
43 / BI, iff y e Y is reachable via exposed node e Y
44 dmflag[y] = < BX, iff y e Y is reachable via exposed node e X
45 \ BR, iff y e Y - (BI u BX)
46 o on return, vector dmwght is filled with the following values:
47 dmwght[SI] - weight of SI dmwght[BI] - weight of BI
48 dmwght[SX] - weight of SX dmwght[BX] - weight of BX
49 dmwght[SR] - weight of SR dmwght[BR] - weight of BR
50 - void DMviaFlow(gbipart_t *Gbipart, int *flow, int *rc, int *dmflag,
51 int *dmwght);
52 o vectors dmflag and dmwght are filled as described above
53
54 ******************************************************************************/
55
56 #include <space.h>
57
58 #define FREE -1
59 #define SOURCE -2
60 #define SINK -3
61
62
63 /*****************************************************************************
64 ******************************************************************************/
65 gbipart_t*
newBipartiteGraph(int nX,int nY,int nedges)66 newBipartiteGraph(int nX, int nY, int nedges)
67 { gbipart_t *Gbipart;
68
69 mymalloc(Gbipart, 1, gbipart_t);
70 Gbipart->G = newGraph(nX+nY, nedges);
71 Gbipart->nX = nX;
72 Gbipart->nY = nY;
73
74 return(Gbipart);
75 }
76
77
78 /*****************************************************************************
79 ******************************************************************************/
80 void
freeBipartiteGraph(gbipart_t * Gbipart)81 freeBipartiteGraph(gbipart_t *Gbipart)
82 {
83 freeGraph(Gbipart->G);
84 free(Gbipart);
85 }
86
87
88 /*****************************************************************************
89 ******************************************************************************/
90 void
printGbipart(gbipart_t * Gbipart)91 printGbipart(gbipart_t *Gbipart)
92 { graph_t *G;
93 int count, u, i, istart, istop;
94
95 G = Gbipart->G;
96 printf("\n#vertices %d (nX %d, nY %d), #edges %d, type %d, totvwght %d\n",
97 G->nvtx, Gbipart->nX, Gbipart->nY, G->nedges >> 1, G->type,
98 G->totvwght);
99 for (u = 0; u < G->nvtx; u++)
100 { count = 0;
101 printf("--- adjacency list of vertex %d (weight %d):\n", u, G->vwght[u]);
102 istart = G->xadj[u];
103 istop = G->xadj[u+1];
104 for (i = istart; i < istop; i++)
105 { printf("%5d", G->adjncy[i]);
106 if ((++count % 16) == 0)
107 printf("\n");
108 }
109 if ((count % 16) != 0)
110 printf("\n");
111 }
112 }
113
114
115 /*****************************************************************************
116 ******************************************************************************/
117 gbipart_t*
setupBipartiteGraph(graph_t * G,int * bipartvertex,int nX,int nY,int * vtxmap)118 setupBipartiteGraph(graph_t *G, int *bipartvertex, int nX, int nY, int *vtxmap)
119 { gbipart_t *Gbipart;
120 int *xadj, *adjncy, *vwght, *xadjGb, *adjncyGb, *vwghtGb;
121 int nvtx, nedgesGb, totvwght, u, x, y, i, j, jstart, jstop, ptr;
122
123 nvtx = G->nvtx;
124 xadj = G->xadj;
125 adjncy = G->adjncy;
126 vwght = G->vwght;
127
128 /* ----------------------------------------------------------------
129 compute number of edges and local indices of vertices in Gbipart
130 ---------------------------------------------------------------- */
131 nedgesGb = 0;
132 for (i = 0; i < nX+nY; i++)
133 { u = bipartvertex[i];
134 if ((u < 0) || (u >= nvtx))
135 { fprintf(stderr, "\nError in function setupBipartiteGraph\n"
136 " node %d does not belong to graph\n", u);
137 quit();
138 }
139 jstart = xadj[u];
140 jstop = xadj[u+1];
141 for (j = jstart; j < jstop; j++)
142 vtxmap[adjncy[j]] = -1;
143 nedgesGb += (jstop - jstart);
144 }
145 for (i = 0; i < nX+nY; i++)
146 { u = bipartvertex[i];
147 vtxmap[u] = i;
148 }
149
150 Gbipart = newBipartiteGraph(nX, nY, nedgesGb);
151 xadjGb = Gbipart->G->xadj;
152 adjncyGb = Gbipart->G->adjncy;
153 vwghtGb = Gbipart->G->vwght;
154
155 /* ---------------------------------
156 build the induced bipartite graph
157 --------------------------------- */
158 totvwght = 0; ptr = 0;
159 for (i = 0; i < nX; i++)
160 { x = bipartvertex[i];
161 xadjGb[i] = ptr;
162 vwghtGb[i] = vwght[x];
163 totvwght += vwght[x];
164 jstart = xadj[x];
165 jstop = xadj[x+1];
166 for (j = jstart; j < jstop; j++)
167 { y = adjncy[j];
168 if (vtxmap[y] >= nX)
169 adjncyGb[ptr++] = vtxmap[y];
170 }
171 }
172 for (i = nX; i < nX+nY; i++)
173 { y = bipartvertex[i];
174 xadjGb[i] = ptr;
175 vwghtGb[i] = vwght[y];
176 totvwght += vwght[y];
177 jstart = xadj[y];
178 jstop = xadj[y+1];
179 for (j = jstart; j < jstop; j++)
180 { x = adjncy[j];
181 if ((vtxmap[x] >= 0) && (vtxmap[x] < nX))
182 adjncyGb[ptr++] = vtxmap[x];
183 }
184 }
185 xadjGb[nX+nY] = ptr;
186 Gbipart->G->type = G->type;
187 Gbipart->G->totvwght = totvwght;
188 return(Gbipart);
189 }
190
191
192 /*****************************************************************************
193 ******************************************************************************/
194 void
maximumMatching(gbipart_t * Gbipart,int * matching)195 maximumMatching(gbipart_t *Gbipart, int *matching)
196 { int *xadj, *adjncy, *level, *marker, *queue, *stack;
197 int top, top2, u, x, x2, y, y2, nX, nY, i, istart, istop;
198 int qhead, qtail, max_level;
199
200 xadj = Gbipart->G->xadj;
201 adjncy = Gbipart->G->adjncy;
202 nX = Gbipart->nX;
203 nY = Gbipart->nY;
204
205 mymalloc(level, (nX+nY), int);
206 mymalloc(marker, (nX+nY), int);
207 mymalloc(queue, nX, int);
208 mymalloc(stack, nY, int);
209
210 /* -------------------
211 initialize matching
212 ------------------- */
213 for (u = 0; u < nX+nY; u++)
214 matching[u] = FREE;
215
216 /* ---------------------------------------------------
217 construct maximal matching in bipartite graph (X,Y)
218 --------------------------------------------------- */
219 for (x = 0; x < nX; x++)
220 { istart = xadj[x];
221 istop = xadj[x+1];
222 for (i = istart; i < istop; i++)
223 { y = adjncy[i];
224 if (matching[y] == FREE)
225 { matching[x] = y;
226 matching[y] = x;
227 break;
228 }
229 }
230 }
231
232 /* --------------------------------------------------------------------
233 construct maximum matching in bipartite graph (X,Y) (Hopcroft, Karp)
234 -------------------------------------------------------------------- */
235 while (TRUE)
236 { for (u = 0; u < nX+nY; u++)
237 level[u] = marker[u] = -1;
238 qhead = qtail = 0; /* fill queue with free X nodes */
239 for (x = 0; x < nX; x++)
240 if (matching[x] == FREE)
241 { queue[qtail++] = x;
242 level[x] = 0;
243 }
244
245 /* --------------------------------------------------------------
246 breadth first search to construct layer network containing all
247 vertex disjoint augmenting paths of minimal length
248 -------------------------------------------------------------- */
249 top = 0;
250 max_level = MAX_INT;
251 while (qhead != qtail)
252 { x = queue[qhead++]; /* note: queue contains only */
253 if (level[x] < max_level) /* nodes from X */
254 { istart = xadj[x];
255 istop = xadj[x+1];
256 for (i = istart; i < istop; i++)
257 { y = adjncy[i];
258 if (level[y] == -1)
259 { level[y] = level[x] + 1;
260 if (matching[y] == FREE)
261 { max_level = level[y]; /* note: stack contains only */
262 stack[top++] = y; /* nodes form Y */
263 }
264 else if (level[y] < max_level)
265 { x2 = matching[y];
266 level[x2] = level[y] + 1;
267 queue[qtail++] = x2;
268 }
269 }
270 }
271 }
272 }
273 if (top == 0) break; /* no augmenting path found */
274
275 /* ------------------------------------------------------------
276 restricted depth first search to construct maximal number of
277 vertex disjoint augmenting paths in layer network
278 ------------------------------------------------------------ */
279 while (top > 0)
280 { top2 = top--;
281 y = stack[top2-1]; /* get the next exposed node in Y */
282 marker[y] = xadj[y]; /* points to next neighbor of y */
283
284 while (top2 > top)
285 { y = stack[top2-1];
286 i = marker[y]++;
287 if (i < xadj[y+1]) /* not all neighbors of y visited */
288 { x = adjncy[i];
289 if ((marker[x] == -1) && (level[x] == level[y]-1))
290 { marker[x] = 0;
291 if (level[x] == 0) /* augmenting path found */
292 while (top2 > top) /* pop stack */
293 { y2 = stack[--top2];
294 x2 = matching[y2]; /* / o == o */
295 matching[x] = y2; /* / */
296 matching[y2] = x; /* x -- y2 == x2 -- y */
297 x = x2; /* \ */
298 } /* \ o == o */
299 else
300 { y2 = matching[x];
301 stack[top2++] = y2;
302 marker[y2] = xadj[y2];
303 }
304 }
305 }
306 else top2--;
307 }
308 }
309 }
310
311 /* -------------------------------
312 free working storage and return
313 ------------------------------- */
314 free(level); free(marker);
315 free(queue); free(stack);
316 }
317
318
319 /*****************************************************************************
320 ******************************************************************************/
321 void
maximumFlow(gbipart_t * Gbipart,int * flow,int * rc)322 maximumFlow(gbipart_t *Gbipart, int *flow, int *rc)
323 { int *xadj, *adjncy, *vwght, *parent, *marker, *queue;
324 int nedges, u, v, x, y, nX, nY, j, i, istart, istop;
325 int qhead, qtail, capacity;
326
327 nedges = Gbipart->G->nedges;
328 xadj = Gbipart->G->xadj;
329 adjncy = Gbipart->G->adjncy;
330 vwght = Gbipart->G->vwght;
331 nX = Gbipart->nX;
332 nY = Gbipart->nY;
333
334 mymalloc(parent, (nX+nY), int);
335 mymalloc(marker, (nX+nY), int);
336 mymalloc(queue, (nX+nY), int);
337
338 /* -------------------------------------
339 initialize flow and residual capacity
340 ------------------------------------- */
341 for (u = 0; u < nX+nY; u++)
342 rc[u] = vwght[u];
343 for (i = 0; i < nedges; i++)
344 flow[i] = 0;
345
346 /* --------------------------------------------------
347 determine an initial flow in the bipartite network
348 -------------------------------------------------- */
349 for (x = 0; x < nX; x++)
350 { istart = xadj[x];
351 istop = xadj[x+1];
352 for (i = istart; i < istop; i++)
353 { y = adjncy[i];
354 capacity = min(rc[x], rc[y]);
355 if (capacity > 0)
356 { rc[x] -= capacity;
357 rc[y] -= capacity;
358 flow[i] = capacity;
359 for (j = xadj[y]; adjncy[j] != x; j++);
360 flow[j] = -capacity;
361 }
362 if (rc[x] == 0) break;
363 }
364 }
365
366 /* -----------------------------------------------------------
367 construct maximum flow in bipartite network (Edmonds, Karp)
368 ----------------------------------------------------------- */
369 while (TRUE)
370 { for (u = 0; u < nX+nY; u++)
371 parent[u] = marker[u] = -1;
372 qhead = qtail = 0; /* fill queue with free X nodes */
373 for (x = 0; x < nX; x++)
374 if (rc[x] > 0)
375 { queue[qtail++] = x;
376 parent[x] = x;
377 }
378
379 /* ---------------------------------------------------------
380 breadth first search to find the shortest augmenting path
381 --------------------------------------------------------- */
382 capacity = 0;
383 while (qhead != qtail)
384 { u = queue[qhead++];
385 istart = xadj[u];
386 istop = xadj[u+1];
387 for (i = istart; i < istop; i++)
388 { v = adjncy[i];
389 if ((parent[v] == -1) && ((v >= nX) || (flow[i] < 0)))
390 /* v >= nX => u->v is a forward edge having infty capacity */
391 /* otherwise u<-v is a backward edge and (v,u) must have */
392 /* positive capacity (i.e. (u,v) has neg. capacity) */
393 { parent[v] = u;
394 marker[v] = i;
395 queue[qtail++] = v;
396 if ((v >= nX) && (rc[v] > 0)) /* found it! */
397 { u = v; /* (v,sink) is below capacity */
398 capacity = rc[u];
399 while (parent[u] != u) /* get minimal residual capa. */
400 { i = marker[u];
401 u = parent[u];
402 if (u >= nX)
403 capacity = min(capacity, -flow[i]);
404 }
405 capacity = min(capacity, rc[u]);
406 rc[v] -= capacity; /* augment flow by min. rc */
407 while (parent[v] != v)
408 { i = marker[v];
409 u = parent[v];
410 flow[i] += capacity;
411 for (j = xadj[v]; adjncy[j] != u; j++);
412 flow[j] = -flow[i];
413 v = u;
414 }
415 rc[v] -= capacity;
416 qhead = qtail; /* escape inner while loop */
417 break;
418 }
419 }
420 }
421 }
422
423 if (capacity == 0)
424 break;
425 }
426
427 free(parent); free(marker);
428 free(queue);
429 }
430
431
432 /*****************************************************************************
433 ******************************************************************************/
434 void
DMviaMatching(gbipart_t * Gbipart,int * matching,int * dmflag,int * dmwght)435 DMviaMatching(gbipart_t *Gbipart, int *matching, int *dmflag, int *dmwght)
436 { int *xadj, *adjncy, *vwght, *queue, qhead, qtail;
437 int u, x, nX, y, nY, i, istart, istop;
438
439 xadj = Gbipart->G->xadj;
440 adjncy = Gbipart->G->adjncy;
441 vwght = Gbipart->G->vwght;
442 nX = Gbipart->nX;
443 nY = Gbipart->nY;
444
445 mymalloc(queue, (nX+nY), int);
446
447 /* ----------------------------------------------------------------------
448 mark all exposed nodes of X with SI and all exposed nodes of Y with BI
449 ---------------------------------------------------------------------- */
450 qhead = qtail = 0;
451 for (x = 0; x < nX; x++)
452 if (matching[x] == FREE)
453 { queue[qtail++] = x;
454 dmflag[x] = SI;
455 }
456 else dmflag[x] = SR;
457 for (y = nX; y < nX+nY; y++)
458 if (matching[y] == FREE)
459 { queue[qtail++] = y;
460 dmflag[y] = BI;
461 }
462 else dmflag[y] = BR;
463
464 /* ------------------------------------------------------------------
465 construct Dulmage-Mendelsohn decomp. starting with SI and BI nodes
466 ------------------------------------------------------------------ */
467 while (qhead != qtail)
468 { u = queue[qhead++];
469 istart = xadj[u];
470 istop = xadj[u+1];
471 switch(dmflag[u])
472 { case SI:
473 for (i = istart; i < istop; i++)
474 { y = adjncy[i];
475 if (dmflag[y] == BR)
476 { queue[qtail++] = y;
477 dmflag[y] = BX;
478 }
479 }
480 break;
481 case BX:
482 x = matching[u];
483 dmflag[x] = SI;
484 queue[qtail++] = x;
485 break;
486 case BI:
487 for (i = istart; i < istop; i++)
488 { x = adjncy[i];
489 if (dmflag[x] == SR)
490 { queue[qtail++] = x;
491 dmflag[x] = SX;
492 }
493 }
494 break;
495 case SX:
496 y = matching[u];
497 dmflag[y] = BI;
498 queue[qtail++] = y;
499 break;
500 }
501 }
502
503 /* ----------------------
504 fill the dmwght vector
505 ---------------------- */
506 dmwght[SI] = dmwght[SX] = dmwght[SR] = 0;
507 for (x = 0; x < nX; x++)
508 switch(dmflag[x])
509 { case SI: dmwght[SI] += vwght[x]; break;
510 case SX: dmwght[SX] += vwght[x]; break;
511 case SR: dmwght[SR] += vwght[x]; break;
512 }
513 dmwght[BI] = dmwght[BX] = dmwght[BR] = 0;
514 for (y = nX; y < nX+nY; y++)
515 switch(dmflag[y])
516 { case BI: dmwght[BI] += vwght[y]; break;
517 case BX: dmwght[BX] += vwght[y]; break;
518 case BR: dmwght[BR] += vwght[y]; break;
519 }
520
521 free(queue);
522 }
523
524
525 /*****************************************************************************
526 ******************************************************************************/
527 void
DMviaFlow(gbipart_t * Gbipart,int * flow,int * rc,int * dmflag,int * dmwght)528 DMviaFlow(gbipart_t *Gbipart, int *flow, int *rc, int *dmflag, int *dmwght)
529 { int *xadj, *adjncy, *vwght, *queue, qhead, qtail;
530 int u, v, x, nX, y, nY, i, istart, istop;
531
532 xadj = Gbipart->G->xadj;
533 adjncy = Gbipart->G->adjncy;
534 vwght = Gbipart->G->vwght;
535 nX = Gbipart->nX;
536 nY = Gbipart->nY;
537
538 mymalloc(queue, (nX+nY), int);
539
540 /* ----------------------------------------------------------
541 mark all nodes reachable from source/sink with SOURCE/SINK
542 ---------------------------------------------------------- */
543 qhead = qtail = 0;
544 for (x = 0; x < nX; x++)
545 if (rc[x] > 0)
546 { queue[qtail++] = x;
547 dmflag[x] = SOURCE;
548 }
549 else dmflag[x] = FREE;
550 for (y = nX; y < nX+nY; y++)
551 if (rc[y] > 0)
552 { queue[qtail++] = y;
553 dmflag[y] = SINK;
554 }
555 else dmflag[y] = FREE;
556
557 /* --------------------------------------------------------------------
558 construct Dulmage-Mendelsohn decomp. starting with SOURCE/SINK nodes
559 -------------------------------------------------------------------- */
560 while (qhead != qtail)
561 { u = queue[qhead++];
562 istart = xadj[u];
563 istop = xadj[u+1];
564 switch(dmflag[u])
565 { case SOURCE:
566 for (i = istart; i < istop; i++)
567 { v = adjncy[i];
568 if ((dmflag[v] == FREE) && ((v >= nX) || (flow[i] < 0)))
569 { queue[qtail++] = v;
570 dmflag[v] = SOURCE; /* v reachable via forward edge u->v */
571 } /* or via backward edge u<-v */
572 }
573 break;
574 case SINK:
575 for (i = istart; i < istop; i++)
576 { v = adjncy[i];
577 if ((dmflag[v] == FREE) && ((v < nX) || (flow[i] > 0)))
578 { queue[qtail++] = v;
579 dmflag[v] = SINK; /* u reachable via forward edge v->u */
580 } /* or via backward edge v<-u */
581 }
582 break;
583 }
584 }
585
586 /* -----------------------------------------------------
587 all nodes x in X with dmflag[x] = SOURCE belong to SI
588 all nodes x in X with dmflag[x] = SINK belong to SX
589 all nodes x in X with dmflag[x] = FREE belong to SR
590 ----------------------------------------------------- */
591 dmwght[SI] = dmwght[SX] = dmwght[SR] = 0;
592 for (x = 0; x < nX; x++)
593 switch(dmflag[x])
594 { case SOURCE: dmflag[x] = SI; dmwght[SI] += vwght[x]; break;
595 case SINK: dmflag[x] = SX; dmwght[SX] += vwght[x]; break;
596 default: dmflag[x] = SR; dmwght[SR] += vwght[x];
597 }
598
599 /* -----------------------------------------------------
600 all nodes y in Y with dmflag[y] = SOURCE belong to BX
601 all nodes y in Y with dmflag[y] = SINK belong to BI
602 all nodes y in Y with dmflag[y] = FREE belong to BR
603 ----------------------------------------------------- */
604 dmwght[BI] = dmwght[BX] = dmwght[BR] = 0;
605 for (y = nX; y < nX+nY; y++)
606 switch(dmflag[y])
607 { case SOURCE: dmflag[y] = BX; dmwght[BX] += vwght[y]; break;
608 case SINK: dmflag[y] = BI; dmwght[BI] += vwght[y]; break;
609 default: dmflag[y] = BR; dmwght[BR] += vwght[y];
610 }
611
612 free(queue);
613 }
614
615