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 <stdio.h>
18 #include <stdlib.h>
19 /*#include <sys/param.h>*/
20 #include <memory.h>
21 #include <math.h>
22
23 /* SYMPHONY include files */
24 #include "sym_constants.h"
25 #include "sym_macros.h"
26 #include "sym_proccomm.h"
27 #include "sym_messages.h"
28 #include "sym_timemeas.h"
29 #include "sym_lp_u.h"
30 #include "sym_dg_params.h"
31 #ifdef MULTI_CRITERIA
32 #include "sym_cg_u.h"
33 #endif
34
35 /* CNRP include files */
36 #include "cnrp_lp.h"
37 #include "cnrp_dg_functions.h"
38 #include "cnrp_macros.h"
39 #include "cnrp_const.h"
40
41 /*===========================================================================*/
42
43 /*===========================================================================*\
44 * This file contains the user-written functions for the LP process.
45 \*===========================================================================*/
46
47 /*===========================================================================*\
48 * Here is where the user must receive all of the data sent from
49 * user_send_lp_data() and set up data structures. Note that this function is
50 * only called if one of COMPILE_IN_LP or COMPILE_IN_TM is FALSE.
51 \*===========================================================================*/
52
user_receive_lp_data(void ** user)53 int user_receive_lp_data(void **user)
54 {
55 cnrp_spec *cnrp;
56 int vertnum;
57 int i, j, k, l;
58 int total_edgenum;
59 int zero_varnum, *zero_vars = NULL;
60
61 cnrp = (cnrp_spec *) calloc (1, sizeof(cnrp_spec));
62
63 *user = (void *)cnrp;
64
65 receive_char_array((char *)(&cnrp->par), sizeof(cnrp_lp_params));
66 receive_int_array(&cnrp->window, 1);
67 receive_int_array(&cnrp->numroutes, 1);
68 receive_int_array(&cnrp->vertnum, 1);
69 vertnum = cnrp->vertnum;
70 cnrp->demand = (double *) calloc (vertnum, sizeof(double));
71 receive_dbl_array(cnrp->demand, vertnum);
72 receive_dbl_array(&cnrp->capacity, 1);
73 total_edgenum = vertnum*(vertnum-1)/2;
74 cnrp->costs = (int *) calloc (total_edgenum, sizeof(int));
75 receive_int_array(cnrp->costs, total_edgenum);
76 receive_int_array(&zero_varnum, 1);
77 if (zero_varnum){
78 zero_vars = (int *) malloc (zero_varnum*sizeof(int));
79 receive_int_array(zero_vars, zero_varnum);
80 }
81 receive_dbl_array(&cnrp->utopia_fixed, 1);
82 receive_dbl_array(&cnrp->utopia_variable, 1);
83 receive_dbl_array(&cnrp->ub, 1);
84
85 /* The one additional edge allocated is for the extra variable when finding
86 nondominated solutions for multi-criteria problems */
87 cnrp->edges = (int *) calloc (2*(total_edgenum+1), sizeof(int));
88
89 /* Create the edge list (we assume a complete graph) The edge is set to
90 (0,0) in the edge list if it was eliminated in preprocessing*/
91 for (i = 1, k = 0, l = 0; i < vertnum; i++){
92 for (j = 0; j < i; j++){
93 if (l < zero_varnum && k == zero_vars[l]){
94 /*This is one of the zero edges*/
95 cnrp->edges[2*k] = cnrp->edges[2*k+1] = 0;
96 l++;
97 k++;
98 continue;
99 }
100 cnrp->edges[2*k] = j;
101 cnrp->edges[2*k+1] = i;
102 k++;
103 }
104 }
105 cnrp->edges[2*total_edgenum] = cnrp->edges[2*total_edgenum + 1] = 0;
106 FREE(zero_vars);
107
108 if (cnrp->par.prob_type == VRP || cnrp->par.prob_type == TSP ||
109 cnrp->par.prob_type == BPP){
110 cnrp->cur_sol = (_node *) calloc (cnrp->vertnum, sizeof(_node));
111 }else{
112 cnrp->cur_sol_tree = (int *) calloc (cnrp->vertnum - 1, ISIZE);
113 }
114
115 cnrp->variable_cost = cnrp->fixed_cost = MAXDOUBLE;
116
117 /*__BEGIN_EXPERIMENTAL_SECTION__*/
118 if (cnrp->window){
119 copy_node_set(cnrp->window, TRUE, (char *)"Weighted solution");
120 copy_node_set(cnrp->window, TRUE, (char *)"Flow solution");
121 }
122
123 /*___END_EXPERIMENTAL_SECTION___*/
124 return(USER_SUCCESS);
125 }
126
127 /*===========================================================================*/
128
129 /*===========================================================================*\
130 * Free all the user data structures
131 \*===========================================================================*/
132
user_free_lp(void ** user)133 int user_free_lp(void **user)
134 {
135 cnrp_spec *cnrp = (cnrp_spec *)(*user);
136
137 #if !(defined(COMPILE_IN_TM) && defined(COMPILE_IN_LP))
138 FREE(cnrp->demand);
139 FREE(cnrp->costs);
140 FREE(cnrp->edges);
141 #endif
142 FREE(cnrp->cur_sol);
143 FREE(cnrp->cur_sol_tree);
144 FREE(cnrp);
145 return(USER_SUCCESS);
146 }
147
148 /*===========================================================================*/
149
150 /*===========================================================================*\
151 * Here is where the user must create the initial LP relaxation for
152 * each search node. See the comments below.
153 \*===========================================================================*/
154
user_create_subproblem(void * user,int * indices,MIPdesc * mip,int * maxn,int * maxm,int * maxnz)155 int user_create_subproblem(void *user, int *indices, MIPdesc *mip,
156 int *maxn, int *maxm, int *maxnz)
157 {
158 return(USER_DEFAULT);
159
160 cnrp_spec *cnrp = (cnrp_spec *)user;
161 int *costs = cnrp->costs;
162 int *edges = cnrp->edges;
163 int i, j, maxvars = 0;
164 char resize = FALSE;
165 int vertnum = cnrp->vertnum;
166 int total_edgenum = vertnum*(vertnum-1)/2;
167 char prob_type = cnrp->par.prob_type, od_const = FALSE, d_x_vars = FALSE;
168 int v0, v1;
169 double flow_capacity;
170 #ifdef ADD_CAP_CUTS
171 int basecutnum = (2 + od_const)*vertnum - 1 + 2*total_edgenum;
172 #elif defined(ADD_FLOW_VARS)
173 int basecutnum = (2 + od_const)*vertnum - 1;
174 #else
175 int basecutnum = (1 + od_const)*vertnum;
176 #endif
177 #ifdef ADD_X_CUTS
178 basecutnum += total_edgenum;
179 #endif
180 #if defined(DIRECTED_X_VARS) && !defined(ADD_FLOW_VARS)
181 #ifdef FIND_NONDOMINATED_SOLUTIONS
182 int edgenum = (mip->n - 1)/2;
183 #else
184 int edgenum = (mip->n)/2;
185 #endif
186 #elif defined(ADD_FLOW_VARS)
187 #ifdef DIRECTED_X_VARS
188 #ifdef FIND_NONDOMINATED_SOLUTIONS
189 int edgenum = (mip->n - 1)/4;
190 #else
191 int edgenum = (mip->n)/4;
192 #endif
193
194 flow_capacity = cnrp->capacity;
195 #else
196 #ifdef FIND_NONDOMINATED_SOLUTIONS
197 int edgenum = (mip->n-1)/3;
198 #else
199 int edgenum = (mip->n)/3;
200 #endif
201
202 if (cnrp->par.prob_type == CSTP || cnrp->par.prob_type == CTP)
203 flow_capacity = cnrp->capacity;
204 else
205 flow_capacity = cnrp->capacity/2;
206 #endif
207 #else
208 #ifdef FIND_NONDOMINATED_SOLUTIONS
209 int edgenum = mip->n - 1;
210 #else
211 int edgenum = mip->n;
212 #endif
213 #endif
214
215 /* set up the inital LP data */
216
217 /*Estimate the number of nonzeros*/
218 #ifdef ADD_CAP_CUTS
219 mip->nz = 12*edgenum;
220 #elif defined(ADD_FLOW_VARS)
221 mip->nz = 8*edgenum;
222 #else
223 mip->nz = 3*edgenum;
224 #endif
225 #ifdef ADD_X_CUTS
226 mip->nz += 2*edgenum;
227 #endif
228 #ifdef FIND_NONDOMINATED_SOLUTIONS
229 mip->nz += mip->n + 1;
230 #endif
231 *maxm = MAX(100, 3 * mip->m);
232 #ifdef ADD_FLOW_VARS
233 *maxn = 3*total_edgenum;
234 #else
235 *maxn = total_edgenum;
236 #endif
237 #ifdef DIRECTED_X_VARS
238 *maxn += total_edgenum;
239 #endif
240 *maxnz = mip->nz + ((*maxm) * (*maxn) / 10);
241
242 /* Allocate the arrays. These are owned by SYMPHONY after returning. */
243 mip->matbeg = (int *) malloc((mip->n + 1) * ISIZE);
244 mip->matind = (int *) malloc((mip->nz) * ISIZE);
245 mip->matval = (double *) malloc((mip->nz) * DSIZE);
246 mip->obj = (double *) malloc(mip->n * DSIZE);
247 mip->ub = (double *) malloc(mip->n * DSIZE);
248 mip->lb = (double *) calloc(mip->n, DSIZE); /* zero lower bounds */
249 mip->rhs = (double *) malloc(mip->m * DSIZE);
250 mip->sense = (char *) malloc(mip->m * CSIZE);
251 mip->rngval = (double *) calloc(mip->m, DSIZE);
252 mip->is_int = (char *) calloc(mip->n, CSIZE);
253
254 #ifdef DIRECTED_X_VARS
255 /*whether or not we will have out-degree constraints*/
256 od_const = (prob_type == VRP || prob_type == TSP || prob_type == BPP);
257 d_x_vars = TRUE;
258 #endif
259
260 for (i = 0, j = 0; i < mip->n; i++){
261 #ifdef FIND_NONDOMINATED_SOLUTIONS
262 if (indices[i] == mip->n - 1){
263 mip->is_int[i] = FALSE;
264 mip->ub[i] = MAXINT;
265 mip->matbeg[i] = j;
266 mip->obj[i] = 1.0;
267 mip->matval[j] = -1.0;
268 mip->matind[j++] = basecutnum;
269 mip->matval[j] = -1.0;
270 mip->matind[j++] = basecutnum + 1;
271 continue;
272 }
273 #endif
274 if (indices[i] < total_edgenum){
275 mip->is_int[i] = TRUE;
276 mip->ub[i] = 1.0;
277 mip->matbeg[i] = j;
278 #ifdef FIND_NONDOMINATED_SOLUTIONS
279 mip->obj[i] = cnrp->par.rho*((double) costs[indices[i]]);
280 mip->matval[j] = cnrp->par.gamma*((double) costs[indices[i]]);
281 mip->matind[j++] = basecutnum;
282 #else
283 mip->obj[i] = cnrp->par.gamma*((double) costs[indices[i]]);
284 #endif
285 if (prob_type == CSTP || prob_type == CTP){
286 /*cardinality constraint*/
287 mip->matind[j] = 0;
288 mip->matval[j++] = 1.0;
289 }
290 /*in-degree constraint*/
291 mip->matval[j] = 1.0;
292 mip->matind[j++] = edges[2*indices[i]+1];
293 #ifdef DIRECTED_X_VARS
294 /*out-degree constraint*/
295 if (od_const){
296 mip->matval[j] = 1.0;
297 mip->matind[j++] = vertnum + edges[2*indices[i]];
298 }
299 #else
300 if (prob_type == VRP || prob_type == TSP ||
301 prob_type == BPP || edges[2*indices[i]]){
302 mip->matval[j] = 1.0;
303 mip->matind[j++] = edges[2*indices[i]];
304 }
305 #endif
306 #ifdef ADD_CAP_CUTS
307 v0 = edges[2*indices[i]];
308 mip->matval[j] = -flow_capacity + (v0 ? cnrp->demand[v0] : 0);
309 mip->matind[j++] = (2 + od_const)*vertnum - 1 + indices[i];
310 #ifndef DIRECTED_X_VARS
311 mip->matval[j] = -flow_capacity +
312 cnrp->demand[edges[2*indices[i] + 1]];
313 mip->matind[j++] = 2*cnrp->vertnum - 1 + total_edgenum +
314 indices[i];
315 #endif
316 #endif
317 #ifdef ADD_X_CUTS
318 mip->matval[j] = 1.0;
319 mip->matind[j++] = (2 + od_const)*vertnum-1 + 2*total_edgenum +
320 indices[i];
321 #endif
322 #ifdef DIRECTED_X_VARS
323 }else if (indices[i] < 2*total_edgenum){
324 mip->is_int[i] = TRUE;
325 mip->ub[i] = 1.0;
326 mip->matbeg[i] = j;
327 #ifdef FIND_NONDOMINATED_SOLUTIONS
328 mip->obj[i] = cnrp->par.rho*((double) costs[indices[i] -
329 total_edgenum]);
330 mip->matval[j] = cnrp->par.gamma*((double) costs[indices[i] -
331 total_edgenum]);
332 mip->matind[j++] = basecutnum;
333 #else
334 mip->obj[i] = cnrp->par.gamma*((double)costs[indices[i] -
335 total_edgenum]);
336 #endif
337 if (prob_type == CSTP || prob_type == CTP){
338 /*cardinality constraint*/
339 mip->matind[j] = 0;
340 mip->matval[j++] = 1.0;
341 }
342 /*in-degree constraint*/
343 if (od_const || edges[2*(indices[i] - total_edgenum)]){
344 mip->matval[j] = 1.0;
345 mip->matind[j++] = edges[2*(indices[i] - total_edgenum)];
346 }
347 /*out-degree constraint*/
348 if (od_const){
349 mip->matval[j] = 1.0;
350 mip->matind[j++] = vertnum + edges[2*(indices[i] -
351 total_edgenum)+1];
352 }
353 #ifdef ADD_CAP_CUTS
354 mip->matval[j] = -flow_capacity +
355 cnrp->demand[edges[2*(indices[i] - total_edgenum) + 1]];
356 mip->matind[j++] = (2 + od_const)*vertnum - 1 + indices[i];
357 #endif
358 #ifdef ADD_X_CUTS
359 mip->matval[j] = 1.0;
360 mip->matind[j++] = (2 + od_const)*vertnum-1 + 2*total_edgenum +
361 indices[i] - total_edgenum;
362 #endif
363 #endif
364 }else if (indices[i] < (2+d_x_vars)*total_edgenum){
365 mip->is_int[i] = FALSE;
366 v0 = edges[2*(indices[i]-(1+d_x_vars)*total_edgenum)];
367 mip->ub[i] = flow_capacity - (v0 ? cnrp->demand[v0] : 0);
368 mip->matbeg[i] = j;
369 #ifdef FIND_NONDOMINATED_SOLUTIONS
370 mip->obj[i] = cnrp->par.rho*((double) costs[indices[i] -
371 (1+d_x_vars)*total_edgenum]);
372 mip->matval[j] = cnrp->par.tau*((double) costs[indices[i]-
373 (1+d_x_vars)*total_edgenum]);
374 mip->matind[j++] = basecutnum + 1;
375 #else
376 mip->obj[i] =
377 cnrp->par.tau*((double) costs[indices[i]-
378 (1+d_x_vars)*total_edgenum]);
379 #endif
380 #ifdef ADD_CAP_CUTS
381 mip->matval[j] = 1.0;
382 mip->matval[j+1] = 1.0;
383 if (edges[2*(indices[i]-(1+d_x_vars)*total_edgenum)])
384 mip->matval[j+2] = -1.0;
385 mip->matind[j++] = (2 + od_const)*vertnum - 1 + indices[i] -
386 (1+d_x_vars)*total_edgenum;
387 mip->matind[j++] = (1+od_const)*vertnum + edges[2*(indices[i] -
388 (1+d_x_vars)*total_edgenum) + 1] - 1;
389 if (edges[2*(indices[i] - (1+d_x_vars)*total_edgenum)])
390 mip->matind[j++] = (1+od_const)*vertnum + edges[2*(indices[i] -
391 (1+d_x_vars)*total_edgenum)] - 1;
392 #else
393 mip->matval[j] = 1.0;
394 if (edges[2*(indices[i]-(1+d_x_vars)*total_edgenum)])
395 mip->matval[j+1] = -1.0;
396 mip->matind[j++] = (1+od_const)*vertnum + edges[2*(indices[i] -
397 (1+d_x_vars)*total_edgenum) + 1] - 1;
398 if (edges[2*(indices[i] - (1+d_x_vars)*total_edgenum)])
399 mip->matind[j++] = (1+od_const)*vertnum + edges[2*(indices[i] -
400 (1+d_x_vars)*total_edgenum)] - 1;
401 #endif
402 }else{
403 mip->is_int[i] = FALSE;
404 v1 = edges[2*(indices[i]-(2+d_x_vars)*total_edgenum) + 1];
405 mip->ub[i] = flow_capacity - cnrp->demand[v1];
406 mip->matbeg[i] = j;
407 #ifdef FIND_NONDOMINATED_SOLUTIONS
408 mip->obj[i] = cnrp->par.rho*((double) costs[indices[i] -
409 (2+d_x_vars)*total_edgenum]);
410 mip->matval[j] = cnrp->par.tau*((double) costs[indices[i]-
411 (2+d_x_vars)*total_edgenum]);
412 mip->matind[j++] = basecutnum + 1;
413 #else
414 mip->obj[i] =
415 cnrp->par.tau*((double) costs[indices[i]-
416 (2+d_x_vars)*total_edgenum]);
417 #endif
418 #ifdef ADD_CAP_CUTS
419 mip->matval[j] = 1.0;
420 mip->matval[j+1] = -1.0;
421 if (edges[2*(indices[i] - (2+d_x_vars)*total_edgenum)])
422 mip->matval[j+2] = 1.0;
423 mip->matind[j++] = (2+od_const)*vertnum - 1 + indices[i] -
424 (1+d_x_vars)*total_edgenum;
425 mip->matind[j++] = (1+od_const)*vertnum + edges[2*(indices[i] -
426 (2+d_x_vars)*total_edgenum)+1] - 1;
427 if (edges[2*(indices[i] - (2+d_x_vars)*total_edgenum)])
428 mip->matind[j++] = (1+od_const)*vertnum + edges[2*(indices[i] -
429 (2+d_x_vars)*total_edgenum)] - 1;
430 #else
431 mip->matval[j] = -1.0;
432 if (edges[2*(indices[i] - (2+d_x_vars)*total_edgenum)])
433 mip->matval[j+1] = 1.0;
434 mip->matind[j++] = (1+od_const)*vertnum + edges[2*(indices[i] -
435 (2+d_x_vars)*total_edgenum)+1] - 1;
436 if (edges[2*(indices[i] - (2+d_x_vars)*total_edgenum)])
437 mip->matind[j++] = (1+od_const)*vertnum + edges[2*(indices[i] -
438 (2+d_x_vars)*total_edgenum)] - 1;
439 #endif
440 }
441 }
442 mip->matbeg[i] = j;
443
444 /* set the initial right hand side */
445 if (od_const){
446 /*degree constraints for the depot*/
447 #if 0
448 mip->rhs[0] = cnrp->numroutes;
449 mip->sense[0] = 'E';
450 mip->rhs[vertnum] = cnrp->numroutes;
451 mip->sense[vertnum] = 'E';
452 #else
453 mip->rhs[0] = 1.0;
454 mip->sense[0] = 'G';
455 mip->rhs[vertnum] = 1.0;
456 mip->sense[vertnum] = 'G';
457 #endif
458 }else if (prob_type == VRP || prob_type == TSP || prob_type == BPP){
459 (mip->rhs[0]) = 2*cnrp->numroutes;
460 mip->sense[0] = 'E';
461 }else{
462 /*cardinality constraint*/
463 mip->rhs[0] = vertnum - 1;
464 mip->sense[0] = 'E';
465 }
466 for (i = vertnum - 1; i > 0; i--){
467 switch (prob_type){
468 case VRP:
469 case TSP:
470 case BPP:
471 if (od_const){
472 mip->rhs[i] = 1.0;
473 mip->sense[i] = 'E';
474 mip->rhs[i+vertnum] = 1.0;
475 mip->sense[i+vertnum] = 'E';
476 }else{
477 mip->rhs[i] = 2.0;
478 mip->sense[i] = 'E';
479 }
480 break;
481 case CSTP:
482 case CTP:
483 mip->rhs[i] = 1.0;
484 #ifdef DIRECTED_X_VARS
485 mip->sense[i] = 'E';
486 #else
487 mip->sense[i] = 'G';
488 #endif
489 break;
490 }
491 #ifdef ADD_FLOW_VARS
492 mip->rhs[(1+od_const)*vertnum + i - 1] = cnrp->demand[i];
493 mip->sense[(1+od_const)*vertnum + i - 1] = 'E';
494 #endif
495 }
496 #ifdef ADD_CAP_CUTS
497 for (i = (2+od_const)*vertnum - 1;
498 i < (2+od_const)*vertnum - 1 + 2*total_edgenum; i++){
499 mip->rhs[i] = 0.0;
500 mip->sense[i] = 'L';
501 }
502 #endif
503 #ifdef ADD_X_CUTS
504 for (i = (2+od_const)*vertnum-1+2*total_edgenum;
505 i < (2+od_const)*vertnum-1+3*total_edgenum; i++){
506 mip->rhs[i] = 1;
507 mip->sense[i] = 'L';
508 }
509 #endif
510 #ifdef FIND_NONDOMINATED_SOLUTIONS
511 mip->rhs[basecutnum] = cnrp->par.gamma*cnrp->utopia_fixed;
512 mip->sense[basecutnum] = 'L';
513 mip->rhs[basecutnum+1] = cnrp->par.tau*cnrp->utopia_variable;
514 mip->sense[basecutnum+1] = 'L';
515 #endif
516 return(USER_SUCCESS);
517 }
518
519
520 /*===========================================================================*/
521
522 /*===========================================================================*\
523 * This function takes an LP solution and checks it for feasibility.
524 * In our case, that means (1) it is integral (2) it is connected, and
525 * (3) the routes obey the capacity constraints.
526 \*===========================================================================*/
527
user_is_feasible(void * user,double lpetol,int varnum,int * indices,double * values,int * feasible,double * true_objval,char branching,double * heur_solution)528 int user_is_feasible(void *user, double lpetol, int varnum, int *indices,
529 double *values, int *feasible, double *true_objval,
530 char branching, double *heur_solution)
531 {
532 cnrp_spec *cnrp = (cnrp_spec *)user;
533 vertex *verts;
534 double *demand = cnrp->demand, capacity = cnrp->capacity, *compdemands;
535 int rcnt, *compnodes;
536 int vertnum = cnrp->vertnum, i, x_varnum;
537 network *n;
538 double *compcuts;
539 int total_edgenum = vertnum*(vertnum - 1)/2;
540 double fixed_cost = 0.0, variable_cost = 0.0;
541 #ifdef DIRECTED_X_VARS
542 char d_x_vars = TRUE;
543 #else
544 char d_x_vars = FALSE;
545 #endif
546
547 #ifdef ADD_FLOW_VARS
548 int tmp = varnum;
549 edge* edge1;
550 double flow_value, real_demand;
551
552 #ifndef ADD_CAP_CUTS
553 n = create_flow_net(indices, values, varnum, lpetol, cnrp->edges, demand,
554 vertnum);
555 #else
556 #ifdef DIRECTED_X_VARS
557 for (x_varnum = 0; x_varnum < varnum && indices[x_varnum] < 2*total_edgenum;
558 x_varnum++);
559 #else
560 for (x_varnum = 0; x_varnum < varnum && indices[x_varnum] < total_edgenum;
561 x_varnum++);
562 #endif
563
564 n = create_net(indices, values, x_varnum, lpetol, cnrp->edges, demand,
565 vertnum);
566 #endif
567 #else
568 n = create_net(indices, values, x_varnum, lpetol, cnrp->edges, demand,
569 vertnum);
570 #endif
571
572 if (!n->is_integral){
573 *feasible = IP_INFEASIBLE;
574 free_net(n);
575 return(USER_SUCCESS);
576 }
577
578 #if defined(ADD_FLOW_VARS) && !defined(ADD_CAP_CUTS)
579 #ifdef DIRECTED_X_VARS
580 for (i = 0, edge1 = n->edges; i < n->edgenum; i++, edge1++){
581 if ((flow_value = edge1->flow1) > lpetol){
582 real_demand = edge1->v0 ? demand[edge1->v0] : 0;
583 if ((capacity - real_demand)*edge1->weight1 < edge1->flow1 -
584 lpetol){
585 *feasible = IP_INFEASIBLE;
586 free_net(n);
587 return(USER_SUCCESS);
588 }
589 }
590 if ((flow_value = edge1->flow2) > lpetol){
591 if ((capacity-demand[edge1->v1])*edge1->weight2<edge1->flow2 -
592 lpetol){
593 *feasible = IP_INFEASIBLE;
594 free_net(n);
595 return(USER_SUCCESS);
596 }
597 }
598 }
599 #else
600 for (i = 0, edge1 = n->edges; i < n->edgenum; i++, edge1++){
601 if (capacity*edge1->weight < edge1->flow1 + edge1->flow2 - lpetol){
602 *feasible = IP_INFEASIBLE;
603 free_net(n);
604 return(USER_SUCCESS);
605 }
606 }
607 #endif
608 #endif
609
610 verts = n->verts;
611 compnodes = (int *) calloc (vertnum + 1, sizeof(int));
612 compdemands = (double *) calloc (vertnum + 1, sizeof(double));
613 compcuts = (double *) calloc (vertnum + 1, sizeof(double));
614 /*get the components of the solution graph without the depot to check if the
615 graph is connected or not*/
616 #if defined(ADD_FLOW_VARS) && !defined(ADD_CAP_CUTS)
617 rcnt = flow_connected(n, compnodes, compdemands, NULL, compcuts, NULL, lpetol);
618 #else
619 rcnt = connected(n, compnodes, compdemands, NULL, compcuts, NULL);
620 #endif
621
622 /*------------------------------------------------------------------------*\
623 * For each component check to see if the cut it induces is nonzero.
624 * Depending on the problem type, each component's cut value
625 * must be either 0, 1, or 2 since we have integrality
626 \*------------------------------------------------------------------------*/
627
628 for (i = 0; i < rcnt; i++){
629 if (compcuts[i+1] < lpetol || compdemands[i+1] > capacity){
630 *feasible = IP_INFEASIBLE;
631 FREE(compnodes);
632 FREE(compdemands);
633 FREE(compcuts);
634 free_net(n);
635 return(USER_SUCCESS);
636 }
637 }
638
639 FREE(compnodes);
640 FREE(compdemands);
641 FREE(compcuts);
642
643 if (cnrp->par.verbosity > 5){
644 display_support_graph(cnrp->window, FALSE, (char *)"Weighted solution",
645 x_varnum, indices, values, .000001, total_edgenum,
646 FALSE);
647 #ifdef ADD_FLOW_VARS
648 display_support_graph_flow(cnrp->window, FALSE, (char *)"Flow solution",
649 tmp, x_varnum, indices, values, .000001,
650 total_edgenum,CTOI_WAIT_FOR_CLICK_AND_REPORT);
651 #endif
652 }
653
654 #if defined(MULTI_CRITERIA) && defined(FIND_NONDOMINATED_SOLUTIONS) && 0
655 if (construct_feasible_solution(cnrp, n, true_objval, lpetol,
656 branching) > 0){
657 *feasible = IP_FEASIBLE_BUT_CONTINUE;
658 }else{
659 *feasible = IP_FEASIBLE;
660 }
661 #else
662 *feasible = IP_FEASIBLE;
663 #endif
664
665 free_net(n);
666
667 return (USER_SUCCESS);
668 }
669
670 /*===========================================================================*/
671
672 /*===========================================================================*\
673 * In my case, a feasible solution is specified most compactly by
674 * essentially a permutation of the customers along with routes numbers,
675 * specifying the order of the customers on their routes. This is just
676 * sent as a character array which then gets cast to an array of
677 * structures, one for each customers specifying the route number and
678 * the next customer on the route.
679 \*===========================================================================*/
680
user_send_feasible_solution(void * user,double lpetol,int varnum,int * indices,double * values)681 int user_send_feasible_solution(void *user, double lpetol, int varnum,
682 int *indices, double *values)
683 {
684 cnrp_spec *cnrp = (cnrp_spec *)user;
685
686 if (cnrp->par.prob_type == TSP || cnrp->par.prob_type == VRP ||
687 cnrp->par.prob_type == BPP)
688 send_char_array((char *)cnrp->cur_sol, cnrp->vertnum*sizeof(_node));
689 else
690 send_int_array(cnrp->cur_sol_tree, (cnrp->vertnum-1) * ISIZE);
691
692 return(USER_SUCCESS);
693 }
694
695
696 /*===========================================================================*/
697
698 /*===========================================================================*\
699 * This function graphically displays the current fractional solution
700 * This is done using the Interactie Graph Drawing program.
701 \*===========================================================================*/
702
user_display_lp_solution(void * user,int which_sol,int varnum,int * indices,double * values)703 int user_display_lp_solution(void *user, int which_sol, int varnum,
704 int *indices, double *values)
705 {
706 cnrp_spec *cnrp = (cnrp_spec *)user;
707 int i, total_edgenum = cnrp->vertnum*(cnrp->vertnum -1)/2;
708
709 #ifdef ADD_FLOW_VARS
710 for (i = 0; i < varnum && indices[i] < 2*total_edgenum; i++);
711 #endif
712
713 if (cnrp->par.verbosity > 10 ||
714 (cnrp->par.verbosity > 8 && (which_sol == DISP_FINAL_RELAXED_SOLUTION))
715 || (cnrp->par.verbosity > 6 && (which_sol == DISP_FEAS_SOLUTION))){
716 display_support_graph(cnrp->window, FALSE, (char *)"Weighted solution",
717 i, indices, values, .000001, total_edgenum, FALSE);
718 display_support_graph_flow(cnrp->window, FALSE, (char *)"Flow solution",
719 varnum, i, indices, values, .000001,
720 total_edgenum,CTOI_WAIT_FOR_CLICK_AND_REPORT);
721 }
722
723 if (which_sol == DISP_FINAL_RELAXED_SOLUTION){
724 return(DISP_NZ_INT);
725 }else{
726 return(USER_SUCCESS);
727 }
728 }
729
730 /*===========================================================================*/
731
732 /*===========================================================================*\
733 * You can add whatever information you want about a node to help you
734 * recreate it. I don't have a use for it, but maybe you will.
735 \*===========================================================================*/
736
user_add_to_desc(void * user,int * desc_size,char ** desc)737 int user_add_to_desc(void *user, int *desc_size, char **desc)
738 {
739 return(USER_DEFAULT);
740 }
741
742 /*===========================================================================*/
743
744 /*===========================================================================*\
745 * Compare cuts to see if they are the same. We use the default, which
746 * is just comparing byte by byte.
747 \*===========================================================================*/
748
user_same_cuts(void * user,cut_data * cut1,cut_data * cut2,int * same_cuts)749 int user_same_cuts(void *user, cut_data *cut1, cut_data *cut2, int *same_cuts)
750 {
751 /*for now, we just compare byte by byte, as in the previous version of the
752 code. Later, we might want to change this to be more efficient*/
753 return(USER_DEFAULT);
754 }
755
756 /*===========================================================================*/
757
758 /*===========================================================================*\
759 * This function receives a cut, unpacks it, and adds it to the set of
760 * rows to be added to the LP.
761 \*===========================================================================*/
762
user_unpack_cuts(void * user,int from,int type,int varnum,var_desc ** vars,int cutnum,cut_data ** cuts,int * new_row_num,waiting_row *** new_rows)763 int user_unpack_cuts(void *user, int from, int type, int varnum,
764 var_desc **vars, int cutnum, cut_data **cuts,
765 int *new_row_num, waiting_row ***new_rows)
766 {
767 int i, j, k, nzcnt = 0, nzcnt_side = 0, nzcnt_across = 0;
768 cnrp_spec *cnrp = (cnrp_spec *)user;
769 int index, v0, v1, *edges = cnrp->edges;
770 double demand;
771 waiting_row **row_list = NULL;
772 int *matind = NULL, *matind_across, *matind_side;
773 cut_data *cut;
774 char *coef;
775 double *matval = NULL;
776 int total_edgenum = cnrp->vertnum*(cnrp->vertnum - 1)/2;
777 int size, vertnum = ((cnrp_spec *)user)->vertnum;
778 int cliquecount = 0, val, edgeind;
779 char *clique_array, first_coeff_found, second_coeff_found, third_coeff_found;
780 #ifdef DIRECTED_X_VARS
781 char d_x_vars = TRUE;
782 #else
783 char d_x_vars = FALSE;
784 #endif
785 #if defined(ADD_FLOW_VARS) && defined(DIRECTED_X_VARS)
786 int numarcs, *arcs;
787 char *coef2;
788 #endif
789 *new_row_num = cutnum;
790 if (cutnum > 0)
791 *new_rows = row_list = (waiting_row **) calloc (cutnum,
792 sizeof(waiting_row *));
793
794 for (j = 0; j < cutnum; j++){
795 coef = (cut = cuts[j])->coef;
796 cuts[j] = NULL;
797 (row_list[j] = (waiting_row *) malloc(sizeof(waiting_row)))->cut = cut;
798 switch (cut->type){
799 /*-------------------------------------------------------------------*\
800 * The subtour elimination constraints are stored as a vector of
801 * bits indicating which side of the cut each customer is on
802 \*-------------------------------------------------------------------*/
803
804 #if 0
805 case SUBTOUR_ELIM:
806 matind_side = (int *) malloc(varnum * ISIZE);
807 matind_across = (int *) malloc(varnum * ISIZE);
808 for (i = 0, nzcnt = 0; i < varnum; i++){
809 #ifdef ADD_FLOW_VARS
810 #ifdef DIRECTED_X_VARS
811 if (vars[i]->userind < 2*total_edgenum){
812 if (vars[i]->userind >= total_edgenum){
813 edgeind = vars[i]->userind - total_edgenum;
814 }else{
815 edgeind = vars[i]->userind;
816 }
817 #else
818 if ((edgeind = vars[i]->userind) < total_edgenum){
819 #endif
820 #else
821 #ifdef DIRECTED_X_VARS
822 {
823 if (vars[i]->userind >= total_edgenum){
824 edgeind = vars[i]->userind - total_edgenum;
825 }else{
826 edgeind = vars[i]->userind;
827 }
828 #else
829 {
830 edgeind = vars[i]->userind;
831 #endif
832 #endif
833 v0 = edges[edgeind << 1];
834 v1 = edges[(edgeind << 1) + 1];
835 if (coef[v0 >> DELETE_POWER] &
836 (1 << (v0 & DELETE_AND)) &&
837 (coef[v1 >> DELETE_POWER]) &
838 (1 << (v1 & DELETE_AND))){
839 matind_side[nzcnt_side++] = i;
840 }else if ((coef[v0 >> DELETE_POWER] >>
841 (v0 & DELETE_AND) & 1) ^
842 (coef[v1 >> DELETE_POWER] >>
843 (v1 & DELETE_AND) & 1)){
844 matind_across[nzcnt_across++] = i;
845 }
846 }
847 }
848 cut->type = nzcnt_side < nzcnt_across ? SUBTOUR_ELIM_SIDE :
849 SUBTOUR_ELIM_ACROSS;
850 cut->deletable = TRUE;
851 switch (cut->type){
852 case SUBTOUR_ELIM_SIDE:
853 row_list[j]->nzcnt = nzcnt_side;
854 row_list[j]->matind = matind_side;
855 cut->rhs = 0; /*RHS(compnodes[i+1],compdemands[i+1], capacity)*/
856 cut->sense = 'L';
857 FREE(matind_across);
858 break;
859
860 case SUBTOUR_ELIM_ACROSS:
861 row_list[j]->nzcnt = nzcnt_across;
862 row_list[j]->matind = matind_across;
863 cut->rhs = 0; /*2*BINS(compdemands[i+1], capacity)*/
864 cut->sense = 'G';
865 FREE(matind_side);
866 break;
867 }
868
869 break;
870 #endif
871 case SUBTOUR_ELIM:
872 case SUBTOUR_ELIM_SIDE:
873 matind = (int *) malloc(varnum * ISIZE);
874 for (i = 0, nzcnt = 0; i < varnum; i++){
875 #ifdef DIRECTED_X_VARS
876 if (vars[i]->userind < 2*total_edgenum){
877 if (vars[i]->userind >= total_edgenum){
878 edgeind = vars[i]->userind - total_edgenum;
879 }else{
880 edgeind = vars[i]->userind;
881 }
882 #else
883 if ((edgeind = vars[i]->userind) < total_edgenum){
884 #endif
885 v0 = edges[edgeind << 1];
886 v1 = edges[(edgeind << 1) + 1];
887 if (coef[v0 >> DELETE_POWER] & (1 << (v0 & DELETE_AND)) &&
888 (coef[v1 >> DELETE_POWER]) & (1 << (v1 & DELETE_AND))){
889 matind[nzcnt++] = i;
890 }
891 }
892 }
893 cut->sense = 'L';
894 cut->deletable = TRUE;
895 cut->branch = DO_NOT_BRANCH_ON_THIS_ROW;
896 break;
897
898 case SUBTOUR_ELIM_ACROSS:
899 matind = (int *) malloc(varnum * ISIZE);
900 for (i = 0, nzcnt = 0; i < varnum; i++){
901 #ifdef DIRECTED_X_VARS
902 if (vars[i]->userind < 2*total_edgenum){
903 if (vars[i]->userind >= total_edgenum){
904 edgeind = vars[i]->userind - total_edgenum;
905 v1 = edges[edgeind << 1];
906 v0 = edges[(edgeind << 1) + 1];
907 }else{
908 edgeind = vars[i]->userind;
909 v0 = edges[edgeind << 1];
910 v1 = edges[(edgeind << 1) + 1];
911 }
912 if ((coef[v1 >> DELETE_POWER] >> (v1 & DELETE_AND) & 1) &&
913 !(coef[v0 >> DELETE_POWER] >> (v0 & DELETE_AND) & 1)){
914 matind[nzcnt++] = i;
915 }
916 }
917 #else
918 if (vars[i]->userind < total_edgenum){
919 edgeind = vars[i]->userind;
920 v0 = edges[edgeind << 1];
921 v1 = edges[(edgeind << 1) + 1];
922 if ((coef[v1 >> DELETE_POWER] >> (v1 & DELETE_AND) & 1) ^
923 (coef[v0 >> DELETE_POWER] >> (v0 & DELETE_AND) & 1)){
924 matind[nzcnt++] = i;
925 }
926 }
927 #endif
928 }
929 cut->sense = 'G';
930 cut->deletable = TRUE;
931 cut->branch = DO_NOT_BRANCH_ON_THIS_ROW;
932 break;
933
934 #if defined(ADD_FLOW_VARS) && defined(DIRECTED_X_VARS)
935 case MIXED_DICUT:
936 matind = (int *) malloc(varnum * ISIZE);
937 matval = (double *) malloc(varnum*DSIZE);
938 demand = ((double *)coef)[0];
939 numarcs = ((int *)(coef + DSIZE))[0];
940 /* Array of the nodes in the set S */
941 coef2 = coef + DSIZE + ISIZE;
942 /* Array of the arcs in the set C */
943 arcs = (int *) (coef + DSIZE + ISIZE + (vertnum >> DELETE_POWER)+1);
944 for (i = 0, nzcnt = 0; i < varnum; i++){
945 if (vars[i]->userind < 2*total_edgenum){
946 if (vars[i]->userind >= total_edgenum){
947 edgeind = vars[i]->userind - total_edgenum;
948 v1 = edges[edgeind << 1];
949 v0 = edges[(edgeind << 1) + 1];
950 }else{
951 edgeind = vars[i]->userind;
952 v0 = edges[edgeind << 1];
953 v1 = edges[(edgeind << 1) + 1];
954 }
955 if ((coef2[v1 >> DELETE_POWER] >> (v1 & DELETE_AND) & 1) &&
956 !(coef2[v0 >> DELETE_POWER] >> (v0 & DELETE_AND) & 1)){
957 for (k = 0; k < numarcs; k++){
958 if (v0 == arcs[k << 1] && v1 == arcs[(k << 1) + 1])
959 break;
960 }
961 if (k == numarcs){
962 matind[nzcnt] = i;
963 matval[nzcnt++] = MIN(cnrp->capacity, demand);
964 }
965 }
966 }else{
967 if (vars[i]->userind < 3*total_edgenum){
968 edgeind = vars[i]->userind - 2*total_edgenum;
969 v0 = edges[edgeind << 1];
970 v1 = edges[(edgeind << 1) + 1];
971 }else{
972 edgeind = vars[i]->userind - 3*total_edgenum;
973 v1 = edges[edgeind << 1];
974 v0 = edges[(edgeind << 1) + 1];
975 }
976 if ((coef2[v1 >> DELETE_POWER] >> (v1 & DELETE_AND) & 1) &&
977 !(coef2[v0 >> DELETE_POWER] >> (v0 & DELETE_AND) & 1)){
978 for (k = 0; k < numarcs; k++){
979 if (v0 == arcs[k << 1] && v1 == arcs[(k << 1) + 1])
980 break;
981 }
982 if (k < numarcs){
983 matind[nzcnt] = i;
984 matval[nzcnt++] = 1.0;
985 }
986 }
987 }
988 }
989 cut->sense = 'G';
990 cut->deletable = TRUE;
991 cut->branch = DO_NOT_BRANCH_ON_THIS_ROW;
992 break;
993 #endif
994
995 #ifdef ADD_FLOW_VARS
996 case FLOW_CAP:
997 matind = (int *) malloc(3 * ISIZE);
998 matval = (double *) malloc(3 * DSIZE);
999
1000 index = ((int *)coef)[0];
1001 v0 = index < total_edgenum ? edges[index << 1] :
1002 edges[(index - total_edgenum) << 1];
1003 v1 = index < total_edgenum ? edges[(index << 1) + 1] :
1004 edges[((index - total_edgenum) << 1) + 1];
1005 if (v0){
1006 demand =
1007 index < total_edgenum ? cnrp->demand[v0] : cnrp->demand[v1];
1008 }else{
1009 demand = 0;
1010 }
1011
1012 first_coeff_found = second_coeff_found = third_coeff_found = FALSE;
1013 for (i = 0; i < varnum && (!first_coeff_found ||
1014 !second_coeff_found); i++){
1015 if (vars[i]->userind == index){
1016 matind[0] = i;
1017 first_coeff_found = TRUE;
1018 }
1019 if (vars[i]->userind == (index + 2 * total_edgenum)){
1020 matind[1] = i;
1021 second_coeff_found = TRUE;
1022 }
1023 }
1024 #ifndef DIRECTED_X_VARS
1025 for (i = 0; i < varnum && !third_coeff_found; i++){
1026 if (vars[i]->userind == (index + total_edgenum)){
1027 matind[2] = i;
1028 third_coeff_found = TRUE;
1029 }
1030 }
1031 #endif
1032 if (first_coeff_found){
1033 #ifdef DIRECTED_X_VARS
1034 matval[0] = -(cnrp->capacity - demand);
1035 #else
1036 if (cnrp->par.prob_type == CSTP || cnrp->par.prob_type == CTP){
1037 matval[0] = -cnrp->capacity;
1038 }else{
1039 matval[0] = -cnrp->capacity/2;
1040 }
1041 #endif
1042 if (second_coeff_found){
1043 if (third_coeff_found){
1044 matval[1] = matval[2] = 1.0;
1045 nzcnt = 3;
1046 }else{
1047 matval[1] = 1.0;
1048 nzcnt = 2;
1049 }
1050 }else if (third_coeff_found){
1051 matind[1] = matind[2];
1052 matval[1] = 1.0;
1053 nzcnt = 2;
1054 }else{
1055 nzcnt = 0;
1056 }
1057 }else if (second_coeff_found){
1058 matind[0] = matind[1];
1059 matval[0] = 1.0;
1060 if (third_coeff_found){
1061 matind[1] = matind[2];
1062 matval[1] = 1.0;
1063 nzcnt = 2;
1064 }else{
1065 nzcnt = 1;
1066 }
1067 }else if (third_coeff_found){
1068 matind[0] = matind[2];
1069 matval[0] = 1.0;
1070 nzcnt = 1;
1071 }else{
1072 nzcnt = 0;
1073 }
1074 cut->sense = 'L';
1075 cut->deletable = FALSE;
1076 cut->branch = DO_NOT_BRANCH_ON_THIS_ROW;
1077 break;
1078
1079 case TIGHT_FLOW:
1080
1081 matind = (int *) malloc((vertnum + 1) * ISIZE);
1082 matval = (double *) malloc((vertnum + 1) * DSIZE);
1083
1084 if ((index = ((int *)coef)[0]) < total_edgenum){
1085 v0 = edges[index << 1];
1086 v1 = edges[(index << 1) + 1];
1087 }else{
1088 v1 = edges[(index - total_edgenum) << 1];
1089 v0 = edges[((index - total_edgenum) << 1) + 1];
1090 }
1091
1092 #ifdef DIRECTED_X_VARS
1093 for (nzcnt = 0, k = 0; k < varnum; k++){
1094 if (vars[k]->userind == index){
1095 matind[nzcnt] = k;
1096 matval[nzcnt++] = v1 ? -cnrp->demand[v1] : 0;
1097 break;
1098 }
1099 }
1100 #else
1101 for (nzcnt = 0, k = 0; k < varnum; k++){
1102 if (vars[k]->userind == (index < total_edgenum ? index :
1103 index - total_edgenum)){
1104 matind[nzcnt] = k;
1105 matval[nzcnt++] = v1 ? -cnrp->demand[v1] : 0;
1106 break;
1107 }
1108 }
1109 #endif
1110 for (k = 0; k < varnum; k++){
1111 if (vars[k]->userind == index + (1 + d_x_vars) * total_edgenum){
1112 matind[nzcnt] = k;
1113 matval[nzcnt++] = 1.0;
1114 break;
1115 }
1116 }
1117 /* This loop is done very inefficiently and should be rewritten */
1118 for (i = 0; i < v1; i++){
1119 index = INDEX(i, v1) + (2 + d_x_vars) * total_edgenum;
1120 for (k = 0; k < varnum; k++){
1121 if (vars[k]->userind == index){
1122 matind[nzcnt] = k;
1123 matval[nzcnt++] = -1.0;
1124 break;
1125 }
1126 }
1127 }
1128 for (i = v1 + 1; i < vertnum; i++){
1129 index = INDEX(i, v1) + (1 + d_x_vars) * total_edgenum;
1130 for (k = 0; k < varnum; k++){
1131 if (vars[k]->userind == index){
1132 matind[nzcnt] = k;
1133 matval[nzcnt++] = -1.0;
1134 break;
1135 }
1136 }
1137 }
1138 cut->sense = 'L';
1139 cut->deletable = FALSE;
1140 cut->branch = DO_NOT_BRANCH_ON_THIS_ROW;
1141
1142 break;
1143 #endif
1144
1145 #ifdef DIRECTED_X_VARS
1146 case X_CUT:
1147 matind = (int *) malloc(2 * ISIZE);
1148 matval = (double *) malloc(2 * DSIZE);
1149 first_coeff_found = second_coeff_found = FALSE;
1150 for (i = 0, nzcnt = 0; i < varnum && (!first_coeff_found ||
1151 !second_coeff_found); i++){
1152 if (vars[i]->userind == ((int *)coef)[0]){
1153 matind[0] = i;
1154 first_coeff_found = TRUE;
1155 }
1156 if (vars[i]->userind == ((int *)coef)[0]+total_edgenum){
1157 matind[1] = i;
1158 second_coeff_found = TRUE;
1159 }
1160 }
1161 if (!first_coeff_found || !second_coeff_found){
1162 printf("ERROR constructing X Cut!!\n\n");
1163 nzcnt = 0;
1164 }else{
1165 matval[0] = matval[1] = 1.0;
1166 nzcnt = 2;
1167 }
1168 cut->sense = 'L';
1169 cut->deletable = FALSE;
1170 cut->branch = DO_NOT_BRANCH_ON_THIS_ROW;
1171 break;
1172 #endif
1173
1174 case OPTIMALITY_CUT_FIXED:
1175 matind = (int *) malloc (varnum * ISIZE);
1176 matval = (double *) malloc (varnum * DSIZE);
1177 for (nzcnt = 0, i = 0; i < varnum; i++){
1178 if (vars[i]->userind < total_edgenum){
1179 matind[nzcnt] = i;
1180 matval[nzcnt++] = cnrp->costs[vars[i]->userind];
1181 }
1182 #ifdef DIRECTED_X_VARS
1183 else if (vars[i]->userind < 2*total_edgenum){
1184 matind[nzcnt] = i;
1185 matval[nzcnt++] = cnrp->costs[vars[i]->userind - total_edgenum];
1186 }
1187 #endif
1188 }
1189 cut->sense = 'L';
1190 cut->deletable = FALSE;
1191 cut->branch = DO_NOT_BRANCH_ON_THIS_ROW;
1192 break;
1193
1194 case OPTIMALITY_CUT_VARIABLE:
1195 matind = (int *) malloc (varnum * ISIZE);
1196 matval = (double *) malloc (varnum * DSIZE);
1197 for (nzcnt = 0, i = 0; i < varnum; i++){
1198 if (vars[i]->userind >= (1 + d_x_vars) * total_edgenum &&
1199 vars[i]->userind < (3 + d_x_vars) * total_edgenum){
1200 matind[nzcnt] = i;
1201 if (vars[i]->userind < (2 + d_x_vars) * total_edgenum){
1202 matval[nzcnt++] =
1203 cnrp->costs[vars[i]->userind-(1+d_x_vars) * total_edgenum];
1204 }else{
1205 matval[nzcnt++] =
1206 cnrp->costs[vars[i]->userind-(2+d_x_vars) * total_edgenum];
1207 }
1208 }
1209 }
1210 cut->sense = 'L';
1211 cut->deletable = FALSE;
1212 cut->branch = DO_NOT_BRANCH_ON_THIS_ROW;
1213 break;
1214
1215 case CLIQUE:
1216 size = (vertnum >> DELETE_POWER) + 1;
1217 memcpy(&cliquecount, coef, ISIZE);
1218 matind = (int *) malloc(cliquecount*varnum*ISIZE);
1219 matval = (double *) malloc(cliquecount*varnum*DSIZE);
1220 coef += ISIZE;
1221 for (nzcnt = 0, i = 0; i < varnum; i++){
1222 #ifdef ADD_FLOW_VARS
1223 #ifdef DIRECTED_X_VARS
1224 if (vars[i]->userind < 2*total_edgenum){
1225 if (vars[i]->userind >= total_edgenum){
1226 edgeind = vars[i]->userind - total_edgenum;
1227 }else{
1228 edgeind = vars[i]->userind;
1229 }
1230 #else
1231 if ((edgeind = vars[i]->userind) < total_edgenum){
1232 #endif
1233 #else
1234 #ifdef DIRECTED_X_VARS
1235 {
1236 if (vars[i]->userind >= total_edgenum){
1237 edgeind = vars[i]->userind - total_edgenum;
1238 }else{
1239 edgeind = vars[i]->userind;
1240 }
1241 #else
1242 {
1243 edgeind = vars[i]->userind;
1244 #endif
1245 #endif
1246 v0 = edges[edgeind << 1];
1247 v1 = edges[(edgeind << 1) + 1];
1248 val = 0;
1249 for (k = 0; k < cliquecount; k++){
1250 clique_array = coef + size * k;
1251 if (clique_array[v0 >> DELETE_POWER] &
1252 (1 << (v0 & DELETE_AND)) &&
1253 (clique_array[v1 >> DELETE_POWER]) &
1254 (1 << (v1 & DELETE_AND))){
1255 val += 1;
1256 }
1257 }
1258 if (val){
1259 matind[nzcnt] = i;
1260 matval[nzcnt++] = val;
1261 }
1262 }
1263 }
1264 cut->branch = DO_NOT_BRANCH_ON_THIS_ROW;
1265 cut->deletable = TRUE;
1266 break;
1267
1268 default:
1269 printf("Unrecognized cut type %i!\n", cut->type);
1270 }
1271
1272 row_list[j]->matind = matind =
1273 (int *) realloc((char *)matind, nzcnt*ISIZE);
1274 row_list[j]->nzcnt = nzcnt;
1275 if (cut->type == SUBTOUR_ELIM || cut->type == SUBTOUR_ELIM_ACROSS ||
1276 cut->type == SUBTOUR_ELIM_SIDE){
1277 row_list[j]->matval = matval = (double *) malloc(nzcnt * DSIZE);
1278 for (i = nzcnt-1; i >= 0; i--)
1279 matval[i] = 1;
1280 cut->branch = ALLOWED_TO_BRANCH_ON;
1281 }else{
1282 row_list[j]->matval=(double *) realloc((char *)matval, nzcnt * DSIZE);
1283 }
1284 }
1285
1286 return(USER_SUCCESS);
1287 }
1288
1289 /*===========================================================================*/
1290
1291 int user_send_lp_solution(void *user, int varnum, var_desc **vars, double *x,
1292 int where)
1293 {
1294 return(SEND_NONZEROS);
1295 }
1296
1297 /*===========================================================================*/
1298
1299 /*===========================================================================*\
1300 * This routine does logical fixing of variables
1301 \*===========================================================================*/
1302
1303 int user_logical_fixing(void *user, int varnum, var_desc **vars, double *x,
1304 char *status, int *num_fixed)
1305 {
1306 cnrp_spec *cnrp = (cnrp_spec *)user;
1307 lp_net *lp_net;
1308 double *compdemands, capacity = cnrp->capacity;
1309 int numchains = 0, v0, v1;
1310 lp_net_node *verts;
1311 int i;
1312 int *edges = cnrp->edges;
1313 int fixed_num = 0;
1314 #ifdef ADD_FLOW_VARS
1315 int total_edgenum = cnrp->vertnum*(cnrp->vertnum - 1)/2;
1316 #endif
1317
1318 return(USER_SUCCESS); /*for now, don't do logical fixing*/
1319
1320 /*This routine could possibly be sped up by using pointers directly
1321 as in the min_cut routine */
1322
1323 /*set up the graph induced by the edges fixed to one*/
1324 lp_net = create_lp_net(cnrp, status, varnum, vars);
1325
1326 verts = lp_net->verts;
1327
1328 compdemands = (double *) calloc (varnum + 1, sizeof(double));
1329
1330 /*get the connected components of the 1-edge graph*/
1331 numchains = cnrp_lp_connected(lp_net, compdemands);
1332
1333 #ifdef ADD_FLOW_VARS
1334 for (i = 0; i < varnum && vars[i]->userind < total_edgenum; i++){
1335 #else
1336 for (i = 0; i < varnum; i++){
1337 #endif
1338 if (!(status[i] & NOT_FIXED) || (status[i] & VARIABLE_BRANCHED_ON))
1339 continue;
1340 v0 = edges[(vars[i]->userind) << 1];
1341 v1 = edges[((vars[i]->userind) << 1) + 1];
1342 if (!v0){
1343 if (verts[0].degree == 2*(cnrp->numroutes)){
1344 /* if the depot has 2*numroutes edges adjacent to it fixed to one,
1345 then we can eliminate all other edges adjacent to the depot
1346 from the problem*/
1347 status[i] = PERM_FIXED_TO_LB;
1348 fixed_num++;
1349 }
1350 }else if ((verts[v0].degree == 2) || (verts[v1].degree == 2)){
1351 /* if a particular node has to fixed-to-one edges adjacent to it,
1352 then we can eliminate all other edges adjacent to that node from
1353 the problem*/
1354 fixed_num++;
1355 status[i] = PERM_FIXED_TO_LB;
1356 }else if (verts[v0].comp == verts[v1].comp){
1357 /*if two vertices are in the same component in the 1-edge graph, then
1358 the edge between them can be eliminated from the problem*/
1359 fixed_num++;
1360 status[i] = PERM_FIXED_TO_LB;
1361 }else if (compdemands[verts[v0].comp] + compdemands[verts[v1].comp]
1362 > capacity){
1363 /*if the sum of the demands in two components of the 1-edge graph is
1364 greater than the capacity of a truck, then these two components
1365 cannot be linked and so any edge that goes bewtween them can be
1366 eliminated from the problem*/
1367 fixed_num++;
1368 status[i] = PERM_FIXED_TO_LB;
1369 }
1370 }
1371
1372 *num_fixed = fixed_num;
1373
1374 free_lp_net(lp_net);
1375
1376 FREE(compdemands);
1377
1378 return(USER_SUCCESS);
1379 }
1380
1381 /*===========================================================================*/
1382
1383 /*===========================================================================*\
1384 * This function generates the 'next' column
1385 \*===========================================================================*/
1386
1387 int user_generate_column(void *user, int generate_what, int cutnum,
1388 cut_data **cuts, int prevind, int nextind,
1389 int *real_nextind, double *colval, int *colind,
1390 int *collen, double *obj, double *lb, double *ub)
1391 {
1392 cnrp_spec *cnrp = (cnrp_spec *)user;
1393 int vh, vl, i;
1394 int total_edgenum = cnrp->vertnum*(cnrp->vertnum-1)/2;
1395
1396 switch (generate_what){
1397 case GENERATE_NEXTIND:
1398 /* Here we just have to generate the specified column. First, we
1399 determine the endpoints */
1400 BOTH_ENDS(nextind, &vh, &vl);
1401 *real_nextind = nextind;
1402 break;
1403 case GENERATE_REAL_NEXTIND:
1404 /* In this case, we have to determine what the "real" next edge is*/
1405 *real_nextind = nextind;
1406 if (prevind >= total_edgenum-1){
1407 *real_nextind = -1;
1408 return(USER_SUCCESS);
1409 }else{
1410 if (nextind == -1) nextind = total_edgenum;
1411 /*first, cycle through the edges that were eliminated in the root*/
1412 for (i = prevind + 1; i < nextind && !cnrp->edges[(i<<1)+1]; i++);
1413 /*now we should have the next nonzero edge*/
1414 vl = cnrp->edges[i << 1];
1415 vh = cnrp->edges[(i << 1) + 1];
1416 }
1417 if (i == nextind)
1418 return(USER_SUCCESS);
1419
1420 *real_nextind = i;
1421 break;
1422 }
1423
1424 /* Now we just have to generate the column corresponding to (vh, vl) */
1425
1426 {
1427 int nzcnt = 0, vertnum = cnrp->vertnum;
1428 char *coef;
1429 cut_data *cut;
1430 int j, size;
1431 int cliquecount = 0, val;
1432 char *clique_array;
1433
1434 colval[0] = 1;
1435 colind[0] = vl; /* supposes vl < vh !!!!!!**********/
1436 colval[1] = 1;
1437 colind[1] = vh;
1438 nzcnt = 2;
1439
1440 /* The coefficient for each row depends on what kind of cut it
1441 is */
1442 for (i = 0; i < cutnum; i++){
1443 coef = (cut = cuts[i])->coef;
1444 switch(cut->type){
1445
1446 case SUBTOUR_ELIM_SIDE:
1447 if (isset(coef, vh) && isset(coef, vl)){
1448 colval[nzcnt] = 1;
1449 colind[nzcnt++] = vertnum + i;
1450 }
1451 break;
1452
1453 case SUBTOUR_ELIM_ACROSS:
1454 /* It's important to have isclr here!!!!! see the macros */
1455 if (isclr(coef, vh) ^ isclr(coef, vl)){
1456 colval[nzcnt] = 1;
1457 colind[nzcnt++] = vertnum + i;
1458 }
1459 break;
1460
1461 case CLIQUE:
1462 size = (vertnum >> DELETE_POWER) + 1;
1463 memcpy(&cliquecount, coef, ISIZE);
1464 coef += ISIZE;
1465 val = 0;
1466 for (j = 0; j < cliquecount; j++){
1467 clique_array = coef + size * j;
1468 if (isset(clique_array, vh) && isset(clique_array, vl))
1469 val += 1;
1470 }
1471 if (val){
1472 colval[nzcnt] = val;
1473 colind[nzcnt++] = vertnum + i;
1474 }
1475 break;
1476
1477 default:
1478 printf("Unrecognized cut type %i!\n", cut->type);
1479 }
1480 }
1481 *collen = nzcnt;
1482 *obj = cnrp->costs[*real_nextind];
1483 }
1484
1485 return(USER_SUCCESS);
1486 }
1487
1488 /*===========================================================================*/
1489
1490 /*===========================================================================*\
1491 * You might want to print some statistics on the types and quantities
1492 * of cuts or something like that.
1493 \*===========================================================================*/
1494
1495 int user_print_stat_on_cuts_added(void *user, int rownum, waiting_row **rows)
1496 {
1497 return(USER_DEFAULT);
1498 }
1499
1500 /*===========================================================================*/
1501
1502 /*===========================================================================*\
1503 * You might want to eliminate rows from the local pool based on
1504 * knowledge of problem structure.
1505 \*===========================================================================*/
1506
1507 int user_purge_waiting_rows(void *user, int rownum, waiting_row **rows,
1508 char *delete_rows)
1509 {
1510 return(USER_DEFAULT);
1511 }
1512
1513 /*===========================================================================*/
1514
1515 /*===========================================================================*\
1516 * The user might want to generate cuts in the LP using information
1517 * about the current tableau, etc. This is for advanced users only.
1518 \*===========================================================================*/
1519
1520 int user_generate_cuts_in_lp(void *user, LPdata *lp_data, int varnum,
1521 var_desc **vars, double *x,
1522 int *new_row_num, cut_data ***cuts)
1523 {
1524 return(GENERATE_CGL_CUTS);
1525 }
1526
1527 /*===========================================================================*/
1528
1529 /*===========================================================================*\
1530 * This function creates a the network of fixed edges that is used in the
1531 * logical fixing routine
1532 \*===========================================================================*/
1533
1534 lp_net *create_lp_net(cnrp_spec *cnrp, char *status, int edgenum,
1535 var_desc **vars)
1536 {
1537 lp_net *n;
1538 lp_net_node *verts;
1539 int nv0 = 0, nv1 = 0;
1540 lp_net_edge *adjlist;
1541 int vertnum = cnrp->vertnum, i;
1542 double *demand = cnrp->demand;
1543 int *edges = cnrp->edges;
1544 #ifdef ADD_FLOW_VARS
1545 int total_edgenum = vertnum*(vertnum-1)/2;
1546 #endif
1547
1548 n = (lp_net *) calloc (1, sizeof(lp_net));
1549 n->vertnum = vertnum;
1550 n->edgenum = vertnum*(vertnum-1)/2;
1551 n->verts = (lp_net_node *) calloc (n->vertnum, sizeof(lp_net_node));
1552 n->adjlist = (lp_net_edge *) calloc (2*(n->edgenum), sizeof(lp_net_edge));
1553 verts = n->verts;
1554 adjlist = n->adjlist;
1555
1556 #ifdef ADD_FLOW_VARS
1557 for (i = 0; i < edgenum && vars[i]->userind < total_edgenum; i++, status++){
1558 #else
1559 for (i = 0; i < edgenum; i++, status++){
1560 #endif
1561 if (*status != PERM_FIXED_TO_UB)
1562 continue;
1563 nv0 = edges[vars[i]->userind << 1];
1564 nv1 = edges[(vars[i]->userind << 1) +1];
1565 if (!verts[nv0].first){
1566 verts[nv0].first = adjlist;
1567 verts[nv0].degree += (int) vars[i]->ub;
1568 }
1569 else{
1570 adjlist->next = verts[nv0].first;
1571 verts[nv0].first = adjlist;
1572 verts[nv0].degree += (int) vars[i]->ub;
1573 }
1574 adjlist->other_end = nv1;
1575 adjlist++;
1576 if (!verts[nv1].first){
1577 verts[nv1].first = adjlist;
1578 verts[nv1].degree += (int) vars[i]->ub;
1579 }
1580 else{
1581 adjlist->next = verts[nv1].first;
1582 verts[nv1].first = adjlist;
1583 verts[nv1].degree += (int) vars[i]->ub;
1584 }
1585 adjlist->other_end = nv0;
1586 adjlist++;
1587 }
1588
1589 for (i=0; i< vertnum; i++)
1590 verts[i].demand = demand[i];
1591
1592 return(n);
1593 }
1594
1595 /*===========================================================================*/
1596
1597 /*===========================================================================*\
1598 * This function constructs the connected components of the 1-edges graph
1599 * used in the logical fixing routine
1600 \*===========================================================================*/
1601
1602 int cnrp_lp_connected(lp_net *n, double *compdemands)
1603 {
1604 int cur_node = 0, cur_comp = 0;
1605 lp_net_node *verts = n->verts;
1606 int vertnum = n->vertnum;
1607 lp_net_edge *cur_edge;
1608 int *nodes_to_scan, num_nodes_to_scan = 0;
1609
1610 nodes_to_scan = (int *) calloc (vertnum, sizeof(int));
1611
1612 while (TRUE){
1613 for (cur_node = 1; cur_node < vertnum; cur_node++)
1614 if (!verts[cur_node].comp){ /* Look for the first node not already
1615 in a component */
1616 break;
1617 }
1618
1619 if (cur_node == n->vertnum) break; /* we are done */
1620
1621 nodes_to_scan[num_nodes_to_scan++] = cur_node;
1622 /* add the cur_node to the list of nodes to be scanned */
1623 verts[cur_node].comp = ++cur_comp;
1624 /* add the current node to the current component */
1625 compdemands[cur_comp] = verts[cur_node].demand;
1626 while(TRUE){
1627 /* In each iteration of this loop, we take the next node off the
1628 list of nodes to be scanned, add it to the current component, and
1629 then add all its neighbors to the list of nodes to be scanned */
1630 for (cur_node = nodes_to_scan[--num_nodes_to_scan],
1631 verts[cur_node].scanned = TRUE,
1632 cur_edge = verts[cur_node].first,
1633 cur_comp = verts[cur_node].comp;
1634 cur_edge; cur_edge = cur_edge->next){
1635 if (cur_edge->other_end){
1636 if (!verts[cur_edge->other_end].comp){
1637 verts[cur_edge->other_end].comp = cur_comp;
1638 compdemands[cur_comp] += verts[cur_edge->other_end].demand;
1639 nodes_to_scan[num_nodes_to_scan++] = cur_edge->other_end;
1640 }
1641 }
1642 }
1643 if (!num_nodes_to_scan) break;
1644 /* when there are no more nodes to scan, we start a new component */
1645 }
1646 }
1647
1648 free((char *) nodes_to_scan);
1649 return(cur_comp);
1650 }
1651
1652 /*===========================================================================*/
1653
1654 /*===========================================================================*\
1655 * Free the data structures associated with the 1-edges graph
1656 \*===========================================================================*/
1657
1658 void free_lp_net(lp_net *n)
1659 {
1660 if (n){
1661 FREE(n->adjlist);
1662 FREE(n->verts);
1663 FREE(n);
1664 }
1665 }
1666
1667 /*===========================================================================*/
1668
1669 #if 0
1670 char construct_feasible_solution(cnrp_spec *cnrp, network *n,
1671 double *true_objval, double etol,
1672 char branching)
1673 {
1674 _node *tour = cnrp->cur_sol;
1675 int cur_vert = 0, prev_vert = 0, cur_route, i, count;
1676 elist *cur_route_start = NULL;
1677 edge *edge_data;
1678 vertex *verts = n->verts;
1679 double fixed_cost = 0.0, variable_cost = 0.0;
1680 int cuts = 0;
1681 char print_solution = FALSE;
1682 char continue_with_node = FALSE;
1683
1684 #ifdef MULTI_CRITERIA
1685 for (i = 0; i < n->edgenum; i++){
1686 fixed_cost += cnrp->costs[INDEX(n->edges[i].v0, n->edges[i].v1)];
1687 #ifdef ADD_FLOW_VARS
1688 variable_cost += (n->edges[i].flow1+n->edges[i].flow2)*
1689 cnrp->costs[INDEX(n->edges[i].v0, n->edges[i].v1)];
1690 #endif
1691 }
1692
1693 #if 0
1694 *true_objval -= cnrp->par.rho*(fixed_cost+variable_cost);
1695 #endif
1696
1697 if (cnrp->ub > 0 &&
1698 *true_objval-cnrp->par.rho*(fixed_cost+variable_cost) > cnrp->ub+etol){
1699 return(FALSE);
1700 }
1701
1702 if (cnrp->par.gamma == 1.0){
1703 if (fixed_cost < cnrp->fixed_cost + etol){
1704 if (fixed_cost < cnrp->fixed_cost - etol ||
1705 (fixed_cost >= cnrp->fixed_cost - etol
1706 && variable_cost < cnrp->variable_cost - etol)){
1707 printf("\nBetter Solution Found:\n");
1708 #ifdef ADD_FLOW_VARS
1709 printf("Solution Fixed Cost: %.1f\n", fixed_cost);
1710 printf("Solution Variable Cost: %.1f\n", variable_cost);
1711 #else
1712 printf("Solution Cost: %.0f\n", fixed_cost);
1713 #endif
1714 cnrp->variable_cost = variable_cost;
1715 cnrp->fixed_cost = fixed_cost;
1716 cnrp->ub = *true_objval-cnrp->par.rho*(fixed_cost+variable_cost);
1717 print_solution = TRUE;
1718 }
1719 /* Add an optimality cut for the second objective */
1720 if (!branching){
1721 cut_data *new_cut = (cut_data *) calloc(1, sizeof(cut_data));
1722 new_cut->coef = NULL;
1723 new_cut->rhs = (int) (variable_cost + etol) - 1;
1724 new_cut->size = 0;
1725 new_cut->type = OPTIMALITY_CUT_VARIABLE;
1726 new_cut->name = CUT__DO_NOT_SEND_TO_CP;
1727 continue_with_node = cg_send_cut(new_cut);
1728 FREE(new_cut);
1729 }else{
1730 continue_with_node = TRUE;
1731 }
1732 }
1733 }else if (cnrp->par.tau == 1.0){
1734 if (variable_cost < cnrp->variable_cost + etol){
1735 if (variable_cost < cnrp->variable_cost - etol ||
1736 (variable_cost >= cnrp->variable_cost - etol
1737 && fixed_cost < cnrp->fixed_cost - etol)){
1738 printf("\nBetter Solution Found:\n");
1739 #ifdef ADD_FLOW_VARS
1740 printf("Solution Fixed Cost: %.1f\n", fixed_cost);
1741 printf("Solution Variable Cost: %.1f\n", variable_cost);
1742 #else
1743 printf("Solution Cost: %.0f\n", fixed_cost);
1744 #endif
1745 cnrp->variable_cost = variable_cost;
1746 cnrp->fixed_cost = fixed_cost;
1747 cnrp->ub = *true_objval-cnrp->par.rho*(fixed_cost+variable_cost);
1748 print_solution = TRUE;
1749 }
1750 /* Add an optimality cut for the second objective */
1751 if (!branching){
1752 cut_data *new_cut = (cut_data *) calloc(1, sizeof(cut_data));
1753 new_cut->coef = NULL;
1754 new_cut->rhs = (int) (fixed_cost + etol) - 1;
1755 new_cut->size = 0;
1756 new_cut->type = OPTIMALITY_CUT_FIXED;
1757 new_cut->name = CUT__DO_NOT_SEND_TO_CP;
1758 continue_with_node = cg_send_cut(new_cut);
1759 FREE(new_cut);
1760 }else{
1761 continue_with_node = TRUE;
1762 }
1763 }
1764 }else{
1765 if ((*true_objval-cnrp->par.rho*(fixed_cost+variable_cost) <
1766 cnrp->ub - etol) ||
1767 (fixed_cost < cnrp->fixed_cost - etol &&
1768 variable_cost < cnrp->variable_cost + etol) ||
1769 (variable_cost < cnrp->variable_cost - etol &&
1770 fixed_cost < cnrp->fixed_cost + etol)){
1771 printf("\nBetter Solution Found:\n");
1772 #ifdef ADD_FLOW_VARS
1773 printf("Solution Fixed Cost: %.1f\n", fixed_cost);
1774 printf("Solution Variable Cost: %.1f\n", variable_cost);
1775 #else
1776 printf("Solution Cost: %.0f\n", fixed_cost);
1777 #endif
1778 cnrp->variable_cost = variable_cost;
1779 cnrp->fixed_cost = fixed_cost;
1780 cnrp->ub = *true_objval-cnrp->par.rho*(fixed_cost+variable_cost);
1781 print_solution = TRUE;
1782 }
1783 if (!branching){
1784 if (cnrp->par.gamma*(fixed_cost - cnrp->utopia_fixed) >
1785 *true_objval-cnrp->par.rho*(fixed_cost+variable_cost)-etol){
1786 /* Add an optimality cut for the second objective */
1787 cut_data *new_cut = (cut_data *) calloc(1, sizeof(cut_data));
1788 new_cut->coef = NULL;
1789 new_cut->rhs = (int) (variable_cost + etol) - 1;
1790 new_cut->size = 0;
1791 new_cut->type = OPTIMALITY_CUT_VARIABLE;
1792 new_cut->name = CUT__DO_NOT_SEND_TO_CP;
1793 continue_with_node = cg_send_cut(new_cut);
1794 FREE(new_cut);
1795 }else{
1796 /* Add an optimality cut for the second objective */
1797 cut_data *new_cut = (cut_data *) calloc(1, sizeof(cut_data));
1798 new_cut->coef = NULL;
1799 new_cut->rhs = (int) (fixed_cost + etol) - 1;
1800 new_cut->size = 0;
1801 new_cut->type = OPTIMALITY_CUT_FIXED;
1802 new_cut->name = CUT__DO_NOT_SEND_TO_CP;
1803 continue_with_node = cg_send_cut(new_cut);
1804 FREE(new_cut);
1805 }
1806 }else{
1807 continue_with_node = TRUE;
1808 }
1809 }
1810 #endif
1811
1812 if (!print_solution){
1813 return(continue_with_node);
1814 }
1815
1816
1817 if (cnrp->par.prob_type == TSP || cnrp->par.prob_type == VRP ||
1818 cnrp->par.prob_type == BPP){
1819 /*construct the tour corresponding to this solution vector*/
1820 for (cur_route_start = verts[0].first, cur_route = 1,
1821 edge_data = cur_route_start->data; cur_route <= cnrp->numroutes;
1822 cur_route++){
1823 edge_data = cur_route_start->data;
1824 edge_data->scanned = TRUE;
1825 cur_vert = edge_data->v1;
1826 tour[prev_vert].next = cur_vert;
1827 tour[cur_vert].route = cur_route;
1828 prev_vert = 0;
1829 while (cur_vert){
1830 if (verts[cur_vert].first->other_end != prev_vert){
1831 prev_vert = cur_vert;
1832 edge_data = verts[cur_vert].first->data;
1833 cur_vert = verts[cur_vert].first->other_end;
1834 }
1835 else{
1836 prev_vert = cur_vert;
1837 edge_data = verts[cur_vert].last->data; /*This statement
1838 could possibly
1839 be taken out to speed
1840 things up a bit*/
1841 cur_vert = verts[cur_vert].last->other_end;
1842 }
1843 tour[prev_vert].next = cur_vert;
1844 tour[cur_vert].route = cur_route;
1845 }
1846 edge_data->scanned = TRUE;
1847
1848 while (cur_route_start->data->scanned){
1849 if (!(cur_route_start = cur_route_start->next_edge)) break;
1850 }
1851 }
1852
1853 /* Display the solution */
1854
1855 cur_vert = tour[0].next;
1856
1857 if (tour[0].route == 1)
1858 printf("\n0 ");
1859 while (cur_vert != 0){
1860 if (tour[prev_vert].route != tour[cur_vert].route){
1861 printf("\nRoute #%i: ", tour[cur_vert].route);
1862 count = 0;
1863 }
1864 printf("%i ", cur_vert);
1865 count++;
1866 if (count > 15){
1867 printf("\n");
1868 count = 0;
1869 }
1870 prev_vert = cur_vert;
1871 cur_vert = tour[cur_vert].next;
1872 }
1873 printf("\n\n");
1874 }else{
1875 for (i = 0; i < n->edgenum; i++){
1876 cnrp->cur_sol_tree[i] = INDEX(n->edges[i].v0, n->edges[i].v1);
1877 }
1878
1879 /* Display the solution */
1880
1881 for (i = 0; i < n->edgenum; i++){
1882 printf("%i %i\n", n->edges[i].v0, n->edges[i].v1);
1883 }
1884 }
1885
1886 return(continue_with_node);
1887 }
1888 #endif
1889
1890 /*===========================================================================*/
1891
1892 /*__BEGIN_EXPERIMENTAL_SECTION__*/
1893 #ifdef TRACE_PATH
1894
1895 #include "sym_lp.h"
1896
1897 void check_lp(lp_prob *p)
1898 {
1899 LPdata *lp_data = p->lp_data;
1900 int i, j, l;
1901 tm_prob *tm = p->tm;
1902 double *x = (double *) malloc(tm->feas_sol_size * DSIZE), lhs, cost = 0;
1903 double *lhs_totals = (double *) calloc(lp_data->m, DSIZE);
1904 MakeMPS(lp_data, 0, 0);
1905
1906 get_x(lp_data);
1907
1908 printf("Optimal Fractional Solution: %.10f\n", lp_data->lpetol);
1909 for (i = 0; i < lp_data->n; i++){
1910 if (lp_data->x[i] > lp_data->lpetol){
1911 printf("uind: %i colind: %i value: %.10f cost: %f\n",
1912 lp_data->vars[i]->userind, i, lp_data->x[i], lp_data->obj[i]);
1913 }
1914 cost += lp_data->obj[i]*lp_data->x[i];
1915 }
1916 printf("Cost: %f\n", cost);
1917
1918 printf("\nFeasible Integer Solution:\n");
1919 for (cost = 0, i = 0; i < tm->feas_sol_size; i++){
1920 printf("uind: %i ", tm->feas_sol[i]);
1921 for (j = 0; j < lp_data->n; j++){
1922 if (lp_data->vars[j]->userind == tm->feas_sol[i]){
1923 cost += lp_data->obj[j];
1924 printf("colind: %i lb: %f ub: %f obj: %f\n", j, lp_data->lb[j],
1925 lp_data->ub[j], lp_data->obj[j]);
1926 break;
1927 }
1928 }
1929 if (j == lp_data->n)
1930 printf("\n\nERROR!!!!!!!!!!!!!!!!\n\n");
1931 x[i] = 1.0;
1932 }
1933 printf("Cost: %f\n", cost);
1934
1935 printf("\nChecking LP....\n\n");
1936 printf("Number of cuts: %i\n", lp_data->m);
1937 for (i = 0; i < tm->feas_sol_size; i++){
1938 for (j = 0; j < lp_data->n; j++){
1939 if (tm->feas_sol[i] == lp_data->vars[j]->userind)
1940 break;
1941 }
1942 for (l = lp_data->matbeg[j]; l < lp_data->matbeg[j] + lp_data->matcnt[j];
1943 l++){
1944 lhs_totals[lp_data->matind[l]] += 1;
1945 }
1946 }
1947 for (i = 0; i < p->base.cutnum; i++){
1948 printf("Cut %i: %f %c %f\n", i, lhs_totals[i], lp_data->sense[i],
1949 lp_data->rhs[i]);
1950 }
1951 for (; i < lp_data->m; i++){
1952 lhs = compute_lhs(tm->feas_sol_size, tm->feas_sol, x,
1953 lp_data->rows[i].cut, p->base.cutnum);
1954 printf("Cut %i: %f %f %c %f\n", i, lhs_totals[i], lhs, lp_data->sense[i],
1955 lp_data->rhs[i]);
1956 if (lp_data->rows[i].cut->sense == 'G' ?
1957 lhs < lp_data->rows[i].cut->rhs : lhs > lp_data->rows[i].cut->rhs){
1958 printf("LP: ERROR -- row is violated by feasible solution!!!\n");
1959 sleep(600);
1960 exit(1);
1961 }
1962 }
1963 }
1964 #endif
1965 /*___END_EXPERIMENTAL_SECTION___*/
1966