1 #include <memory.h>
2 #include <stdlib.h>
3 
4 #include "BB_constants.h"
5 #include "BB_macros.h"
6 #include "min_cut.h"
7 #include "vrp_macros.h"
8 #include "capforest.h"
9 #include "vrp_const.h"
10 
11 /*---------------------------------------------------------------------------*\
12 | This routine finds small cuts in the network by applying the ideas of       |
13 | Nagamochi and Ibaraki                                                       |
14 \*---------------------------------------------------------------------------*/
15 
min_cut(cg_vrp_spec * vrp,network * n,double etol)16 int min_cut(cg_vrp_spec *vrp, network *n, double etol)
17 {
18   cg_user_params *user_par = &vrp->par;
19   float mincut = 2;
20   register elist  *ne;
21   register vertex *nv = NULL;
22   register vertex *nu = NULL;
23   register vertex *nw = NULL;
24   int i, j;
25   float contr_above = 2;
26   float nodecut = 0;
27   int vertnum = n->vertnum;
28   int edgenum = n->edgenum;
29   char scannedmark = 0;
30   vertex **nen;
31   int itnum =0;
32   int cut_size = (vertnum >> DELETE_POWER) + 1;
33   int total_demand = 0, capacity  = vrp->capacity, num_cuts = 0;
34   cut_data *cut;
35   char change_in_vertnum = TRUE, *coef;
36   float max_q = 0;
37   vertex *verts = n->verts;
38   int v0, v1, cur_edge1_end;
39   elist *cur_edge1, *cur_edge2, *cur_edge3, *prev_edge;
40 
41   /*allocate memory for existing nodes list and the binary tree used by
42      capforest*/
43   nen = n->enodes = (vertex **) calloc (n->vertnum, sizeof(vertex *));
44   n->tnodes = (vertex **) calloc (n->vertnum, sizeof(vertex *));
45 
46   /*set up the data structure for the vertices*/
47   for(i = 0, total_demand = 0; i < vertnum; i++){
48     verts[i].scanned = 0;
49     verts[i].orignodenum = i;
50     verts[i].orig_node_list_size = 1;
51     verts[i].orig_node_list = (char *) calloc (cut_size, sizeof(char));
52     verts[i].orig_node_list[i >> DELETE_POWER] |= (1 << (i & DELETE_AND));
53     total_demand += verts[i].demand;
54   }
55 
56   /*shrink all the 1-edges in the current graph -- this operation preserves
57     all small cuts in the graph and reduces the size of the problem*/
58   num_cuts += shrink_one_edges(vrp, n, &vertnum, &edgenum, capacity, etol);
59 
60   cut = (cut_data *) calloc (1, sizeof(cut_data));
61   cut->size = cut_size;
62   coef = (char *) calloc (cut_size, sizeof(char));
63 
64   while (change_in_vertnum){/*if there hasn't been a change in the number of
65 			      vertices since the last iteration, then we are
66 			      done*/
67     itnum++;
68     change_in_vertnum = FALSE;
69 
70     scannedmark = 1-scannedmark;
71     max_q = capforest(n, vertnum, scannedmark);/*construct the capacitated
72 						 forests
73 						 discussed by N-I*/
74     contr_above = MIN(contr_above, max_q);/*reset the contraction limit*/
75 
76     /*--------------------------------------------------------------------*\
77     | run through the list of "super nodes" and contract any edge whose    |
78     | weight is above the limit "contr_above"                              |
79     \*--------------------------------------------------------------------*/
80 
81     for( i=1; i<vertnum; i++){
82       nv = nen[i];     /*get the next existing node*/
83       ne = nv->first;  /*get the first edge in its adjacency list*/
84       v0 = nv->orignodenum;
85       if (!v0) continue;  /*the depot should not be considered as part of any
86 			    cut*/
87 
88       /*-------------------------------------------------------------------*\
89       | scan the adjacency list of this node and contract all eligible edges|
90       | in its adjacency list                                               |
91       \*-------------------------------------------------------------------*/
92 
93       while (ne != NULL){
94 	if ((ne->data->q >= contr_above - etol) && (ne->data->v0) &&
95 	    (ne->data->v1)){
96 	  v1 = OTHER_END(ne, v0);
97 	  nw = verts+v1; /*nw is a pointer to the other end vertex of the
98 			   current node*/
99 
100 	  /*-----------------------------------------------------------------*\
101 	  | merge super node v0 with super node v1. "orig_node_list" is a bit |
102 	  | array containing the list of the original indices of nodes that   |
103 	  | have been contracted into this super node. the demand of the super|
104 	  | node is the sum of the demands of all the nodes that have been    |
105 	  | merged into it.                                                   |
106 	  \*-----------------------------------------------------------------*/
107 	  for (j = 0; j <cut_size; j++)
108 	    nv->orig_node_list[j] |= nw->orig_node_list[j];
109 	  nv->orig_node_list_size+=nw->orig_node_list_size;
110 	  nv->demand += nw->demand;
111 	  nw->scanned = TRUE;
112 
113 	  /*-----------------------------------------------------------------*\
114 	  | Now update the adjacency lists to reflect the contraction of edge |
115 	  | (v0, v1) To do this, we scan through the adjacency list of node v1|
116 	  | For each of these edges, we look at the other end vertex of the   |
117 	  | edge Then there are two cases. if the other end vertex is not also|
118 	  | adjacent to vertex v0, then we simply change the edge data so that|
119 	  | it originates from v0 instead of v1. If on the other hand, the    |
120 	  | vertex in question is also adjacent to v0, then there already     |
121 	  | exists an edge from v0 to the vertex in question and in this case,|
122 	  | we have to sum the weights of the two edges that are now to       |
123 	  | originate from v0 and merge them into one edge.                   |
124 	  \*-----------------------------------------------------------------*/
125 	  for (cur_edge1 = nw->first; cur_edge1;
126 	       cur_edge1 = cur_edge1->next_edge){
127 	    cur_edge1_end = OTHER_END(cur_edge1, v1);/*get the other end of
128 						       the edge*/
129 
130 	    for (cur_edge2 = nv->first, prev_edge = NULL; cur_edge2;
131 		 prev_edge = cur_edge2, cur_edge2 = cur_edge2->next_edge){
132 
133 	      /*---remove edge (v0, v1) from the adjacency list of v0--------*/
134 	      if ((OTHER_END(cur_edge2, v0)) == v1){
135 		cur_edge2->data->weight = 0;
136 		cur_edge2->data->deleted = TRUE;
137 		if (prev_edge)
138 		  prev_edge->next_edge = cur_edge2->next_edge;
139 		else
140 		  nv->first = cur_edge2->next_edge;
141 		if (!cur_edge2->next_edge) nv->last = prev_edge;
142 		nv->degree--;
143 		if (!nv->first) break;
144 	      }
145 
146 	      /*if cur_edge1_end is also adjacent to v0, then we must merge
147 		two edges*/
148 	      else if ((OTHER_END(cur_edge2, v0)) == cur_edge1_end){
149 		cur_edge2->data->weight += cur_edge1->data->weight;
150 		cur_edge1->data->deleted = TRUE;
151 		cur_edge1->data->weight = 0;
152 		nu = verts+cur_edge1_end;
153 
154 	        /*check whether we have found a violated cut*/
155 		if ((user_par->do_extra_checking) && (cur_edge1_end) &&
156 		    (cur_edge2->data->weight > 2 -
157 		     BINS(nv->demand + nu->demand, capacity))){
158 		  cur_edge3 = nv->first;
159 		  nodecut = 0;
160 		  while (cur_edge3 != NULL){
161 		    if ((cur_edge3->data != cur_edge2->data) &&
162 			(cur_edge3->data != ne->data))
163 		      nodecut += cur_edge3->data->weight;
164 		    cur_edge3 = cur_edge3->next_edge;
165 		  }
166 		  cur_edge3 = nu->first;
167 		  while (cur_edge3 != NULL){
168 		    if ((cur_edge3->data != cur_edge2->data) &&
169 			(cur_edge3->data != cur_edge1->data))
170 		      nodecut += cur_edge3->data->weight;
171 		    cur_edge3 = cur_edge3->next_edge;
172 		  }
173 		  cur_edge3 = nw->first;
174 		  while (cur_edge3 != NULL){
175 		    if ((cur_edge3->data != cur_edge1->data) &&
176 			(cur_edge3->data != ne->data))
177 		      nodecut += cur_edge3->data->weight;
178 		    cur_edge3 = cur_edge3->next_edge;
179 		  }
180 		  /*if the cut is violated, then send it back to the LP*/
181 		  if (nodecut <
182 		      2*(BINS(nv->demand + nu->demand, capacity))- etol){
183 		    cut->type =
184 		      (nv->orig_node_list_size+nu->orig_node_list_size <
185 		       n->vertnum/2 ? SUBTOUR_ELIM_SIDE:SUBTOUR_ELIM_ACROSS);
186 		    cut->rhs  =
187 		      (cut->type == SUBTOUR_ELIM_SIDE ?
188 		       RHS(nv->orig_node_list_size+nu->orig_node_list_size,
189 			   nv->demand+nu->demand, capacity):
190 		       2*(BINS(nv->demand + nu->demand, capacity)));
191 		    for (j = 0; j <cut_size; j++){
192 		      coef[j] = (nv->orig_node_list[j]|nu->orig_node_list[j]);
193 		    }
194 		    cut->coef = coef;
195 		    num_cuts += cg_send_cut(cut);
196 		  }
197 		} /*if ((cur_edge1_end) && (cur_edge2->data->weight > 2 - ...*/
198 
199 	        nu->degree--;
200 
201 		/*----update the adjacency list of cur_edge1_end------------*/
202 		for (prev_edge = NULL, cur_edge3 = nu->first;;
203 		     prev_edge = cur_edge3, cur_edge3 = cur_edge3->next_edge){
204 		  if ((OTHER_END(cur_edge3, cur_edge1_end)) == v1){
205 		    if (prev_edge){
206 		      prev_edge->next_edge = cur_edge3->next_edge;
207 		    }
208 		    else{
209 		      nu->first = cur_edge3->next_edge;
210 		    }
211 		    if (!cur_edge3->next_edge){
212 		      nu->last = prev_edge;
213 		    }
214 		    break;
215 		  }
216 		}
217 		cur_edge1->data->weight = 0;
218 	      } /*--else if ((OTHER_END(cur_edge2, v0)) == cur_edge1_end)*/
219 	    } /*--for (cur_edge2 = nv->first, prev_edge = NULL; cur_edge2;...*/
220 	    if (!nv->first) break;
221 	  } /*--for (cur_edge1 = nw->first; cur_edge1;...*/
222 
223 	  /*-----------------------------------------------------------------*\
224 	  | Now merge the remaining adjacency lists of v0 and v1              |
225 	  \*-----------------------------------------------------------------*/
226 
227 	  for (cur_edge1 = nw->first; cur_edge1;
228 	       cur_edge1 = cur_edge1->next_edge){
229 	    if (cur_edge1){
230 	      if (cur_edge1->data->weight){
231 		if (nv->last){
232 		  nv->last->next_edge = cur_edge1;
233 		  nv->last = cur_edge1;
234 		  nv->degree++;
235 		  if (cur_edge1->data->v0 == v1)
236 		    cur_edge1->data->v0 = v0;
237 		  else
238 		    cur_edge1->data->v1 = v0;
239 		}
240 		else{
241 		  nv->first = nv->last = cur_edge1;
242 		  if (cur_edge1->data->v0 == v1)
243 		    cur_edge1->data->v0 = v0;
244 		  else
245 		    cur_edge1->data->v1 = v0;
246 		}
247 	      }
248 	    }
249 	  }
250 	  nv->last->next_edge = NULL;
251 	  nw->degree = 0;
252 	  nw->first = nw->last = NULL;
253 
254 	  /*-----------------------------------------------------------------*\
255 	  | Now we update the existing node list. We want to delete v1 from   |
256 	  | the list but we have to be careful when we do this. We want to    |
257 	  | make sure that all the still existing nodes that have not been    |
258 	  | scanned yet in this iteration of the algorithm remain in a        |
259 	  | position in the array such that they will be scanned. In other    |
260 	  | words, in one of the positions with index >= i. There are several |
261 	  | cases that must be onsidered depending on where the node to be    |
262 	  | deleted is positioned in the existing nodes array                 |
263 	  \*-----------------------------------------------------------------*/
264 
265 	  if (nw->enodenum == i-1){
266 	    nen[nw->enodenum] = nen[--vertnum];
267 	    nen[vertnum]->enodenum = nw->enodenum;
268 	    i--;
269 	  }
270 	  else if (nw->enodenum < i-1){
271 	    nen[nw->enodenum] = nen[i-1];
272 	    nen[i-1]->enodenum = nw->enodenum;
273 	    nen[i-1] = nen[--vertnum];
274 	    nen[vertnum]->enodenum = i-1;
275 	    i--;
276 	  }
277 	  else{
278 	    nen[nw->enodenum] = nen[--vertnum];
279 	    nen[vertnum]->enodenum = nw->enodenum;
280 	  }
281 	  change_in_vertnum = TRUE;
282 	  edgenum--;
283 
284 	  /*now check to see if the cut induced by super node v0 is violated*/
285 	  cur_edge1 = nv->first;
286 	  nodecut = 0;
287 	  while (cur_edge1 != NULL){
288 	    nodecut += cur_edge1->data->weight;
289 	    cur_edge1 = cur_edge1->next_edge;
290 	  }
291 	  if (nodecut < 2*(BINS(nv->demand, capacity)) - etol){
292 	    cut->type = (nv->orig_node_list_size < n->vertnum/2 ?
293 			 SUBTOUR_ELIM_SIDE:SUBTOUR_ELIM_ACROSS);
294 	    cut->rhs  = (cut->type == SUBTOUR_ELIM_SIDE ?
295 			 RHS(nv->orig_node_list_size, nv->demand, capacity):
296 			 2*(BINS(nv->demand, capacity)));
297 	    cut->coef = nv->orig_node_list;
298 	    num_cuts += cg_send_cut(cut);
299 	  }
300 	  if (nodecut < mincut-etol){
301 	    mincut = nodecut;
302 	    if (user_par->update_contr_above)
303 	      contr_above = mincut;
304 	  }
305 	  ne = ne->next_edge;
306 	} /*---------if ((ne->data->q >= contr_above - etol) ...*/
307 	else{
308 	  ne = ne->next_edge;
309 	} /*---------if ((ne->data->q >= contr_above - etol) ...*/
310       } /*-----------while (ne != NULL) */
311     } /*-------------for( i=1; i<vertnum; i++)*/
312 
313     /*-----------------------------------------------------------------------*\
314     | finally as an extra check, we can scan through all the nodes after each |
315     | iteration is completely finished and see if the cuts induced by any of  |
316     | these nodes are violated but this probably isn't necessary and may      |
317     | result in the same cut being imposed several times                      |
318     \*-----------------------------------------------------------------------*/
319 
320 #if 0
321     if ((vertnum > 1) && (user_par->do_extra_checking)){
322       for ( i = 1; i<vertnum; i++ ){
323 	nodecut = 0;
324 	nv = nen[i];
325 	ne = nv->first;
326 	while (ne != NULL){
327 	  nodecut += ne->data->weight;
328 	  ne = ne->next_edge;
329 	}
330 	if (nodecut < 2*(BINS(nv->demand, capacity)) - etol){
331 	  cut->type = (nv->orig_node_list_size < n->vertnum/2 ?
332 		       SUBTOUR_ELIM_SIDE:SUBTOUR_ELIM_ACROSS);
333 	  cut->rhs  = (cut->type == SUBTOUR_ELIM_SIDE ?
334 		       RHS(nv->orig_node_list_size, nv->demand, capacity):
335 		       2*(BINS(nv->demand, capacity)));
336 	  cut->coef = nv->orig_node_list;
337 	  num_cuts += cg_send_cut(cut);
338 	}
339 	if (nodecut < mincut-etol){
340 	  mincut = nodecut;
341 	  if (user_par->update_contr_above)
342 	    contr_above = mincut;
343 	}
344       }
345     }
346     else
347 #endif
348     contr_above = mincut;
349   }
350 
351   FREE(coef);
352   FREE(cut);
353 
354   n->mincut = mincut;
355 
356   for (i = 0; i < n->vertnum; i++)
357     if (n->verts[i].orig_node_list)
358       free((char *)n->verts[i].orig_node_list);
359   free((char *)n->tnodes);
360   free((char *)n->enodes);
361   return(num_cuts);
362 }
363 
364 /*===========================================================================*/
365 
366 /*---------------------------------------------------------------------------*\
367 | This function shrinks 1-edges in the graph that are not adjacent to the     |
368 | depot                                                                       |
369 \*---------------------------------------------------------------------------*/
370 
shrink_one_edges(cg_vrp_spec * vrp,network * n,int * cur_verts,int * cur_edges,int capacity,double etol)371 int shrink_one_edges(cg_vrp_spec *vrp, network *n, int *cur_verts,
372 		     int *cur_edges, int capacity, double etol)
373 {
374   cg_user_params *user_par = &vrp->par;
375   register vertex *nw;
376   register vertex *nv = NULL;
377   register vertex *nu;
378   register elist *ne;
379   register edge *dat = NULL;
380   vertex **nen = n->enodes;
381   int cut_size = (n->vertnum >> DELETE_POWER) + 1;
382   vertex *verts = n->verts;
383   int num_cuts = 0, itnum = 0;
384   float nodecut;
385 
386   elist *new, *ne2, *nep, *prev;
387   int vertnum = *cur_verts;
388   int edgenum = 0;
389   int i, j;
390   cut_data *cut;
391   char change_in_vertnum = TRUE, *coef;
392   int v0, v1, cur_edge1_end;
393   elist *cur_edge1, *cur_edge2, *cur_edge3, *np, *prev_edge;
394 
395   /*----------------------------------------------------------------------*\
396   | First we contract all chains of 1-edges, called 1-paths. We simply look|
397   | at all nodes of degree 2 and keep following the 1-path in each         |
398   | direction from that node until we reach a node that is not of degree 2 |
399   | Contracting chains is much easier to do than actually contracting all  |
400   | the 1-edges and so we do it first                                      |
401   \*----------------------------------------------------------------------*/
402 
403   cut = (cut_data *) calloc (1, sizeof(cut_data));
404   cut->size = cut_size;
405   coef = (char *) calloc (cut_size, sizeof(char));
406 
407   for (i = 1; i < vertnum - 1; i++){
408     /*-----------check whether we have a degree 2 node-----------*/
409     if ((verts[i].degree != 2) || (verts[i].scanned))
410       continue;
411 
412     nv = verts+i;
413     ne2 = (ne = nv->first)->next_edge;
414     /*-----------------------------------------------------------------------*\
415     | follow the 1-path from vertex i until we hit either the depot or a node |
416     | that is not of degree 2                                                 |
417     \*-----------------------------------------------------------------------*/
418     for (nw = ne->other; nw->degree==2 && nw->orignodenum != 0; nw=ne->other){
419       if (nw == nv){/*if we come back to the same node, that means we have a
420 		      subtour and we can simply impose the appropriate
421 		      constraint */
422 	cut->coef = nv->orig_node_list;
423 	cut->type = (nv->orig_node_list_size < n->vertnum/2 ?
424 		     SUBTOUR_ELIM_SIDE:SUBTOUR_ELIM_ACROSS);
425 	cut->rhs  = (cut->type == SUBTOUR_ELIM_SIDE ?
426 		     RHS(nv->orig_node_list_size, nv->demand, capacity):
427 		     2*(BINS(nv->demand, capacity)));
428 	num_cuts += cg_send_cut(cut);
429 	nv->scanned = TRUE;
430 	nv->first = nv->last = NULL;
431 	break;
432       }
433       dat = ne->data;
434       nep = ne;
435       ne = nw->first;
436       if (ne->data == dat)
437 	ne = ne->next_edge;
438       if (!ne->other_end){
439 	ne = nep;
440 	break;
441       }
442 
443       nw->scanned = TRUE;
444 
445       /*As we go along, we contract all the nodes on the 1-path into node i*/
446       nv->demand += nw->demand;
447       for (j = 0; j <cut_size; j++)
448 	nv->orig_node_list[j] |= nw->orig_node_list[j];
449       nv->orig_node_list_size+=nw->orig_node_list_size;
450       nw->first = nw->last = NULL;
451       if (nv->demand > capacity){
452 	cut->coef = nv->orig_node_list;
453 	cut->type = (nv->orig_node_list_size < n->vertnum/2 ?
454 		     SUBTOUR_ELIM_SIDE:SUBTOUR_ELIM_ACROSS);
455 	cut->rhs  = (cut->type == SUBTOUR_ELIM_SIDE ?
456 		     RHS(nv->orig_node_list_size, nv->demand, capacity):
457 		     2*(BINS(nv->demand, capacity)));
458 	num_cuts += cg_send_cut(cut);
459       }
460     } /*--for (nw = ne->other; nw->degree==2 && nw->orignodenum != 0;
461 	nw=ne->other)..*/
462 
463     if (nv->scanned) continue; /*in this case we had a subtour and so we don't
464 				 need to go on*/
465     if (nw->orignodenum){ /*otherwise, nw is a pointer to one of the endpoints
466 			   of the 1-path. make new be the last edge on the
467 			   path that is adjacent to nw*/
468       dat = ne->data;
469       for (ne=nw->first; ne->data != dat; ne=ne->next_edge);
470       new = ne;
471     }
472     else{ /*this indicates that nw points to the depot which means that node i
473 	    is joined to the depot by a one edge. since we do not want to
474 	    contract one edges that are adjacent to the depot, we ignore this
475 	    edge and consider i to be the endpoint of the 1-path*/
476       if (!((ne = nv->first)->other_end))
477 	ne=ne->next_edge;
478       new = ne;
479     }
480 
481     /*-----------------------------------------------------------------------*\
482     | Now we find the other end point of the 1-path in a similar fashion      |
483     \*-----------------------------------------------------------------------*/
484     ne = ne2;
485     for (nu = ne->other; nu->degree==2 && nu->orignodenum != 0;
486 	 nu = ne->other){
487        dat = ne->data;
488        nep = ne;
489        ne = nu->first;
490        if (ne->data == dat)
491 	  ne = ne->next_edge;
492        if (!ne->other_end){
493 	  ne = nep;
494 	  break;
495        }
496 
497        nu->scanned = TRUE;
498 
499        nv->demand += nu->demand;
500        for (j = 0; j <cut_size; j++)
501 	  nv->orig_node_list[j] |= nu->orig_node_list[j];
502        nv->orig_node_list_size+=nu->orig_node_list_size;
503        nu->first = nu->last = NULL;
504        if (nv->demand > capacity){
505 	  cut->coef = nv->orig_node_list;
506 	  cut->type = (nv->orig_node_list_size < n->vertnum/2 ?
507 		     SUBTOUR_ELIM_SIDE:SUBTOUR_ELIM_ACROSS);
508 	  cut->rhs  = (cut->type == SUBTOUR_ELIM_SIDE ?
509 		       RHS(nv->orig_node_list_size, nv->demand, capacity):
510 		       2*(BINS(nv->demand, capacity)));
511 	  num_cuts += cg_send_cut(cut);
512        }
513     }
514     if(nu->orignodenum){
515        dat = ne->data;
516        for (ne=nu->first; ne->data != dat; ne=ne->next_edge);
517     }
518     else{
519        if (!((ne = nv->first)->other_end))
520 	  ne = ne->next_edge;
521     }
522 
523     /*-------------------------------------------------------------------*\
524     | Now we update the adjacency lists appropriately. Remember that if   |
525     | either nu or nw point to the depot, we must not consider these as   |
526     | endpoints of the 1-path so we have several cases                    |
527     \*-------------------------------------------------------------------*/
528 
529     if(nu->orignodenum && nw->orignodenum){
530       dat = ne->data = new->data;
531       ne->other = nw;
532       new->other = nu;
533       ne->data->v0 = nu->orignodenum;
534       ne->data->v1 = nw->orignodenum;
535       nv->scanned = TRUE;
536       /*------------- contract nv into nu ----------------------------------*/
537       nu->demand += nv->demand;
538       for (j = 0; j <cut_size; j++)
539 	nu->orig_node_list[j] |= nv->orig_node_list[j];
540       nu->orig_node_list_size+=nv->orig_node_list_size;
541       nv->first = nv->last = NULL;
542       /*-------- if nu induces a violated constraint, then impose it -------*/
543       if (nu->demand > capacity){
544 	cut->coef = nu->orig_node_list;
545 	cut->type = (nu->orig_node_list_size < n->vertnum/2 ?
546 		     SUBTOUR_ELIM_SIDE:SUBTOUR_ELIM_ACROSS);
547 	cut->rhs  = (cut->type == SUBTOUR_ELIM_SIDE ?
548 		     RHS(nu->orig_node_list_size, nu->demand, capacity):
549 		     2*(BINS(nu->demand, capacity)));
550 	num_cuts += cg_send_cut(cut);
551       }
552     }
553     else if (!nu->orignodenum){/*in this case, nu points to the depot*/
554       dat = ne->data = new->data;
555       ne->other = nw;
556       new->other = nv;
557       ne->data->v0 = nw->orignodenum;
558       ne->data->v1 = nv->orignodenum;
559       nu = nv;
560     }
561     else if (!nw->orignodenum){/*in this case, nw points to the depot*/
562       dat = ne->data = new->data;
563       ne->other = nv;
564       new->other = nu;
565       ne->data->v0 = nu->orignodenum;
566       ne->data->v1 = nv->orignodenum;
567       nw = nv;
568     }
569 
570     /*--------------------------------------------------------------------*\
571     | finally we check to see if there is another edge connecting nu to nw.|
572     | If so, then there must be a violated constraint since this implies a |
573     | cut of less than two in the graph. Also, this requires a bit of care |
574     | in updating the adjacency lists since then we will get a duplicate   |
575     | edge between nu and nw                                               |
576     \*--------------------------------------------------------------------*/
577 
578     for (prev = NULL, ne = nu->first; ne; prev = ne, ne = ne->next_edge){
579       if ((ne->other == nw) && (ne->data != dat)){
580 	dat->weight += ne->data->weight;
581 	ne->data->weight = 0;
582 	ne->data->deleted = TRUE;
583 	if (prev)
584 	  prev->next_edge = ne->next_edge;
585 	else
586 	  nu->first = ne->next_edge;
587 	if (!ne->next_edge)
588 	  nu->last = prev;
589 	cut->type =
590 	  (nu->orig_node_list_size+nw->orig_node_list_size < n->vertnum/2 ?
591 	   SUBTOUR_ELIM_SIDE:SUBTOUR_ELIM_ACROSS);
592 	cut->rhs  = (cut->type == SUBTOUR_ELIM_SIDE ?
593 		     RHS(nu->orig_node_list_size+nw->orig_node_list_size,
594 			 nu->demand+nw->demand, capacity):
595 		     2*(BINS(nu->demand+nw->demand, capacity)));
596 	for (j = 0; j <cut_size; j++){
597 	  coef[j] = nu->orig_node_list[j] | nw->orig_node_list[j];
598 	}
599 	cut->coef = coef;
600 	num_cuts += cg_send_cut(cut);
601 	for (prev = NULL, ne = nw->first; ne; prev = ne, ne = ne->next_edge){
602 	  if (ne->data->deleted){
603 	    if (prev)
604 	      prev->next_edge = ne->next_edge;
605 	    else
606 	      nw->first = ne->next_edge;
607 	    if(!ne->next_edge)
608 	      nw->last = prev;
609 	    break;
610 	  }
611 	}
612 	break;
613       }
614     }
615   }
616 
617   /*----- here we construct the list of existing nodes ----------------*/
618   for (i=0, j = 0; i<vertnum; j++){
619     if (!verts[j].scanned){
620       nen[i] = verts+j;
621       verts[j].enodenum = i;
622       i++;
623       edgenum += verts[j].degree;
624     }else{
625       vertnum--;
626     }
627   }
628   edgenum /= 2;
629 
630   if (!user_par->shrink_one_edges){
631     *cur_verts = vertnum;
632     *cur_edges = edgenum;
633     return(num_cuts);
634   }
635 
636   /*-----------------------------------------------------------------------*\
637   | Here we can optionally shrink all the 1-edges in the graph. I am not    |
638   | sure whether this is advantageous or not but it makes it easier to spot |
639   | violated cuts. These edges would get contracted anyway in the next phase|
640   | of this algorithm but it might make sense to contract them beforehand.  |
641   | The rest of this function is uncommented because it is exactly the same |
642   | as the main loop in the min_cut routine except that we only consider    |
643   | 1-edges for contraction instead of any edge above the threshold         |
644   \*-----------------------------------------------------------------------*/
645 
646   while (change_in_vertnum){
647     itnum++;
648     change_in_vertnum = FALSE;
649     for( i=1; i<vertnum; i++){
650       nv = nen[i];
651       ne = nv->first;
652       v0 = nv->orignodenum;
653       while (ne != NULL){
654 	if ((ne->data->weight >= 1 - etol) && (ne->data->v0) &&
655 	    (ne->data->v1)){
656 	  v1 = OTHER_END(ne, v0);
657 	  nw = verts+v1;
658 	  for (j = 0; j <cut_size; j++)
659 	    nv->orig_node_list[j] |= nw->orig_node_list[j];
660 	  nv->orig_node_list_size+=nw->orig_node_list_size;
661 	  nv->demand += nw->demand;
662 	  nw->scanned = TRUE;
663 	  if (nv->demand > capacity){
664 	    cut->coef = nv->orig_node_list;
665 	    cut->type = (nv->orig_node_list_size < n->vertnum/2 ?
666 			 SUBTOUR_ELIM_SIDE:SUBTOUR_ELIM_ACROSS);
667 	    cut->rhs  = (cut->type == SUBTOUR_ELIM_SIDE ?
668 			 RHS(nv->orig_node_list_size, nv->demand, capacity):
669 			 2*(BINS(nv->demand, capacity)));
670 	    num_cuts += cg_send_cut(cut);
671 	  }
672 
673 	  for (cur_edge1 = nw->first; cur_edge1;
674 	       cur_edge1 = cur_edge1->next_edge){
675 	    if ((cur_edge1_end = OTHER_END(cur_edge1, v1)) == v0){
676 	      cur_edge1->data->weight = 0;
677 	      cur_edge1->data->deleted = TRUE;
678 	    }
679 
680 	    for (cur_edge2 = nv->first, prev_edge = NULL; cur_edge2;
681 		 prev_edge = cur_edge2, cur_edge2 = cur_edge2->next_edge){
682 	      if ((OTHER_END(cur_edge2, v0)) == v1){
683 		cur_edge2->data->weight = 0;
684 		cur_edge2->data->deleted = TRUE;
685 		if (prev_edge)
686 		  prev_edge->next_edge = cur_edge2->next_edge;
687 		else
688 		  nv->first = cur_edge2->next_edge;
689 		if (!cur_edge2->next_edge) verts[v0].last = prev_edge;
690 		nv->degree--;
691 		if (!nv->first) break;
692 	      }
693 	      else if ((OTHER_END(cur_edge2, v0)) == cur_edge1_end){
694 		cur_edge2->data->weight += cur_edge1->data->weight;
695 		cur_edge1->data->deleted = TRUE;
696 		cur_edge1->data->weight = 0;
697 		nu = verts + cur_edge1_end;
698 		if ((cur_edge1_end) && (cur_edge2->data->weight > 2 -
699 		    BINS(nv->demand + nu->demand, capacity) + etol)){
700 		  cut->type =
701 		    (nv->orig_node_list_size+nu->orig_node_list_size <
702 		     n->vertnum/2 ? SUBTOUR_ELIM_SIDE:SUBTOUR_ELIM_ACROSS);
703 		  cut->rhs  =
704 		    (cut->type == SUBTOUR_ELIM_SIDE ?
705 		     RHS(nv->orig_node_list_size+nu->orig_node_list_size,
706 			 nv->demand+nu->demand, capacity):
707 		     2*(BINS(nv->demand+nu->demand, capacity)));
708 		  for (j = 0; j <cut_size; j++){
709 		    coef[j] = nv->orig_node_list[j] | nu->orig_node_list[j];
710 		  }
711 		  cut->coef = coef;
712 		  num_cuts += cg_send_cut(cut);
713 		}
714 		nu->degree--;
715 		for (prev_edge = NULL, cur_edge3 = nu->first;;
716 		     prev_edge = cur_edge3, cur_edge3 = cur_edge3->next_edge){
717 		  if ((OTHER_END(cur_edge3, cur_edge1_end)) == v1){
718 		    if (prev_edge){
719 		      prev_edge->next_edge = cur_edge3->next_edge;
720 		    }
721 		    else{
722 		      nu->first = cur_edge3->next_edge;
723 		    }
724 		    if (!cur_edge3->next_edge){
725 		      nu->last = prev_edge;
726 		    }
727 		    break;
728 		  }
729 		}
730 		cur_edge1->data->weight = 0;
731 	      }
732 	    }
733 	    if (!nv->first) break;
734 	  }
735 
736 	  for (cur_edge1 = nw->first; cur_edge1;
737 	       cur_edge1 = cur_edge1->next_edge){
738 	    if (cur_edge1){
739 	      if (cur_edge1->data->weight){
740 		if (nv->last){
741 		  nv->last->next_edge = cur_edge1;
742 		  nv->last = cur_edge1;
743 		  nv->degree++;
744 		  if (cur_edge1->data->v0 == v1)
745 		    cur_edge1->data->v0 = v0;
746 		  else
747 		    cur_edge1->data->v1 = v0;
748 		}
749 		else{
750 		  nv->first = nv->last = cur_edge1;
751 		  if (cur_edge1->data->v0 == v1)
752 		    cur_edge1->data->v0 = v0;
753 		  else
754 		    cur_edge1->data->v1 = v0;
755 		}
756 	      }
757 	    }
758 	  }
759 	  nv->last->next_edge = NULL;
760 	  nw->degree = 0;
761 	  nw->first = nw->last = NULL;
762 	  if (nw->enodenum == i-1){
763 	    nen[nw->enodenum] = nen[--vertnum];
764 	    nen[vertnum]->enodenum = nw->enodenum;
765 	    i--;
766 	  }
767 	  else if (nw->enodenum < i-1){
768 	    nen[nw->enodenum] = nen[i-1];
769 	    nen[i-1]->enodenum = nw->enodenum;
770 	    nen[i-1] = nen[--vertnum];
771 	    nen[vertnum]->enodenum = i-1;
772 	    i--;
773 	  }
774 	  else{
775 	    nen[nw->enodenum] = nen[--vertnum];
776 	    nen[vertnum]->enodenum = nw->enodenum;
777 	  }
778 	  change_in_vertnum = TRUE;
779 	  edgenum--;
780 	  np = ne;
781 	  ne = ne->next_edge;
782 	}
783 	else{
784 	  np = ne;
785 	  ne = ne->next_edge;
786 	}
787       }
788     }
789   }
790 
791   if ((vertnum > 1) && (user_par->do_extra_checking)){
792     for ( i = 1; i<vertnum; i++ ){
793       nodecut = 0;
794       nv = nen[i];
795       ne = nv->first;
796       while (ne != NULL){
797 	nodecut += ne->data->weight;
798 	ne = ne->next_edge;
799       }
800       if (nodecut < 2*(BINS(nv->demand, capacity)) - etol){
801 	cut->type = (nv->orig_node_list_size < n->vertnum/2 ?
802 		     SUBTOUR_ELIM_SIDE:SUBTOUR_ELIM_ACROSS);
803 	cut->rhs  = (cut->type == SUBTOUR_ELIM_SIDE ?
804 		     RHS(nv->orig_node_list_size, nv->demand, capacity):
805 		     2*(BINS(nv->demand, capacity)));
806 	cut->coef = nv->orig_node_list;
807 	num_cuts += cg_send_cut(cut);
808       }
809     }
810   }
811 
812   *cur_verts = vertnum;
813   *cur_edges = edgenum;
814 
815   FREE(coef);
816   FREE(cut);
817 
818   return(num_cuts);
819 }
820