Scippy

SCIP

Solving Constraint Integer Programs

graph_mcut.c
Go to the documentation of this file.
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the program and library */
4 /* SCIP --- Solving Constraint Integer Programs */
5 /* */
6 /* Copyright (C) 2002-2022 Konrad-Zuse-Zentrum */
7 /* fuer Informationstechnik Berlin */
8 /* */
9 /* SCIP is distributed under the terms of the ZIB Academic License. */
10 /* */
11 /* You should have received a copy of the ZIB Academic License */
12 /* along with SCIP; see the file COPYING. If not visit scipopt.org. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15 
16 /**@file graph_mcut.c
17  * @brief Minimum cut routine for Steiner problems
18  * @author Gerald Gamrath
19  * @author Thorsten Koch
20  * @author Daniel Rehfeldt
21  *
22  * This file implements a graph minimum cut routine for Steiner problems. For more details see \ref STP_MINCUT page.
23  *
24  * @page STP_MINCUT Graph minimum cut routine
25  *
26  * The implemented algorithm is described in "A Faster Algorithm for Finding the Minimum Cut in a Graph" by Hao and Orlin.
27  *
28  * A list of all interface methods can be found in graph.h.
29  */
30 
31 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
32 
33 
34 #include <stdio.h>
35 #include <assert.h>
36 #include <stdlib.h>
37 #include "portab.h"
38 #include "graph.h"
39 
40 #define DEBUG 0 /* 0 = No, 1 = Validation, 2 = Show flow */
41 #define STATIST 0
42 #define CHECK 0
43 
44 #ifndef RESTRICT
45 #define RESTRICT restrict
46 #endif
47 
48 #ifdef NDEBUG
49 #undef STATIST
50 #undef DEBUG
51 #define STATIST 0
52 #define DEBUG 0
53 #endif
54 
55 #define Q_NULL -1 /* NULL element of queue/list */
56 #define GLOBALRELABEL_ADD 10 /* constant for global relabel heuristic */
57 #define GLOBALRELABEL_MULT 10 /* constant for global relabel heuristic */
58 
59 /** remove first element from active (singly-linked) list */
60 #define listdelete(node,nodedist,next,headactive)\
61 {\
62  assert(node >= 0);\
63  assert(nodedist >= 0);\
64  assert(headactive != NULL);\
65  headactive[nodedist] = next[node];\
66 }
67 
68 /** add element to active (singly-linked) list */
69 #define listinsert(node,nodedist,next,headactive,actmin,actmax)\
70 {\
71  assert(node >= 0);\
72  assert(next != NULL);\
73  assert(headactive != NULL);\
74  assert(nodedist >= 0);\
75  next[node] = headactive[nodedist];\
76  headactive[nodedist] = node;\
77  if( nodedist < (actmin) )\
78  actmin = nodedist;\
79  if( nodedist > (actmax) )\
80  actmax = nodedist;\
81 }
82 #if 0
83 /** add element to active (singly-linked) list */
84 #define listinsert(node,nodedist,next,headactive,actmin,actmax)\
85 {\
86  assert(node >= 0);\
87  assert(next != NULL);\
88  assert(headactive != NULL);\
89  assert(nodedist >= 0);\
90  next[node] = headactive[nodedist];\
91  headactive[nodedist] = node;\
92  if( nodedist < (actmin) )\
93  actmin = nodedist;\
94  if( nodedist > (actmax) )\
95  actmax = nodedist;\
96 }
97 #endif
98 
99 
100 /** remove first element from active (singly-linked) list */
101 #define activedelete(node,nodedist,next,headactive)\
102 {\
103  assert(node >= 0);\
104  assert(nodedist >= 0);\
105  assert(headactive != NULL);\
106  headactive[nodedist] = next[node];\
107 }
108 
109 /** add element to active (singly-linked) list */
110 #define activeinsert(node,nodedist,next,headactive,actmin,actmax,glbmax)\
111 {\
112  assert(node >= 0);\
113  assert(next != NULL);\
114  assert(headactive != NULL);\
115  assert(nodedist >= 0);\
116  next[node] = headactive[nodedist];\
117  headactive[nodedist] = node;\
118  if( nodedist < (actmin) )\
119  actmin = nodedist;\
120  if( nodedist > (actmax) )\
121  actmax = nodedist;\
122  if( (actmax) > (glbmax) )\
123  glbmax = (actmax);\
124 }
125 
126 /** remove element from inactive (doubly-linked) list */
127 #define inactivedelete(node,nodedist,next,prev,headinactive,nextnode,prevnode)\
128 {\
129  assert(node >= 0);\
130  assert(nodedist >= 0);\
131  assert(next != NULL);\
132  assert(prev != NULL);\
133  assert(headinactive != NULL);\
134  nextnode = next[node];\
135  if( headinactive[nodedist] == node )\
136  {\
137  headinactive[nodedist] = nextnode;\
138  if( nextnode >= 0 )\
139  prev[nextnode] = Q_NULL;\
140  }\
141  else\
142  {\
143  prevnode = prev[node];\
144  assert(prevnode >= 0);\
145  assert(next[prevnode] == node);\
146  next[prevnode] = nextnode;\
147  if( nextnode >= 0 )\
148  prev[nextnode] = prevnode;\
149  }\
150 }
151 
152 /** add element to inactive (doubly-linked) list */
153 #define inactiveinsert(node,nodedist,next,prev,headinactive,nextnode)\
154 {\
155  assert(node >= 0);\
156  assert(nodedist >= 0);\
157  assert(next != NULL);\
158  assert(prev != NULL);\
159  assert(headinactive != NULL);\
160  nextnode = headinactive[nodedist];\
161  next[node] = nextnode;\
162  prev[node] = Q_NULL;\
163  if( nextnode >= 0 )\
164  prev[nextnode] = node;\
165  headinactive[nodedist] = node;\
166 }
167 
168 
169 #if DEBUG
170 static int is_valid(const GRAPH*, const int, const int, const int*, const int *);
171 static void show_flow(const GRAPH*, const int*, const int*);
172 #endif
173 #if STATIST
174 static void cut_statist(void);
175 static void cut_sum(const GRAPH*, const int*, const int*);
176 #endif
177 
178 #if STATIST
179 static int s_pushes = 0;
180 static int n_pushes = 0;
181 static int m_pushes = 0;
182 static int x_pushes = 0;
183 static int relabels = 0;
184 static int s_sleeps = 0;
185 static int m_sleeps = 0;
186 static int searches = 0;
187 static int cutsums = 0;
188 #endif
189 
190 
191 #if 1
192 
193 #if CHECK
194 /** checker */
195 static int is_valid_arr(
196  const GRAPH* p,
197  const int s,
198  const int t,
199  const int* startarr,
200  const int* headarr,
201  const int* revedgearr,
202  const int* resarr,
203  const int* capa,
204  const int* w)
205 {
206  int* e;
207  int* r;
208  int* dist;
209 
210  int j;
211  int k;
212  int end;
213  int head;
214 
215  assert(p != NULL);
216  assert(p->mincut_e != NULL);
217  assert(p->mincut_r != NULL);
218 
219  e = p->mincut_e;
220  r = p->mincut_r;
221  dist = p->mincut_dist;
222 
223  for(j = 0; j < p->knots; j++)
224  {
225  if( e[j] < 0 )
226  return ((void) fprintf(stderr, "Negativ Execess Node %d\n", j), FALSE);
227 
228  if( dist[j] >= p->knots )
229  return ((void) fprintf(stderr, "Distance too big Node %d dist: %d\n", j, dist[j]), FALSE);
230 
231  /* extended dormacy property */
232  for( k = startarr[j], end = startarr[j + 1]; k != end; k++ )
233  {
234  head = headarr[k];
235 
236  if( r[k] > 0 )
237  {
238  if( dist[j] > dist[head] + 1 && w[j] == 0 && w[head] == 0 )
239  {
240  printf("distance fail! %d->%d (%d)\n", j, head, k);
241  return FALSE;
242  }
243 
244  if( (w[j] && (w[j] < w[head])) )
245  {
246  printf("Extended Dormacy Violation %d %d \n", j, head);
247  }
248 
249  if( w[j] && !w[head] )
250  {
251  printf("Simple Dormacy Violation %d %d \n", j, head);
252  }
253 
254  if( (w[j] && !w[head]) || (w[j] && (w[j] < w[head])) )
255  {
256  (void) printf("e=%d r[e]=%d head=%d tail=%d w[h]=%d w[t]=%d\n",
257  k, r[k], head, p->tail[k], w[head],
258  w[p->tail[k]]);
259 
260  return ((void) fprintf(stderr,
261  "Dormacy Violation Knot %d\n", j), FALSE);
262  }
263  }
264 
265  if( r[k] < 0 )
266  return ((void) fprintf(stderr, "Negativ Residual Edge %d\n", k), FALSE);
267  }
268  }
269 
270  return(TRUE);
271 }
272 #endif
273 
274 #if 0
275 /** backwards bfs from t along arcs of cap > 0 and set distance labels according to distance from t */
276 static int bfs(
277  const GRAPH* p,
278  const int s,
279  const int t,
280  int* dist,
281  int* temp,
282  int* w,
283  const int* capa,
284  const int* edgestart,
285  const int* edgearr,
286  const int* headarr
287  )
288 {
289  int i;
290  int j;
291  int k;
292  int l;
293  int end;
294  int nnodes;
295  int dplus1;
296  int visited = 0;
297 
298  assert(temp != NULL);
299  assert(dist != NULL);
300  assert(w != NULL);
301  assert(s >= 0);
302  assert(s < p->knots);
303  assert(t >= 0);
304  assert(t < p->knots);
305 
306  /* initialize */
307  temp[visited++] = t;
308  dist[t] = 0;
309 
310  nnodes = p->knots;
311 
312  /* bfs loop */
313  for( j = 0; j < visited; j++ )
314  {
315  assert(visited <= nnodes);
316 
317  i = temp[j];
318  dplus1 = dist[i] + 1;
319 
320  assert(i >= 0);
321  assert(i < nnodes);
322  assert(dist[i] >= 0);
323  assert(dist[i] < visited);
324  assert(w[i] == 0);
325 
326  for( l = edgestart[i], end = edgestart[i + 1]; l != end; l++ )
327  {
328  k = headarr[l];
329 
330  /* not visited yet? */
331  if( dist[k] < 0 && w[k] == 0 )
332  {
333  assert(w[k] == 0);
334 
335  dist[k] = dplus1;
336  temp[visited++] = k;
337  assert(dist[k] < nnodes);
338  }
339  }
340  } /* bfs loop */
341 
342  return (visited);
343 }
344 #endif
345 
346 
347 /** initializes minimum cut arrays */
348 static
349 SCIP_RETCODE mincutInit(
350  SCIP* scip, /**< SCIP data structure */
351  int nnodes, /**< number of nodes */
352  int nedges, /**< number of edges */
353  GRAPH* p /**< graph data structure */
354 )
355 {
356  assert(nnodes > 0 && nedges > 0);
357 
358  p->mincut_nnodes = nnodes;
359  p->mincut_nedges = nedges;
360  SCIP_CALL( SCIPallocMemoryArray(scip, &(p->mincut_dist), nnodes + 1) );
361  SCIP_CALL( SCIPallocMemoryArray(scip, &(p->mincut_head), nnodes + 1) );
362  SCIP_CALL( SCIPallocMemoryArray(scip, &(p->mincut_head_inact), nnodes + 1) );
363  SCIP_CALL( SCIPallocMemoryArray(scip, &(p->mincut_numb), nnodes + 1) );
364  SCIP_CALL( SCIPallocMemoryArray(scip, &(p->mincut_prev), nnodes + 1) );
365  SCIP_CALL( SCIPallocMemoryArray(scip, &(p->mincut_next), nnodes + 1) );
366  SCIP_CALL( SCIPallocMemoryArray(scip, &(p->mincut_temp), nnodes + 1) );
367  SCIP_CALL( SCIPallocMemoryArray(scip, &(p->mincut_e), nnodes + 1) );
368  SCIP_CALL( SCIPallocMemoryArray(scip, &(p->mincut_r), nedges) );
369 
370  return SCIP_OKAY;
371 }
372 
373 /** global relabel heuristic that sets distance of sink to zero and relabels all other nodes using backward bfs on residual
374  * graph, starting from the sink. */
375 static void globalrelabel(
376  const int s,
377  const int t,
378  const int nnodesreal,
379  int* RESTRICT dist,
380  int* RESTRICT headactive,
381  int* RESTRICT headinactive,
382  int* RESTRICT edgecurr,
383  int* RESTRICT next,
384  int* RESTRICT prev,
385  int* RESTRICT excess,
386  int* RESTRICT residual,
387  int* RESTRICT w,
388  const int* edgestart,
389  const int* edgearr,
390  const int* headarr,
391  int* actmin,
392  int* actmax,
393  int* glbmax,
394  int* dormmax
395  )
396 {
397  int currdist;
398  int nextnode;
399  int actmaxloc;
400  int actminloc;
401  int glbmaxloc;
402  int dormmaxnext;
403  const int nnodes = nnodesreal;
404  const int end = *glbmax;
405  SCIP_Bool hit;
406 
407  assert(s >= 0);
408  assert(s < nnodesreal);
409  assert(t >= 0);
410  assert(t < nnodesreal);
411  assert(s != t);
412  assert(w != NULL);
413  assert(headactive != NULL);
414  assert(headinactive != NULL);
415  assert(edgecurr != NULL);
416  assert(prev != NULL);
417  assert(next != NULL);
418  assert(residual != NULL);
419  assert(excess != NULL);
420  assert(edgearr != NULL);
421  assert(headarr != NULL);
422  assert(edgecurr != NULL);
423  assert(edgestart != NULL);
424  assert(headactive != NULL);
425  assert(headinactive != NULL);
426 
427  assert(end < nnodes);
428  assert(w[t] == 0);
429 
430  for( int i = 0; i <= end; i++ )
431  {
432  headactive[i] = Q_NULL;
433  headinactive[i] = Q_NULL;
434  if( w[i] == 0 )
435  w[i] = -1;
436  }
437 
438 #ifndef NDEBUG
439  for( int i = 0; i < nnodes; i++ )
440  assert(headactive[i] == Q_NULL && headinactive[i] == Q_NULL);
441 #endif
442 
443  for( int i = end + 1; i < nnodes; i++ )
444  {
445  if( w[i] == 0 )
446  w[i] = -1;
447  }
448 
449  assert(w[t] == -1);
450 
451  w[t] = 0;
452  dist[t] = 0;
453 
454  glbmaxloc = 0;
455  actmaxloc = 0;
456  actminloc = nnodes;
457 
458  /* add t to inactive list */
459  inactiveinsert(t, 0, next, prev, headinactive, nextnode);
460 
461  /* bfs loop */
462  for( currdist = 0; TRUE; currdist++ )
463  {
464  const int nextdist = currdist + 1;
465 
466  /* no more nodes with current distance? */
467  if( (headactive[currdist] < 0) && (headinactive[currdist] < 0) )
468  break;
469 
470  /* check inactive nodes */
471  for( int i = headinactive[currdist]; i >= 0; i = next[i] )
472  {
473  const int edges_start = edgestart[i];
474  const int edges_end = edgestart[i + 1];
475  for( int q = edges_start; q != edges_end; q++ )
476  {
477  if( residual[edgearr[q]] != 0 )
478  {
479  if( w[headarr[q]] < 0 )
480  {
481  const int k = headarr[q];
482 
483  w[k] = 0;
484 
485  dist[k] = nextdist;
486  edgecurr[k] = edgestart[k];
487 
488  if( excess[k] > 0 )
489  {
490  activeinsert(k, nextdist, next, headactive, actminloc, actmaxloc, glbmaxloc);
491  }
492  else
493  {
494  inactiveinsert(k, nextdist, next, prev, headinactive, nextnode);
495  }
496  }
497  }
498  }
499  }
500 
501  /* check active nodes */
502  for( int i = headactive[currdist]; i >= 0; i = next[i] )
503  {
504  const int edges_start = edgestart[i];
505  const int edges_end = edgestart[i + 1];
506  for( int q = edges_start; q != edges_end; q++ )
507  {
508  if( residual[edgearr[q]] != 0 )
509  {
510  const int k = headarr[q];
511 
512  if( w[k] < 0 )
513  {
514  w[k] = 0;
515 
516  dist[k] = nextdist;
517  edgecurr[k] = edgestart[k];
518 
519  if( excess[k] > 0 )
520  {
521  activeinsert(k, nextdist, next, headactive, actminloc, actmaxloc, glbmaxloc);
522  }
523  else
524  {
525  inactiveinsert(k, nextdist, next, prev, headinactive, nextnode);
526  }
527  }
528  }
529  }
530  }
531  }
532 
533  assert((headactive[currdist] == Q_NULL) && (headinactive[currdist] == Q_NULL));
534  assert(currdist == 0 || (headactive[currdist - 1] != Q_NULL) || (headinactive[currdist - 1] != Q_NULL));
535 
536  /* set global maximum label */
537  if( (currdist - 1) > (glbmaxloc) )
538  glbmaxloc = currdist - 1;
539 
540  hit = FALSE;
541 
542  assert(w[t] == 0);
543 
544  dormmaxnext = (*dormmax) + 1;
545  for( int i = 0; i < nnodes; i++ )
546  {
547  if( w[i] < 0 )
548  {
549  w[i] = dormmaxnext;
550  hit = TRUE;
551  }
552  }
553  if( hit )
554  *dormmax = dormmaxnext;
555 
556  *glbmax = glbmaxloc;
557  *actmin = actminloc;
558  *actmax = actmaxloc;
559 
560  return;
561 }
562 
563 
564 /** initialize data for the first max-flow run */
565 static
566 void initialise(
567  const int s,
568  const int t,
569  const int nnodesreal,
570  const int rootcutsize,
571  const int* rootcut,
572  const int* capa,
573  int* RESTRICT dist,
574  int* RESTRICT headactive,
575  int* RESTRICT headinactive,
576  int* RESTRICT edgecurr,
577  int* RESTRICT next,
578  int* RESTRICT prev,
579  int* RESTRICT temp,
580  int* RESTRICT excess,
581  int* RESTRICT residual,
582  int* RESTRICT w,
583  const int* edgestart,
584  const int* edgearr,
585  const int* headarr,
586  int* dormmax,
587  int* actmin,
588  int* actmax,
589  int* glbmax
590  )
591 {
592 #if 0
593  int k;
594  int e;
595  int q;
596  int end;
597  int head;
598  int nedges;
599 #endif
600  const int nnodes = nnodesreal;
601  int nextnode;
602  int actmaxloc;
603  int actminloc;
604  int glbmaxloc;
605 
606  assert(s >= 0);
607  assert(s < nnodes);
608  assert(t >= 0);
609  assert(t < nnodes);
610  assert(s != t);
611  assert(w != NULL);
612  assert(headactive != NULL);
613  assert(headinactive != NULL);
614  assert(prev != NULL);
615  assert(next != NULL);
616  assert(temp != NULL);
617  assert(residual != NULL);
618  assert(excess != NULL);
619  assert(edgearr != NULL);
620  assert(headarr != NULL);
621  assert(edgecurr != NULL);
622  assert(edgestart != NULL);
623  assert(headactive != NULL);
624  assert(headinactive != NULL);
625 
626 #if 0
627  nedges = p->edges;
628 #endif
629 
630  actmaxloc = 0;
631  actminloc = nnodes;
632  glbmaxloc = 0;
633 
634  *dormmax = 1;
635 
636  assert(w[s] > 0);
637  assert(w[t] == 0);
638 
639  /* set distance labels, and add nodes to lists */
640  for( int i = nnodes - 1; i >= 0; i-- )
641  {
642  dist[i] = 1;
643 
644  /* in root cut? */
645  if( w[i] > 0 )
646  continue;
647 
648  /* sink? */
649  if( i == t )
650  {
651  dist[i] = 0;
652  inactiveinsert(i, 0, next, prev, headinactive, nextnode);
653  continue;
654  }
655 
656  assert(w[i] == 0);
657 
658  if( excess[i] > 0 )
659  {
660  /* add to active list */
661  activeinsert(i, 1, next, headactive, actminloc, actmaxloc, glbmaxloc);
662  }
663  else
664  {
665  /* add to inactive list */
666  inactiveinsert(i, 1, next, prev, headinactive, nextnode);
667  }
668  }
669 
670  glbmaxloc = 1;
671 
672 #if 0 /* both */
673  assert(nnodes <= p->edges);
674 
675  for( i = 0; i < nnodes; i++ )
676  {
677  dist[i] = -1;
678  excess[i] = 0;
679  edgecurr[i] = edgestart[i];
680  headactive[i] = Q_NULL;
681  headinactive[i] = Q_NULL;
682 
683  residual[i] = capa[i];
684  }
685 
686  for( e = nnodes; e < nedges; e++ )
687  residual[e] = capa[e];
688 #endif
689 
690 
691 #if 0 /* simple */
692 
693  /* push from source */
694  for( q = edgestart[s], end = edgestart[s + 1]; q != end; q++ )
695  {
696  e = edgearr[q];
697 
698  if( residual[e] == 0 )
699  continue;
700 
701  i = headarr[q];
702 
703  assert(w[i] == 0);
704  assert(residual[e] == capa[e]);
705 
706  residual[flipedge(e)] += residual[e];
707  excess[i] += residual[e];
708  residual[e] = 0; /* -= residual[e] */
709  }
710 
711  /* set distance labels, and add nodes to lists */
712  for( i = 0; i < nnodes; i++ )
713  {
714  /* sink? */
715  if( i == t )
716  {
717  dist[i] = 0;
718  inactiveinsert(i, 0, next, prev, headinactive, nextnode);
719  continue;
720  }
721 
722  dist[i] = 1;
723 
724  /* source? */
725  if( i == s )
726  {
727  assert(excess[s] == 0);
728 
729  /* put s (source) into lowest dormant set */
730  w[i] = 1; /* *dormmax */
731  continue;
732  }
733 
734  assert(w[i] == 0);
735 
736  if( excess[i] > 0 )
737  {
738  /* add to active list */
739  activeinsert(i, 1, next, headactive, actminloc, actmaxloc, glbmaxloc);
740  }
741  else
742  {
743  /* add to inactive list */
744  inactiveinsert(i, 1, next, prev, headinactive, nextnode);
745  }
746  }
747  glbmaxloc = 1;
748 #endif
749 
750 #if 0 /* with initial global relabel */
751  /* backward bfs from t */
752  (void)bfs(p, s, t, dist, temp, w, capa, edgestart, edgearr, headarr);
753 
754  /* put unreachable nodes to sleep */
755  for( i = 0; i < nnodes; i++ )
756  if( dist[i] < 0 )
757  {
758  if( (*dormmax) == 1 )
759  *dormmax = 2;
760  w[i] = 2;
761  }
762 
763  /* put s (source) into lowest dormant set */
764  w[s] = 1; /* *dormmax */
765 
766  /* s not reached? */
767  if( dist[s] < 0 )
768  return;
769 
770  assert(w[s] > 0);
771  assert(dist[s] > 0);
772  assert(w[t] == 0);
773  assert(dist[t] == 0);
774 
775  /* push from source */
776  for( q = edgestart[s], end = edgestart[s + 1]; q != end; q++ )
777  {
778  e = edgearr[q];
779 
780  if( residual[e] == 0 )
781  continue;
782 
783  i = headarr[q];
784 
785  assert(residual[e] == capa[e]);
786 
787  residual[flipedge(e)] += residual[e];
788  excess[i] += residual[e];
789  residual[e] = 0; /* -= residual[e] */
790  }
791 
792  w[t] = 1;
793  inactiveinsert(t, 0, next, prev, headinactive, nextnode);
794 
795  for( i = 0; i < nnodes; i++ )
796  {
797  /* including t and s */
798  if( w[i] )
799  continue;
800 
801  assert(w[i] == 0);
802  assert(s != i);
803 
804  if( excess[i] > 0 )
805  {
806  /* add to active list */
807  activeinsert(i, dist[i], next, headactive, actminloc, actmaxloc, glbmaxloc);
808  }
809  else
810  {
811  if( dist[i] > (glbmaxloc) )
812  glbmaxloc = dist[i];
813 
814  /* add to inactive list */
815  inactiveinsert(i, dist[i], next, prev, headinactive, nextnode);
816  }
817  }
818 
819  w[t] = 0;
820 #endif
821 
822  *glbmax = glbmaxloc;
823  *actmax = actmaxloc;
824  *actmin = actminloc;
825 
826 #if DEBUG > 0
827  assert(is_valid(p, s, t, capa, w));
828 #endif
829 }
830 
831 /** initialize data for the repeated max-flow run */
832 static
833 void reinitialise(
834  const int s,
835  const int t,
836  const int nnodesreal,
837  const int rootcutsize,
838  const int* rootcut,
839  const int* capa,
840  int* RESTRICT dist,
841  int* RESTRICT headactive,
842  int* RESTRICT headinactive,
843  int* RESTRICT edgecurr,
844  int* RESTRICT next,
845  int* RESTRICT prev,
846  int* RESTRICT temp,
847  int* RESTRICT excess,
848  int* RESTRICT residual,
849  int* RESTRICT w,
850  const int* edgestart,
851  const int* edgearr,
852  const int* headarr,
853  int* dormmax,
854  int* actmin,
855  int* actmax,
856  int* glbmax
857  )
858 {
859  const int nnodes = nnodesreal;
860  int visited;
861 #if 0
862  int a;
863  int headnode;
864 #endif
865  int nextnode;
866  int actmaxloc;
867  int actminloc;
868  int glbmaxloc;
869  int dormmaxloc;
870  int dormmaxlocp1;
871  SCIP_Bool hit;
872 
873  assert(s >= 0);
874  assert(s < nnodes);
875  assert(t >= 0);
876  assert(t < nnodes);
877  assert(s != t);
878  assert(w != NULL);
879  assert(dist != NULL);
880  assert(prev != NULL);
881  assert(next != NULL);
882  assert(temp != NULL);
883  assert(residual != NULL);
884  assert(excess != NULL);
885  assert(edgearr != NULL);
886  assert(headarr != NULL);
887  assert(edgecurr != NULL);
888  assert(edgestart != NULL);
889  assert(headactive != NULL);
890  assert(headinactive != NULL);
891 
892  /* initialize */
893 
894  assert(w[s] == 1);
895 
896  dormmaxloc = 1;
897  actmaxloc = 0;
898  actminloc = nnodes;
899  glbmaxloc = 0;
900 
901  /* t already awake? */
902  if( w[t] == 0 )
903  {
904  for( int i = nnodes - 1; i >= 0; i-- )
905  {
906  headactive[i] = Q_NULL;
907  headinactive[i] = Q_NULL;
908 
909  /* reset distance label for awake nodes and adapt incident edges */
910  if( w[i] == 0 )
911  {
912  dist[i] = -1;
913  edgecurr[i] = edgestart[i];
914  }
915  else if( w[i] > (dormmaxloc) )
916  {
917  dormmaxloc = w[i];
918  }
919  }
920  }
921  else
922  {
923  const int wt = w[t];
924  for( int i = nnodes - 1; i >= 0; i-- )
925  {
926  headactive[i] = Q_NULL;
927  headinactive[i] = Q_NULL;
928 
929  /* wake up nodes in higher or equal dormant layer */
930  if( (w[i] == 0) || (w[i] >= wt) )
931  {
932  w[i] = 0;
933  dist[i] = -1;
934  edgecurr[i] = edgestart[i];
935  }
936  else if( w[i] > (dormmaxloc) )
937  {
938  dormmaxloc = w[i];
939  }
940  }
941  }
942 
943  /* backward bfs from t */
944  visited = 0;
945 
946  inactiveinsert(t, 0, next, prev, headinactive, nextnode);
947 
948  /* initialize */
949  temp[visited++] = t;
950  dist[t] = 0;
951 
952  /* bfs loop */
953  for( int j = 0; j < visited; j++ )
954  {
955  const int i = temp[j];
956 
957  assert(visited <= nnodes);
958  assert(i >= 0);
959  assert(i < nnodes);
960  assert(dist[i] >= 0);
961  assert(dist[i] < visited);
962  assert(w[i] == 0);
963 
964  for( int outedge = edgestart[i], end = edgestart[i + 1]; outedge != end; outedge++ )
965  {
966  const int k = headarr[outedge];
967 
968  /* not visited yet? */
969  if( dist[k] < 0 && w[k] == 0 )
970  {
971  if( residual[edgearr[outedge]] > 0 )
972  {
973  dist[k] = dist[i] + 1;
974  temp[visited++] = k;
975  assert(dist[k] < nnodes);
976 
977  if( excess[k] > 0 )
978  {
979  /* add to active list */
980  activeinsert(k, dist[k], next, headactive, actminloc, actmaxloc, glbmaxloc);
981  }
982  else
983  {
984  if( dist[k] > (glbmaxloc) )
985  glbmaxloc = dist[k];
986 
987  /* add to inactive list */
988  inactiveinsert(k, dist[k], next, prev, headinactive, nextnode);
989  }
990  }
991  }
992  }
993  } /* bfs loop */
994 
995  hit = FALSE;
996  dormmaxlocp1 = dormmaxloc + 1;
997  for( int i = nnodes - 1; i >= 0; i-- )
998  {
999  /* unreachable non-dormant node? */
1000  if( dist[i] < 0 && w[i] == 0 )
1001  {
1002  dist[i] = -1;
1003  w[i] = dormmaxlocp1;
1004  hit = TRUE;
1005  }
1006  }
1007 
1008  if( hit )
1009  dormmaxloc++;
1010 #if 0
1011  for( i = nnodes - 1; i >= 0; i-- )
1012  {
1013  /* is the node awake? */
1014  if( w[i] == 0 )
1015  {
1016  /*
1017  * update excess and residual weights
1018  * */
1019 
1020  excess[i] = 0;
1021 
1022  /* iterate over outgoing arcs of i */
1023  for( j = edgestart[i], end = edgestart[i + 1]; j != end; j++ )
1024  {
1025  headnode = headarr[j];
1026  a = edgearr[j];
1027 
1028  if( w[headnode] )
1029  {
1030  assert(residual[flipedge(a)] == 0);
1031  excess[i] += capa[flipedge(a)];
1032  }
1033  else
1034  {
1035  residual[a] = capa[a];
1036  }
1037  }
1038  }
1039  }
1040 
1041  assert(w[t] == 0);
1042  assert(dist[t] == 0);
1043  assert(temp[0] == t);
1044 
1045 
1046  for( i = 1; i < visited && 0; i++ )
1047  {
1048  k = temp[i];
1049  assert(w[k] == 0);
1050  assert(k != s);
1051  assert(k != t);
1052 
1053  if( excess[k] > 0 )
1054  {
1055  /* add to active list */
1056  activeinsert(k, dist[k], next, headactive, actminloc, actmaxloc, glbmaxloc);
1057  }
1058  else
1059  {
1060  if( dist[k] > (glbmaxloc) )
1061  glbmaxloc = dist[k];
1062 
1063  /* add to inactive list */
1064  inactiveinsert(k, dist[k], next, prev, headinactive, nextnode);
1065  }
1066  }
1067 #endif
1068  *dormmax = dormmaxloc;
1069  *actmin = actminloc;
1070  *actmax = actmaxloc;
1071  *glbmax = glbmaxloc;
1072 
1073 #if DEBUG > 0
1074  assert(is_valid(p, s, t, capa, w));
1075 #endif
1076 }
1077 
1078 
1079 
1080 /** sets default values*/
1081 void graph_mincut_setDefaultVals(
1082  GRAPH* g /**< graph data structure */
1083  )
1084 {
1085  const int nnodes = g->mincut_nnodes;
1086  int* RESTRICT excess = g->mincut_e;
1087  int* RESTRICT headactive = g->mincut_head;
1088  int* RESTRICT headinactive = g->mincut_head_inact;
1089 
1090  assert(g && excess && headactive && headinactive);
1091  assert(g->mincut_nnodes >= g->knots);
1092 
1093  for( int k = 0; k < nnodes; k++ )
1094  headactive[k] = Q_NULL;
1095 
1096  for( int k = 0; k < nnodes; k++ )
1097  headinactive[k] = Q_NULL;
1098 
1099  for( int k = 0; k < nnodes; k++ )
1100  excess[k] = 0;
1101 }
1102 
1103 
1104 /** initialize min cut arrays */
1105 SCIP_RETCODE graph_mincut_init(
1106  SCIP* scip, /**< SCIP data structure */
1107  GRAPH* p /**< graph data structure */
1108  )
1109 {
1110  assert(p != NULL);
1111  assert(!graph_mincut_isInitialized(p));
1112 
1113  SCIP_CALL( mincutInit(scip, p->knots, p->edges, p) );
1114 
1115  return SCIP_OKAY;
1116 }
1117 
1118 
1119 /** reinitializes minimum cut arrays */
1120 SCIP_RETCODE graph_mincut_reInit(
1121  SCIP* scip, /**< SCIP data structure */
1122  int nnodes, /**< number of nodes */
1123  int nedges, /**< number of edges */
1124  GRAPH* p /**< graph data structure */
1125  )
1126 {
1127  if( graph_mincut_isInitialized(p) )
1128  {
1129  graph_mincut_exit(scip, p);
1130  }
1131 
1132  SCIP_CALL( mincutInit(scip, nnodes, nedges, p) );
1133 
1134  return SCIP_OKAY;
1135 }
1136 
1137 
1138 /** is the minimum data initialized? */
1139 SCIP_Bool graph_mincut_isInitialized(
1140  const GRAPH* p /**< graph data structure */
1141  )
1142 {
1143  assert(p);
1144 
1145  if( p->mincut_dist )
1146  {
1147  assert(p->mincut_dist != NULL);
1148  assert(p->mincut_head != NULL);
1149  assert(p->mincut_numb != NULL);
1150  assert(p->mincut_prev != NULL);
1151  assert(p->mincut_next != NULL);
1152  assert(p->mincut_temp != NULL);
1153  assert(p->mincut_e != NULL);
1154  assert(p->mincut_nnodes > 0);
1155 
1156  return TRUE;
1157  }
1158 
1159  assert(p->mincut_nnodes == 0);
1160  assert(p->mincut_dist == NULL);
1161  assert(p->mincut_head == NULL);
1162  assert(p->mincut_numb == NULL);
1163  assert(p->mincut_prev == NULL);
1164  assert(p->mincut_next == NULL);
1165  assert(p->mincut_temp == NULL);
1166  assert(p->mincut_e == NULL);
1167  assert(p->mincut_x == NULL);
1168  assert(p->mincut_r == NULL);
1169 
1170  return FALSE;
1171 }
1172 
1173 
1174 /** frees min cut arrays */
1175 void graph_mincut_exit(
1176  SCIP* scip, /**< SCIP data structure */
1177  GRAPH* p /**< graph data structure */
1178  )
1179 {
1180  assert(scip && p);
1181  assert(graph_mincut_isInitialized(p));
1182  assert(p->mincut_nnodes > 0);
1183 
1184  p->mincut_nnodes = 0;
1185  p->mincut_nedges = 0;
1186  SCIPfreeMemoryArray(scip, &(p->mincut_r));
1187  SCIPfreeMemoryArray(scip, &(p->mincut_e));
1188  SCIPfreeMemoryArray(scip, &(p->mincut_temp));
1189  SCIPfreeMemoryArray(scip, &(p->mincut_next));
1190  SCIPfreeMemoryArray(scip, &(p->mincut_prev));
1191  SCIPfreeMemoryArray(scip, &(p->mincut_numb));
1192  SCIPfreeMemoryArray(scip, &(p->mincut_head_inact));
1193  SCIPfreeMemoryArray(scip, &(p->mincut_head));
1194  SCIPfreeMemoryArray(scip, &(p->mincut_dist));
1195 
1196 #if STATIST
1197  cut_statist();
1198 #endif
1199 }
1200 
1201 /** finds a minimum s-t cut */
1202 void graph_mincut_exec(
1203  const GRAPH* p,
1204  const int s,
1205  const int t,
1206  const int nnodesreal,
1207  const int nedgesreal,
1208  const int rootcutsize,
1209  const int* rootcut,
1210  const int* capa,
1211  int* RESTRICT w,
1212  const int* edgestart,
1213  const int* edgearr,
1214  const int* headarr,
1215  const SCIP_Bool isRerun
1216  )
1217 {
1218  int l;
1219  int actmax;
1220  int actmin;
1221  int glbmax;
1222  int dminus1;
1223  int dormmax;
1224  int mindist;
1225  int minnode;
1226  int minedgestart;
1227  int* RESTRICT e;
1228  int* RESTRICT r;
1229  int* RESTRICT dist;
1230  int* RESTRICT prev;
1231  int* RESTRICT next;
1232  int* RESTRICT temp;
1233  int* RESTRICT edgecurr;
1234  int* RESTRICT headactive;
1235  int* RESTRICT headinactive;
1236  int j;
1237  int end;
1238  int nnodes;
1239  int headnode;
1240  int prevnode;
1241  int nextnode;
1242  int maxresval;
1243  int maxresedge;
1244  int maxresnode;
1245  int relabeltrigger;
1246  int relabelupdatebnd;
1247 
1248  assert(p != NULL);
1249  assert(s >= 0);
1250  assert(s < p->knots);
1251  assert(t >= 0);
1252  assert(t < p->knots);
1253  assert(s != t);
1254  assert(w != NULL);
1255  assert(p->mincut_dist != NULL);
1256  assert(p->mincut_numb != NULL);
1257  assert(p->mincut_head != NULL);
1258  assert(p->mincut_prev != NULL);
1259  assert(p->mincut_next != NULL);
1260  assert(p->mincut_temp != NULL);
1261  assert(p->mincut_e != NULL);
1262  assert(p->mincut_r != NULL);
1263 
1264  dist = p->mincut_dist;
1265  headactive = p->mincut_head;
1266  headinactive = p->mincut_head_inact;
1267  edgecurr = p->mincut_numb;
1268  prev = p->mincut_prev;
1269  next = p->mincut_next;
1270  temp = p->mincut_temp;
1271  e = p->mincut_e;
1272  r = p->mincut_r;
1273 
1274  nnodes = p->knots;
1275  minnode = -1;
1276  maxresedge = -1;
1277 
1278  relabelupdatebnd = (GLOBALRELABEL_MULT * nnodesreal) + nedgesreal;
1279  relabeltrigger = 0;
1280 
1281  if( !isRerun )
1282  {
1283  initialise(s, t, nnodesreal, rootcutsize, rootcut, capa, dist, headactive, headinactive, edgecurr, next, prev, temp, e, r, w, edgestart,
1284  edgearr, headarr, &dormmax, &actmin, &actmax, &glbmax);
1285  }
1286  else
1287  {
1288  reinitialise(s, t, nnodesreal, rootcutsize, rootcut, capa, dist, headactive, headinactive, edgecurr, next, prev, temp, e, r, w, edgestart,
1289  edgearr, headarr, &dormmax, &actmin, &actmax, &glbmax);
1290  }
1291 
1292  /* main loop: get highest label node */
1293  while( actmax >= actmin )
1294  {
1295  const int i = headactive[actmax];
1296 
1297  /* no active vertices with distance label == actmax? */
1298  if( i < 0 )
1299  {
1300  actmax--;
1301  continue;
1302  }
1303 
1304 
1305  assert(actmax <= glbmax);
1306  assert(dist[i] == actmax);
1307  assert(w[i] == 0);
1308  assert(i != t);
1309  assert(i != s);
1310  assert(e[i] > 0);
1311 
1312  /* delete i from active list */
1313  activedelete(i,actmax,next,headactive);
1314 
1315  end = edgestart[i + 1];
1316 
1317  assert(end > edgestart[i]);
1318 
1319  /* try to discharge i (repeatedly) */
1320  for( ;; )
1321  {
1322  /* iterate over outgoing arcs of i */
1323  for( j = edgecurr[i]; j != end; j++ )
1324  {
1325  /* non-saturated edge? */
1326  if( r[j] != 0 )
1327  {
1328  /* non-dormant head node? */
1329  if( w[headarr[j]] == 0 )
1330  {
1331  headnode = headarr[j];
1332 
1333  assert(w[headnode] == 0);
1334  assert(headnode != s);
1335 
1336  dminus1 = dist[i] - 1;
1337 
1338  /* admissible arc? */
1339  if( dist[headnode] == dminus1 )
1340  {
1341  assert(Min(e[i], r[j]) > 0);
1342 
1343  /* is head node now active? */
1344  if( e[headnode] == 0 )
1345  {
1346  if( headnode != t )
1347  {
1348  /* remove head node from inactive list */
1349  inactivedelete(headnode, dminus1, next, prev, headinactive, nextnode, prevnode);
1350 
1351  /* add head node to active list */
1352  activeinsert(headnode, dminus1, next, headactive, actmin, actmax, glbmax);
1353  }
1354  }
1355 
1356  /* not more residual capacity than excess? */
1357  if( r[j] < e[i] )
1358  {
1359  e[i] -= r[j];
1360  e[headnode] += r[j];
1361  r[edgearr[j]] += r[j];
1362  r[j] = 0; /* -= r[a] */
1363 
1364  assert(e[i] > 0);
1365  }
1366  /* more residual capacity than excess */
1367  else
1368  {
1369  r[edgearr[j]] += e[i];
1370  r[j] -= e[i];
1371  e[headnode] += e[i];
1372  e[i] = 0; /* -= e[i] */
1373 
1374  /* excess vanished, so stop discharging */
1375  break;
1376  }
1377  } /* admissible arc */
1378  } /* non-dormant head node */
1379  } /* non-saturated edge */
1380  } /* all outgoing arcs */
1381 
1382  /* is there still excess on i? */
1383  if( j == end )
1384  {
1385  assert(e[i] > 0);
1386 
1387  /*
1388  * relabel i
1389  * */
1390 
1391  /* i only node of distance dist[i]? */
1392  if( (headactive[dist[i]] < 0) && (headinactive[dist[i]] < 0) )
1393  {
1394  /* put i into new dormant set */
1395  w[i] = ++dormmax;
1396 
1397  assert(dormmax <= nnodes);
1398  assert(dist[i] < nnodes);
1399 
1400  /* remove nodes of distance label >= dist[i] */
1401  for( j = dist[i] + 1; j <= glbmax; j++ )
1402  {
1403  assert(headactive[j] == Q_NULL);
1404 
1405  for( l = headinactive[j]; l >= 0; l = next[l] )
1406  {
1407  if( w[l] == 0 )
1408  w[l] = dormmax;
1409  }
1410  headinactive[j] = Q_NULL;
1411  }
1412 
1413  actmax = dist[i] - 1;
1414  glbmax = actmax;
1415 
1416  break;
1417  }
1418 
1419  j = edgestart[i];
1420 
1421  assert(end > j);
1422  assert(end == edgestart[i + 1]);
1423 
1424  relabeltrigger += GLOBALRELABEL_ADD + (end - j);
1425 
1426  mindist = nnodes;
1427  maxresval = 0;
1428 
1429  /* todo only for debugging, could be deleted */
1430  maxresnode = -1;
1431  minedgestart = -1;
1432 
1433  /* find the first (!) minimum */
1434  for( ; j != end; j++ )
1435  {
1436  /* useable edge? */
1437  if( r[j] != 0 )
1438  {
1439  /* non-dormant node? */
1440  if( w[headarr[j]] == 0 )
1441  {
1442 #if 0
1443  if( dist[headarr[j]] < mindist )
1444  {
1445  mindist = dist[headarr[j]];
1446  minedgestart = j;
1447  }
1448 #endif
1449  if( dist[headarr[j]] <= mindist )
1450  {
1451  if( dist[headarr[j]] < mindist )
1452  {
1453  mindist = dist[headarr[j]];
1454  minedgestart = j;
1455  maxresval = r[j];
1456  maxresedge = j;
1457  minnode = headarr[j];
1458  maxresnode = headarr[j];
1459  }
1460  else if( r[j] > maxresval )
1461  {
1462  maxresval = r[j];
1463  maxresedge = j;
1464  maxresnode = headarr[j];
1465  }
1466  }
1467  }
1468  }
1469  }
1470 
1471  if( (++mindist) < nnodes )
1472  {
1473  assert(minedgestart >= 0);
1474  assert(r[minedgestart] > 0);
1475 
1476  dist[i] = mindist;
1477  //edgecurr[i] = minedgestart; // iff: #if 1
1478 
1479  if( glbmax < mindist )
1480  glbmax = mindist;
1481 
1482  assert(maxresnode >= 0);
1483  assert(maxresnode == headarr[maxresedge]);
1484 #if 1
1485  /* can we completely discharge i already? */
1486  if( r[minedgestart] >= e[i] )
1487  {
1488  /* is node now active? */
1489  if( e[minnode] == 0 )
1490  {
1491  if( minnode != t )
1492  {
1493  inactivedelete(minnode, dist[minnode], next, prev, headinactive, nextnode, prevnode);
1494  activeinsert(minnode, dist[minnode], next, headactive, actmin, actmax, glbmax);
1495  }
1496  }
1497 
1498  r[edgearr[minedgestart]] += e[i];
1499  r[minedgestart] -= e[i];
1500  e[minnode] += e[i];
1501  e[i] = 0; /* -= e[i] */
1502 
1503  inactiveinsert(i, mindist, next, prev, headinactive, nextnode);
1504 
1505  edgecurr[i] = minedgestart;
1506 
1507  /* excess vanished, so stop discharging */
1508  break;
1509  }
1510 
1511  /* is node now active? */
1512  if( e[minnode] == 0 )
1513  {
1514  if( minnode != t )
1515  {
1516  inactivedelete(minnode, dist[minnode], next, prev, headinactive, nextnode, prevnode);
1517  activeinsert(minnode, dist[minnode], next, headactive, actmin, actmax, glbmax);
1518  }
1519  }
1520 
1521  e[i] -= r[minedgestart];
1522  e[minnode] += r[minedgestart];
1523  r[edgearr[minedgestart]] += r[minedgestart];
1524  r[minedgestart] = 0; /* -= r[minedgestart] */
1525 
1526  assert(e[i] > 0);
1527 
1528  if( maxresval >= e[i] && minedgestart != maxresedge ) // todo
1529  {
1530  /* is node now active? */
1531  if( e[maxresnode] == 0 )
1532  {
1533  if( maxresnode != t )
1534  {
1535  /* remove head node from inactive list */
1536  inactivedelete(maxresnode, dist[maxresnode], next, prev, headinactive, nextnode, prevnode);
1537 
1538  /* add head node to active list */
1539  activeinsert(maxresnode, dist[maxresnode], next, headactive, actmin, actmax, glbmax);
1540  }
1541  }
1542 
1543  r[edgearr[maxresedge]] += e[i];
1544  r[maxresedge] -= e[i];
1545  e[maxresnode] += e[i];
1546  e[i] = 0; /* -= e[i] */
1547 
1548  inactiveinsert(i, mindist, next, prev, headinactive, nextnode);
1549 
1550  edgecurr[i] = minedgestart;
1551 
1552  /* excess vanished, so stop discharging */
1553  break;
1554  }
1555  edgecurr[i] = minedgestart + 1;
1556 #endif
1557 
1558  continue;
1559  }
1560  /* could not relabel i */
1561  else
1562  {
1563  /* put i into new dormant set */
1564  w[i] = ++dormmax;
1565 
1566  assert(dormmax <= nnodes);
1567 
1568  assert(!((headactive[dist[i]] == Q_NULL) && (headinactive[dist[i]] == Q_NULL)));
1569 
1570  if( relabeltrigger > relabelupdatebnd )
1571  {
1572  if( actmax < actmin )
1573  break;
1574 
1575  /* execute global relabel heuristic */
1576  globalrelabel(s, t, nnodesreal, dist, headactive, headinactive, edgecurr, next, prev,
1577  e, r, w, edgestart, edgearr, headarr, &actmin, &actmax, &glbmax, &dormmax);
1578 
1579  relabeltrigger = 0;
1580  }
1581 
1582  break;
1583  }
1584  }
1585  /* no excess on i */
1586  else
1587  {
1588  assert(e[i] == 0);
1589 
1590  edgecurr[i] = j;
1591 
1592  inactiveinsert(i, dist[i], next, prev, headinactive, nextnode);
1593 
1594  break;
1595  }
1596  }
1597 
1598  if( relabeltrigger > relabelupdatebnd )
1599  {
1600  if( actmax < actmin )
1601  break;
1602 
1603  /* execute global relabel heuristic */
1604  globalrelabel(s, t, nnodesreal, dist, headactive, headinactive, edgecurr, next, prev,
1605  e, r, w, edgestart, edgearr, headarr, &actmin, &actmax, &glbmax, &dormmax);
1606 
1607  relabeltrigger = 0;
1608  }
1609  }
1610  assert(w[s]);
1611  assert(!w[t]);
1612 
1613 #if CHECK
1614  assert(capa);
1615  if( !is_valid_arr(p, s, t, edgestart, headarr, edgearr, r, capa, w) )
1616  printf("flow is not valid \n");
1617  assert(is_valid_arr(p, s, t, edgestart, headarr, edgearr, r, capa, w));
1618 
1619 #endif
1620 
1621 
1622 #if STATIST
1623  cut_sum(p, capa, w);
1624 #endif
1625 
1626 #if DEBUG > 0
1627  assert(is_valid(p, s, t, capa, w));
1628 #endif
1629 }
1630 
1631 #if STATIST
1632 static void cut_statist(void)
1633 {
1634  (void)printf("Mincut Statistics:\n");
1635  (void)printf("Node-Searches=%d, Cut Sums=%d\n",
1636  searches, cutsums);
1637  (void)printf("S-Pushes=%d, N-Pushes=%d, X-Pushes=%d, M-Pushes=%d\n",
1638  s_pushes, n_pushes, x_pushes, m_pushes);
1639  (void)printf("Relabels=%d, S-Sleeps=%d, M-Sleeps=%d\n\n",
1640  relabels, s_sleeps, m_sleeps);
1641 
1642  s_pushes = 0;
1643  n_pushes = 0;
1644  x_pushes = 0;
1645  m_pushes = 0;
1646  relabels = 0;
1647  s_sleeps = 0;
1648  m_sleeps = 0;
1649  searches = 0;
1650  cutsums = 0;
1651 }
1652 
1653 static void cut_sum(
1654  const GRAPH* p,
1655  const int* capa,
1656  const int* w)
1657 {
1658  int sum = 0;
1659  int i;
1660  int j;
1661  int k;
1662 
1663  assert(p != NULL);
1664  assert(capa != NULL);
1665  assert(w != NULL);
1666 
1667  for(k = 0; k < p->edges; k++)
1668  {
1669  i = p->head[k];
1670  j = p->tail[k];
1671 
1672  if ((w[i] && !w[j]) || (!w[i] && w[j]))
1673  sum += capa[k];
1674  }
1675 #if DEBUG > 0
1676  (void)printf("Cut Sum=%d\n", sum);
1677 #endif
1678  cutsums += sum;
1679 }
1680 #endif
1681 
1682 
1683 
1684 
1685 #if DEBUG > 0
1686 static int is_valid(
1687  const GRAPH* p,
1688  const int s,
1689  const int t,
1690  const int* capa,
1691  const int* w)
1692 {
1693  int* e;
1694  int* r;
1695  int* dist;
1696 #if 0
1697  int* x;
1698 #endif
1699  int j;
1700  int k;
1701 
1702  assert(p != NULL);
1703  assert(p->mincut_e != NULL);
1704  assert(p->mincut_r != NULL);
1705  assert(p->mincut_x != NULL);
1706 
1707  e = p->mincut_e;
1708  r = p->mincut_r;
1709  dist = p->mincut_dist;
1710 #if 0
1711  x = p->mincut_x;
1712 #endif
1713 
1714 #if 0
1715  for (j = 0; j < p->knots; j++)
1716  {
1717  for (k = p->outbeg[j]; k != EAT_LAST; k = p->oeat[k])
1718  {
1719  if( r[k] > 0 )
1720  {
1721  if( dist[j] > dist[p->head[k]] + 1 )
1722  {
1723  printf("distance fail! %d->%d (%d)\n", j, p->head[k], k);
1724  return FALSE;
1725  }
1726  }
1727  }
1728  }
1729 #endif
1730 
1731  for(j = 0; j < p->knots; j++)
1732  {
1733 #if 0
1734  if ((q[j] >= 0) && (a[q[j]] != j))
1735  return((void)fprintf(stderr, "Queue Error 1 Knot %d\n", j), FALSE);
1736 
1737  if (!w[j] && (q[j] < 0) && (e[j] > 0) && (j != t))
1738  return((void)fprintf(stderr, "Queue Error 2 Knot %d\n", j), FALSE);
1739 
1740  if (!w[j] && (q[j] >= 0) && (e[j] == 0))
1741  return((void)fprintf(stderr, "Queue Error 3 Knot %d\n", j), FALSE);
1742 
1743  if (w[j] && (q[j] >= 0))
1744  return((void)fprintf(stderr, "Queue Error 4 Knot %d\n", j), FALSE);
1745 #endif
1746  if (e[j] < 0)
1747  return((void)fprintf(stderr, "Negativ Execess Knot %d\n", j), FALSE);
1748 #if 0
1749  if (p->mincut_dist[j] >= p->knots)
1750  return((void)fprintf(stderr, "Distance too big Knot %d\n", j), FALSE);
1751 #endif
1752  /* Extended Dormacy Property
1753  */
1754  for(k = p->outbeg[j]; k != EAT_LAST; k = p->oeat[k])
1755  {
1756  if (r[k] > 0)
1757  {
1758 #if 1
1759  if( dist[j] > dist[p->head[k]] + 1 && w[j] == 0 && w[p->head[k]] == 0 )
1760  {
1761  printf("distance fail! %d->%d (%d)\n", j, p->head[k], k);
1762  return FALSE;
1763  }
1764 #endif
1765 
1766  if((w[j] && (w[j] < w[p->head[k]])))
1767  {
1768  printf("fail %d %d\n", j, p->head[k]);
1769  }
1770 
1771  if( w[j] && !w[p->head[k]] )
1772  {
1773  printf("fail2 %d -> %d\n", j, p->head[k]);
1774  printf("w: %d -> w: %d \n", w[j], w[p->head[k]]);
1775  printf("edge %d \n", k);
1776  }
1777 
1778  if ((w[j] && !w[p->head[k]]) || (w[j] && (w[j] < w[p->head[k]])))
1779  {
1780  (void)printf("k=%d r[k]=%d head=%d tail=%d w[h]=%d w[t]=%d\n",
1781  k, r[k], p->head[k], p->tail[k], w[p->head[k]], w[p->tail[k]]);
1782 
1783  return((void)fprintf(stderr, "Extended Dormacy Violation Knot %d\n", j), FALSE);
1784  }
1785  }
1786  }
1787  }
1788  for(j = 0; j < p->edges; j++)
1789  {
1790  if (r[j] < 0)
1791  return((void)fprintf(stderr, "Negativ Residual Edge %d\n", j), FALSE);
1792 
1793  if (r[j] + r[Edge_anti(j)] != capa[j] + capa[Edge_anti(j)] && w[j] == 0 && w[p->head[k]] == 0 )
1794  return((void)fprintf(stderr, "Wrong Capacity Equation Edge %d\n", j), FALSE);
1795 #if 0
1796  if (x[j] < 0)
1797  return((void)fprintf(stderr, "Negativ Flow Edge %d\n", j), FALSE);
1798 
1799  if (r[j] != capa[j] - x[j] + x[Edge_anti(j)])
1800  return((void)fprintf(stderr, "Wrong Residual Equation Edge %d\n", j), FALSE);
1801 #endif
1802  }
1803  return(TRUE);
1804 }
1805 
1806 static void show_flow(
1807  const GRAPH* p,
1808  const int* capa,
1809  const int* w)
1810 {
1811  int j;
1812  int* head;
1813  int* numb;
1814  int* prev;
1815  int* next;
1816  int* e;
1817  int* x;
1818  int* r;
1819 
1820  assert(p != NULL);
1821  assert(w != NULL);
1822  assert(p->mincut_numb != NULL);
1823  assert(p->mincut_head != NULL);
1824  assert(p->mincut_prev != NULL);
1825  assert(p->mincut_next != NULL);
1826  assert(p->mincut_e != NULL);
1827  assert(p->mincut_x != NULL);
1828  assert(p->mincut_r != NULL);
1829 
1830  head = p->mincut_head;
1831  numb = p->mincut_numb;
1832  prev = p->mincut_prev;
1833  next = p->mincut_next;
1834  e = p->mincut_e;
1835  x = p->mincut_x;
1836  r = p->mincut_r;
1837 
1838 
1839 
1840  (void)printf(" ");
1841  for(j = 0; j < p->edges; j++)
1842  (void)printf("%6d ", j);
1843  (void)fputc('\n', stdout);
1844 
1845  (void)printf("ta:");
1846  for(j = 0; j < p->edges; j++)
1847  (void)printf("%6d ", p->tail[j]);
1848  (void)fputc('\n', stdout);
1849 
1850  (void)printf("he:");
1851  for(j = 0; j < p->edges; j++)
1852  (void)printf("%6d ", p->head[j]);
1853  (void)fputc('\n', stdout);
1854 
1855  (void)printf("x: ");
1856  for(j = 0; j < p->edges; j++)
1857  (void)printf("%6d ", x[j]);
1858  (void)fputc('\n', stdout);
1859 
1860  (void)printf("r: ");
1861  for(j = 0; j < p->edges; j++)
1862  (void)printf("%6d ", r[j]);
1863  (void)fputc('\n', stdout);
1864 
1865  (void)printf("ca:");
1866  for(j = 0; j < p->edges; j++)
1867  (void)printf("%6d ", capa[j]);
1868  (void)fputc('\n', stdout);
1869 
1870  (void)printf("w: ");
1871  for(j = 0; j < p->knots; j++)
1872  (void)printf("%2d ", w[j]);
1873  (void)fputc('\n', stdout);
1874 
1875  (void)printf("d: ");
1876  for(j = 0; j < p->knots; j++)
1877  (void)printf("%2d ", p->mincut_dist[j]);
1878  (void)fputc('\n', stdout);
1879 
1880  (void)printf("n: ");
1881  for(j = 0; j < p->knots; j++)
1882  (void)printf("%2d ", numb[j]);
1883  (void)fputc('\n', stdout);
1884 
1885  (void)printf("h: ");
1886  for(j = 0; j < p->knots; j++)
1887  (void)printf("%2d ", head[j]);
1888  (void)fputc('\n', stdout);
1889 
1890  (void)printf("p: ");
1891  for(j = 0; j < p->knots; j++)
1892  (void)printf("%2d ", prev[j]);
1893  (void)fputc('\n', stdout);
1894 
1895  (void)printf("n: ");
1896  for(j = 0; j < p->knots; j++)
1897  (void)printf("%2d ", next[j]);
1898  (void)fputc('\n', stdout);
1899 
1900  (void)printf("e: ");
1901  for(j = 0; j < p->knots; j++)
1902  (void)printf("%2d ", e[j]);
1903  (void)fputc('\n', stdout);
1904 }
1905 
1906 #endif /* DEBUG > 0 */
1907 
1908 
1909 
1910 
1911 
1912 
1913 #else
1914 /*---------------------------------------------------------------------------*/
1915 /*--- Name : GRAPH MINimumCUT INITialise ---*/
1916 /*--- Function : Holt den Speicher fuer die Hilfsarrays die wir brauchen. ---*/
1917 /*--- Parameter: Graph ---*/
1918 /*--- Returns : Nichts ---*/
1919 /*---------------------------------------------------------------------------*/
1920 SCIP_RETCODE graph_mincut_init(
1921  SCIP* scip, /**< SCIP data structure */
1922  GRAPH* p /**< graph data structure */
1923  )
1924 {
1925  assert(p != NULL);
1926  assert(p->mincut_dist == NULL);
1927  assert(p->mincut_head == NULL);
1928  assert(p->mincut_numb == NULL);
1929  assert(p->mincut_prev == NULL);
1930  assert(p->mincut_next == NULL);
1931  assert(p->mincut_temp == NULL);
1932  assert(p->mincut_e == NULL);
1933  assert(p->mincut_x == NULL);
1934  assert(p->mincut_r == NULL);
1935 
1936 #if BLK
1937  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(p->mincut_dist), p->knots) );
1938  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(p->mincut_head), p->knots) );
1939  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(p->mincut_numb), p->knots) );
1940  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(p->mincut_prev), p->knots) );
1941  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(p->mincut_next), p->knots) );
1942  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(p->mincut_temp), p->knots) );
1943  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(p->mincut_e), p->knots) );
1944  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(p->mincut_x), p->edges) );
1945  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(p->mincut_r), p->edges) );
1946 #else
1947  SCIP_CALL( SCIPallocMemoryArray(scip, &(p->mincut_dist), p->knots) );
1948  SCIP_CALL( SCIPallocMemoryArray(scip, &(p->mincut_head), p->knots) );
1949  SCIP_CALL( SCIPallocMemoryArray(scip, &(p->mincut_numb), p->knots) );
1950  SCIP_CALL( SCIPallocMemoryArray(scip, &(p->mincut_prev), p->knots) );
1951  SCIP_CALL( SCIPallocMemoryArray(scip, &(p->mincut_next), p->knots) );
1952  SCIP_CALL( SCIPallocMemoryArray(scip, &(p->mincut_temp), p->knots) );
1953  SCIP_CALL( SCIPallocMemoryArray(scip, &(p->mincut_e), p->knots) );
1954  SCIP_CALL( SCIPallocMemoryArray(scip, &(p->mincut_x), p->edges) );
1955  SCIP_CALL( SCIPallocMemoryArray(scip, &(p->mincut_r), p->edges) );
1956 #endif
1957 
1958  return SCIP_OKAY;
1959 }
1960 
1961 /*---------------------------------------------------------------------------*/
1962 /*--- Name : GRAPH MINimumCUT EXIT ---*/
1963 /*--- Function : Gibt den Speicher fuer die Hilfsarrays wieder frei. ---*/
1964 /*--- Parameter: Keine ---*/
1965 /*--- Returns : Nichts ---*/
1966 /*---------------------------------------------------------------------------*/
1967 void graph_mincut_exit(
1968  SCIP* scip, /**< SCIP data structure */
1969  GRAPH* p /**< graph data structure */
1970  )
1971 {
1972  assert(p->mincut_dist != NULL);
1973  assert(p->mincut_head != NULL);
1974  assert(p->mincut_numb != NULL);
1975  assert(p->mincut_prev != NULL);
1976  assert(p->mincut_next != NULL);
1977  assert(p->mincut_temp != NULL);
1978  assert(p->mincut_e != NULL);
1979  assert(p->mincut_x != NULL);
1980  assert(p->mincut_r != NULL);
1981 
1982 #if BLK
1983  SCIPfreeBlockMemoryArray(scip, &(p->mincut_r), p->edges);
1984  SCIPfreeBlockMemoryArray(scip, &(p->mincut_x), p->edges);
1985  SCIPfreeBlockMemoryArray(scip, &(p->mincut_e), p->knots);
1986  SCIPfreeBlockMemoryArray(scip, &(p->mincut_temp), p->knots);
1987  SCIPfreeBlockMemoryArray(scip, &(p->mincut_next), p->knots);
1988  SCIPfreeBlockMemoryArray(scip, &(p->mincut_prev), p->knots);
1989  SCIPfreeBlockMemoryArray(scip, &(p->mincut_numb), p->knots);
1990  SCIPfreeBlockMemoryArray(scip, &(p->mincut_head), p->knots);
1991  SCIPfreeBlockMemoryArray(scip, &(p->mincut_dist), p->knots);
1992 #else
1993  SCIPfreeMemoryArray(scip, &(p->mincut_r));
1994  SCIPfreeMemoryArray(scip, &(p->mincut_x));
1995  SCIPfreeMemoryArray(scip, &(p->mincut_e));
1996  SCIPfreeMemoryArray(scip, &(p->mincut_temp));
1997  SCIPfreeMemoryArray(scip, &(p->mincut_next));
1998  SCIPfreeMemoryArray(scip, &(p->mincut_prev));
1999  SCIPfreeMemoryArray(scip, &(p->mincut_numb));
2000  SCIPfreeMemoryArray(scip, &(p->mincut_head));
2001  SCIPfreeMemoryArray(scip, &(p->mincut_dist));
2002 #endif
2003 
2004 #if STATIST
2005  cut_statist();
2006 #endif
2007 }
2008 
2009 #if 0
2010 inline static void delete(
2011  const int knot,
2012  int* q_dist,
2013  int* q_head,
2014  int* q_prev,
2015  int* q_next)
2016 {
2017  assert(knot >= 0);
2018  assert(q_dist != NULL);
2019  assert(q_head != NULL);
2020  assert(q_prev != NULL);
2021  assert(q_next != NULL);
2022  assert(q_dist[knot] > 0);
2023 
2024  if (q_next[knot] != Q_NM)
2025  {
2026  /* Etwa Erster ?
2027  */
2028  if (q_prev[knot] == Q_NULL)
2029  {
2030  assert(q_dist[knot] >= 0);
2031  assert(q_head[q_dist[knot]] == knot);
2032 
2033  q_head[q_dist[knot]] = q_next[knot];
2034  }
2035  else
2036  {
2037  assert(q_prev[knot] >= 0);
2038  assert(q_next[q_prev[knot]] == knot);
2039 
2040  q_next[q_prev[knot]] = q_next[knot];
2041  }
2042 
2043  /* Sind wir auch nicht letzter ?
2044  */
2045  if (q_next[knot] != Q_NULL)
2046  {
2047  assert(q_next[knot] >= 0);
2048  assert(q_prev[q_next[knot]] == knot);
2049 
2050  q_prev[q_next[knot]] = q_prev[knot];
2051  }
2052  q_next[knot] = Q_NM;
2053  q_prev[knot] = Q_NM;
2054  }
2055  assert(q_next[knot] == Q_NM);
2056  assert(q_prev[knot] == Q_NM);
2057 }
2058 
2059 inline static int insert(
2060  const int knot,
2061  int* q_dist,
2062  int* q_head,
2063  int* q_prev,
2064  int* q_next,
2065  int dmin,
2066  int* dlmax)
2067 {
2068  assert(knot >= 0);
2069  assert(q_dist != NULL);
2070  assert(q_head != NULL);
2071  assert(q_prev != NULL);
2072  assert(q_next != NULL);
2073  assert(q_dist[knot] > 0);
2074  assert(dmin >= 1);
2075 
2076  if( q_prev[knot] == Q_NM )
2077  {
2078  q_prev[knot] = Q_NULL;
2079  q_next[knot] = q_head[q_dist[knot]];
2080 
2081  if( q_next[knot] != Q_NULL )
2082  q_prev[q_next[knot]] = knot;
2083 
2084  q_head[q_dist[knot]] = knot;
2085 
2086  if( q_dist[knot] < dmin )
2087  dmin = q_dist[knot];
2088  if( q_dist[knot] > (*dlmax) )
2089  *dlmax = q_dist[knot];
2090  }
2091  assert(q_next[knot] != Q_NM);
2092  assert(q_prev[knot] != Q_NM);
2093  assert(q_dist[knot] >= dmin);
2094 
2095  return(dmin);
2096 }
2097 #endif
2098 
2099 static int bfs(
2100  const GRAPH* p,
2101  const int s,
2102  const int t,
2103  int* q_dist,
2104  int* q_numb,
2105  int* q_temp,
2106  int* w,
2107  const int* edgestart,
2108  const int* edgearr,
2109  const int* headarr)
2110 {
2111  int i;
2112  int j;
2113  int k;
2114  int l;
2115  int end;
2116  int visited = 0;
2117 
2118  assert(q_temp != NULL);
2119  assert(q_dist != NULL);
2120  assert(q_numb != NULL);
2121  assert(w != NULL);
2122  assert(s >= 0);
2123  assert(s < p->knots);
2124  assert(t >= 0);
2125  assert(t < p->knots);
2126 
2127  /* Beginnen wir bei der Senke */
2128  q_temp[visited++] = t;
2129  q_dist[t] = 0;
2130 
2131  /* Solange noch schon besuchte Knoten da sind, von denen aus nicht
2132  * versucht wurde weiter zu kommen:
2133  */
2134  for(j = 0; (j < visited) && (visited < p->knots); j++)
2135  {
2136  assert(visited < p->knots);
2137  assert(j < visited);
2138 
2139  i = q_temp[j];
2140 
2141  assert(i >= 0);
2142  assert(i < p->knots);
2143  assert(q_dist[i] >= 0);
2144  assert(q_dist[i] < visited);
2145  assert(w[i] == 0);
2146 
2147  assert((j == 0) || (q_dist[i] >= q_dist[q_temp[j - 1]]));
2148  assert((j == visited - 1) || (q_dist[i] <= q_dist[q_temp[j + 1]]));
2149 
2150  /* Wo koennen wir den ueberall hin:
2151  */
2152  for( k = edgestart[i], end = edgestart[i + 1]; k != end; k++ )
2153  {
2154  l = headarr[k];
2155 
2156  /* Waren wir da noch nicht ?
2157  */
2158  assert(!w[l] || (q_dist[l] >= 0));
2159 
2160  if( q_dist[l] < 0 )
2161  {
2162  q_dist[l] = q_dist[i] + 1;
2163  q_temp[visited++] = l;
2164  q_numb[q_dist[l]]++;
2165 
2166  assert(q_dist[l] < p->knots);
2167  }
2168  }
2169  }
2170  return(visited);
2171 }
2172 
2173 
2174 /** global relabel heuristic that sets distance of sink to zero and relabels all other nodes using backward bfs on residual
2175  * graph, starting from the sink. */
2176 static void globalrelabel(
2177  const GRAPH* p,
2178  const int s,
2179  const int t,
2180  int* q_dist,
2181  int* q_numb,
2182  int* q_head,
2183  int* q_prev,
2184  int* q_next,
2185  int* q_temp,
2186  int* excess,
2187  int* transx,
2188  int* residual,
2189  int* w,
2190  int* dormmax,
2191  int* dlmax,
2192  int* dmin)
2193 {
2194  int i;
2195  int k;
2196  int l;
2197  int j;
2198  int hit;
2199  int dist1;
2200  int nnodes;
2201  int visited;
2202 
2203  assert(p != NULL);
2204  assert(s >= 0);
2205  assert(s < p->knots);
2206  assert(t >= 0);
2207  assert(t < p->knots);
2208  assert(s != t);
2209  assert(w != NULL);
2210  assert(q_head != NULL);
2211  assert(q_dist != NULL);
2212  assert(q_numb != NULL);
2213  assert(q_prev != NULL);
2214  assert(q_next != NULL);
2215  assert(q_temp != NULL);
2216  assert(excess != NULL);
2217  assert(transx != NULL);
2218  assert(residual != NULL);
2219 
2220  nnodes = p->knots;
2221 
2222  assert(w[s]);
2223 
2224  /* backwards bfs starting from sink */
2225 
2226  for( i = 0; i < nnodes; i++ )
2227  {
2228  q_next[i] = Q_NULL;
2229  q_head[i] = Q_NULL;
2230  q_numb[i] = 0;
2231 
2232  if( w[i] == 0 )
2233  w[i] = -1;
2234  }
2235 
2236  *dlmax = 0;
2237  (*dmin) = nnodes;
2238 
2239  visited = 0;
2240 
2241  q_temp[visited++] = t;
2242  q_dist[t] = 0;
2243  w[t] = 0;
2244 
2245  for( j = 0; (j < visited) && (visited < nnodes); j++ )
2246  {
2247  assert(visited < nnodes);
2248  assert(j < visited);
2249 
2250  i = q_temp[j];
2251 
2252  assert(i >= 0);
2253  assert(i < nnodes);
2254  assert(q_dist[i] < nnodes);
2255  assert(w[i] == 0);
2256 
2257  dist1 = q_dist[i] + 1;
2258 
2259  assert(dist1 < nnodes);
2260 
2261  /* check all neighbors */
2262  for( k = p->outbeg[i]; k != EAT_LAST; k = p->oeat[k] )
2263  {
2264  l = p->head[k];
2265 
2266  /* is not node awake, allows for positive flow along edge, and has not been visited so far? */
2267  if( (w[l] < 0) && (residual[flipedge(k)] > 0) )
2268  {
2269  assert(l != s);
2270  w[l] = 0;
2271  q_temp[visited++] = l;
2272 
2273  q_dist[l] = dist1;
2274  q_numb[q_dist[l]]++;
2275 
2276  if( excess[l] > 0 )
2277  {
2278  /* renew the distance label */
2279  listinsert(l, dist1, q_next, q_head, (*dmin), (*dlmax));
2280 
2281  assert(q_dist[l] < nnodes);
2282  }
2283  }
2284  }
2285  }
2286 
2287  hit = FALSE;
2288  for( i = 0; i < nnodes; i++ )
2289  if( w[i] < 0 )
2290  {
2291  hit = TRUE;
2292  w[i] = *dormmax + 1;
2293  }
2294 
2295  if( hit )
2296  *dormmax = (*dormmax) + 1;
2297 
2298  return;
2299 }
2300 
2301 
2302 /*---------------------------------------------------------------------------*/
2303 /*--- Name : INITialise LABELS ---*/
2304 /*--- Function : Fuehrt eine BFS durch um die Distanzen der Knoten von ---*/
2305 /*--- der Senke zu ermitten. Kanten ohne Kapazitaet werden ---*/
2306 /*--- nicht beruecksichtigt, ebenso Knoten die nur ueber die ---*/
2307 /*--- Quelle zu erreichen sind. ---*/
2308 /*--- Parameter: Graph, Quelle, Senke, Kantenkapazitaeten. ---*/
2309 /*--- Returns : Anzahl der Aktiven (Erreichbaren) Knoten, ---*/
2310 /*--- Fuellt a[] mit den Nummern dieser Knoten und d[] mit den ---*/
2311 /*--- Entfernungen der zugehoerigen Knoten zur Senke. ---*/
2312 /*--- a[] ist aufsteigend sortiert. ---*/
2313 /*---------------------------------------------------------------------------*/
2314 static void initialise(
2315  const GRAPH* p,
2316  const int s,
2317  const int t,
2318  const int* capa,
2319  int* q_dist,
2320  int* q_numb,
2321  int* q_head,
2322  int* q_prev,
2323  int* q_next,
2324  int* q_temp,
2325  int* excess,
2326  int* transx,
2327  int* residual,
2328  int* w,
2329  const int* edgestart,
2330  const int* edgearr,
2331  const int* headarr,
2332  int* dormmax,
2333  int* dlmax)
2334 {
2335 
2336  int i;
2337  int j;
2338  int k;
2339  int q;
2340  int end;
2341  int actminloc;
2342 
2343  assert(p != NULL);
2344  assert(s >= 0);
2345  assert(s < p->knots);
2346  assert(t >= 0);
2347  assert(t < p->knots);
2348  assert(s != t);
2349  assert(capa != NULL);
2350  assert(w != NULL);
2351  assert(q_head != NULL);
2352  assert(q_dist != NULL);
2353  assert(q_numb != NULL);
2354  assert(q_prev != NULL);
2355  assert(q_next != NULL);
2356  assert(q_temp != NULL);
2357  assert(excess != NULL);
2358  assert(transx != NULL);
2359  assert(residual != NULL);
2360  assert(p->mincut_r != NULL);
2361  assert(p->mincut_x != NULL);
2362 
2363  /* Knotenarrays initialisieren
2364  */
2365  *dormmax = 1;
2366  *dlmax = -1;
2367  actminloc = p->knots;
2368 
2369  for(i = 0; i < p->knots; i++)
2370  {
2371  excess[i] = 0;
2372  w [i] = 0;
2373  q_next[i] = Q_NULL;
2374  q_head[i] = Q_NULL;
2375  q_numb[i] = 0;
2376  q_dist[i] = -1;
2377  }
2378  /* Jetzt die Kantenarrays.
2379  */
2380  for(k = 0; k < p->edges; k++)
2381  {
2382  transx[k] = 0;
2383  residual[k] = capa[k];
2384  }
2385  /* Jetzt noch dist und numb.
2386  */
2387  (void)bfs(p, s, t, q_dist, q_numb, q_temp, w, edgestart, edgearr, headarr);
2388 
2389  /* Alles was wir nicht erreichen konnten schlafen legen.
2390  */
2391  for(i = 0; i < p->knots; i++)
2392  if (q_dist[i] < 0)
2393  w[i] = *dormmax + 1;
2394 
2395  /* put sink into lowest dormant set */
2396  w[s] = 1; /* dormmax */
2397 
2398  /* Falls wir die Quelle s nicht erreichen konnten sind wir fertig.
2399  */
2400  if (q_dist[s] < 0)
2401  return;
2402 
2403  assert(w[s] > 0);
2404  assert(q_dist[s] > 0);
2405  assert(w[t] == 0);
2406  assert(q_dist[t] == 0);
2407 
2408  /* Label der Quelle abziehen
2409  */
2410  q_numb[q_dist[s]]--;
2411 
2412  /* Von der Quelle alles wegschieben wofuer Kapazitaeten da sind.
2413  */
2414  for( q = edgestart[s], end = edgestart[s + 1]; q != end; q++ )
2415  {
2416  k = edgearr[q];
2417 
2418  if (capa[k] == 0)
2419  continue;
2420 
2421  j = headarr[q];
2422 
2423  transx[k] = capa[k];
2424  residual[k] = 0; /* -= transx[k] */
2425 
2426  residual[Edge_anti(k)] += transx[k]; /* Ueberfluessig weil w[s] == 1 */
2427  excess[j] += transx[k];
2428 
2429  if (j != t)
2430  {
2431  listinsert(j, q_dist[j], q_next, q_head, actminloc, (*dlmax));
2432  }
2433 
2434  assert(w[j] == 0);
2435  assert(excess[j] > 0);
2436 
2437  assert((p->mincut_r)[k] + (p->mincut_r)[Edge_anti(k)] == capa[k] + capa[Edge_anti(k)]);
2438  assert((p->mincut_x)[k] >= 0);
2439  assert((p->mincut_x)[k] <= capa[k]);
2440  assert((p->mincut_r)[k] == capa[k] - (p->mincut_x)[k] + (p->mincut_x)[Edge_anti(k)]);
2441  }
2442 #if DEBUG > 1
2443  show_flow(p, capa, w);
2444 #endif
2445 #if DEBUG > 0
2446  assert(is_valid(p, s, t, capa, w));
2447 #endif
2448 }
2449 
2450 static void reinitialise(
2451  const GRAPH* p,
2452  const int s,
2453  const int t,
2454  const int* capa,
2455  int* q_dist,
2456  int* q_numb,
2457  int* q_head,
2458  int* q_prev,
2459  int* q_next,
2460  int* q_temp,
2461  int* excess,
2462  int* transx,
2463  int* residual,
2464  int* w,
2465  const int* edgestart,
2466  const int* edgearr,
2467  const int* headarr,
2468  int* dormmax,
2469  int* dlmax)
2470 {
2471  int wt;
2472  int i;
2473  int j;
2474  int k;
2475  int visited;
2476  int actminloc;
2477 
2478  assert(p != NULL);
2479  assert(s >= 0);
2480  assert(s < p->knots);
2481  assert(t >= 0);
2482  assert(t < p->knots);
2483  assert(s != t);
2484  assert(capa != NULL);
2485  assert(w != NULL);
2486  assert(q_head != NULL);
2487  assert(q_dist != NULL);
2488  assert(q_numb != NULL);
2489  assert(q_prev != NULL);
2490  assert(q_next != NULL);
2491  assert(q_temp != NULL);
2492  assert(excess != NULL);
2493  assert(transx != NULL);
2494  assert(residual != NULL);
2495 
2496  /* initialize */
2497  assert(w[s]);
2498 
2499  wt = (w[t] == 0) ? p->knots + 1 : w[t];
2500  actminloc = p->knots;
2501  *dormmax = 1;
2502  *dlmax = 1;
2503 
2504  for( i = 0; i < p->knots; i++ )
2505  {
2506  q_numb[i] = 0;
2507  q_head[i] = Q_NULL;
2508  q_next[i] = Q_NULL;
2509 
2510  /* wake up nodes in higher or equal dormant layer */
2511  if( (w[i] == 0) || (w[i] >= wt) )
2512  {
2513  w[i] = 0;
2514  q_dist[i] = -1;
2515  }
2516  else if( w[i] > *dormmax )
2517  *dormmax = w[i];
2518  }
2519  /* Jetzt noch dist und numb.
2520  */
2521  visited = bfs(p, s, t, q_dist, q_numb, q_temp, w, edgestart, edgearr, headarr);
2522 
2523  /* Alles was wir nicht erreichen konnten schlafen legen.
2524  */
2525  for( i = 0; i < p->knots; i++ )
2526  {
2527  if( q_dist[i] < 0 )
2528  w[i] = *dormmax + 1;
2529  else if( (w[i] == 0) && (q_dist[i] > *dlmax) )
2530  *dlmax = q_dist[i];
2531  }
2532 
2533  /* Jetzt die Kantenarrays und ggf. e updaten.
2534  */
2535  for( k = 0; k < p->edges; k += 2 )
2536  {
2537  i = p->head[k];
2538  j = p->tail[k];
2539 
2540  /* both tail and head awake? */
2541  if (!w[i] && !w[j])
2542  {
2543  assert(w[s]);
2544 
2545  excess[i] += transx[k + 1] - transx[k];
2546  excess[j] += transx[k] - transx[k + 1];
2547  transx[k] = 0;
2548  residual[k] = capa[k];
2549  transx[k + 1] = 0;
2550  residual[k + 1] = capa[k + 1];
2551  }
2552  }
2553 
2554  assert(w[t] == 0);
2555  assert(q_dist[t] == 0);
2556  assert(q_temp[0] == t);
2557 
2558  /* Jetzt noch die mit Excess einsortieren.
2559  */
2560  for( i = 1; i < visited; i++ )
2561  {
2562  assert(w[q_temp[i]] == 0);
2563  assert(q_temp[i] != s);
2564  assert(q_temp[i] != t);
2565 
2566  if (excess[q_temp[i]] > 0)
2567  {
2568  listinsert(q_temp[i], q_dist[q_temp[i]], q_next, q_head, actminloc, *dlmax);
2569  }
2570  }
2571 #if DEBUG > 1
2572  show_flow(p, capa, w);
2573 #endif
2574 #if DEBUG > 0
2575  assert(is_valid(p, s, t, capa, w));
2576 #endif
2577 }
2578 
2579 /*---------------------------------------------------------------------------*/
2580 /*--- Name : GRAPH MINimumCUT EXECute ---*/
2581 /*--- Function : Fuehrt den Mincut Algorithmus durch und findet ---*/
2582 /*--- (hoffentlich) einen Minimalen (s,t) Schnitt im Graphen. ---*/
2583 /*--- Parameter: Graph, Quelle, Senke, Kantenkapazitaeten, Zustandsarray, ---*/
2584 /*--- Flag um vorhandenen Fluss zu belassen. ---*/
2585 /*--- Returns : Nichts, fuellt aber w[] mit nicht Nulleintraegen fuer ---*/
2586 /*--- die Knoten, die auf der Quellenseite des Schnittes ---*/
2587 /*--- liegen und Null fuer die auf der Senkenseite. ---*/
2588 /*---------------------------------------------------------------------------*/
2589 void graph_mincut_exec(
2590  GRAPH* p,
2591  int s,
2592  int t,
2593  int nedges,
2594  const int rootcutsize,
2595  const int* rootcut,
2596  const int* capa,
2597  int* w,
2598  const int* edgestart,
2599  const int* edgearr,
2600  const int* headarr,
2601  int rerun)
2602 {
2603  int min_dist;
2604  int min_capa;
2605  int min_knot;
2606  int min_arc;
2607  int dormmax;
2608  int dlmax;
2609  int i;
2610  int k;
2611  int di1;
2612  int glbadd;
2613  int nnodes;
2614  int startind;
2615  int dmin = 1;
2616  int* dist;
2617  int* head;
2618  int* numb;
2619  int* prev;
2620  int* next;
2621  int* temp;
2622  int* e;
2623  int* x;
2624  int* r;
2625  int j;
2626  int end;
2627  int headnode;
2628  int relabelupdatebnd;
2629  int relabeltrigger;
2630 
2631  assert(p != NULL);
2632  assert(s >= 0);
2633  assert(s < p->knots);
2634  assert(t >= 0);
2635  assert(t < p->knots);
2636  assert(s != t);
2637  assert(capa != NULL);
2638  assert(w != NULL);
2639  assert(p->mincut_dist != NULL);
2640  assert(p->mincut_numb != NULL);
2641  assert(p->mincut_head != NULL);
2642  assert(p->mincut_prev != NULL);
2643  assert(p->mincut_next != NULL);
2644  assert(p->mincut_temp != NULL);
2645  assert(p->mincut_e != NULL);
2646  assert(p->mincut_x != NULL);
2647  assert(p->mincut_r != NULL);
2648 
2649  dist = p->mincut_dist;
2650  head = p->mincut_head;
2651  numb = p->mincut_numb;
2652  prev = p->mincut_prev;
2653  next = p->mincut_next;
2654  temp = p->mincut_temp;
2655  e = p->mincut_e;
2656  x = p->mincut_x;
2657  r = p->mincut_r;
2658  nnodes = p->knots;
2659 
2660  relabelupdatebnd = (GLOBALRELABEL_MULT * nnodes) + nedges;
2661  relabeltrigger = 0;
2662 
2663  if( !rerun )
2664  initialise(p, s, t, capa, dist, numb, head, prev, next, temp, e, x, r, w, edgestart,
2665  edgearr, headarr, &dormmax, &dlmax);
2666  else
2667  reinitialise(p, s, t, capa, dist, numb, head, prev, next, temp, e, x, r, w, edgestart,
2668  edgearr, headarr, &dormmax, &dlmax);
2669 
2670  /* main loop */
2671  while( dlmax >= dmin )
2672  {
2673  if( (head[dlmax] == Q_NULL) )
2674  {
2675  dlmax--;
2676  continue;
2677  }
2678 
2679  /* get highest label node */
2680  i = head[dlmax];
2681 
2682  listdelete(i, dlmax, next, head);
2683 
2684  assert(dist[i] == dlmax);
2685 
2686  startind = edgestart[i];
2687  end = edgestart[i + 1];
2688  glbadd = (end - startind);
2689 
2690  /* discharge i */
2691  for( ;; )
2692  {
2693  assert(w[i] == 0);
2694  assert(i != t);
2695  assert(i != s);
2696  assert(e[i] > 0);
2697 
2698  min_knot = -1;
2699  min_dist = nnodes;
2700  min_capa = 0;
2701  min_arc = -1;
2702 
2703  di1 = dist[i] - 1;
2704 
2705  /* traverse incident edges */
2706  for( j = startind; j != end; j++ )
2707  {
2708  /* res. cap. > 0? */
2709  if( r[edgearr[j]] > 0 )
2710  {
2711  /* head node awake? */
2712  if( w[headarr[j]] == 0 )
2713  {
2714  k = edgearr[j];
2715  headnode = headarr[j];
2716 
2717  /* non-admissible edge? */
2718  if( di1 != dist[headnode] )
2719  {
2720  assert(dist[i] <= dist[headnode]);
2721 
2722  if( (dist[headnode] < min_dist) || ((dist[headnode] == min_dist) && (r[k] > min_capa)) )
2723  {
2724  min_knot = headnode;
2725  min_dist = dist[min_knot];
2726  min_capa = r[k];
2727  min_arc = k;
2728  }
2729  }
2730  else /* admissible edge */
2731  {
2732  assert(Min(e[i], r[k]) > 0);
2733 
2734  assert(headnode == p->head[k]);
2735 
2736  if( e[headnode] == 0 && headnode != t )
2737  {
2738  listinsert(headnode, dist[headnode], next, head, dmin, (dlmax));
2739  }
2740 
2741  if( e[i] <= r[k] )
2742  {
2743  #if STATIST
2744  (e[i] == r[k]) ? x_pushes++ : n_pushes++;
2745  #endif
2746  x[k] += e[i];
2747  r[k] -= e[i];
2748  r[Edge_anti(k)] += e[i];
2749  e[headnode] += e[i];
2750  e[i] = 0; /* -= e[i] */
2751 
2752  assert(e[headnode] > 0);
2753  assert(w[headnode] == 0);
2754 
2755  assert(r[k] + r[Edge_anti(k)] == capa[k] + capa[Edge_anti(k)]);
2756  assert(r[k] == capa[k] - x[k] + x[Edge_anti(k)]);
2757 
2758  break;
2759  }
2760 
2761  /* e[i] > r[k] */
2762  r[Edge_anti(k)] += r[k];
2763  e[headnode] += r[k];
2764  e[i] -= r[k];
2765  x[k] += r[k];
2766  r[k] = 0; /* -= r[k] */
2767 
2768  assert(r[k] + r[Edge_anti(k)] == capa[k] + capa[Edge_anti(k)]);
2769  assert(r[k] == capa[k] - x[k] + x[Edge_anti(k)]);
2770  assert(e[i] > 0);
2771 
2772  } /* admissible edge */
2773  } /* head node awake? */
2774  } /* res. cap. > 0? */
2775  #if STATIST
2776  s_pushes++;
2777  #endif
2778  } /* traverse incident edges */
2779 
2780  /* e[i] == 0? */
2781  if( j != end )
2782  {
2783  assert(e[i] == 0);
2784  break;
2785  }
2786 
2787  assert(e[i] > 0);
2788 
2789  /*
2790  * relabel
2791  */
2792  assert(numb[dist[i]] > 0);
2793 
2794  if( numb[dist[i]] == 1 )
2795  {
2796  w[i] = ++dormmax;
2797 
2798  numb[dist[i]] = 0;
2799 
2800 
2801  assert(dormmax <= nnodes);
2802 
2803  for( k = 0; k < nnodes; k++ )
2804  {
2805  /* put all nodes with distance >= dist[i] (including i) into new dormant set */
2806  if (!w[k] && (dist[i] < dist[k]))
2807  {
2808  numb[dist[k]]--;
2809  w[k] = dormmax;
2810  }
2811  }
2812  #if STATIST
2813  m_sleeps++;
2814  #endif
2815  break;
2816  }
2817 
2818  /* no nodes reachable from i? */
2819  if( min_knot < 0 )
2820  {
2821  /* put i into new dormant set */
2822  numb[dist[i]]--;
2823  w[i] = ++dormmax;
2824 
2825  assert(dormmax <= nnodes);
2826 
2827  #if STATIST
2828  s_sleeps++;
2829  #endif
2830  break;
2831  }
2832  else
2833  {
2834  /* set distance label of i to minimum (distance + 1) among reachable, active adjacent nodes */
2835  assert(min_dist < nnodes);
2836  assert(min_capa > 0);
2837  assert(min_knot >= 0);
2838  assert(min_arc >= 0);
2839 
2840  relabeltrigger += GLOBALRELABEL_ADD + glbadd;
2841 
2842  numb[dist[i]]--;
2843 
2844  dist[i] = min_dist + 1;
2845 
2846  numb[dist[i]]++;
2847 
2848  assert(dist[i] < nnodes);
2849 
2850  assert(min_capa > 0);
2851  assert(min_capa == r[min_arc]);
2852  assert(p->head[min_arc] == min_knot);
2853  assert(p->tail[min_arc] == i);
2854  assert(dist[i] == dist[min_knot] + 1);
2855  assert(w[min_knot] == 0);
2856 
2857  if (e[i] <= min_capa)
2858  {
2859  if( (e[p->head[min_arc]] == 0) && (p->head[min_arc] != t) )
2860  {
2861  listinsert(p->head[min_arc], dist[p->head[min_arc]], next, head, dmin, (dlmax));
2862  }
2863 
2864  x[min_arc] += e[i];
2865  r[min_arc] -= e[i];
2866  r[Edge_anti(min_arc)] += e[i];
2867  e[min_knot] += e[i];
2868  e[i] = 0; /* -= e[i] */
2869 
2870  assert(r[min_arc] + r[Edge_anti(min_arc)] == capa[min_arc] + capa[Edge_anti(min_arc)]);
2871  assert(r[min_arc] >= 0);
2872  assert(r[min_arc] == capa[min_arc] - x[min_arc] + x[Edge_anti(min_arc)]);
2873  #if STATIST
2874  m_pushes++;
2875  #endif
2876 
2877  if( relabeltrigger > relabelupdatebnd )
2878  {
2879  /* start global relabel heuristic */
2880  globalrelabel(p, s, t, dist, numb, head, prev, next, temp, e, x, r, w, &dormmax, &dlmax, &dmin);
2881  relabeltrigger = 0;
2882  }
2883 
2884  break;
2885  }
2886  #if STATIST
2887  relabels++;
2888  #endif
2889  if( relabeltrigger > relabelupdatebnd )
2890  {
2891  /* start global relabel heuristic */
2892  globalrelabel(p, s, t, dist, numb, head, prev, next, temp, e, x, r, w, &dormmax, &dlmax, &dmin);
2893  relabeltrigger = 0;
2894  break;
2895  }
2896  }
2897  } /* discharge i */
2898  } /* main loop */
2899 
2900 #if 0
2901  FILE *fptr;
2902 
2903  fptr = fopen("f2", "a");
2904 
2905  fprintf(fptr, "t: %d flow %d \n", t, e[t]);
2906 
2907  fclose(fptr);
2908 #endif
2909 
2910  assert(w[s]);
2911  assert(!w[t]);
2912 
2913 #if STATIST
2914  cut_sum(p, capa, w);
2915 #endif
2916 #if DEBUG > 1
2917  show_flow(p, capa, w);
2918 #endif
2919 #if DEBUG > 0
2920  assert(is_valid(p, s, t, capa, w));
2921 #endif
2922 }
2923 
2924 #if STATIST
2925 static void cut_statist(void)
2926 {
2927  (void)printf("Mincut Statistics:\n");
2928  (void)printf("Node-Searches=%d, Cut Sums=%d\n",
2929  searches, cutsums);
2930  (void)printf("S-Pushes=%d, N-Pushes=%d, X-Pushes=%d, M-Pushes=%d\n",
2931  s_pushes, n_pushes, x_pushes, m_pushes);
2932  (void)printf("Relabels=%d, S-Sleeps=%d, M-Sleeps=%d\n\n",
2933  relabels, s_sleeps, m_sleeps);
2934 
2935  s_pushes = 0;
2936  n_pushes = 0;
2937  x_pushes = 0;
2938  m_pushes = 0;
2939  relabels = 0;
2940  s_sleeps = 0;
2941  m_sleeps = 0;
2942  searches = 0;
2943  cutsums = 0;
2944 }
2945 
2946 static void cut_sum(
2947  const GRAPH* p,
2948  const int* capa,
2949  const int* w)
2950 {
2951  int sum = 0;
2952  int i;
2953  int j;
2954  int k;
2955 
2956  assert(p != NULL);
2957  assert(capa != NULL);
2958  assert(w != NULL);
2959 
2960  for(k = 0; k < p->edges; k++)
2961  {
2962  i = p->head[k];
2963  j = p->tail[k];
2964 
2965  if ((w[i] && !w[j]) || (!w[i] && w[j]))
2966  sum += capa[k];
2967  }
2968 #if DEBUG > 0
2969  (void)printf("Cut Sum=%d\n", sum);
2970 #endif
2971  cutsums += sum;
2972 }
2973 #endif
2974 
2975 #if DEBUG > 0
2976 static int is_valid(
2977  const GRAPH* p,
2978  const int s,
2979  const int t,
2980  const int* capa,
2981  const int* w)
2982 {
2983  int* e;
2984  int* r;
2985  int* x;
2986  int j;
2987  int k;
2988 
2989  assert(p != NULL);
2990  assert(p->mincut_e != NULL);
2991  assert(p->mincut_r != NULL);
2992  assert(p->mincut_x != NULL);
2993 
2994  e = p->mincut_e;
2995  r = p->mincut_r;
2996  x = p->mincut_x;
2997 
2998  for(j = 0; j < p->knots; j++)
2999  {
3000 #if 0
3001  if ((q[j] >= 0) && (a[q[j]] != j))
3002  return((void)fprintf(stderr, "Queue Error 1 Knot %d\n", j), FALSE);
3003 
3004  if (!w[j] && (q[j] < 0) && (e[j] > 0) && (j != t))
3005  return((void)fprintf(stderr, "Queue Error 2 Knot %d\n", j), FALSE);
3006 
3007  if (!w[j] && (q[j] >= 0) && (e[j] == 0))
3008  return((void)fprintf(stderr, "Queue Error 3 Knot %d\n", j), FALSE);
3009 
3010  if (w[j] && (q[j] >= 0))
3011  return((void)fprintf(stderr, "Queue Error 4 Knot %d\n", j), FALSE);
3012 #endif
3013  if (e[j] < 0)
3014  return((void)fprintf(stderr, "Negativ Execess Knot %d\n", j), FALSE);
3015 
3016  if (p->mincut_dist[j] >= p->knots)
3017  return((void)fprintf(stderr, "Distance too big Knot %d\n", j), FALSE);
3018 
3019  /* Extended Dormacy Property
3020  */
3021  for(k = p->outbeg[j]; k != EAT_LAST; k = p->oeat[k])
3022  {
3023  if (r[k] > 0)
3024  {
3025  if ((w[j] && !w[p->head[k]]) || (w[j] && (w[j] < w[p->head[k]])))
3026  {
3027  (void)printf("k=%d r[k]=%d head=%d tail=%d w[h]=%d w[t]=%d\n",
3028  k, r[k], p->head[k], p->tail[k], w[p->head[k]], w[p->tail[k]]);
3029 
3030  return((void)fprintf(stderr, "Extended Dormacy Violation Knot %d\n", j), FALSE);
3031  }
3032  }
3033  }
3034  }
3035  for(j = 0; j < p->edges; j++)
3036  {
3037  if (r[j] < 0)
3038  return((void)fprintf(stderr, "Negativ Residual Edge %d\n", j), FALSE);
3039 
3040  if (x[j] < 0)
3041  return((void)fprintf(stderr, "Negativ Flow Edge %d\n", j), FALSE);
3042 
3043  if (r[j] + r[Edge_anti(j)] != capa[j] + capa[Edge_anti(j)])
3044  return((void)fprintf(stderr, "Wrong Capacity Equation Edge %d\n", j), FALSE);
3045 
3046  if (r[j] != capa[j] - x[j] + x[Edge_anti(j)])
3047  return((void)fprintf(stderr, "Wrong Residual Equation Edge %d\n", j), FALSE);
3048  }
3049  return(TRUE);
3050 }
3051 
3052 static void show_flow(
3053  const GRAPH* p,
3054  const int* capa,
3055  const int* w)
3056 {
3057  int j;
3058  int* head;
3059  int* numb;
3060  int* prev;
3061  int* next;
3062  int* e;
3063  int* x;
3064  int* r;
3065 
3066  assert(p != NULL);
3067  assert(w != NULL);
3068  assert(p->mincut_numb != NULL);
3069  assert(p->mincut_head != NULL);
3070  assert(p->mincut_prev != NULL);
3071  assert(p->mincut_next != NULL);
3072  assert(p->mincut_e != NULL);
3073  assert(p->mincut_x != NULL);
3074  assert(p->mincut_r != NULL);
3075 
3076  head = p->mincut_head;
3077  numb = p->mincut_numb;
3078  prev = p->mincut_prev;
3079  next = p->mincut_next;
3080  e = p->mincut_e;
3081  x = p->mincut_x;
3082  r = p->mincut_r;
3083 
3084 
3085 
3086  (void)printf(" ");
3087  for(j = 0; j < p->edges; j++)
3088  (void)printf("%6d ", j);
3089  (void)fputc('\n', stdout);
3090 
3091  (void)printf("ta:");
3092  for(j = 0; j < p->edges; j++)
3093  (void)printf("%6d ", p->tail[j]);
3094  (void)fputc('\n', stdout);
3095 
3096  (void)printf("he:");
3097  for(j = 0; j < p->edges; j++)
3098  (void)printf("%6d ", p->head[j]);
3099  (void)fputc('\n', stdout);
3100 
3101  (void)printf("x: ");
3102  for(j = 0; j < p->edges; j++)
3103  (void)printf("%6d ", x[j]);
3104  (void)fputc('\n', stdout);
3105 
3106  (void)printf("r: ");
3107  for(j = 0; j < p->edges; j++)
3108  (void)printf("%6d ", r[j]);
3109  (void)fputc('\n', stdout);
3110 
3111  (void)printf("ca:");
3112  for(j = 0; j < p->edges; j++)
3113  (void)printf("%6d ", capa[j]);
3114  (void)fputc('\n', stdout);
3115 
3116  (void)printf("w: ");
3117  for(j = 0; j < p->knots; j++)
3118  (void)printf("%2d ", w[j]);
3119  (void)fputc('\n', stdout);
3120 
3121  (void)printf("d: ");
3122  for(j = 0; j < p->knots; j++)
3123  (void)printf("%2d ", p->mincut_dist[j]);
3124  (void)fputc('\n', stdout);
3125 
3126  (void)printf("n: ");
3127  for(j = 0; j < p->knots; j++)
3128  (void)printf("%2d ", numb[j]);
3129  (void)fputc('\n', stdout);
3130 
3131  (void)printf("h: ");
3132  for(j = 0; j < p->knots; j++)
3133  (void)printf("%2d ", head[j]);
3134  (void)fputc('\n', stdout);
3135 
3136  (void)printf("p: ");
3137  for(j = 0; j < p->knots; j++)
3138  (void)printf("%2d ", prev[j]);
3139  (void)fputc('\n', stdout);
3140 
3141  (void)printf("n: ");
3142  for(j = 0; j < p->knots; j++)
3143  (void)printf("%2d ", next[j]);
3144  (void)fputc('\n', stdout);
3145 
3146  (void)printf("e: ");
3147  for(j = 0; j < p->knots; j++)
3148  (void)printf("%2d ", e[j]);
3149  (void)fputc('\n', stdout);
3150 }
3151 
3152 #endif /* DEBUG > 0 */
3153 #endif
SCIPfreeBlockMemoryArray
#define SCIPfreeBlockMemoryArray(scip, ptr, num)
Definition: scip_mem.h:101
GRAPH::mincut_e
int *RESTRICT mincut_e
Definition: graphdefs.h:250
SCIPallocBlockMemoryArray
#define SCIPallocBlockMemoryArray(scip, ptr, num)
Definition: scip_mem.h:84
GRAPH::head
int *RESTRICT head
Definition: graphdefs.h:224
GRAPH::mincut_x
int *RESTRICT mincut_x
Definition: graphdefs.h:251
SCIPfreeMemoryArray
#define SCIPfreeMemoryArray(scip, ptr)
Definition: scip_mem.h:71
graph_mincut_init
SCIP_RETCODE graph_mincut_init(SCIP *scip, GRAPH *p)
Definition: graph_mcut.c:1105
Q_NULL
#define Q_NULL
Definition: graph_mcut.c:55
initialise
static void initialise(const int s, const int t, const int nnodesreal, const int rootcutsize, const int *rootcut, const int *capa, int *RESTRICT dist, int *RESTRICT headactive, int *RESTRICT headinactive, int *RESTRICT edgecurr, int *RESTRICT next, int *RESTRICT prev, int *RESTRICT temp, int *RESTRICT excess, int *RESTRICT residual, int *RESTRICT w, const int *edgestart, const int *edgearr, const int *headarr, int *dormmax, int *actmin, int *actmax, int *glbmax)
Definition: graph_mcut.c:566
inactivedelete
#define inactivedelete(node, nodedist, next, prev, headinactive, nextnode, prevnode)
Definition: graph_mcut.c:127
SCIPallocMemoryArray
#define SCIPallocMemoryArray(scip, ptr, num)
Definition: scip_mem.h:55
GRAPH::mincut_temp
int *RESTRICT mincut_temp
Definition: graphdefs.h:249
FALSE
#define FALSE
Definition: def.h:87
TRUE
#define TRUE
Definition: def.h:86
SCIP_RETCODE
enum SCIP_Retcode SCIP_RETCODE
Definition: type_retcode.h:54
graph.h
includes various files containing graph methods used for Steiner tree problems
mincutInit
static SCIP_RETCODE mincutInit(SCIP *scip, int nnodes, int nedges, GRAPH *p)
Definition: graph_mcut.c:349
EAT_LAST
#define EAT_LAST
Definition: graphdefs.h:58
activedelete
#define activedelete(node, nodedist, next, headactive)
Definition: graph_mcut.c:101
x
SCIP_VAR ** x
Definition: circlepacking.c:54
w
SCIP_VAR * w
Definition: circlepacking.c:58
GRAPH::oeat
int *RESTRICT oeat
Definition: graphdefs.h:231
graph_mincut_reInit
SCIP_RETCODE graph_mincut_reInit(SCIP *scip, int nnodes, int nedges, GRAPH *p)
Definition: graph_mcut.c:1120
GRAPH::mincut_dist
int *RESTRICT mincut_dist
Definition: graphdefs.h:243
graph_mincut_setDefaultVals
void graph_mincut_setDefaultVals(GRAPH *g)
Definition: graph_mcut.c:1081
Min
#define Min(a, b)
Definition: portab.h:55
reinitialise
static void reinitialise(const int s, const int t, const int nnodesreal, const int rootcutsize, const int *rootcut, const int *capa, int *RESTRICT dist, int *RESTRICT headactive, int *RESTRICT headinactive, int *RESTRICT edgecurr, int *RESTRICT next, int *RESTRICT prev, int *RESTRICT temp, int *RESTRICT excess, int *RESTRICT residual, int *RESTRICT w, const int *edgestart, const int *edgearr, const int *headarr, int *dormmax, int *actmin, int *actmax, int *glbmax)
Definition: graph_mcut.c:833
GLOBALRELABEL_MULT
#define GLOBALRELABEL_MULT
Definition: graph_mcut.c:57
GRAPH::mincut_head_inact
int *RESTRICT mincut_head_inact
Definition: graphdefs.h:245
NULL
#define NULL
Definition: lpi_spx1.cpp:155
listdelete
#define listdelete(node, nodedist, next, headactive)
Definition: graph_mcut.c:60
SCIP_OKAY
Definition: type_retcode.h:33
Edge_anti
#define Edge_anti(a)
Definition: graphdefs.h:106
GRAPH::knots
int knots
Definition: graphdefs.h:190
SCIP_CALL
#define SCIP_CALL(x)
Definition: def.h:384
GRAPH::mincut_nedges
int mincut_nedges
Definition: graphdefs.h:242
listinsert
#define listinsert(node, nodedist, next, headactive, actmin, actmax)
Definition: graph_mcut.c:69
globalrelabel
static void globalrelabel(const int s, const int t, const int nnodesreal, int *RESTRICT dist, int *RESTRICT headactive, int *RESTRICT headinactive, int *RESTRICT edgecurr, int *RESTRICT next, int *RESTRICT prev, int *RESTRICT excess, int *RESTRICT residual, int *RESTRICT w, const int *edgestart, const int *edgearr, const int *headarr, int *actmin, int *actmax, int *glbmax, int *dormmax)
Definition: graph_mcut.c:375
GRAPH::mincut_head
int *RESTRICT mincut_head
Definition: graphdefs.h:244
SCIP_Bool
#define SCIP_Bool
Definition: def.h:84
GRAPH::tail
int *RESTRICT tail
Definition: graphdefs.h:223
flipedge
#define flipedge(edge)
Definition: graphdefs.h:84
GRAPH::mincut_nnodes
int mincut_nnodes
Definition: graphdefs.h:241
GRAPH::mincut_prev
int *RESTRICT mincut_prev
Definition: graphdefs.h:247
portab.h
Portable definitions.
GRAPH::mincut_numb
int *RESTRICT mincut_numb
Definition: graphdefs.h:246
r
SCIP_Real * r
Definition: circlepacking.c:50
a
SCIP_VAR * a
Definition: circlepacking.c:57
activeinsert
#define activeinsert(node, nodedist, next, headactive, actmin, actmax, glbmax)
Definition: graph_mcut.c:110
GRAPH::outbeg
int *RESTRICT outbeg
Definition: graphdefs.h:204
GRAPH::edges
int edges
Definition: graphdefs.h:219
GRAPH::mincut_r
int *RESTRICT mincut_r
Definition: graphdefs.h:252
nnodes
#define nnodes
Definition: gastrans.c:65
graph_mincut_exit
void graph_mincut_exit(SCIP *scip, GRAPH *p)
Definition: graph_mcut.c:1175
inactiveinsert
#define inactiveinsert(node, nodedist, next, prev, headinactive, nextnode)
Definition: graph_mcut.c:153
graph_mincut_isInitialized
SCIP_Bool graph_mincut_isInitialized(const GRAPH *p)
Definition: graph_mcut.c:1139
graph_mincut_exec
void graph_mincut_exec(const GRAPH *p, const int s, const int t, const int nnodesreal, const int nedgesreal, const int rootcutsize, const int *rootcut, const int *capa, int *RESTRICT w, const int *edgestart, const int *edgearr, const int *headarr, const SCIP_Bool isRerun)
Definition: graph_mcut.c:1202
GLOBALRELABEL_ADD
#define GLOBALRELABEL_ADD
Definition: graph_mcut.c:56
GRAPH::mincut_next
int *RESTRICT mincut_next
Definition: graphdefs.h:248