# 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 /* */
7 /* fuer Informationstechnik Berlin */
8 /* */
10 /* */
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 */
61 {\
62  assert(node >= 0);\
63  assert(nodedist >= 0);\
66 }
67
70 {\
71  assert(node >= 0);\
72  assert(next != NULL);\
74  assert(nodedist >= 0);\
77  if( nodedist < (actmin) )\
78  actmin = nodedist;\
79  if( nodedist > (actmax) )\
80  actmax = nodedist;\
81 }
82 #if 0
85 {\
86  assert(node >= 0);\
87  assert(next != NULL);\
89  assert(nodedist >= 0);\
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 */
102 {\
103  assert(node >= 0);\
104  assert(nodedist >= 0);\
107 }
108
111 {\
112  assert(node >= 0);\
113  assert(next != NULL);\
115  assert(nodedist >= 0);\
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 */
128 {\
129  assert(node >= 0);\
130  assert(nodedist >= 0);\
131  assert(next != NULL);\
132  assert(prev != NULL);\
134  nextnode = next[node];\
135  if( headinactive[nodedist] == node )\
136  {\
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
154 {\
155  assert(node >= 0);\
156  assert(nodedist >= 0);\
157  assert(next != NULL);\
158  assert(prev != NULL);\
161  next[node] = nextnode;\
162  prev[node] = Q_NULL;\
163  if( nextnode >= 0 )\
164  prev[nextnode] = 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,
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;
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  {
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",
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,
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  {
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
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,
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,
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);
415  assert(edgecurr != NULL);
416  assert(prev != NULL);
417  assert(next != NULL);
418  assert(residual != NULL);
419  assert(excess != NULL);
420  assert(edgearr != NULL);
422  assert(edgecurr != NULL);
423  assert(edgestart != NULL);
426
427  assert(end < nnodes);
428  assert(w[t] == 0);
429
430  for( int i = 0; i <= end; i++ )
431  {
434  if( w[i] == 0 )
435  w[i] = -1;
436  }
437
438 #ifndef NDEBUG
439  for( int i = 0; i < nnodes; i++ )
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? */
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
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
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,
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,
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;
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);
614  assert(prev != NULL);
615  assert(next != NULL);
616  assert(temp != NULL);
617  assert(residual != NULL);
618  assert(excess != NULL);
619  assert(edgearr != NULL);
621  assert(edgecurr != NULL);
622  assert(edgestart != 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];
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
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
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
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,
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,
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;
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);
887  assert(edgecurr != NULL);
888  assert(edgestart != 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  {
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  {
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  {
1026  a = edgearr[j];
1027
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*/
1082  GRAPH* g /**< graph data structure */
1083  )
1084 {
1085  const int nnodes = g->mincut_nnodes;
1086  int* RESTRICT excess = g->mincut_e;
1089
1091  assert(g->mincut_nnodes >= g->knots);
1092
1093  for( int k = 0; k < nnodes; k++ )
1095
1096  for( int k = 0; k < nnodes; k++ )
1098
1099  for( int k = 0; k < nnodes; k++ )
1100  excess[k] = 0;
1101 }
1102
1103
1104 /** initialize min cut arrays */
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 */
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 {
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? */
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);
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);
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 */
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));
1194  SCIPfreeMemoryArray(scip, &(p->mincut_dist));
1195
1196 #if STATIST
1197  cut_statist();
1198 #endif
1199 }
1200
1201 /** finds a minimum s-t cut */
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,
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;
1236  int j;
1237  int end;
1238  int nnodes;
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);
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;
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 */
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  {
1332
1335
1336  dminus1 = dist[i] - 1;
1337
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 */
1350
1353  }
1354  }
1355
1356  /* not more residual capacity than excess? */
1357  if( r[j] < e[i] )
1358  {
1359  e[i] -= 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];
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]? */
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  {
1404
1405  for( l = headinactive[j]; l >= 0; l = next[l] )
1406  {
1407  if( w[l] == 0 )
1408  w[l] = dormmax;
1409  }
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  {
1446  minedgestart = j;
1447  }
1448 #endif
1449  if( dist[headarr[j]] <= mindist )
1450  {
1451  if( dist[headarr[j]] < mindist )
1452  {
1454  minedgestart = j;
1455  maxresval = r[j];
1456  maxresedge = j;
1459  }
1460  else if( r[j] > maxresval )
1461  {
1462  maxresval = r[j];
1463  maxresedge = 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);
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
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
1569
1570  if( relabeltrigger > relabelupdatebnd )
1571  {
1572  if( actmax < actmin )
1573  break;
1574
1575  /* execute global relabel heuristic */
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 */
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  {
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
1779  {
1780  (void)printf("k=%d r[k]=%d head=%d tail=%d w[h]=%d w[t]=%d\n",
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;
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);
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
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++)
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++)
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 /*---------------------------------------------------------------------------*/
1921  SCIP* scip, /**< SCIP data structure */
1922  GRAPH* p /**< graph data structure */
1923  )
1924 {
1925  assert(p != NULL);
1926  assert(p->mincut_dist == 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
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);
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);
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));
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,
2014  int* q_prev,
2015  int* q_next)
2016 {
2017  assert(knot >= 0);
2018  assert(q_dist != 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);
2032
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,
2063  int* q_prev,
2064  int* q_next,
2065  int dmin,
2066  int* dlmax)
2067 {
2068  assert(knot >= 0);
2069  assert(q_dist != 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;
2080
2081  if( q_next[knot] != Q_NULL )
2082  q_prev[q_next[knot]] = knot;
2083
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,
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  {
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,
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);
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;
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  {
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,
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,
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);
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;
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
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,
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,
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);
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;
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  {
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,
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;
2613  int nnodes;
2614  int startind;
2615  int dmin = 1;
2616  int* dist;
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;
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);
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;
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,
2666  else
2667  reinitialise(p, s, t, capa, dist, numb, head, prev, next, temp, e, x, r, w, edgestart,
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 */
2681
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];
2716
2718  if( di1 != dist[headnode] )
2719  {
2721
2722  if( (dist[headnode] < min_dist) || ((dist[headnode] == min_dist) && (r[k] > min_capa)) )
2723  {
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
2735
2737  {
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];
2750  e[i] = 0; /* -= e[i] */
2751
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];
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
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]);
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  {
2860  {
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  {
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  {
3026  {
3027  (void)printf("k=%d r[k]=%d head=%d tail=%d w[h]=%d w[t]=%d\n",
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;
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);
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
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++)
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++)
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
#define SCIPfreeBlockMemoryArray(scip, ptr, num)
Definition: scip_mem.h:101
int *RESTRICT mincut_e
Definition: graphdefs.h:250
#define SCIPallocBlockMemoryArray(scip, ptr, num)
Definition: scip_mem.h:84
Definition: graphdefs.h:224
int *RESTRICT mincut_x
Definition: graphdefs.h:251
#define SCIPfreeMemoryArray(scip, ptr)
Definition: scip_mem.h:71
SCIP_RETCODE graph_mincut_init(SCIP *scip, GRAPH *p)
Definition: graph_mcut.c:1105
#define Q_NULL
Definition: graph_mcut.c:55
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
#define inactivedelete(node, nodedist, next, prev, headinactive, nextnode, prevnode)
Definition: graph_mcut.c:127
#define SCIPallocMemoryArray(scip, ptr, num)
Definition: scip_mem.h:55
int *RESTRICT mincut_temp
Definition: graphdefs.h:249
#define FALSE
Definition: def.h:87
#define TRUE
Definition: def.h:86
enum SCIP_Retcode SCIP_RETCODE
Definition: type_retcode.h:54
includes various files containing graph methods used for Steiner tree problems
static SCIP_RETCODE mincutInit(SCIP *scip, int nnodes, int nedges, GRAPH *p)
Definition: graph_mcut.c:349
#define EAT_LAST
Definition: graphdefs.h:58
Definition: graph_mcut.c:101
SCIP_VAR ** x
Definition: circlepacking.c:54
SCIP_VAR * w
Definition: circlepacking.c:58
int *RESTRICT oeat
Definition: graphdefs.h:231
SCIP_RETCODE graph_mincut_reInit(SCIP *scip, int nnodes, int nedges, GRAPH *p)
Definition: graph_mcut.c:1120
int *RESTRICT mincut_dist
Definition: graphdefs.h:243
void graph_mincut_setDefaultVals(GRAPH *g)
Definition: graph_mcut.c:1081
#define Min(a, b)
Definition: portab.h:55
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
#define GLOBALRELABEL_MULT
Definition: graph_mcut.c:57
Definition: graphdefs.h:245
#define NULL
Definition: lpi_spx1.cpp:155
Definition: graph_mcut.c:60
#define Edge_anti(a)
Definition: graphdefs.h:106
int knots
Definition: graphdefs.h:190
#define SCIP_CALL(x)
Definition: def.h:384
int mincut_nedges
Definition: graphdefs.h:242
#define listinsert(node, nodedist, next, headactive, actmin, actmax)
Definition: graph_mcut.c:69
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
Definition: graphdefs.h:244
#define SCIP_Bool
Definition: def.h:84
int *RESTRICT tail
Definition: graphdefs.h:223
#define flipedge(edge)
Definition: graphdefs.h:84
int mincut_nnodes
Definition: graphdefs.h:241
int *RESTRICT mincut_prev
Definition: graphdefs.h:247
Portable definitions.
int *RESTRICT mincut_numb
Definition: graphdefs.h:246
SCIP_Real * r
Definition: circlepacking.c:50
SCIP_VAR * a
Definition: circlepacking.c:57
#define activeinsert(node, nodedist, next, headactive, actmin, actmax, glbmax)
Definition: graph_mcut.c:110
int *RESTRICT outbeg
Definition: graphdefs.h:204
int edges
Definition: graphdefs.h:219
int *RESTRICT mincut_r
Definition: graphdefs.h:252
#define nnodes
Definition: gastrans.c:65
void graph_mincut_exit(SCIP *scip, GRAPH *p)
Definition: graph_mcut.c:1175
#define inactiveinsert(node, nodedist, next, prev, headinactive, nextnode)
Definition: graph_mcut.c:153
SCIP_Bool graph_mincut_isInitialized(const GRAPH *p)
Definition: graph_mcut.c:1139
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