1 /*===========================================================================*/
2 /*                                                                           */
3 /* This file is part of a demonstration application for use with the         */
4 /* SYMPHONY Branch, Cut, and Price Library. This application is a solver for */
5 /* Capacitated Network Routing Problems.                                     */
6 /*                                                                           */
7 /* (c) Copyright 2000-2007 Ted Ralphs. All Rights Reserved.                  */
8 /*                                                                           */
9 /* This application was developed by Ted Ralphs (ted@lehigh.edu)             */
10 /*                                                                           */
11 /* This software is licensed under the Eclipse Public License. Please see    */
12 /* accompanying file for terms.                                              */
13 /*                                                                           */
14 /*===========================================================================*/
15 
16 /* system include files */
17 #include <stdlib.h>
18 #include <math.h>
19 #include <string.h>           /* memset() is defined here in LINUX */
20 
21 /* SYMPHONY include files */
22 #include "sym_constants.h"
23 #include "sym_macros.h"
24 #include "sym_messages.h"
25 #include "sym_proccomm.h"
26 #include "sym_cg.h"
27 
28 /* CNRP include files */
29 #include "cnrp_cg.h"
30 #include "cnrp_macros.h"
31 #include "network.h"
32 
33 /*__BEGIN_EXPERIMENTAL_SECTION__*/
34 #if 0
35 #include "util.h"
36 extern CCrandstate rand_state;
37 #endif
38 /*___END_EXPERIMENTAL_SECTION___*/
39 
40 /*===========================================================================*/
41 
42 #define SEND_DIR_SUBTOUR_CONSTRAINT(num_nodes, total_demand)                 \
43 new_cut->type = (num_nodes < vertnum/2 ?                                     \
44 		SUBTOUR_ELIM_SIDE:SUBTOUR_ELIM_ACROSS);                      \
45 new_cut->rhs = (new_cut->type == SUBTOUR_ELIM_SIDE ?                         \
46                RHS(num_nodes, total_demand, capacity) :                      \
47 	       BINS(total_demand, capacity));                                \
48 cuts_found += cg_send_cut(new_cut, num_cuts, alloc_cuts, cuts);              \
49 
50 #define SEND_SUBTOUR_CONSTRAINT(num_nodes, total_demand)                     \
51 if (mult - 1){                                                               \
52    new_cut->type = (num_nodes < vertnum/2 ?                                  \
53 		    SUBTOUR_ELIM_SIDE:SUBTOUR_ELIM_ACROSS);                  \
54    new_cut->rhs = (new_cut->type == SUBTOUR_ELIM_SIDE ?                      \
55 		   RHS(num_nodes, total_demand, capacity) :                  \
56 		   mult*BINS(total_demand, capacity));                       \
57    cuts_found += cg_send_cut(new_cut, num_cuts, alloc_cuts, cuts);           \
58 }else{                                                                       \
59    new_cut->type = SUBTOUR_ELIM_ACROSS;                                      \
60    new_cut->rhs = mult*BINS(total_demand, capacity);                         \
61    cuts_found += cg_send_cut(new_cut, num_cuts, alloc_cuts, cuts);           \
62 }                                                                            \
63 
64 /*===========================================================================*/
65 
66 /*===========================================================================*\
67  * This file implements the greedy shrinking algorithm of Augerat, et al.
68  * The original implementation was done by Leonid Kopman.
69 \*===========================================================================*/
70 
reduce_graph(network * n,double etol,double * demand,double capacity,int mult,cut_data * new_cut,int * num_cuts,int * alloc_cuts,cut_data *** cuts)71 int reduce_graph(network *n, double etol, double *demand, double capacity,
72 		 int mult, cut_data *new_cut, int *num_cuts, int *alloc_cuts,
73 		 cut_data ***cuts)
74 {
75    elist *e1, *e2, *e3;
76    edge *cur_edge;
77    int v1, v2, deg, count, i, k, vertnum = n->vertnum;
78    edge *edges = n->edges;
79    int num_edges = n->edgenum, cuts_found = 0;
80    int edges_deleted = 0;
81    vertex *verts = n->verts;
82    vertex *v2_pt, *third_node;
83    char *coef;
84 
85    new_cut->size = (vertnum >> DELETE_POWER) + 1;
86    new_cut->coef = coef = (char *) (calloc(new_cut->size, CSIZE));
87 
88    while(TRUE){
89       num_edges = n->edgenum;
90       edges_deleted = 0;
91       for (i = 0; i < num_edges; i++){
92 	 cur_edge = edges + i;
93 	 if (cur_edge->weight >= 1 - etol && cur_edge->v0 &&
94 	     cur_edge->v1 && !(cur_edge->deleted)){
95 	    cur_edge->deleted = TRUE;
96 	    n->edgenum--;
97 	    edges_deleted++;
98 	    v1 = (verts[cur_edge->v0].degree ==
99 		  MIN(verts[cur_edge->v0].degree, verts[cur_edge->v1].degree))?
100 	          cur_edge->v0 : cur_edge->v1;
101 	    v2 = (v1 == cur_edge->v0) ? cur_edge->v1 : cur_edge->v0;
102 	    if (cur_edge->weight + verts[v1].weight + verts[v2].weight >
103 		RHS(verts[v1].orig_node_list_size +
104 		    verts[v2].orig_node_list_size + 2,
105 		    demand[v1]+demand[v2], capacity)){
106 	       memset(coef, 0, new_cut->size*CSIZE);
107 	       for (k = 0; k < verts[v1].orig_node_list_size; k++)
108 		  (coef[(verts[v1].orig_node_list)[k] >>
109 		       DELETE_POWER]) |=
110 		     (1 << ((verts[v1].orig_node_list)[k] &
111 			    DELETE_AND));
112 	       (coef[v1 >> DELETE_POWER]) |=
113 			(1 << (v1 & DELETE_AND));
114 	       for (k = 0; k < verts[v2].orig_node_list_size; k++)
115 		  (coef[(verts[v2].orig_node_list)[k] >>
116 		       DELETE_POWER]) |=
117 		     (1 << ((verts[v2].orig_node_list)[k] &
118 			    DELETE_AND));
119 	       (coef[v2 >> DELETE_POWER]) |=
120 			(1 << (v2 & DELETE_AND));
121 #if 0
122 	       new_cut->type = SUBTOUR_ELIM_SIDE;
123 	       new_cut->rhs =  RHS(verts[v1].orig_node_list_size +
124 				   verts[v2].orig_node_list_size + 2,
125 				   demand[v1]+demand[v2], capacity);
126 	       cuts_found += cg_send_cut(new_cut);
127 #endif
128 #ifdef DIRECTED_X_VARS
129 	       SEND_DIR_SUBTOUR_CONSTRAINT(verts[v1].orig_node_list_size +
130 					   verts[v2].orig_node_list_size + 2,
131 					   demand[v1]+demand[v2]);
132 #else
133 	       SEND_SUBTOUR_CONSTRAINT(verts[v1].orig_node_list_size +
134 				       verts[v2].orig_node_list_size + 2,
135 				       demand[v1]+demand[v2]);
136 #endif
137 	    }
138 	    verts[v1].deleted = TRUE;
139 	    demand[v2] += demand[v1];
140 	    demand[v1] = 0;
141 	    verts[v2].weight += cur_edge->weight + verts[v1].weight;
142 	    v2_pt = verts + v2;
143 	    v2_pt->degree--;
144 	    if (v2_pt->first->other_end == v1){
145 	       v2_pt->first = v2_pt->first->next_edge;
146 	    }else{
147 	       for (e3 = v2_pt->first; e3 && e3->next_edge; e3 = e3->next_edge)
148 		  if (e3->next_edge->other_end == v1){
149 		     e3->next_edge = e3->next_edge->next_edge;
150 		     if (e3->next_edge == NULL) v2_pt->last = e3;
151 		     break;
152 		  }
153 	    }
154 
155 	    if (!(v2_pt->orig_node_list_size))
156 	       v2_pt->orig_node_list = (int *) malloc(vertnum*ISIZE);
157 	    (v2_pt->orig_node_list)[(v2_pt->orig_node_list_size)++] = v1;
158 
159 	    for (k = 0; k < verts[v1].orig_node_list_size; k++){
160 	       (v2_pt->orig_node_list)[(v2_pt->orig_node_list_size)++] =
161 		  (verts[v1].orig_node_list)[k];
162 	    }
163 	    deg = verts[v1].degree;
164 
165 	    for (e1 = verts[v1].first, count=0; e1 && (count < deg); count++ ){
166 	       third_node = e1->other;
167 	       if (third_node->orignodenum == v2){
168 		  e1 = e1->next_edge;
169 		  continue;
170 	       }
171 	       for (e2 = v2_pt->first; e2; e2 = e2->next_edge){
172 		  if (e2->other_end == e1->other_end ){
173 		     e2->data->weight += e1->data->weight;
174 		     e1->data->deleted = TRUE;
175 		     edges_deleted++;
176 		     (third_node->degree)--;
177 		     if (third_node->first->other_end == v1){
178 			third_node->first=third_node->first->next_edge;
179 		     }else{
180 			for (e3 = third_node->first; e3 && e3->next_edge;
181 			     e3 = e3->next_edge)
182 			   if (e3->next_edge->other_end == v1){
183 			      e3->next_edge = e3->next_edge->next_edge;
184 			      if (e3->next_edge == NULL) third_node->last = e3;
185 			      break;
186 			   }
187 		     }
188 		     break;
189 		  }
190 	       }
191 	       if (e2){
192 		  e1 = e1->next_edge;
193 		  continue;
194 	       }
195 	       /* ok, so e1->other_node is not incident to v2 */
196 	       for (e3 = third_node->first; e3 ; e3 = e3->next_edge){
197 		  if (e3->other_end == v1){
198 		     e3->other = v2_pt;
199 		     e3->other_end = v2;
200 		     e3->data->v0 = MIN(v2, third_node->orignodenum);
201 		     e3->data->v1 = MAX(v2, third_node->orignodenum);
202 		     break;
203 		  }
204 	       }
205 	       v2_pt->last->next_edge = e1;
206 	       v2_pt->last = e1;
207 	       v2_pt->degree++;
208 	       e1=e1->next_edge;
209 	       v2_pt->last->next_edge = NULL;
210 	    }
211 	 }
212       }
213       if (!edges_deleted) break;
214    }
215    FREE(new_cut->coef);
216    return(cuts_found);
217 }
218 
219 /*===========================================================================*/
220 
greedy_shrinking1(network * n,double capacity,double etol,int max_num_cuts,cut_data * new_cut,int * compnodes,int * compmembers,int compnum,char * in_set,double * cut_val,int * ref,char * cut_list,double * demand,int mult,int * num_cuts,int * alloc_cuts,cut_data *** cuts)221 int greedy_shrinking1(network *n, double capacity, double etol,
222 		      int max_num_cuts, cut_data *new_cut,
223 		      int *compnodes, int *compmembers, int compnum,
224 		      char *in_set, double *cut_val, int *ref, char *cut_list,
225 		      double *demand, int mult, int *num_cuts, int *alloc_cuts,
226 		      cut_data ***cuts)
227 {
228    double set_weight, set_demand;
229    vertex *verts = n->verts;
230    elist *e;
231    int cuts_found = 0, i, j, k;
232    char *pt, *cutpt;
233    int *ipt;
234    double  *dpt;
235    int vertnum = n->vertnum;
236 
237    int max_vert = 0, set_size, begin = 1, cur_comp, end = 1, other_end;
238    char *coef;
239    double maxval, weight;
240    vertex *cur_nodept;
241 
242    new_cut->size = (vertnum >> DELETE_POWER) + 1;
243    new_cut->coef = coef = (char *) (calloc(new_cut->size, sizeof(char)));
244    memset(cut_list, 0, new_cut->size * (max_num_cuts + 1));
245 
246    *in_set = 0;
247 
248    for (i = 1; i < vertnum;  i++){
249       if (verts[compmembers[i]].deleted) compmembers[i] = 0;
250       ref[compmembers[i]] = i;
251    }
252    *ref = 0;
253    /* ref is a reference array for compmembers: gives a place
254       in which a vertex is listed in  compmembers */
255 
256    for (cur_comp = 1; cur_comp <= compnum;
257 	begin += compnodes[cur_comp], cur_comp++)  /* for every component */
258       for (i = begin, end = begin + compnodes[cur_comp]; i < end; i++){
259 	 if (compmembers[i] == 0) continue;
260 	 /* for every node as a starting one */
261 	 /* initialize the data structures */
262 	 memset(in_set  + begin, 0, compnodes[cur_comp] * sizeof(char));
263 	 memset(cut_val + begin, 0, compnodes[cur_comp] * sizeof(double));
264 	 in_set[i] = 1;
265 	 set_size = 1 + verts[compmembers[i]].orig_node_list_size;
266 	 set_weight = verts[compmembers[i]].weight;
267 	 for (e = verts[compmembers[i]].first; e; e = e->next_edge){
268 	    if (e->other_end)
269 	       cut_val[ref[e->other_end]] = e->data->weight;
270 	 }
271 	 set_demand = demand[compmembers[i]];
272 
273 	 while(TRUE){
274 	    if (set_weight > RHS(set_size, set_demand, capacity) + etol &&
275 		set_size > 2){
276 	       memset(coef, 0, new_cut->size*sizeof(char));
277 	       for (j = begin, ipt = compmembers + begin; j < end; j++, ipt++){
278 		  if (in_set[j]){
279 		     cur_nodept = verts + (*ipt);
280 		     if (cur_nodept->orig_node_list_size)
281 			for (k = 0; k < cur_nodept->orig_node_list_size; k++)
282 			   (coef[(cur_nodept->orig_node_list)[k] >>
283 				DELETE_POWER]) |=
284 			      (1 << ((cur_nodept->orig_node_list)[k] &
285 				     DELETE_AND));
286 		     (coef[(*ipt) >> DELETE_POWER]) |=
287 			(1 << ((*ipt) & DELETE_AND));
288 		  }
289 	       }
290 	       for (k = 0, cutpt = cut_list; k < cuts_found; k++,
291 		       cutpt += new_cut->size)
292 		  if (!memcmp(coef, cutpt, new_cut->size*sizeof(char)))
293 		     break;/* same cuts */
294 	       if (k >= cuts_found){
295 #if 0
296 		  new_cut->type = SUBTOUR_ELIM_SIDE;
297 		  new_cut->rhs = RHS(set_size, set_demand, capacity);
298 		  cuts_found += cg_send_cut(new_cut);
299 #endif
300 #ifdef DIRECTED_X_VARS
301 		  SEND_DIR_SUBTOUR_CONSTRAINT(set_size, set_demand);
302 #else
303 		  SEND_SUBTOUR_CONSTRAINT(set_size, set_demand);
304 #endif
305 		  memcpy(cutpt, coef, new_cut->size);
306 	       }
307 	       if (cuts_found > max_num_cuts){
308 		  FREE(new_cut->coef);
309 		  return(cuts_found);
310 	       }
311 	    }
312 	    for (maxval = -1, pt = in_set + begin, dpt = cut_val + begin,
313 		    j = begin; j < end; pt++, dpt++, j++){
314 	       if (!(*pt) && *dpt > maxval){
315 		  maxval = cut_val[j];
316 		  max_vert = j;
317 	       }
318 	    }
319 	    if (maxval > 0){    /* add the vertex to the set */
320 	       in_set[max_vert]=1;
321 	       set_size += 1+ verts[compmembers[max_vert]].orig_node_list_size;
322 	       set_demand += demand[compmembers[max_vert]];
323 	       set_weight += verts[compmembers[max_vert]].weight;
324 	       cut_val[max_vert] = 0;
325 	       for (e=verts[compmembers[max_vert]].first; e; e = e->next_edge){
326 		  other_end = ref[e->other_end];
327 		  weight = e->data->weight;
328 		  set_weight += (in_set[other_end]) ? weight : -weight;
329 		  cut_val[other_end] += (in_set[other_end]) ? 0 : weight;
330 	       }
331 	    }
332 	    else{ /* can't add anything to the set */
333 	       break;
334 	    }
335 	 }
336       }
337    FREE(new_cut->coef);
338    return(cuts_found);
339 }
340 
341 /*===========================================================================*/
342 
343 #if defined(ADD_FLOW_VARS) && defined(DIRECTED_X_VARS)
greedy_shrinking1_dicut(network * n,double capacity,double etol,int max_num_cuts,cut_data * new_cut,int * compnodes,int * compmembers,int compnum,char * in_set,double * cut_val,int * ref,char * cut_list,double * demand,int mult,int * num_cuts,int * alloc_cuts,cut_data *** cuts)344 int greedy_shrinking1_dicut(network *n, double capacity, double etol,
345 			    int max_num_cuts, cut_data *new_cut,
346 			    int *compnodes, int *compmembers, int compnum,
347 			    char *in_set, double *cut_val, int *ref,
348 			    char *cut_list, double *demand, int mult,
349 			    int *num_cuts, int *alloc_cuts,
350 			    cut_data ***cuts)
351 {
352    double set_demand, set_cut_val, tmp_cut_val;
353    vertex *verts = n->verts;
354    elist *e;
355    int cuts_found = 0, i, j;
356    char *pt;
357    int vertnum = n->vertnum;
358 
359    int min_vert = 0, set_size, begin = 1, cur_comp, end = 1;
360    char *coef;
361    double min_val;
362    int max_size, numarcs, *arcs;
363 
364    max_size = DSIZE + ISIZE + (vertnum >> DELETE_POWER) + 1
365       + 2*vertnum*vertnum*ISIZE;
366    new_cut->coef = (char *) calloc(max_size, sizeof(char));
367    coef = new_cut->coef + DSIZE + ISIZE;
368    arcs = (int *) (coef + (vertnum >> DELETE_POWER) + 1);
369 
370    *in_set = 0;
371 
372    for (i = 1; i < vertnum;  i++){
373       if (verts[compmembers[i]].deleted) compmembers[i] = 0;
374       ref[compmembers[i]] = i;
375    }
376    *ref = 0;
377    /* ref is a reference array for compmembers: gives a place
378       in which a vertex is listed in  compmembers */
379 
380    for (cur_comp = 1; cur_comp <= compnum; begin += compnodes[cur_comp],
381 	   cur_comp++){  /* for every component */
382       for (i = begin, end = begin + compnodes[cur_comp]; i < end; i++){
383 	 if (compmembers[i] == 0) continue;
384 	 /* for every node as a starting one */
385 	 /* initialize the data structures */
386 	 memset(in_set  + begin, 0, compnodes[cur_comp] * sizeof(char));
387 	 in_set[i] = 1;
388 	 set_size = 1 + verts[compmembers[i]].orig_node_list_size;
389 	 set_demand = demand[compmembers[i]];
390 	 set_cut_val = 0;
391 	 for (e = verts[compmembers[i]].first; e; e = e->next_edge){
392 	    if (e->other_end < compmembers[i]){
393 	       set_cut_val += MIN(e->data->flow1,
394 				  MIN(set_demand, capacity)*e->data->weight1);
395 	    }else{
396 	       set_cut_val += MIN(e->data->flow2,
397 				  MIN(set_demand, capacity)*e->data->weight2);
398 	    }
399 	 }
400 
401 	 while(TRUE){
402 	    if (set_cut_val + etol < set_demand){
403 	       memset(coef, 0, ((vertnum >> DELETE_POWER) + 1)*sizeof(char));
404 	       numarcs = 0;
405 	       for (j = begin; j < end; j++){
406 		  if (in_set[j]){
407 		     (coef[(compmembers[j]) >> DELETE_POWER]) |=
408 			(1 << ((compmembers[j]) & DELETE_AND));
409 		     for (e = verts[compmembers[j]].first; e; e=e->next_edge){
410 			if (e->other_end < compmembers[j]){
411 			   if (!in_set[ref[e->other_end]] &&
412 			       set_demand*e->data->weight1 >
413 			       e->data->flow1+etol){
414 			      arcs[numarcs << 1] = e->other_end;
415 			      arcs[(numarcs << 1) + 1] = compmembers[j];
416 			      numarcs++;
417 			   }
418 			}else{
419 			   if (!in_set[ref[e->other_end]] &&
420 			       set_demand*e->data->weight2 >
421 			       e->data->flow2+etol){
422 			      arcs[numarcs << 1] = e->other_end;
423 			      arcs[(numarcs << 1) + 1] = compmembers[j];
424 			      numarcs++;
425 			   }
426 			}
427 		     }
428 		  }
429 	       }
430 	       ((double *)(new_cut->coef))[0] = set_demand;
431 	       ((int *)(new_cut->coef + DSIZE))[0] = numarcs;
432 	       new_cut->size = DSIZE + ISIZE + (vertnum >> DELETE_POWER)
433 		  + 1 + 2 * numarcs * ISIZE;
434 	       new_cut->type = MIXED_DICUT;
435 	       new_cut->rhs = set_demand;
436 	       new_cut->name  = CUT__SEND_TO_CP;
437 	       cuts_found += cg_send_cut(new_cut, num_cuts, alloc_cuts, cuts);
438 	       if (cuts_found > max_num_cuts){
439 		  FREE(new_cut->coef);
440 		  return(cuts_found);
441 	       }
442 	    }
443 	    for (min_val = 0, pt = in_set + begin, j = begin; j < end;
444 		 pt++, j++){
445 	       if (!in_set[j]){
446 		  for (tmp_cut_val = 0, e = verts[compmembers[j]].first;
447 		       e; e = e->next_edge){
448 		     if (in_set[ref[e->other_end]]){
449 			if (compmembers[j] < e->other_end){
450 			   tmp_cut_val -=
451 			      MIN(e->data->flow1,
452 				  MIN(set_demand, capacity)*e->data->weight1);
453 			}else{
454 			   tmp_cut_val -=
455 			      MIN(e->data->flow2,
456 				  MIN(set_demand, capacity)*e->data->weight2);
457 			}
458 		     }else{
459 			if (compmembers[j] < e->other_end){
460 			   tmp_cut_val +=
461 			      MIN(e->data->flow2,
462 				  MIN(set_demand, capacity)*e->data->weight2);
463 			}else{
464 			   tmp_cut_val +=
465 			      MIN(e->data->flow1,
466 				  MIN(set_demand, capacity)*e->data->weight1);
467 			}
468 		     }
469 		  }
470 		  if (tmp_cut_val < min_val - etol){
471 		     min_vert = j;
472 		     min_val = tmp_cut_val;
473 		  }
474 	       }
475 	    }
476 	    if (min_val < 0){    /* add the vertex to the set */
477 	       in_set[min_vert] = 1;
478 	       set_size += 1 +
479 		  verts[compmembers[min_vert]].orig_node_list_size;
480 	       set_demand += demand[compmembers[min_vert]];
481 	       set_cut_val += min_val;
482 	    }
483 	    else{ /* can't add anything to the set */
484 	       break;
485 	    }
486 	 }
487       }
488    }
489    FREE(new_cut->coef);
490    return(cuts_found);
491 }
492 #endif
493 
494 /*===========================================================================*/
495 
greedy_shrinking6(network * n,double capacity,double etol,cut_data * new_cut,int * compnodes,int * compmembers,int compnum,char * in_set,double * cut_val,int * ref,char * cut_list,int max_num_cuts,double * demand,int trial_num,double prob,int mult,int * num_cuts,int * alloc_cuts,cut_data *** cuts)496 int greedy_shrinking6(network *n, double capacity, double etol,
497 		      cut_data *new_cut, int *compnodes,
498 		      int *compmembers, int compnum,char *in_set,
499 		      double *cut_val, int *ref, char *cut_list,
500 		      int max_num_cuts, double *demand, int trial_num,
501 		      double prob, int mult, int *num_cuts, int *alloc_cuts,
502 		      cut_data ***cuts)
503 {
504    double set_weight, set_demand;
505    vertex  *verts = n->verts;
506    elist *e;
507    int i, j, k, cuts_found = 0;
508    char *pt, *cutpt;
509    double *dpt;
510    int vertnum = n->vertnum;
511 
512    int max_vert = 0, set_size, begin = 1, cur_comp, end = 1, num_trials;
513    char *coef;
514    double maxval;
515    double denominator=pow(2.0,31.0)-1.0;
516    double r, q;
517 
518    int other_end;
519    double weight;
520    int *ipt;
521    vertex *cur_nodept;
522 
523    new_cut->size = (vertnum >> DELETE_POWER) + 1;
524    new_cut->coef =coef= (char *) (calloc(new_cut->size,sizeof(char)));
525    memset(cut_list, 0, new_cut->size * (max_num_cuts +1));
526 
527 
528    *in_set=0;
529 
530    for(i = 1; i < vertnum; i++){
531       if (verts[compmembers[i]].deleted) compmembers[i] = 0;
532       ref[compmembers[i]] = i;
533    }
534    *ref = 0;
535 
536    /* ref is a reference array for compmembers: gives a place
537       in which a vertex is listed in  compmembers */
538 
539    for (cur_comp = 1; cur_comp <= compnum; begin += compnodes[cur_comp],
540 	   cur_comp++){
541       /* for every component */
542       if (compnodes[cur_comp] <= 7) continue;
543 
544       for (num_trials = 0; num_trials < trial_num * compnodes[cur_comp];
545 	num_trials++){
546 	 end = begin + compnodes[cur_comp];
547 	 /*initialize the data structures */
548 	 memset(in_set + begin, 0, compnodes[cur_comp] * sizeof(char));
549 	 memset(cut_val+ begin, 0, compnodes[cur_comp] * sizeof(double));
550 	 set_size = 0;
551 	 set_demand = 0;
552 	 set_weight = 0;
553          for (i = begin; i < end; i++ ){
554 	    if (compmembers[i] == 0) continue;
555 /*__BEGIN_EXPERIMENTAL_SECTION__*/
556 #if 0
557 	    r  = CCutil_lprand(&rand_state)/CC_PRANDMAX;
558 #endif
559 /*___END_EXPERIMENTAL_SECTION___*/
560 	    r = RANDOM()/denominator;
561 	    q = (prob/compnodes[cur_comp]);
562 	    if (r < q){
563 	       in_set[i] = 1;
564 	       set_size += 1 + verts[compmembers[i]].orig_node_list_size;
565 	       set_demand += demand[compmembers[i]];
566 	       set_weight += verts[compmembers[i]].weight;
567 	       for (e = verts[compmembers[i]].first; e; e = e-> next_edge){
568 		  other_end = ref[e->other_end];
569 		  weight = e->data->weight;
570 		  set_weight += (in_set[other_end]) ? weight : -weight;
571 		  cut_val[other_end] += (in_set[other_end]) ? 0 : weight;
572 	       }
573 	    }
574 	 }
575 	 while(set_size){
576 	    if (set_weight > RHS(set_size, set_demand, capacity) + etol &&
577 		set_size > 2){
578 	       memset(coef, 0, new_cut->size*sizeof(char));
579 	       for (j = begin, ipt = compmembers + begin; j < end; j++, ipt++){
580 		  if (in_set[j]){
581 		     cur_nodept = verts + (*ipt);
582 		     if (cur_nodept->orig_node_list_size)
583 			for (k = 0; k < cur_nodept->orig_node_list_size; k++)
584 			   (coef[(cur_nodept->orig_node_list)[k] >>
585 				DELETE_POWER]) |=
586 			      (1 << ((cur_nodept->orig_node_list)[k] &
587 				     DELETE_AND));
588 		     (coef[(*ipt) >> DELETE_POWER]) |= (1 << ((*ipt) &
589 							      DELETE_AND));
590 		  }
591 	       }
592 	       for (k = 0, cutpt = cut_list; k < cuts_found; k++,
593 		       cutpt += new_cut->size)
594 		  if (!memcmp(coef, cutpt, new_cut->size*sizeof(char))) break;
595 	       if ( k >= cuts_found){
596 #if 0
597 		  new_cut->type = SUBTOUR_ELIM_SIDE;
598 		  new_cut->rhs =  RHS(set_size, set_demand, capacity);
599 		  cuts_found += cg_send_cut(new_cut);
600 #endif
601 #ifdef DIRECTED_X_VARS
602 		  SEND_DIR_SUBTOUR_CONSTRAINT(set_size, set_demand);
603 #else
604 		  SEND_SUBTOUR_CONSTRAINT(set_size, set_demand);
605 #endif
606 		  memcpy(cutpt, coef, new_cut->size);
607 	       }
608 
609 	       if ( cuts_found > max_num_cuts){
610 		  FREE(new_cut->coef);
611 		  return(cuts_found);
612 	       }
613 	    }
614 	    for (maxval = -1, pt = in_set+begin, dpt = cut_val+begin,
615 		    j = begin; j < end; pt++, dpt++, j++){
616 	       if (!(*pt) && *dpt > maxval){
617 		  maxval = cut_val[j];
618 		  max_vert = j;
619 	       }
620 	    }
621 	    if (maxval > 0){    /* add the vertex to the set */
622 	       in_set[max_vert]=1;
623 	       set_size+=1+ verts[compmembers[max_vert]].orig_node_list_size;
624 	       set_demand += demand[compmembers[max_vert]];
625 	       set_weight += verts[compmembers[max_vert]].weight;
626 	       cut_val[max_vert]=0;
627 	       for (e = verts[compmembers[max_vert]].first; e;
628 		    e = e->next_edge){
629 		  other_end = ref[e->other_end];
630 		  weight = e->data->weight;
631 		  set_weight += (in_set[other_end]) ? weight : -weight;
632 		  cut_val[other_end]+=(in_set[other_end]) ? 0 : weight;
633 	       }
634 	    }
635 	    else{ /* can't add anything to the set */
636 	       break;
637 	    }
638 	 }
639       }
640    }
641 
642    FREE(new_cut->coef);
643    return(cuts_found);
644 }
645 
646 /*===========================================================================*/
647 
648 #if defined(ADD_FLOW_VARS) && defined(DIRECTED_X_VARS)
greedy_shrinking6_dicut(network * n,double capacity,double etol,cut_data * new_cut,int * compnodes,int * compmembers,int compnum,char * in_set,double * cut_val,int * ref,char * cut_list,int max_num_cuts,double * demand,int trial_num,double prob,int mult,int * num_cuts,int * alloc_cuts,cut_data *** cuts)649 int greedy_shrinking6_dicut(network *n, double capacity, double etol,
650 			    cut_data *new_cut, int *compnodes,
651 			    int *compmembers, int compnum,char *in_set,
652 			    double *cut_val, int *ref, char *cut_list,
653 			    int max_num_cuts, double *demand, int trial_num,
654 			    double prob, int mult, int *num_cuts,
655 			    int *alloc_cuts, cut_data ***cuts)
656 {
657    double set_demand, set_cut_val, tmp_cut_val;
658    vertex  *verts = n->verts;
659    elist *e;
660    int cuts_found = 0, i, j;
661    char *pt, *cutpt;
662    double *dpt;
663    int vertnum = n->vertnum;
664 
665    int min_vert = 0, set_size, begin = 1, cur_comp, end = 1, num_trials;
666    char *coef;
667    double min_val;
668    int max_size, numarcs, *arcs;
669    double denominator=pow(2.0,31.0)-1.0;
670    double r, q;
671 
672    max_size = DSIZE + ISIZE + (vertnum >> DELETE_POWER) + 1
673       + 2*vertnum*vertnum*ISIZE;
674    new_cut->coef = (char *) calloc(max_size, sizeof(char));
675    coef = new_cut->coef + DSIZE + ISIZE;
676    arcs = (int *) (coef + (vertnum >> DELETE_POWER) + 1);
677 
678    *in_set=0;
679 
680    for(i = 1; i < vertnum; i++){
681       if (verts[compmembers[i]].deleted) compmembers[i] = 0;
682       ref[compmembers[i]] = i;
683    }
684    *ref = 0;
685    /* ref is a reference array for compmembers: gives a place
686       in which a vertex is listed in  compmembers */
687 
688    for (cur_comp = 1; cur_comp <= compnum; begin += compnodes[cur_comp],
689 	   cur_comp++){
690       /* for every component */
691       if (compnodes[cur_comp] <= 7) continue;
692 
693       for (num_trials = 0; num_trials < trial_num * compnodes[cur_comp];
694 	num_trials++){
695 	 end = begin + compnodes[cur_comp];
696 	 /*initialize the data structures */
697 	 memset(in_set + begin, 0, compnodes[cur_comp] * sizeof(char));
698 	 set_size = 0;
699 	 set_demand = 0;
700 	 set_cut_val = 0;
701          for (i = begin; i < end; i++ ){
702 	    if (compmembers[i] == 0) continue;
703 /*__BEGIN_EXPERIMENTAL_SECTION__*/
704 #if 0
705 	    r  = CCutil_lprand(&rand_state)/CC_PRANDMAX;
706 #endif
707 /*___END_EXPERIMENTAL_SECTION___*/
708 	    r = RANDOM()/denominator;
709 	    q = (prob/compnodes[cur_comp]);
710 	    if (r < q){
711 	       in_set[i] = 1;
712 	       set_size += 1 + verts[compmembers[i]].orig_node_list_size;
713 	       set_demand += demand[compmembers[i]];
714 	       for (e = verts[compmembers[i]].first; e; e = e-> next_edge){
715 		  if (in_set[ref[e->other_end]]){
716 		     if (compmembers[i] < e->other_end){
717 			set_cut_val -=
718 			   MIN(e->data->flow1,
719 			       MIN(set_demand, capacity)*e->data->weight1);
720 		     }else{
721 			set_cut_val -=
722 			   MIN(e->data->flow2,
723 			       MIN(set_demand, capacity)*e->data->weight2);
724 		     }
725 		  }else{
726 		     if (compmembers[i] < e->other_end){
727 			set_cut_val +=
728 			   MIN(e->data->flow2,
729 			       MIN(set_demand, capacity)*e->data->weight2);
730 		     }else{
731 			set_cut_val +=
732 			   MIN(e->data->flow1,
733 			       MIN(set_demand, capacity)*e->data->weight1);
734 		     }
735 		  }
736 	       }
737 	    }
738 	 }
739 	 while(TRUE){
740 	    if (set_cut_val + etol < set_demand){
741 	       memset(coef, 0, ((vertnum >> DELETE_POWER) + 1)*sizeof(char));
742 	       numarcs = 0;
743 	       for (j = begin; j < end; j++){
744 		  if (in_set[j]){
745 		     (coef[(compmembers[j]) >> DELETE_POWER]) |=
746 			(1 << ((compmembers[j]) & DELETE_AND));
747 		     for (e = verts[compmembers[j]].first; e; e=e->next_edge){
748 			if (e->other_end < compmembers[j]){
749 			   if (!in_set[ref[e->other_end]] &&
750 			       set_demand*e->data->weight1 >
751 			       e->data->flow1+etol){
752 			      arcs[numarcs << 1] = e->other_end;
753 			      arcs[(numarcs << 1) + 1] = compmembers[j];
754 			      numarcs++;
755 			   }
756 			}else{
757 			   if (!in_set[ref[e->other_end]] &&
758 			       set_demand*e->data->weight2 >
759 			       e->data->flow2+etol){
760 			      arcs[numarcs << 1] = e->other_end;
761 			      arcs[(numarcs << 1) + 1] = compmembers[j];
762 			      numarcs++;
763 			   }
764 			}
765 		     }
766 		  }
767 	       }
768 	       ((double *)(new_cut->coef))[0] = set_demand;
769 	       ((int *)(new_cut->coef + DSIZE))[0] = numarcs;
770 	       new_cut->size = DSIZE + ISIZE + (vertnum >> DELETE_POWER)
771 		  + 1 + 2 * numarcs * ISIZE;
772 	       new_cut->type = MIXED_DICUT;
773 	       new_cut->rhs = set_demand;
774 	       new_cut->name  = CUT__SEND_TO_CP;
775 	       cuts_found += cg_send_cut(new_cut, num_cuts, alloc_cuts, cuts);
776 	       if (cuts_found > max_num_cuts){
777 		  FREE(new_cut->coef);
778 		  return(cuts_found);
779 	       }
780 	    }
781 	    for (min_val = 0, pt = in_set + begin, j = begin; j < end;
782 		 pt++, j++){
783 	       if (!in_set[j]){
784 		  for (tmp_cut_val = 0, e = verts[compmembers[j]].first;
785 		       e; e = e->next_edge){
786 		     if (in_set[ref[e->other_end]]){
787 			if (compmembers[j] < e->other_end){
788 			   tmp_cut_val -=
789 			      MIN(e->data->flow1,
790 				  MIN(set_demand, capacity)*e->data->weight1);
791 			}else{
792 			   tmp_cut_val -=
793 			      MIN(e->data->flow2,
794 				  MIN(set_demand, capacity)*e->data->weight2);
795 			}
796 		     }else{
797 			if (compmembers[j] < e->other_end){
798 			   tmp_cut_val +=
799 			      MIN(e->data->flow2,
800 				  MIN(set_demand, capacity)*e->data->weight2);
801 			}else{
802 			   tmp_cut_val +=
803 			      MIN(e->data->flow1,
804 				  MIN(set_demand, capacity)*e->data->weight1);
805 			}
806 		     }
807 		  }
808 		  if (tmp_cut_val < min_val - etol){
809 		     min_vert = j;
810 		     min_val = tmp_cut_val;
811 		  }
812 	       }
813 	    }
814 	    if (min_val < 0){    /* add the vertex to the set */
815 	       in_set[min_vert] = 1;
816 	       set_size += 1 +
817 		  verts[compmembers[min_vert]].orig_node_list_size;
818 	       set_demand += demand[compmembers[min_vert]];
819 	       set_cut_val += min_val;
820 	    }
821 	    else{ /* can't add anything to the set */
822 	       break;
823 	    }
824 	 }
825       }
826    }
827    FREE(new_cut->coef);
828    return(cuts_found);
829 }
830 #endif
831 
832 /*===========================================================================*/
833 
greedy_shrinking1_one(network * n,double capacity,double etol,int max_num_cuts,cut_data * new_cut,char * in_set,double * cut_val,char * cut_list,int num_routes,double * demand,int mult,int * num_cuts,int * alloc_cuts,cut_data *** cuts)834 int greedy_shrinking1_one(network *n, double capacity, double etol,
835 			  int max_num_cuts, cut_data *new_cut,char *in_set,
836 			  double *cut_val, char *cut_list, int num_routes,
837 			  double *demand, int mult, int *num_cuts,
838 			  int *alloc_cuts, cut_data ***cuts)
839 {
840 
841    double set_weight, set_cut_val, set_demand;
842    vertex  *verts = n->verts;
843    elist *e;
844    int i, j, k, cuts_found = 0;
845    char *pt, *cutpt;
846    double  *dpt;
847    int vertnum = n->vertnum;
848    int max_vert = 0;
849    int set_size;
850    /* int flag=0; */
851 
852    double complement_demand, total_demand = verts[0].demand;
853    double complement_cut_val;
854    int complement_size;
855    char *coef;
856    double maxval;
857    int other_end;
858    double weight;
859    vertex *cur_nodept;
860 
861    new_cut->size = (vertnum >> DELETE_POWER) + 1;
862    new_cut->coef = coef = (char *) (calloc(new_cut->size,sizeof(char)));
863    memset(cut_list, 0, new_cut->size * (max_num_cuts + 1));
864 
865    for (i = 1; i < vertnum; i++ ){
866       if (verts[i].deleted) continue;/* for every node as a starting one */
867       /*initialize the data structures */
868       memset(in_set, 0, vertnum*sizeof(char));
869       memset(cut_val, 0,vertnum* sizeof(double));
870       in_set[i] = 1;
871       set_size = 1 + verts[i].orig_node_list_size;
872       set_cut_val = 0;
873       set_weight = verts[i].weight;
874       for (e= verts[i].first; e; e = e-> next_edge){
875 	 cut_val[e->other_end] = e->data->weight;
876 	 set_cut_val += e->data->weight;
877       }
878       set_demand = demand[i];
879 
880       while(TRUE){
881 	 if (set_weight > RHS(set_size, set_demand, capacity) + etol &&
882 	     set_size > 2){
883 	    memset(coef, 0, new_cut->size*sizeof(char));
884 	    /* printf("%d :", i); */
885 	    /*  printf("%d ", j); */
886 	    for (j = 1; j < vertnum; j++)
887 	       if (in_set[j]){
888 		  cur_nodept = verts + j;
889 		  if (cur_nodept->orig_node_list_size)
890 		     for (k = 0; k < cur_nodept->orig_node_list_size; k++)
891 			(coef[(cur_nodept->orig_node_list)[k] >>
892 			     DELETE_POWER]) |=
893 			   (1 << ((cur_nodept->orig_node_list)[k] &
894 				  DELETE_AND));
895 		  (coef[j>> DELETE_POWER]) |= (1 << ( j & DELETE_AND));
896 	       }
897 	    /*  printf("%f ", set_demand);
898 	    printf("%f \n",set_cut_val);*/
899 	    for (k = 0, cutpt = cut_list; k < cuts_found; k++,
900 		    cutpt += new_cut->size)
901 		  if (!memcmp(coef, cutpt, new_cut->size*sizeof(char)))
902 		     break; /* same cuts */
903 	    if ( k >= cuts_found){
904 #if 0
905 	       new_cut->type = SUBTOUR_ELIM_SIDE;
906 	       new_cut->rhs =  RHS(set_size, set_demand, capacity);
907 	       cuts_found += cg_send_cut(new_cut);
908 #endif
909 #ifdef DIRECTED_X_VARS
910 	       SEND_DIR_SUBTOUR_CONSTRAINT(set_size, set_demand);
911 #else
912 	       SEND_SUBTOUR_CONSTRAINT(set_size, set_demand);
913 #endif
914 	       memcpy(cutpt, coef, new_cut->size);
915 	    }
916 
917 	    if ( cuts_found > max_num_cuts){
918 	       FREE(new_cut->coef);
919 	       return(cuts_found);
920 	    }
921 	 }
922 	 /* check the complement */
923 
924 	 complement_demand = total_demand - set_demand;
925 	 complement_cut_val = set_cut_val- 2*(*cut_val) + 2*num_routes;
926 	 complement_size = vertnum - 1 - set_size;
927 	 if (complement_cut_val< mult*(ceil(complement_demand/capacity))-etol
928 	     && complement_size > 2){
929 	    memset(coef, 0, new_cut->size*sizeof(char));
930 	    for (j = 1; j < vertnum; j++){
931 	       if (!(in_set[j]) && !(verts[j].deleted)){
932 		  cur_nodept = verts + j;
933 		  if (cur_nodept->orig_node_list_size)
934 		     for (k = 0; k < cur_nodept->orig_node_list_size; k++)
935 			(coef[(cur_nodept->orig_node_list)[k] >>
936 			     DELETE_POWER]) |=
937 			   (1 << ((cur_nodept->orig_node_list)[k] &
938 				  DELETE_AND));
939 		  (coef[j>> DELETE_POWER]) |= (1 << ( j & DELETE_AND));
940 	       }
941 	    }
942 	    for (k=0, cutpt = cut_list; k < cuts_found; k++,
943 		    cutpt += new_cut->size)
944 		  if (!memcmp(coef, cutpt, new_cut->size*sizeof(char))) break;
945 	    if ( k >= cuts_found){
946 #if 0
947 	       new_cut->type = SUBTOUR_ELIM_SIDE;
948 	       new_cut->rhs =  RHS(complement_size, complement_demand,
949 				   capacity);
950 	       cuts_found += cg_send_cut(new_cut);
951 #endif
952 #ifdef DIRECTED_X_VARS
953 	       SEND_DIR_SUBTOUR_CONSTRAINT(complement_size, complement_demand);
954 #else
955 	       SEND_SUBTOUR_CONSTRAINT(complement_size, complement_demand);
956 #endif
957 	       memcpy(cutpt, coef, new_cut->size);
958 	    }
959 
960 	    if (cuts_found > max_num_cuts){
961 	       FREE(new_cut->coef);
962 	       return(cuts_found);
963 	    }
964 	 }
965 
966 	 for (maxval = -1, pt = in_set, dpt = cut_val,pt++, dpt++,
967 		 j = 1; j < vertnum; pt++, dpt++, j++){
968 	    if (!(*pt) && *dpt > maxval){
969 	       maxval = cut_val[j];
970 	       max_vert = j;
971 	    }
972 	 }
973 	 if (maxval > 0){    /* add the vertex to the set */
974 	    in_set[max_vert] = 1;
975 	    set_size += 1 + verts[max_vert].orig_node_list_size ;
976 	    set_demand += demand[max_vert];
977 	    set_weight += verts[max_vert].weight;
978 	    cut_val[max_vert] = 0;
979 	    for (e = verts[max_vert].first; e; e = e-> next_edge){
980 	       other_end = e->other_end;
981 	       weight = e->data->weight;
982 	       set_weight += (in_set[other_end]) ? weight : -weight;
983 	       set_cut_val += (in_set[other_end]) ? (-weight) : weight;
984 	       cut_val[other_end] += weight;
985 	    }
986 	 }
987 	 else{ /* can't add anything to the set */
988 	    break;
989 	 }
990       }
991    }
992    FREE(new_cut->coef);
993    return(cuts_found);
994 }
995 
996 /*===========================================================================*/
997 
greedy_shrinking6_one(network * n,double capacity,double etol,cut_data * new_cut,char * in_set,double * cut_val,int num_routes,char * cut_list,int max_num_cuts,double * demand,int trial_num,double prob,int mult,int * num_cuts,int * alloc_cuts,cut_data *** cuts)998 int greedy_shrinking6_one(network *n, double capacity,
999 			  double etol, cut_data *new_cut,
1000 			  char *in_set, double *cut_val, int num_routes,
1001 			  char *cut_list, int max_num_cuts, double *demand,
1002 			  int trial_num, double prob, int mult, int *num_cuts,
1003 			  int *alloc_cuts, cut_data ***cuts)
1004 {
1005 
1006    double set_weight, set_cut_val, set_demand;
1007    vertex  *verts=n->verts;
1008    elist *e;
1009    int i, j, k, cuts_found = 0;
1010    char *pt, *cutpt;
1011    double  *dpt;
1012    int vertnum = n->vertnum;
1013 
1014    int max_vert = 0, set_size, begin = 1, end = 1, num_trials;
1015    char *coef;
1016    double maxval, r, q;
1017    double denominator = pow(2.0, 31.0) - 1.0;
1018 
1019    int other_end;
1020    double weight;
1021 
1022    double complement_demand, total_demand = verts[0].demand;
1023    double complement_cut_val;
1024    int complement_size;
1025    vertex *cur_nodept;
1026    /* int flag=0;*/
1027 
1028    new_cut->size = (vertnum >> DELETE_POWER) + 1;
1029    new_cut->coef = coef = (char *) (calloc(new_cut->size,sizeof(char)));
1030      memset(cut_list, 0, new_cut->size * (max_num_cuts +1));
1031 
1032    *in_set = 0;
1033 
1034    for (num_trials = 0; num_trials < trial_num*vertnum ; num_trials++){
1035 
1036       /*initialize the data structures */
1037       memset(in_set, 0, vertnum*sizeof(char));
1038       memset(cut_val, 0,vertnum* sizeof(double));
1039 
1040       set_cut_val = 0;
1041       set_size = 0;
1042       set_demand = 0;
1043       set_weight = 0;
1044       for (i = 1 ; i < vertnum; i++ ){
1045 	 if (verts[i].deleted) continue;
1046 /*__BEGIN_EXPERIMENTAL_SECTION__*/
1047 #if 0
1048 	 r  = CCutil_lprand(&rand_state)/CC_PRANDMAX;
1049 #endif
1050 /*___END_EXPERIMENTAL_SECTION___*/
1051 	 r = RANDOM()/denominator;
1052 	 q = (prob/vertnum);
1053 	 if (r < q){
1054 	    in_set[i] = 1;
1055 	    set_size += 1 + verts[i].orig_node_list_size;
1056 	    set_demand += demand[i];
1057 	    set_weight += verts[i].weight;
1058 	    for (e = verts[i].first; e; e = e-> next_edge){
1059 		other_end = e->other_end;
1060 		weight  = e->data->weight;
1061 		set_weight += (in_set[other_end]) ? weight : -weight;
1062 		set_cut_val += (in_set[other_end]) ? (-weight) : weight;
1063 		cut_val[other_end] += (in_set[other_end]) ? 0 : weight;
1064 	    }
1065 	 }
1066       }
1067       while(set_size){
1068 	 if (set_weight > RHS(set_size, set_demand, capacity) + etol &&
1069 	     set_size > 2){
1070 	    memset(coef, 0, new_cut->size*sizeof(char));
1071 	    for (j = 1; j < vertnum; j++ ){
1072 	       if (in_set[j]){
1073 		  cur_nodept = verts + j;
1074 		  if (cur_nodept->orig_node_list_size)
1075 		     for (k = 0; k < cur_nodept->orig_node_list_size; k++)
1076 			(coef[(cur_nodept->orig_node_list)[k] >>
1077 			     DELETE_POWER]) |=
1078 			   (1 << ((cur_nodept->orig_node_list)[k] &
1079 				  DELETE_AND));
1080 		  (coef[j>> DELETE_POWER]) |= (1 << ( j & DELETE_AND));
1081 	       }
1082 	    }
1083 	    for (k = 0, cutpt = cut_list; k < cuts_found; k++,
1084 		    cutpt += new_cut->size)
1085 	       if (!memcmp(coef, cutpt, new_cut->size*sizeof(char))) break;
1086 	    if ( k >= cuts_found){
1087 #if 0
1088 	       new_cut->type = SUBTOUR_ELIM_SIDE;
1089 	       new_cut->rhs =  RHS(set_size, set_demand,
1090 				   capacity);
1091 	       cuts_found += cg_send_cut(new_cut);
1092 #endif
1093 #ifdef DIRECTED_X_VARS
1094 	       SEND_DIR_SUBTOUR_CONSTRAINT(set_size, set_demand);
1095 #else
1096 	       SEND_SUBTOUR_CONSTRAINT(set_size, set_demand);
1097 #endif
1098 	       memcpy(cutpt, coef, new_cut->size);
1099 	    }
1100 	    if (cuts_found > max_num_cuts){
1101 	       FREE(new_cut->coef);
1102 	       return(cuts_found);
1103 	    }
1104 	 }
1105 
1106 	 /* check the complement */
1107 
1108 	 complement_demand = total_demand - set_demand;
1109 	 complement_cut_val = set_cut_val - 2*(*cut_val) + 2*num_routes;
1110 	 complement_size = vertnum - 1 - set_size;
1111 	 if (complement_cut_val< mult*(ceil(complement_demand/capacity))-
1112 	     etol && complement_size > 2){
1113 	    memset(coef, 0, new_cut->size*sizeof(char));
1114 	    for (j = 1; j < vertnum; j++){
1115 	       if (!(in_set[j])&& !(verts[j].deleted)){
1116 		  cur_nodept = verts + j;
1117 		  if (cur_nodept->orig_node_list_size)
1118 		     for (k = 0; k < cur_nodept->orig_node_list_size; k++)
1119 			(coef[(cur_nodept->orig_node_list)[k] >>
1120 			     DELETE_POWER]) |=
1121 			   (1 << ((cur_nodept->orig_node_list)[k] &
1122 				  DELETE_AND));
1123 		  (coef[j>> DELETE_POWER]) |= (1 << ( j & DELETE_AND));
1124 	       }
1125 	    }
1126 	    for (k = 0, cutpt = cut_list; k < cuts_found; k++,
1127 		    cutpt += new_cut->size)
1128 	       if (!memcmp(coef, cutpt, new_cut->size*sizeof(char))) break;
1129 	    if ( k >= cuts_found){
1130 #if 0
1131 	       new_cut->type = SUBTOUR_ELIM_SIDE;
1132 	       new_cut->rhs =  RHS(complement_size, complement_demand,
1133 				   capacity);
1134 	       cuts_found += cg_send_cut(new_cut);
1135 #endif
1136 #ifdef DIRECTED_X_VARS
1137 	       SEND_DIR_SUBTOUR_CONSTRAINT(complement_size, complement_demand);
1138 #else
1139 	       SEND_SUBTOUR_CONSTRAINT(complement_size, complement_demand);
1140 #endif
1141 	       memcpy(cutpt, coef, new_cut->size);
1142 	    }
1143 
1144 	    if (cuts_found > max_num_cuts){
1145 	       FREE(new_cut->coef);
1146 	       return(cuts_found);
1147 	    }
1148 	 }
1149 
1150 	 for (maxval = -1, pt = in_set + begin, dpt = cut_val + begin,
1151 		 j = begin; j < end; pt++, dpt++, j++){
1152 	    if (!(*pt) && *dpt > maxval){
1153 	       maxval = cut_val[j];
1154 		  max_vert = j;
1155 	    }
1156 	 }
1157 	 if (maxval > 0){    /* add the vertex to the set */
1158 	    in_set[max_vert] = 1;
1159 	    set_size += 1 + verts[max_vert].orig_node_list_size ;
1160 	    set_demand += demand[max_vert];
1161 	    set_weight += verts[max_vert].weight;
1162 	    cut_val[max_vert] = 0;
1163 	    for (e = verts[max_vert].first; e; e = e-> next_edge){
1164 	       other_end = e->other_end;
1165 	       weight  = e->data->weight;
1166 	       set_weight += (in_set[other_end]) ? weight : -weight;
1167 	       set_cut_val += (in_set[other_end]) ? (-weight) : weight;
1168 	       cut_val[other_end] += weight;
1169 	    }
1170 	 }
1171 	 else{ /* can't add anything to the set */
1172 	    break;
1173 	 }
1174       }
1175    }
1176 
1177    FREE(new_cut->coef);
1178    return(cuts_found);
1179 }
1180 
1181 /*===========================================================================*/
1182 
greedy_shrinking2_one(network * n,double capacity,double etol,cut_data * new_cut,char * in_set,double * cut_val,int num_routes,double * demand,int mult,int * num_cuts,int * alloc_cuts,cut_data *** cuts)1183 int greedy_shrinking2_one(network *n, double capacity,
1184 			  double etol, cut_data *new_cut,
1185 			  char *in_set, double *cut_val, int num_routes,
1186 			  double *demand, int mult, int *num_cuts,
1187 			  int *alloc_cuts, cut_data ***cuts)
1188 {
1189 
1190    double set_cut_val, set_demand;
1191    vertex *verts = n->verts;
1192    elist *e, *cur_edge1, *cur_edge2;
1193    int j, k, cuts_found = 0;
1194    char *pt;
1195    double  *dpt;
1196    int vertnum = n->vertnum;
1197 
1198    int max_vert = 0, set_size, begin = 1, end = 1;
1199    char *coef;
1200    double maxval;
1201 
1202    int other_end;
1203    double weight;
1204 
1205    double complement_demand, total_demand = verts[0].demand;
1206    double complement_cut_val;
1207    int complement_size;
1208    vertex *cur_nodept;
1209 
1210    new_cut->size = (vertnum >> DELETE_POWER) + 1;
1211    new_cut->coef =coef= (char *) (calloc(new_cut->size,sizeof(char)));
1212 
1213    *in_set=0;
1214 
1215    for (cur_edge1 = verts[0].first; cur_edge1;
1216 	cur_edge1 = cur_edge1->next_edge){
1217       for (cur_edge2 = cur_edge1->next_edge; cur_edge2;
1218 	   cur_edge2 = cur_edge2->next_edge){
1219 
1220 	 /*initialize the data structures */
1221 	 memset(in_set, 0, vertnum*sizeof(char));
1222 	 memset(cut_val, 0,vertnum* sizeof(double));
1223 
1224 	 set_cut_val = 2;
1225 	 set_size = 2 + cur_edge1->other->orig_node_list_size +
1226 	    cur_edge2->other->orig_node_list_size;
1227 	 set_demand = demand[cur_edge1->other_end] +
1228 	    demand[cur_edge2->other_end];
1229 	 in_set[cur_edge1->other_end] = 1;
1230 
1231 	 for (e = verts[cur_edge1->other_end].first; e; e = e-> next_edge){
1232 	    cut_val[e->other_end] += e->data->weight;
1233 	 }
1234 
1235 	 in_set[cur_edge2->other_end] = 1;
1236 	 for (e = verts[cur_edge2->other_end].first; e; e = e-> next_edge){
1237 	    other_end = e->other_end;
1238 	    weight = e->data->weight;
1239 	    set_cut_val += (in_set[other_end]) ? (-weight) : weight;
1240 	    cut_val[other_end] += (in_set[other_end]) ? 0 : weight;
1241 	 }
1242 	 while(set_size){
1243 	    if (set_cut_val < mult*(ceil(set_demand/capacity)) - etol &&
1244 		set_size > 2){
1245 	    memset(coef, 0, new_cut->size*sizeof(char));
1246 	    for (j = 1; j < vertnum; j++ ){
1247 	       if (in_set[j]){
1248 		  cur_nodept = verts + j;
1249 		  if (cur_nodept->orig_node_list_size)
1250 		     for (k = 0; k < cur_nodept->orig_node_list_size; k++)
1251 			(coef[(cur_nodept->orig_node_list)[k] >>
1252 			     DELETE_POWER]) |=
1253 			   (1 << ((cur_nodept->orig_node_list)[k] &
1254 				  DELETE_AND));
1255 		  (coef[j >> DELETE_POWER]) |= (1 << (j & DELETE_AND));
1256 	       }
1257 	    }
1258 #if 0
1259 	    new_cut->type = SUBTOUR_ELIM_SIDE;
1260 	    new_cut->rhs =  RHS(set_size, set_demand, capacity);
1261 	    cuts_found += cg_send_cut(new_cut);
1262 #endif
1263 #ifdef DIRECTED_X_VARS
1264 	    SEND_DIR_SUBTOUR_CONSTRAINT(set_size, set_demand);
1265 #else
1266 	    SEND_SUBTOUR_CONSTRAINT(set_size, set_demand);
1267 #endif
1268 	 }
1269 
1270 	 /* check the complement */
1271 
1272 	 complement_demand = total_demand - set_demand;
1273 	 complement_cut_val = set_cut_val - 2*(*cut_val) + 2*num_routes;
1274 	 complement_size = vertnum - 1 - set_size;
1275 	 if (complement_cut_val< mult*(ceil(complement_demand/capacity))-etol
1276 	     && complement_size > 2){
1277 	    memset(coef, 0, new_cut->size*sizeof(char));
1278 	    for (j = 1; j < vertnum; j++){
1279 	       if (!in_set[j]){
1280 	       cur_nodept=verts + j;
1281 		  if (cur_nodept->orig_node_list_size)
1282 		     for (k = 0; k < cur_nodept->orig_node_list_size; k++)
1283 			(coef[(cur_nodept->orig_node_list)[k] >>
1284 			     DELETE_POWER]) |=
1285 			   (1 << ((cur_nodept->orig_node_list)[k] &
1286 				  DELETE_AND));
1287 		  (coef[j>> DELETE_POWER]) |= (1 << ( j & DELETE_AND));
1288 	       }
1289 	    }
1290 #if 0
1291 	    new_cut->type = SUBTOUR_ELIM_SIDE;
1292 	    new_cut->rhs =  RHS(complement_size, complement_demand, capacity);
1293 	    cuts_found += cg_send_cut(new_cut);
1294 #endif
1295 #ifdef DIRECTED_X_VARS
1296 	    SEND_DIR_SUBTOUR_CONSTRAINT(complement_size, complement_demand);
1297 #else
1298 	    SEND_SUBTOUR_CONSTRAINT(complement_size, complement_demand);
1299 #endif
1300 	    SEND_SUBTOUR_CONSTRAINT(complement_size, complement_demand);
1301 	 }
1302 
1303 	 for (maxval = -1, pt = in_set+begin, dpt = cut_val+begin,
1304 		 j = begin; j < end; pt++, dpt++, j++){
1305 	    if (!(*pt) && *dpt > maxval){
1306 	       maxval = cut_val[j];
1307 		  max_vert = j;
1308 	    }
1309 	 }
1310 	 if (maxval > 0){    /* add the vertex to the set */
1311 	    in_set[max_vert] = 1;
1312 	    set_size += 1 + verts[max_vert].orig_node_list_size ;
1313 	    set_demand += demand[max_vert];
1314 	    cut_val[max_vert] = 0;
1315 	    for (e = verts[max_vert].first; e; e = e-> next_edge){
1316 	       other_end = e->other_end;
1317 	       weight  = e->data->weight;
1318 	       set_cut_val += (in_set[other_end]) ? (-weight) : weight;
1319 	       cut_val[other_end] += weight;
1320 	    }
1321 	 }
1322 	 else{ /* can't add anything to the set */
1323 	    break;
1324 	 }
1325       }
1326    }
1327    }
1328    FREE(new_cut->coef);
1329    return(cuts_found);
1330 }
1331