# 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 /* */
7 /* fuer Informationstechnik Berlin */
8 /* */
10 /* */
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 */
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 /** initialize min cut arrays */
276  SCIP* scip, /**< SCIP data structure */
277  GRAPH* p /**< graph data structure */
278  )
279 {
280  assert(p != NULL);
281  assert(p->mincut_dist == 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 */
306  SCIP* scip, /**< SCIP data structure */
307  GRAPH* p /**< graph data structure */
308  )
309 {
310  assert(p->mincut_dist != 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));
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,
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  {
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,
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,
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);
452  assert(edgecurr != NULL);
453  assert(prev != NULL);
454  assert(next != NULL);
455  assert(residual != NULL);
456  assert(excess != NULL);
457  assert(edgearr != NULL);
459  assert(edgecurr != NULL);
460  assert(edgestart != 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  {
476  if( w[i] == 0 )
477  w[i] = -1;
478  }
479
480  for( i = 0; i < nnodes; i++ )
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? */
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  {
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  {
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
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,
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,
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;
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);
650  assert(prev != NULL);
651  assert(next != NULL);
652  assert(temp != NULL);
653  assert(residual != NULL);
654  assert(excess != NULL);
655  assert(edgearr != NULL);
657  assert(edgecurr != NULL);
658  assert(edgestart != 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];
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
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
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,
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,
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;
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);
930  assert(edgecurr != NULL);
931  assert(edgestart != 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  {
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  {
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  {
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  {
1072  a = edgearr[j];
1073
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 */
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,
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;
1160  int j;
1161  int end;
1162  int nnodes;
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);
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;
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
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 */
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  {
1252
1255
1256  dminus1 = dist[i] - 1;
1257
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 */
1270
1273  }
1274  }
1275
1276  /* not more residual capacity than excess? */
1277  if( r[j] < e[i] )
1278  {
1279  e[i] -= 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];
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]? */
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  {
1324
1325  for( l = headinactive[j]; l >= 0; l = next[l] )
1326  {
1327  if( w[l] == 0 )
1328  w[l] = dormmax;
1329  }
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  {
1366  minedgestart = j;
1367  }
1368 #endif
1369  if( dist[headarr[j]] <= mindist )
1370  {
1371  if( dist[headarr[j]] < mindist )
1372  {
1374  minedgestart = j;
1375  maxresval = r[j];
1376  maxresedge = j;
1379  }
1380  else if( r[j] > maxresval )
1381  {
1382  maxresval = r[j];
1383  maxresedge = 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);
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
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
1489
1490  if( relabeltrigger > relabelupdatebnd )
1491  {
1492  if( actmax < actmin )
1493  break;
1494
1495  /* execute global relabel heuristic */
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 */
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  {
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
1698  {
1699  (void)printf("k=%d r[k]=%d head=%d tail=%d w[h]=%d w[t]=%d\n",
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;
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);
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
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++)
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++)
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 /*---------------------------------------------------------------------------*/
1840  SCIP* scip, /**< SCIP data structure */
1841  GRAPH* p /**< graph data structure */
1842  )
1843 {
1844  assert(p != NULL);
1845  assert(p->mincut_dist == 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
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);
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);
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));
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,
1933  int* q_prev,
1934  int* q_next)
1935 {
1936  assert(knot >= 0);
1937  assert(q_dist != 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);
1951
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,
1982  int* q_prev,
1983  int* q_next,
1984  int dmin,
1985  int* dlmax)
1986 {
1987  assert(knot >= 0);
1988  assert(q_dist != 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;
1999
2000  if( q_next[knot] != Q_NULL )
2001  q_prev[q_next[knot]] = knot;
2002
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,
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  {
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,
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);
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;
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  {
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,
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,
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);
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;
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
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,
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,
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);
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;
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  {
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,
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;
2532  int nnodes;
2533  int startind;
2534  int dmin = 1;
2535  int* dist;
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;
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);
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;
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,
2585  else
2586  reinitialise(p, s, t, capa, dist, numb, head, prev, next, temp, e, x, r, w, edgestart,
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 */
2600
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];
2635
2637  if( di1 != dist[headnode] )
2638  {
2640
2641  if( (dist[headnode] < min_dist) || ((dist[headnode] == min_dist) && (r[k] > min_capa)) )
2642  {
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
2654
2656  {
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];
2669  e[i] = 0; /* -= e[i] */
2670
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];
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
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]);
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  {
2779  {
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  {
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  {
2945  {
2946  (void)printf("k=%d r[k]=%d head=%d tail=%d w[h]=%d w[t]=%d\n",
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;
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);
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
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++)
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++)
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
#define SCIPfreeBlockMemoryArray(scip, ptr, num)
Definition: scip_mem.h:97
int *RESTRICT mincut_e
Definition: grph.h:113
#define SCIPallocBlockMemoryArray(scip, ptr, num)
Definition: scip_mem.h:80
Definition: grph.h:96
int *RESTRICT mincut_x
Definition: grph.h:114
Definition: grph.h:57
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
#define SCIPfreeMemoryArray(scip, ptr)
Definition: scip_mem.h:69
#define SCIPallocMemoryArray(scip, ptr, num)
Definition: scip_mem.h:53
#define activeinsert(node, nodedist, next, headactive, actmin, actmax, glbmax)
Definition: grphmcut.c:110
#define EAT_LAST
Definition: grph.h:31
#define Edge_anti(a)
Definition: grph.h:171
int *RESTRICT mincut_temp
Definition: grph.h:112
#define FALSE
Definition: def.h:73
#define TRUE
Definition: def.h:72
enum SCIP_Retcode SCIP_RETCODE
Definition: type_retcode.h:54
Definition: grphmcut.c:101
SCIP_VAR ** x
Definition: circlepacking.c:54
#define inactivedelete(node, nodedist, next, prev, headinactive, nextnode, prevnode)
Definition: grphmcut.c:127
Definition: grphmcut.c:56
SCIP_VAR * w
Definition: circlepacking.c:58
int *RESTRICT oeat
Definition: grph.h:103
int *RESTRICT mincut_dist
Definition: grph.h:106
#define Min(a, b)
Definition: portab.h:51
#define Q_NULL
Definition: grphmcut.c:55
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
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
Definition: grph.h:108
#define NULL
Definition: lpi_spx1.cpp:155
int knots
Definition: grph.h:62
#define SCIP_CALL(x)
Definition: def.h:370
SCIP_RETCODE graph_mincut_init(SCIP *scip, GRAPH *p)
Definition: grphmcut.c:275
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
Definition: grph.h:107
#define GLOBALRELABEL_MULT
Definition: grphmcut.c:57
#define SCIP_Bool
Definition: def.h:70
int *RESTRICT tail
Definition: grph.h:95
int *RESTRICT mincut_prev
Definition: grph.h:110
includes various files containing graph methods used for Steiner tree problems
Portable defintions.
int *RESTRICT mincut_numb
Definition: grph.h:109
SCIP_Real * r
Definition: circlepacking.c:50
void graph_mincut_exit(SCIP *scip, GRAPH *p)
Definition: grphmcut.c:305
SCIP_VAR * a
Definition: circlepacking.c:57
int *RESTRICT outbeg
Definition: grph.h:76
#define listinsert(node, nodedist, next, headactive, actmin, actmax)
Definition: grphmcut.c:69
int edges
Definition: grph.h:92
#define flipedge(edge)
Definition: grph.h:150
int *RESTRICT mincut_r
Definition: grph.h:115
#define nnodes
Definition: gastrans.c:65