Scippy

SCIP

Solving Constraint Integer Programs

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