1 /*===========================================================================*/
2 /* */
3 /* This file is part of the SYMPHONY MILP Solver Framework. */
4 /* */
5 /* SYMPHONY was jointly developed by Ted Ralphs (ted@lehigh.edu) and */
6 /* Laci Ladanyi (ladanyi@us.ibm.com). */
7 /* */
8 /* (c) Copyright 2000-2019 Ted Ralphs. All Rights Reserved. */
9 /* */
10 /* The OSI interface in this file was written by Menal Guzelsoy. */
11 /* The OSL interface was written by Ondrej Medek. */
12 /* */
13 /* This software is licensed under the Eclipse Public License. Please see */
14 /* accompanying file for terms. */
15 /* */
16 /*===========================================================================*/
17
18 #include <stdlib.h> /* free() is here on AIX ... */
19 #include <math.h>
20 //#include <string.h>
21
22 #include "sym_lp_solver.h"
23 #include "sym_constants.h"
24 #include "sym_macros.h"
25 #include "sym_qsort.h"
26
27 #ifdef PRINT
28 #undef PRINT
29 #endif
30 #define PRINT(a, b, c) \
31 if ((a) > (b)) printf c
32
33 /*===========================================================================*/
34
35 /*===========================================================================*\
36 * This file contains the interface with the LP Solver.
37 * The first few routines are independent of what LP solver is being used.
38 \*===========================================================================*/
39
dot_product(double * val,int * ind,int collen,double * col)40 double dot_product(double *val, int *ind, int collen, double *col)
41 {
42 const int* lastind = ind + collen;
43 double prod = 0;
44 while (ind != lastind)
45 prod += (*val++) * col[*ind++];
46 return(prod);
47 }
48
49 /*===========================================================================*/
50
free_lp_arrays(LPdata * lp_data)51 void free_lp_arrays(LPdata *lp_data)
52 {
53 FREE(lp_data->not_fixed);
54 FREE(lp_data->status);
55 FREE(lp_data->x);
56 FREE(lp_data->dj);
57 FREE(lp_data->dualsol);
58 FREE(lp_data->slacks);
59 FREE(lp_data->random_hash);
60 FREE(lp_data->hashes);
61 FREE(lp_data->accepted_ind);
62 FREE(lp_data->heur_solution);
63 FREE(lp_data->col_solution);
64 #ifdef __CPLEX__
65 FREE(lp_data->lb);
66 FREE(lp_data->ub);
67 #endif
68 FREE(lp_data->vars);
69 FREE(lp_data->tmp.c);
70 FREE(lp_data->tmp.i1);
71 FREE(lp_data->tmp.i2);
72 FREE(lp_data->tmp.d);
73 FREE(lp_data->tmp.p1);
74 FREE(lp_data->tmp.p2);
75 FREE(lp_data->tmp.cv);
76 FREE(lp_data->tmp.iv);
77 FREE(lp_data->tmp.dv);
78
79 FREE(lp_data->tmp1.i1);
80 FREE(lp_data->tmp1.d);
81 FREE(lp_data->tmp1.c);
82 FREE(lp_data->tmp2.i1);
83 FREE(lp_data->tmp2.d);
84 FREE(lp_data->tmp2.c);
85 }
86
87 /*===========================================================================*/
88
free_mip_desc(MIPdesc * mip)89 void free_mip_desc(MIPdesc *mip)
90 {
91 int j, n = 0;
92
93 FREE(mip->matbeg);
94 FREE(mip->matind);
95 FREE(mip->matval);
96 FREE(mip->col_lengths);
97 FREE(mip->row_matbeg);
98 FREE(mip->row_matind);
99 FREE(mip->row_matval);
100 FREE(mip->row_lengths);
101 FREE(mip->orig_sense);
102 FREE(mip->orig_ind);
103
104 FREE(mip->obj);
105 FREE(mip->obj1);
106 FREE(mip->obj2);
107 FREE(mip->rhs);
108 FREE(mip->rngval);
109 FREE(mip->sense);
110 FREE(mip->lb);
111 FREE(mip->ub);
112 FREE(mip->is_int);
113 if (mip->colname){
114 n = mip->n;
115 if(mip->alloc_n > n){
116 n = mip->alloc_n;
117 }
118 for (j = 0; j < n; j++){
119 FREE(mip->colname[j]);
120 }
121 FREE(mip->colname);
122 }
123 #if 0
124 if(mip->fixed_name){
125 for(j = 0; j < mip->fixed_n; j++){
126 FREE(mip->fixed_name[j]);
127 }
128 FREE(mip->fixed_name);
129 }
130 #endif
131
132 if(mip->fixed_n){
133 FREE(mip->fixed_val);
134 FREE(mip->fixed_ind);
135 }
136
137 if(mip->aggr_n){
138 FREE(mip->aggr_ind);
139 FREE(mip->aggr_to_ind);
140 }
141
142 if(mip->subs_n){
143 FREE(mip->subs_ind);
144 FREE(mip->subs_aval);
145 FREE(mip->subs_rhs);
146 FREE(mip->subs_rbeg);
147 FREE(mip->subs_rind);
148 FREE(mip->subs_rval);
149 }
150
151 if(mip->cru_vars_num){
152 FREE(mip->cru_vars);
153 }
154
155 if(mip->mip_inf){
156 FREE(mip->mip_inf->c_ind);
157 FREE(mip->mip_inf->c_val);
158 FREE(mip->mip_inf->c_beg);
159 FREE(mip->mip_inf->c_sense);
160 FREE(mip->mip_inf->c_rhs);
161 FREE(mip->mip_inf->c_tmp);
162 FREE(mip->mip_inf->rows);
163 FREE(mip->mip_inf->cols);
164 FREE(mip->mip_inf);
165 }
166 }
167
168 /*===========================================================================*/
169
size_lp_arrays(LPdata * lp_data,char do_realloc,char set_max,int row_num,int col_num,int nzcnt)170 void size_lp_arrays(LPdata *lp_data, char do_realloc, char set_max,
171 int row_num, int col_num, int nzcnt)
172 {
173 char resize_m = FALSE;
174 char resize_n = FALSE;
175 int maxm, maxn, maxnz, maxmax;
176
177 if (set_max){
178 maxm = row_num;
179 maxn = col_num;
180 maxnz = nzcnt;
181 }else{
182 maxm = lp_data->m + row_num;
183 maxn = lp_data->n + col_num;
184 maxnz = lp_data->nz + nzcnt;
185 }
186
187 if (maxm > lp_data->maxm){
188 resize_m = TRUE;
189 lp_data->maxm = maxm + (set_max ? 0 : (int)(BB_BUNCH));
190 if (! do_realloc){
191 FREE(lp_data->dualsol);
192 lp_data->dualsol = (double *) malloc(lp_data->maxm * DSIZE);
193 FREE(lp_data->slacks);
194 lp_data->slacks = (double *) malloc(lp_data->maxm * DSIZE);
195 }else{
196 lp_data->dualsol = (double *) realloc((char *)lp_data->dualsol,
197 lp_data->maxm * DSIZE);
198 lp_data->slacks = (double *) realloc((void *)lp_data->slacks,
199 lp_data->maxm * DSIZE);
200 }
201 /* rows is realloc'd in either case just to keep the base constr */
202 lp_data->rows = (row_data *) realloc((char *)lp_data->rows,
203 lp_data->maxm*sizeof(row_data));
204 }
205 if (maxn > lp_data->maxn){
206 // int oldmaxn = MAX(lp_data->maxn, lp_data->n);
207 resize_n = TRUE;
208 lp_data->maxn = maxn + (set_max ? 0 : 5 * (int)(BB_BUNCH));
209 if (! do_realloc){
210 FREE(lp_data->x);
211 lp_data->x = (double *) malloc(lp_data->maxn * DSIZE);
212 FREE(lp_data->dj);
213 lp_data->dj = (double *) malloc(lp_data->maxn * DSIZE);
214 FREE(lp_data->status);
215 lp_data->status = (char *) malloc(lp_data->maxn * CSIZE);
216 FREE(lp_data->random_hash);
217 lp_data->random_hash = (double *) malloc(lp_data->maxn * DSIZE);
218 FREE(lp_data->heur_solution);
219 lp_data->heur_solution = (double *) malloc(lp_data->maxn * DSIZE);
220 FREE(lp_data->col_solution);
221 lp_data->col_solution = (double *) malloc(lp_data->maxn * DSIZE);
222 #ifdef __CPLEX__
223 FREE(lp_data->lb);
224 lp_data->lb = (double *) malloc(lp_data->maxn * DSIZE);
225 FREE(lp_data->ub);
226 lp_data->ub = (double *) malloc(lp_data->maxn * DSIZE);
227 #endif
228 }else{
229 lp_data->x = (double *) realloc((char *)lp_data->x,
230 lp_data->maxn * DSIZE);
231 lp_data->dj = (double *) realloc((char *)lp_data->dj,
232 lp_data->maxn * DSIZE);
233 lp_data->status = (char *) realloc((char *)lp_data->status,
234 lp_data->maxn * CSIZE);
235 lp_data->random_hash = (double *) realloc((char *)lp_data->random_hash,
236 lp_data->maxn * DSIZE);
237 lp_data->heur_solution = (double *) realloc((char *)
238 lp_data->heur_solution, lp_data->maxn * DSIZE);
239 lp_data->col_solution = (double *) realloc((char *)
240 lp_data->col_solution, lp_data->maxn * DSIZE);
241 #ifdef __CPLEX__
242 lp_data->lb = (double *) realloc((char *)lp_data->lb,
243 lp_data->maxn * DSIZE);
244 lp_data->ub = (double *) realloc((char *)lp_data->ub,
245 lp_data->maxn * DSIZE);
246 #endif
247 }
248 }
249 if (maxnz > lp_data->maxnz){
250 lp_data->maxnz = maxnz + (set_max ? 0 : 20 * (int)(BB_BUNCH));
251 }
252
253 /* re(m)alloc the tmp arrays */
254 if (resize_m || resize_n){
255 temporary *tmp = &lp_data->tmp;
256 maxm = lp_data->maxm;
257 maxn = lp_data->maxn;
258 maxmax = MAX(maxm, maxn);
259 /* anything with maxm and maxn in it has to be resized */
260 FREE(tmp->c);
261 FREE(tmp->i1);
262 FREE(tmp->d);
263 tmp->c = (char *) malloc(CSIZE * 4 * maxmax);
264 tmp->i1 = (int *) malloc(ISIZE * MAX(4*maxm, 4*maxn + 1));
265 tmp->d = (double *) malloc(DSIZE * 4 * maxmax);
266 /* These have to be resized only if maxm changes */
267 if (resize_m){
268 FREE(tmp->i2);
269 FREE(tmp->p1);
270 FREE(tmp->p2);
271 tmp->i2 = (int *) malloc(2*maxmax * ISIZE);
272 tmp->p1 = (void **) malloc(maxm * sizeof(void *));
273 tmp->p2 = (void **) malloc(maxm * sizeof(void *));
274 }
275 }
276 }
277
278 #ifdef USE_GLPMPL
279
280 /*This function reads in the GNU MathProg model file and returns either 1 if
281 it is succeded or 0 otherwise.*/
282
read_gmpl(MIPdesc * mip,char * modelfile,char * datafile,char * probname)283 int read_gmpl(MIPdesc *mip, char *modelfile, char *datafile, char *probname)
284 {
285 glp_tran *tran;
286 glp_prob *prob;
287 int errors;
288 int i, j, k, length, type, nonzeros;
289 double *matval;
290 int *matind;
291 int * matbeg;
292
293 double *row_lb;
294 double *row_ub;
295
296 int *indices;
297 double *values;
298 double inf = MAXDOUBLE;//SYM_INFINITY;
299
300 tran = glp_mpl_alloc_wksp(); /* initialize the translator */
301
302 /*if there are some errors in reading the file(s): then errors != 0*/
303 if (glp_mpl_read_model(tran, modelfile, FALSE)){
304 printf("\nError in reading the model (or data) file!");
305 glp_mpl_free_wksp(tran); /* free all the mpl related stuff */
306 return(0);
307 }
308
309 /*if the data is not in the model file and will be given seperately,
310 then errors=1! */
311 if (glp_mpl_read_data(tran, datafile)){
312 printf("\nError in reading the model (or data) file!");
313 glp_mpl_free_wksp(tran); /* free all the mpl related stuff */
314 return(0);
315 }
316
317 /*Generate the model variables, constraints, objective by storing in the
318 translator database.It is possible to capture the messages in a file by
319 passing the filename instead of NULL.*/
320
321 if (glp_mpl_generate(tran, NULL)){
322 printf("\nError in generating the model!");
323 glp_mpl_free_wksp(tran); /* free all the mpl related stuff */
324 return(0);
325 }
326
327 prob = glp_create_prob();
328
329 glp_mpl_build_prob(tran, prob);
330
331 strncpy(probname, glp_get_prob_name(prob), 80); /* name the problem */
332
333 /* get num of rows and cols */
334 mip->m = glp_get_num_rows(prob)-1; /* subtract the objective row */
335 mip->n = glp_get_num_cols(prob);
336 mip->nz = 0; /* for now... */
337
338 /*Indices and values of nonzeros will return beginning with indices[1] and
339 values[1]. Also note that row and column indices begin with 1 in glpmpl*/
340
341 /*get mip->nz and mip->obj*/
342 mip->obj = (double *) calloc(DSIZE, mip->n);
343 mip->obj1 = NULL;
344 mip->obj2 = NULL;
345
346 indices = (int *) malloc(ISIZE * (mip->n + 1));
347 values = (double *) malloc(DSIZE * (mip->n + 1));
348
349 if (glp_get_obj_dir(prob) == GLP_MIN){
350 mip->obj_sense = SYM_MINIMIZE;
351 for (int i = 0; i < mip->n; i++){
352 mip->obj[i] = glp_get_obj_coef(prob, i+1);
353 }
354 }else{
355 mip->obj_sense = SYM_MAXIMIZE;
356 for (int i = 0; i < mip->n; i++){
357 mip->obj[i] = -glp_get_obj_coef(prob, i+1);
358 }
359 }
360 mip->obj_offset = glp_get_obj_coef(prob, 0);
361
362 for(i = 1; i < mip->m + 1; i++){
363 mip->nz += glp_get_mat_row(prob, i+1, NULL, NULL);
364 }
365
366 /* Define a row ordered dummy constraint matrix since glpmpl returns the
367 constraint definitions as row ordered, we will change its order later. */
368
369 /* fill the dummy matbeg, matind, matval, row_lb and row_ub arrays */
370 matbeg = (int *) malloc(ISIZE * (mip->m + 1));
371 matind = (int *) malloc(ISIZE * mip->nz);
372 matval = (double *) malloc(DSIZE * mip->nz);
373
374 row_ub = (double *) malloc(DSIZE * mip->m);
375 row_lb = (double *) malloc(DSIZE * mip->m);
376
377 matbeg[0] = 0;
378 nonzeros = 0;
379 for(i = 1, k = 0; i < mip->m + 1; i++){
380 /* read the nonzeros in row i+1 */
381 length = glp_get_mat_row(prob, i+1, indices, values);
382 /* get the row bounds. we use k instead of i since we have the obj
383 row somewhere. */
384 row_lb[k] = glp_get_row_lb(prob, i+1);
385 row_ub[k] = glp_get_row_ub(prob, i+1);
386 type = glp_get_row_type(prob, i+1);
387 switch(type){
388 case GLP_FR: /* free */
389 row_lb[k] = -inf;
390 row_ub[k] = inf;
391 break;
392 case GLP_LO: /* has lower bound */
393 row_ub[k] = inf;
394 break;
395 case GLP_UP: /* has upper bound */
396 row_lb[k] = -inf;
397 break;
398 default: /* is bounded from both sides or is an equality */
399 break;
400 }
401 for (j = 0; j < length; j++){
402 matind[matbeg[k]+j] = indices[j+1] - 1;
403 matval[matbeg[k]+j] = values[j+1];
404 }
405 nonzeros += length;
406 k++;
407 matbeg[k] = nonzeros;
408 }
409
410 /* fill the column related definitions: ub, lb, is_int and colname arrays */
411
412 mip->ub = (double *) malloc(DSIZE * mip->n);
413 mip->lb = (double *) malloc(DSIZE * mip->n);
414 mip->is_int = (char *) calloc(CSIZE, mip->n);
415 mip->colname = (char **) malloc(sizeof(char *) * mip->n);
416
417 for (j = 0; j < mip->n; j++){
418 mip->lb[j] = glp_get_col_lb(prob, j+1);
419 mip->ub[j] = glp_get_col_ub(prob, j+1);
420 type = glp_get_col_type(prob, j+1);
421 switch(type){
422 case GLP_FR: /* free */
423 mip->lb[j] = -inf;
424 mip->ub[j] = inf;
425 break;
426 case GLP_LO: /* has lower bound */
427 mip->ub[j] = inf;
428 break;
429 case GLP_UP: /* has upper bound */
430 mip->lb[j] = -inf;
431 break;
432 default: /* has both lower and upper bound or is a fixed variable */
433 break;
434 }
435
436 type = glp_get_col_kind(prob, j+1);
437 if(type == GLP_IV || type == GLP_BV){
438 mip->is_int[j] = TRUE;
439 }
440 /* bounds for binary variables were probably not assigned.
441 So assign them! */
442 if(type == GLP_BV){
443 mip->ub[j] = 1.0;
444 mip->lb[j] = 0.0;
445
446 }
447
448 mip->colname[j] = (char *) malloc(CSIZE * MAX_NAME_SIZE);
449 strncpy(mip->colname[j], glp_get_col_name(prob, j+1), MAX_NAME_SIZE);
450 mip->colname[j][MAX_NAME_SIZE-1] = 0; /* ??? */
451 }
452
453 /*load the definitions to a CoinPackedMatrix as row ordered and get the
454 column ordered matrix after reversing its order in order to fill the
455 matrix definitons as column ordered*/
456
457 mip->matbeg = (int *) calloc(ISIZE, (mip->n + 1));
458 mip->matval = (double *) malloc(DSIZE * mip->nz);
459 mip->matind = (int *) malloc(ISIZE * mip->nz);
460
461 #if 0
462 //make CoinPackedMatrix help us for now!!!
463 CoinPackedMatrix matrixByCol (false, mip->n,
464 mip->m, mip->nz, matval, matind, matbeg, 0);
465 matrixByCol.reverseOrdering();
466
467
468 memcpy(mip->matbeg, matrixByCol.getVectorStarts(), ISIZE * (mip->n + 1));
469 memcpy(mip->matval, matrixByCol.getElements(), DSIZE * mip->nz);
470 memcpy(mip->matind, matrixByCol.getIndices(), ISIZE * mip->nz);
471 #endif
472
473 /* what if the user doesn't have COIN, is that possible?:) */
474 nonzeros = 0;
475 for(j = 0; j < mip->n; j++){
476 for(i = 0; i < mip->m; i++){
477 for(k = matbeg[i]; k < matbeg[i+1]; k++){
478 if(matind[k] == j){
479 mip->matind[nonzeros] = i;
480 mip->matval[nonzeros] = matval[k];
481 nonzeros++;
482 break;
483 }
484 }
485 }
486 mip->matbeg[j+1] = nonzeros;
487 }
488
489 /*get the other definitions: rhs, sense and rngval from row_lb and row_ub*/
490
491 mip->rhs = (double *) malloc(DSIZE * mip->m);
492 mip->sense = (char *) malloc(CSIZE * mip->m);
493 mip->rngval = (double *) malloc(DSIZE * mip->m);
494
495 /* convertBoundToSense: stolen from COIN :) */
496 for(i = 0; i < mip->m; i++) {
497 mip->rngval[i] = 0.0;
498 if (row_lb[i] > -inf) {
499 if (row_ub[i] < inf) {
500 mip->rhs[i] = row_ub[i];
501 if (row_lb[i] == row_ub[i]) {
502 mip->sense[i] = 'E';
503 } else {
504 mip->sense[i] = 'R';
505 mip->rngval[i] = row_ub[i] - row_lb[i];
506 }
507 }else{
508 mip->sense[i] = 'G';
509 mip->rhs[i] = row_lb[i];
510 }
511 }else{
512 if (row_ub[i] < inf) {
513 mip->sense[i] = 'L';
514 mip->rhs[i] = row_ub[i];
515 }else{
516 mip->sense[i] = 'N';
517 mip->rhs[i] = 0.0;
518 }
519 }
520 }
521
522 FREE(matind);
523 FREE(matval);
524 FREE(matbeg);
525 FREE(row_lb);
526 FREE(row_ub);
527
528 /* if you could reach here by chance, then you are safe anymore:) */
529 return(1);
530 }
531 #endif /* USE_GLPMPL */
532
533 #ifdef __OSL__
534
535 /*****************************************************************************/
536 /*****************************************************************************/
537 /******* *******/
538 /******* routines when OSL is used *******/
539 /******* *******/
540 /******* WARNING! Not well tested. Please, report bugs. *******/
541 /*****************************************************************************/
542 /*****************************************************************************/
543
544 /*============================================================================
545 * - no fastmip is used
546 * - no scaling is used - cannot test it
547 * - possible problems with getting reduced costs
548 * - LPdata->tmp field mostly not used for safe. malloc and free for temporary
549 * fields are used instead => slow down
550 *============================================================================
551 */
552
553 static int osllib_status;
554
555 #include <memory.h>
556
OSL_check_error(const char * erring_func)557 void OSL_check_error(const char *erring_func)
558 {
559 if (osllib_status){
560 printf("!!! OSL status is nonzero !!! [%s, %i]\n",
561 erring_func, osllib_status);
562 }
563 }
564
565 /*===========================================================================*/
open_lp_solver(LPdata * lp_data)566 void open_lp_solver(LPdata *lp_data)
567 {
568 EKKModel *baseModel;
569
570 lp_data->env = ekk_initializeContext();
571 osllib_status = (lp_data->env == NULL);
572 OSL_check_error("open_lp_solver - ekk_initializeContext");
573 baseModel = ekk_baseModel(lp_data->env);
574 osllib_status = (baseModel == NULL);
575 OSL_check_error("open_lp_solver - ekk_baseModel");
576 ekk_setDebug(baseModel, -1, 0);
577 ekk_setIloglevel(baseModel, 2);
578 /* 1 - 2999 informational messsages
579 3000 - 5999 warn
580 6000 - 6999 error, but keep running
581 7000 - 8999 error and stop running */
582 /* osllib_status = ekk_messagesPrintOn(baseModel, 1, 8999); */
583 /* OSL_check_error("open_lp_solver - ekk_messagePrintOn"); */
584 osllib_status = ekk_messagesPrintOff(baseModel, 1, 5999);
585 OSL_check_error("open_lp_solver - ekk_messagePrintOff");
586
587 /* default is to minimize */
588 /* osllib_status = ekk_setMinimize(baseModel); */
589 /* OSL_check_error("open_lp_solver - ekk_setMinimize"); */
590 /* This should be infeasibility tolerance.*/
591 lp_data->lpetol = ekk_getRtoldinf(baseModel);
592
593 /* Speed up for large sparse problems. Test it, if it's faster or not. */
594 osllib_status = ekk_setIuseRowCopy(baseModel, 1);
595 OSL_check_error("open_lp_solver - ekk_setIuseRowCopy");
596 }
597
598 /*===========================================================================*/
close_lp_solver(LPdata * lp_data)599 void close_lp_solver(LPdata *lp_data)
600 {
601 if (lp_data->lp != NULL) {
602 osllib_status = ekk_deleteModel(lp_data->lp);
603 OSL_check_error("close_lp_solver - ekk_deleteModel");
604 lp_data->lp = NULL;
605 }
606 ekk_endContext(lp_data->env);
607 }
608
609 /*===========================================================================*/
610
611 /*===========================================================================*\
612 * This function loads the data of an lp into the lp solver.
613 \*===========================================================================*/
load_lp_prob(LPdata * lp_data,int scaling,int fastmip)614 void load_lp_prob(LPdata *lp_data, int scaling, int fastmip)
615 {
616 int i;
617 double *lr = lp_data->tmp.d, *ur = lp_data->tmp.d + lp_data->n;
618
619 lp_data->lp = ekk_newModel(lp_data->env, NULL);
620 osllib_status = (lp_data->env == NULL);
621 OSL_check_error("load_lp_prob - ekk_newModel");
622
623 for (i = 0; i < lp_data->m; i++) {
624 switch (lp_data->mip->sense[i]) {
625 case 'E': lr[i] = ur[i] = lp_data->mip->rhs[i]; break;
626 case 'L': lr[i] = - OSL_INFINITY; ur[i] = lp_data->mip->rhs[i]; break;
627 case 'G': lr[i] = lp_data->mip->rhs[i]; ur[i] = OSL_INFINITY; break;
628 case 'R':
629 if (lp_data->mip->rngval[i] >= 0) {
630 ur[i] = lp_data->mip->rhs[i];
631 lr[i] = ur[i] - lp_data->mip->rngval[i];
632 } else {
633 ur[i] = lp_data->mip->rhs[i];
634 lr[i] = ur[i] + lp_data->mip->rngval[i];
635 }
636 break;
637 default: /* This should never happen ... */
638 osllib_status = -1;
639 OSL_check_error("load_lp - unknown sense");
640 }
641 }
642 osllib_status =
643 ekk_loadRimModel(lp_data->lp, lp_data->m, lr, ur, lp_data->n,
644 lp_data->mip->obj, lp_data->mip->lb,
645 lp_data->mip->ub);
646 OSL_check_error("load_lp - ekk_loadRimModel");
647 osllib_status =
648 ekk_addColumnElementBlock(lp_data->lp, lp_data->n, lp_data->mip->matind,
649 lp_data->mip->matbeg, lp_data->mip->matval);
650 OSL_check_error("load_lp - ekk_addColumnElementBlock");
651 /* Not sure we need this since there's only one block */
652 osllib_status = ekk_mergeBlocks(lp_data->lp, 1);
653 OSL_check_error("load_lp - ekk_mergeBlocks");
654
655 /* lp_data->scaling = scaling; */
656 }
657
658 /*===========================================================================*/
659
unload_lp_prob(LPdata * lp_data)660 void unload_lp_prob(LPdata *lp_data)
661 {
662 osllib_status = ekk_deleteModel(lp_data->lp);
663 OSL_check_error("unload_lp - ekk_deleteModel");
664 lp_data->lp = NULL;
665
666 lp_data->m = lp_data->n = lp_data->nz = 0;
667 }
668
669 /*===========================================================================*/
670
load_basis(LPdata * lp_data,int * cstat,int * rstat)671 void load_basis(LPdata *lp_data, int *cstat, int *rstat)
672 {
673 int *stat, i;
674
675 if (cstat != NULL) {
676 stat = ekk_getColstat(lp_data->lp);
677 for (i = lp_data->n - 1; i >= 0; i--) {
678 stat[i] &= 0x1fffffff;
679 switch (cstat[i]) {
680 case VAR_BASIC: stat[i] |= 0x80000000; break;
681 case VAR_FREE: stat[i] |= 0x60000000; break;
682 case VAR_AT_UB: stat[i] |= 0x40000000; break;
683 case VAR_AT_LB: stat[i] |= 0x20000000; break;
684 case VAR_FIXED: stat[i] |= 0x00000000; break;
685 default: break; /* should never happen */
686 }
687 }
688 osllib_status = ekk_setColstat(lp_data->lp, stat);
689 OSL_check_error("load_basis - ekk_setColstat");
690 ekk_free(stat);
691 }
692 if (rstat != NULL) {
693 stat = ekk_getRowstat(lp_data->lp);
694 for (i = lp_data->m - 1; i >= 0; i--) {
695 stat[i] &= 0x1fffffff;
696 switch (rstat[i]) {
697 case SLACK_BASIC: stat[i] |= 0x80000000; break;
698 case SLACK_FREE: stat[i] |= 0x60000000; break;
699 case SLACK_AT_UB: stat[i] |= 0x40000000; break;
700 case SLACK_AT_LB: stat[i] |= 0x20000000; break;
701 case SLACK_FIXED: stat[i] |= 0x00000000; break;
702 }
703 }
704 osllib_status = ekk_setRowstat(lp_data->lp, stat);
705 OSL_check_error("load_basis - ekk_setRowstat");
706 ekk_free(stat);
707 }
708 lp_data->lp_is_modified = LP_HAS_NOT_BEEN_MODIFIED;
709 }
710
711 /*===========================================================================*/
712
refactorize(LPdata * lp_data)713 void refactorize(LPdata *lp_data)
714 {
715 fprintf(stderr, "Function not implemented yet.");
716 exit(-1);
717 }
718
719 /*===========================================================================*/
720
add_rows(LPdata * lp_data,int rcnt,int nzcnt,double * rhs,char * sense,int * rmatbeg,int * rmatind,double * rmatval)721 void add_rows(LPdata *lp_data, int rcnt, int nzcnt, double *rhs,
722 char *sense, int *rmatbeg, int *rmatind, double *rmatval)
723 {
724 int i;
725 double *lr, *ur;
726 /* double *lr = lp_data->tmp.d, *ur = lp_data->tmp.d + lp_data->n; */
727
728 lr = (double *) malloc(rcnt * DSIZE);
729 ur = (double *) malloc(rcnt * DSIZE);
730 for (i = rcnt - 1; i >= 0; i--) {
731 switch (sense[i]) {
732 case 'E': lr[i] = ur[i] = rhs[i]; break;
733 case 'L': lr[i] = - OSL_INFINITY; ur[i] = rhs[i]; break;
734 case 'G': lr[i] = rhs[i]; ur[i] = OSL_INFINITY; break;
735 case 'R': lr[i] = ur[i] = lp_data->mip->rhs[i]; break;
736 /* Range will be added later in change_range */
737 default: /*This should never happen ... */
738 osllib_status = -1;
739 OSL_check_error("add_rows - unknown sense");
740 }
741 }
742 osllib_status = ekk_addRows(lp_data->lp, rcnt, lr, ur, rmatbeg, rmatind,
743 rmatval);
744 OSL_check_error("add_rows - ekk_addRows");
745
746 /* Merge block can make comutation faster */
747 osllib_status = ekk_mergeBlocks(lp_data->lp, 1);
748 OSL_check_error("add_rows - ekk_mergeBlocks");
749
750 FREE(lr);
751 FREE(ur);
752
753 lp_data->m += rcnt;
754 lp_data->nz += nzcnt;
755 lp_data->lp_is_modified = LP_HAS_BEEN_MODIFIED;
756 }
757
758 /*===========================================================================*/
759
add_cols(LPdata * lp_data,int ccnt,int nzcnt,double * obj,int * cmatbeg,int * cmatind,double * cmatval,double * lb,double * ub,char * where_to_move)760 void add_cols(LPdata *lp_data, int ccnt, int nzcnt, double *obj,
761 int *cmatbeg, int *cmatind, double *cmatval,
762 double *lb, double *ub, char *where_to_move)
763 {
764 osllib_status = ekk_addColumns(lp_data->lp, ccnt, obj, lb, ub,
765 cmatbeg, cmatind, cmatval);
766 OSL_check_error("add_cols - ekk_addColumns");
767 osllib_status = ekk_mergeBlocks(lp_data->lp, 1);
768 OSL_check_error("add_cols - ekk_mergeBlocks");
769 lp_data->n += ccnt;
770 lp_data->nz += nzcnt;
771 }
772
773 /*===========================================================================*/
774
change_row(LPdata * lp_data,int row_ind,char sense,double rhs,double range)775 void change_row(LPdata *lp_data, int row_ind,
776 char sense, double rhs, double range)
777 {
778 /*can be sped up using ekk_rowlower - direct acces to internal data*/
779 double lr, ur;
780 switch (sense) {
781 case 'E': lr = ur = rhs; break;
782 case 'L': lr = - OSL_INFINITY; ur = rhs; break;
783 case 'G': lr = rhs; ur = OSL_INFINITY; break;
784 case 'R':
785 if (range >= 0) {
786 lr = rhs; ur = lr + range;
787 } else {
788 ur = rhs; lr = ur + range;
789 }
790 break;
791 default: /*This should never happen ... */
792 osllib_status = -1;
793 OSL_check_error("change_row - default");
794 }
795 osllib_status = ekk_copyRowlower(lp_data->lp, &lr, row_ind, row_ind + 1);
796 OSL_check_error("change_row - ekk_copyRowlower");
797 osllib_status = ekk_copyRowupper(lp_data->lp, &ur, row_ind, row_ind + 1);
798 OSL_check_error("change_row - ekk_copyRowupper");
799 lp_data->lp_is_modified = LP_HAS_BEEN_MODIFIED;
800 }
801
802 /*===========================================================================*/
803
change_col(LPdata * lp_data,int col_ind,char sense,double lb,double ub)804 void change_col(LPdata *lp_data, int col_ind,
805 char sense, double lb, double ub)
806 {
807 switch (sense){
808 case 'E': change_lbub(lp_data, col_ind, lb, ub); break;
809 case 'R': change_lbub(lp_data, col_ind, lb, ub); break;
810 case 'G': change_lb(lp_data, col_ind, lb); break;
811 case 'L': change_ub(lp_data, col_ind, ub); break;
812 default: /*This should never happen ... */
813 osllib_status = -1;
814 OSL_check_error("change_col - default");
815 }
816 }
817
818 /*===========================================================================*/
819
820 /*===========================================================================*\
821 * Solve the lp specified in lp_data->lp with dual simplex. The number of
822 * iterations is returned in 'iterd'. The return value of the function is
823 * the termination code of the dual simplex method.
824 \*===========================================================================*/
825 /* Basis head in the end of this function not finished yet */
dual_simplex(LPdata * lp_data,int * iterd)826 int dual_simplex(LPdata *lp_data, int *iterd)
827 {
828 int term;
829
830 ekk_mergeBlocks(lp_data->lp, 1);
831 ekk_setIiternum(lp_data->lp, 0);
832
833 /*PreSolve seems to cause some problems -- not sure exactly why, but we
834 leave it turned off for now. */
835 #if 0
836 if ((osllib_status = ekk_preSolve(lp_data->lp, 3, NULL)) == 1){
837 /* This means infeasibility was detected during preSolve */
838 term = lp_data->termcode = D_UNBOUNDED;
839 *iterd = 0;
840 lp_data->lp_is_modified = LP_HAS_NOT_BEEN_MODIFIED;
841 return(term);
842 }
843 #endif
844
845 if (lp_data->lp_is_modified == LP_HAS_BEEN_ABANDONED) {
846 /* osllib_status = ekk_crash(lp_data->lp, 2); */
847 /* OSL_check_error("dual_simplex - ekk_crash"); */
848 osllib_status = ekk_allSlackBasis(lp_data->lp);
849 OSL_check_error("dual_simplex - ekk_allSlackBasis");
850 }
851
852 #if 1
853 ekk_dualSimplex(lp_data->lp);
854 #else
855 ekk_simplex(lp_data->lp, 256 + 32); // no presolve and no scaling
856 #endif
857
858 term = ekk_getIprobstat(lp_data->lp);
859 /*Turn postSolve back on if we figure out preSolve problem */
860 #if 0
861 osllib_status = ekk_postSolve(lp_data->lp, NULL);*/
862 #endif
863
864 /* We don't need this if we are not using preSolve. Not sure we need it
865 anyway... */
866 #if 0
867 /* Once more without preSolve. Dual simplex is run again
868 only if solution is not optimal */
869 if ((term == 1) || (term == 2)) {
870 term = ekk_dualSimplex(lp_data->lp);
871 term = ekk_primalSimplex(lp_data->lp, 3);
872 }
873 #endif
874
875 #if 0
876 If (term == 2) {
877 /* Dual infeas. This is impossible, so we must have had iteration
878 * limit AND bound shifting AND dual feasibility not restored within
879 * the given iteration limit. */
880 maxiter = ekk_getImaxiter(lp_data->lp);
881 osllib_status = ekk_setImaxiter(lp_data->lp, LP_MAX_ITER);
882 OSL_check_error("dual_simplex - ekk_setImaxiter");
883 term = ekk_dualSimplex(lp_data->lp);
884 osllib_status = ekk_setImaxiter(lp_data->lp, maxiter);
885 OSL_check_error("dual_simplex - ekk_setImaxiter");
886 }
887 #endif
888
889 switch (term) {
890 case 0:
891 term = LP_OPTIMAL;
892 break;
893 case 1:
894 term = LP_D_UNBOUNDED;
895 ekk_infeasibilities(lp_data->lp, 1, 1, NULL, NULL);
896 break;
897 case 2:
898 term = LP_D_INFEASIBLE;
899 break;
900 case 3:
901 term = LP_D_ITLIM;
902 break;
903 case 4:
904 osllib_status = -1;
905 OSL_check_error("osllib_status-ekk_dualSimplex found no solution!");
906 term = LP_ABANDONED;
907 break;
908 case 5:
909 LP_D_OBJLIM;
910 break;
911 case 6:
912 osllib_status = -1;
913 OSL_check_error("osllib_status-ekk_dualSimplex lack of dstorage file"
914 "space!");
915 term = LP_ABANDONED;
916 break;
917 default: term = LP_ABANDONED;
918 break;
919 }
920
921 lp_data->termcode = term;
922
923 if (term != LP_ABANDONED){
924 *iterd = ekk_getIiternum(lp_data->lp);
925 lp_data->objval = ekk_getRobjvalue(lp_data->lp);
926 lp_data->lp_is_modified = LP_HAS_NOT_BEEN_MODIFIED;
927 }else{
928 lp_data->lp_is_modified = LP_HAS_BEEN_ABANDONED;
929 }
930 return(term);
931 }
932
933 /*===========================================================================*/
solve_hotstart(LPdata * lp_data,int * iterd)934 int solve_hotstart(LPdata *lp_data, int *iterd)
935 {
936 return(dual_simplex(lp_data,iterd));
937 }
938 /*===========================================================================*/
unmark_hotstart(LPdata * lp_data)939 int unmark_hotstart(LPdata *lp_data)
940 {
941 /* only when using osi */
942 return (0);
943 }
944
945 /*===========================================================================*/
mark_hotstart(LPdata * lp_data)946 int mark_hotstart(LPdata *lp_data)
947 {
948 /* only when using osi */
949 return (0);
950 }
951
952 /*===========================================================================*/
953
btran(LPdata * lp_data,double * col)954 void btran(LPdata *lp_data, double *col)
955 {
956 osllib_status = ekk_formBInverseTransposeb(lp_data->lp, col);
957 OSL_check_error("btran - ekk_formBInverseTransposeb");
958 }
959
960 /*===========================================================================*/
961 /* This function is not used currently ... */
962
get_binvcol(LPdata * lp_data,int j,double * col)963 void get_binvcol(LPdata *lp_data, int j, double *col)
964 {
965 fprintf(stderr, "Function not implemented yet.");
966 exit(-1);
967 }
968
969 /*===========================================================================*/
970 /* This function is used only together with get_proof_of_infeasibility... */
971
get_binvrow(LPdata * lp_data,int i,double * row)972 void get_binvrow(LPdata *lp_data, int i, double *row)
973 {
974 fprintf(stderr, "Function not implemented yet.");
975 exit(-1);
976 }
977
978 /*===========================================================================*/
979 /* This function is never called either... */
980
get_basis_header(LPdata * lp_data)981 void get_basis_header(LPdata *lp_data)
982 {
983 fprintf(stderr, "Function not implemented yet.");
984 exit(-1);
985 }
986
987 /*===========================================================================*/
988
get_basis(LPdata * lp_data,int * cstat,int * rstat)989 void get_basis(LPdata *lp_data, int *cstat, int *rstat)
990 {
991 int i, temp_stat;
992 const int *stat;
993
994 if (cstat != NULL) {
995 stat = ekk_colstat(lp_data->lp);
996 for (i = lp_data->n - 1; i >= 0; i--) {
997 if ((stat[i] & 0x80000000) != 0) {
998 cstat[i] = VAR_BASIC;
999 } else {
1000 temp_stat = stat[i] & 0x60000000;
1001 switch (temp_stat) {
1002 case 0x60000000: cstat[i] = VAR_FREE; break;
1003 case 0x40000000: cstat[i] = VAR_AT_UB; break;
1004 case 0x20000000: cstat[i] = VAR_AT_LB; break;
1005 case 0x00000000: cstat[i] = VAR_FIXED; break;
1006 }
1007 }
1008 }
1009 }
1010 if (rstat != NULL) {
1011 stat = ekk_rowstat(lp_data->lp);
1012 for (i = lp_data->m - 1; i >= 0; i--) {
1013 if ((stat[i] & 0x80000000) != 0) {
1014 rstat[i] = SLACK_BASIC;
1015 } else {
1016 temp_stat = stat[i] & 0x60000000;
1017 switch (temp_stat) {
1018 case 0x60000000: rstat[i] = SLACK_FREE; break;
1019 case 0x40000000: rstat[i] = SLACK_AT_UB; break;
1020 case 0x20000000: rstat[i] = SLACK_AT_LB; break;
1021 case 0x00000000: rstat[i] = SLACK_FIXED; break;
1022 }
1023 }
1024 }
1025 }
1026 }
1027
1028 /*===========================================================================*/
1029
1030 /*===========================================================================*\
1031 * Set an upper limit on the objective function value.
1032 \*===========================================================================*/
set_obj_upper_lim(LPdata * lp_data,double lim)1033 void set_obj_upper_lim(LPdata *lp_data, double lim)
1034 {
1035 /* Not sure how to do this in OSL */
1036 fprintf(stderr, "Function not implemented yet.");
1037 exit(-1);
1038 }
1039
1040 /*===========================================================================*/
1041
1042 /*===========================================================================*\
1043 * Set an upper limit on the number of iterations. If itlim < 0 then set
1044 * it to the maximum.
1045 \*===========================================================================*/
1046
set_itlim(LPdata * lp_data,int itlim)1047 void set_itlim(LPdata *lp_data, int itlim)
1048 {
1049 if (itlim < 0) itlim = LP_MAX_ITER;
1050 osllib_status = ekk_setImaxiter(lp_data->lp, itlim);
1051 OSL_check_error("set_itlim - ekk_setImaxiter");
1052 }
1053
1054 /*===========================================================================*/
set_itlim_hotstart(LPdata * lp_data,int itlim)1055 void set_itlim_hotstart(LPdata *lp_data, int itlim)
1056 {
1057 /* read being and nothingness -- Jean Paul Sartre */
1058 }
1059
1060 /*===========================================================================*/
1061
get_column(LPdata * lp_data,int j,double * colval,int * colind,int * collen,double * cj)1062 void get_column(LPdata *lp_data, int j,
1063 double *colval, int *colind, int *collen, double *cj)
1064 {
1065 EKKVector vec;
1066 vec = ekk_getColumn(lp_data->lp, j);
1067 *collen = vec.numNonZero;
1068 memcpy(colind, vec.index, *collen * ISIZE);
1069 memcpy(colval, vec.element, *collen * DSIZE);
1070 ekk_freeVector(&vec);
1071 get_objcoef(lp_data, j, cj);
1072 }
1073
1074 /*===========================================================================*/
get_row(LPdata * lp_data,int i,double * rowval,int * rowind,int * rowlen,double * rowub,double * rowlb)1075 void get_row(LPdata *lp_data, int i,
1076 double *rowval, int *rowind, int *rowlen,
1077 double *rowub, double *rowlb)
1078 {
1079 EKKVector vec;
1080 vec = ekk_getRow(lp_data->lp, i);
1081 *rowlen = vec.numNonZero;
1082 memcpy(rowind, vec.index, *rowlen * ISIZE);
1083 memcpy(rowval, vec.element, *rowlen * DSIZE);
1084 ekk_freeVector(&vec);
1085 }
1086
1087 /*===========================================================================*/
1088 /* This routine returns the index of a row which proves the lp to be primal
1089 * infeasible. It is only needed when column generation is used. */
1090 /*===========================================================================*/
1091
get_proof_of_infeas(LPdata * lp_data,int * infind)1092 int get_proof_of_infeas(LPdata *lp_data, int *infind)
1093 {
1094 fprintf(stderr, "Function not implemented yet.");
1095 return(0);
1096 }
1097
1098 /*===========================================================================*\
1099 * Get the solution (values of the structural variables in an optimal
1100 * solution) to the lp (specified by lp_data->lp) into the vector
1101 * lp_data->x.
1102 \*===========================================================================*/
get_x(LPdata * lp_data)1103 void get_x(LPdata *lp_data)
1104 {
1105 memcpy(lp_data->x, ekk_colsol(lp_data->lp), lp_data->n * DSIZE);
1106 }
1107
1108 /*===========================================================================*/
get_dj_pi(LPdata * lp_data)1109 void get_dj_pi(LPdata *lp_data)
1110 {
1111 /*If scaling, fast integer or compress is used, maybe some changes will be
1112 needed */
1113 /* OSL returns changed sign - is it good or not? */
1114 memcpy(lp_data->dualsol, ekk_rowduals(lp_data->lp), lp_data->m * DSIZE);
1115
1116 # if 0
1117 /* changing the sign */
1118 for (i = lp_data->m - 1; i >= 0; i --) {
1119 lp_data->dualsol[i] = - lp_data->dualsol[i];
1120 }
1121 #endif
1122
1123 memcpy(lp_data->dj, ekk_colrcosts(lp_data->lp), lp_data->n * DSIZE);
1124
1125 #if 0
1126 for (i = lp_data->n - 1; i >= 0; i --) {
1127 lp_data->dj[i] = - lp_data->dj[i];
1128 }
1129 #endif
1130 }
1131
1132 /*===========================================================================*/
1133 /* Possible improper implementetion. */
1134
get_slacks(LPdata * lp_data)1135 void get_slacks(LPdata *lp_data)
1136 {
1137 row_data *rows = lp_data->rows;
1138 double *slacks = lp_data->slacks;
1139 const double *racts;
1140 int i, m = lp_data->m;
1141
1142 racts = ekk_rowacts(lp_data->lp);
1143
1144 for (i = m - 1; i >= 0; i--) {
1145 if ((rows[i].cut->sense == 'R') && (rows[i].cut->range < 0) ) {
1146 slacks[i] = - rows[i].cut->rhs + racts[i];
1147 } else {
1148 slacks[i] = rows[i].cut->rhs - racts[i];
1149 }
1150 }
1151 }
1152
1153 /*===========================================================================*/
1154
change_range(LPdata * lp_data,int rowind,double value)1155 void change_range(LPdata *lp_data, int rowind, double value)
1156 {
1157 const double *lrow, *urow;
1158 double lr, ur;
1159 lrow = ekk_rowlower(lp_data->lp);
1160 urow = ekk_rowupper(lp_data->lp);
1161 if (value >= 0) {
1162 lr = urow[rowind] - value;
1163 } else {
1164 lr = lrow[rowind] + value;
1165 }
1166 osllib_status = ekk_copyRowlower(lp_data->lp, &ur, rowind, rowind + 1);
1167 OSL_check_error("change_range - ekk_copyRowupper");
1168 lp_data->lp_is_modified = LP_HAS_BEEN_MODIFIED;
1169 }
1170
1171 /*===========================================================================*/
1172 /* This function is never called ... */
1173
change_rhs(LPdata * lp_data,int rownum,int * rhsind,double * rhsval)1174 void change_rhs(LPdata *lp_data, int rownum, int *rhsind, double *rhsval)
1175 {
1176 fprintf(stderr, "Function not implemented yet.");
1177 exit(-1);
1178 }
1179
1180 /*===========================================================================*/
1181 /* This function is never called ...*/
1182
change_sense(LPdata * lp_data,int cnt,int * index,char * sense)1183 void change_sense(LPdata *lp_data, int cnt, int *index, char *sense)
1184 {
1185 fprintf(stderr, "Function not implemented yet.");
1186 exit(-1);
1187 }
1188
1189 /*===========================================================================*/
1190
change_bounds(LPdata * lp_data,int cnt,int * index,char * lu,double * bd)1191 void change_bounds(LPdata *lp_data, int cnt, int *index, char *lu, double *bd)
1192 {
1193 double *lb, *ub;
1194 int i, j;
1195 lb = ekk_getCollower(lp_data->lp);
1196 ub = ekk_getColupper(lp_data->lp);
1197 for (i = cnt - 1; i >= 0; i--) {
1198 j = index[i];
1199 switch (lu[i]) {
1200 case 'L': lb[j] = bd[i];break;
1201 case 'U': ub[j] = bd[i];break;
1202 case 'B': lb[j] = ub[j] = bd[i];break;
1203 default: /*This should never happen ... */
1204 osllib_status = -1;
1205 OSL_check_error("change_bounds - default");
1206 }
1207 }
1208 osllib_status = ekk_setCollower(lp_data->lp, lb);
1209 OSL_check_error("change_bounds - ekk_setCollower");
1210 ekk_free(lb);
1211 osllib_status = ekk_setColupper(lp_data->lp, ub);
1212 OSL_check_error("change_bounds - ekk_setColupper");
1213 ekk_free(ub);
1214 lp_data->lp_is_modified = LP_HAS_BEEN_MODIFIED;
1215 }
1216
1217 /*===========================================================================*/
1218
change_lbub(LPdata * lp_data,int j,double lb,double ub)1219 void change_lbub(LPdata *lp_data, int j, double lb, double ub)
1220 {
1221 osllib_status = ekk_copyColupper(lp_data->lp, &ub, j, j + 1);
1222 OSL_check_error("change_lbub - ekk_copyColupper");
1223 osllib_status = ekk_copyCollower(lp_data->lp, &lb, j, j + 1);
1224 OSL_check_error("change_lbub - ekk_copyCollower");
1225 lp_data->lp_is_modified = LP_HAS_BEEN_MODIFIED;
1226 }
1227
1228 /*===========================================================================*/
1229
change_ub(LPdata * lp_data,int j,double ub)1230 void change_ub(LPdata *lp_data, int j, double ub)
1231 {
1232 osllib_status = ekk_copyColupper(lp_data->lp, &ub, j, j + 1);
1233 OSL_check_error("change_ub - ekk_copyColupper");
1234 lp_data->lp_is_modified = LP_HAS_BEEN_MODIFIED;
1235 }
1236
1237 /*===========================================================================*/
1238
change_lb(LPdata * lp_data,int j,double lb)1239 void change_lb(LPdata *lp_data, int j, double lb)
1240 {
1241 osllib_status = ekk_copyCollower(lp_data->lp, &lb, j, j + 1);
1242 OSL_check_error("change_lb - ekk_copyCollower");
1243 lp_data->lp_is_modified = LP_HAS_BEEN_MODIFIED;
1244 }
1245
1246 /*===========================================================================*/
1247
get_ub(LPdata * lp_data,int j,double * ub)1248 void get_ub(LPdata *lp_data, int j, double *ub)
1249 {
1250 /* Maybe some range checking could be added ...*/
1251 const double *uc = ekk_colupper(lp_data->lp);
1252 *ub = uc[j];
1253 }
1254
1255 /*===========================================================================*/
1256
get_lb(LPdata * lp_data,int j,double * lb)1257 void get_lb(LPdata *lp_data, int j, double *lb)
1258 {
1259 /* Maybe some range checking could be added ...*/
1260 const double *lc = ekk_collower(lp_data->lp);
1261 *lb = lc[j];
1262 }
1263
1264 /*===========================================================================*/
1265
get_bounds(LPdata * lp_data)1266 void get_bounds(LPdata *lp_data)
1267 {
1268 lp_data->ub = const_cast<double *>(ekk_colupper(lp_data->lp));
1269 lp_data->lb = const_cast<double *>(ekk_collower(lp_data->lp));
1270 }
1271
1272 /*===========================================================================*/
1273
get_objcoef(LPdata * lp_data,int j,double * objcoef)1274 void get_objcoef(LPdata *lp_data, int j, double *objcoef)
1275 {
1276 /* Maybe some range checking could be added ...*/
1277 const double *oc = ekk_objective(lp_data->lp);
1278 *objcoef = oc[j];
1279 }
1280
1281 /*===========================================================================*/
1282
delete_rows(LPdata * lp_data,int deletable,int * free_rows)1283 void delete_rows(LPdata *lp_data, int deletable, int *free_rows)
1284 {
1285 int i, m = lp_data->m;
1286 int *which = lp_data->tmp.i1 + lp_data->m;
1287 int delnum = 0;
1288
1289 /* which = calloc(delnum, ISIZE); */
1290 for (i = m - 1, delnum = 0; i >= 0; i--){
1291 if (free_rows[i]){
1292 which[delnum++] = i;
1293 }
1294 }
1295 osllib_status = ekk_deleteRows(lp_data->lp, delnum, which);
1296 OSL_check_error("delete_rows - ekk_deleteRows");
1297
1298 #if 0
1299 /* Make result as CPLEX does*/
1300 for (i = 0, delnum = 0; i < m; i++){
1301 if (free_rows[i])
1302 free_rows[i] = -1;
1303 else
1304 free_rows[i] = delnum++;
1305 }
1306 #endif
1307
1308 lp_data->m -= delnum;
1309 lp_data->lp_is_modified = LP_HAS_BEEN_MODIFIED;
1310 }
1311
1312 /*===========================================================================*/
1313
delete_cols(LPdata * lp_data,int delnum,int * delstat)1314 int delete_cols(LPdata *lp_data, int delnum, int *delstat)
1315 {
1316 int i, n = lp_data->n;
1317 int *which = (int *) calloc(delnum, ISIZE);
1318 int num_to_delete = 0, num_to_keep = 0;
1319 double *dj = lp_data->dj;
1320 double *x = lp_data->x;
1321 char *status = lp_data->status;
1322
1323 for (i = n - 1, num_to_delete = 0; i >= 0; i--) {
1324 if (delstat[i]) {
1325 which[num_to_delete++] = i;
1326 }
1327 }
1328
1329 if (!num_to_delete) return(0);
1330
1331 osllib_status = ekk_deleteColumns(lp_data->lp, num_to_delete, which);
1332 OSL_check_error("delete_cols - ekk_deleteCols");
1333 FREE(which);
1334
1335 lp_data->nz = ekk_getInumels(lp_data->lp);
1336 OSL_check_error("delete_cols - ekk_getInumels");
1337
1338 for (i = 0, num_to_keep = 0; i < lp_data->n; i++){
1339 if (delstat[i]){
1340 delstat[i] = -1;
1341 }else{
1342 delstat[i] = num_to_keep++;
1343 dj[delstat[i]] = dj[i];
1344 x[delstat[i]] = x[i];
1345 status[delstat[i]] = status[i];
1346 }
1347 }
1348
1349 lp_data->n = num_to_keep;
1350 return (num_to_delete);
1351 }
1352
1353 /*===========================================================================*/
1354 /* Original (CPLEX) implementation is nothing :-) */
1355 /*===========================================================================*/
1356
release_var(LPdata * lp_data,int j,int where_to_move)1357 void release_var(LPdata *lp_data, int j, int where_to_move)
1358 {
1359 #if 0
1360 switch (where_to_move){
1361 case MOVE_TO_UB:
1362 lp_data->lpbas.cstat[j] = 2; break; /* non-basic at its upper bound */
1363 case MOVE_TO_LB:
1364 lp_data->lpbas.cstat[j] = 0; break; /* non-basic at its lower bound */
1365 }
1366 lp_data->lp_is_modified = LP_HAS_BEEN_MODIFIED;
1367 #endif
1368 }
1369
1370 /*===========================================================================*/
1371 /* There were some side effects setting "temp" fields of lp_data. */
1372
free_row_set(LPdata * lp_data,int length,int * index)1373 void free_row_set(LPdata *lp_data, int length, int *index)
1374 {
1375 int i, j;
1376 double *lb = (double *) ekk_getRowlower(lp_data->lp);
1377 double *ub = (double *) ekk_getRowupper(lp_data->lp);
1378
1379 for (i = length - 1; i >= 0; i--) {
1380 j = index[i];
1381 lb[j] = - OSL_INFINITY;
1382 ub[j] = OSL_INFINITY;
1383 }
1384 osllib_status = ekk_setRowlower(lp_data->lp, lb);
1385 OSL_check_error("free_row_set ekk_setRowLower");
1386 ekk_free(lb);
1387 osllib_status = ekk_setRowupper(lp_data->lp, ub);
1388 OSL_check_error("free_row_set ekk_setRowUpper");
1389 ekk_free(ub);
1390 lp_data->lp_is_modified = LP_HAS_BEEN_MODIFIED;
1391 }
1392
1393 /*===========================================================================*/
1394 /* There were some side effects setting "temp" fileds of lp_data. */
1395
constrain_row_set(LPdata * lp_data,int length,int * index)1396 void constrain_row_set(LPdata *lp_data, int length, int *index)
1397 {
1398 int i, j = 0;
1399 double *lb = ekk_getRowlower(lp_data->lp);
1400 double *ub = ekk_getRowupper(lp_data->lp);
1401 row_data *rows = lp_data->rows;
1402 cut_data *cut;
1403
1404 for (i = length - 1; i >= 0; i--) {
1405 j = index[i];
1406 cut = rows[j].cut;
1407 switch (cut->sense){
1408 case 'E': lb[j] = ub[j] = cut->rhs; break;
1409 case 'L': lb[j] = - OSL_INFINITY; ub[j] = cut->rhs; break;
1410 case 'G': lb[j] = cut->rhs; ub[j] = OSL_INFINITY; break;
1411 case 'R':
1412 if (lp_data->mip->rngval[j] >= 0) {
1413 ub[j] = cut->rhs; lb[j] = ub[j] - lp_data->mip->rngval[j];
1414 } else {
1415 ub[j] = cut->rhs; lb[j] = ub[j] + lp_data->mip->rngval[j];
1416 }
1417 break;
1418 default: /*This should never happen ... */
1419 osllib_status = -1;
1420 OSL_check_error("load_lp - unknown type of constraint");
1421 }
1422 }
1423
1424 j = 0;
1425 if (j){
1426 ekk_free(lb);
1427 }
1428
1429 osllib_status = ekk_setRowlower(lp_data->lp, lb);
1430 OSL_check_error("constrain_row_set ekk_setRowLower");
1431 ekk_free(lb);
1432 osllib_status = ekk_setRowupper(lp_data->lp, ub);
1433 OSL_check_error("constrain_row_set ekk_setRowUpper");
1434 ekk_free(ub);
1435 lp_data->lp_is_modified = LP_HAS_BEEN_MODIFIED;
1436 }
1437
1438 /*===========================================================================*/
1439
read_mps(MIPdesc * desc,char * infile,char * probname,int verbosity)1440 int read_mps(MIPdesc *desc, char *infile, char *probname, int verbosity)
1441 {
1442 printf("\nMps-format file can be read only through OSI interface.\n");
1443
1444 return(1);
1445 }
1446
1447 /*===========================================================================*/
1448
read_lp(MIPdesc * desc,char * infile,char * probname,int verbosity)1449 int read_lp(MIPdesc *desc, char *infile, char *probname, int verbosity)
1450 {
1451 printf("\nLP-format file can be read only through OSI interface.\n");
1452
1453 return(1);
1454 }
1455
1456 /*===========================================================================*/
1457
write_mip_desc_mps(MIPdesc * mip,char * fname)1458 void write_mip_desc_mps(MIPdesc *mip, char *fname)
1459 {
1460 fprintf(stderr, "Function not implemented yet.");
1461 exit(-1);
1462 }
1463
1464
1465 /*===========================================================================*/
1466
write_mip_desc_lp(MIPdesc * mip,char * fname)1467 void write_mip_desc_lp(MIPdesc *mip, char *fname)
1468 {
1469 fprintf(stderr, "Function not implemented yet.");
1470 exit(-1);
1471 }
1472
1473 /*===========================================================================*/
1474
write_mps(LPdata * lp_data,char * fname)1475 void write_mps(LPdata *lp_data, char *fname)
1476 {
1477 osllib_status = ekk_exportModel(lp_data->lp, fname, 1, 2);
1478 OSL_check_error("write_mps");
1479 }
1480
1481 /*===========================================================================*/
1482
write_sav(LPdata * lp_data,char * fname)1483 void write_sav(LPdata *lp_data, char *fname)
1484 {
1485 osllib_status = ekk_saveModel(lp_data->lp, fname);
1486 OSL_check_error("write_sav");
1487 }
1488
1489 /*===========================================================================*/
1490
1491 #ifdef USE_CGL_CUTS
generate_cgl_cuts(LPdata * lp_data,int * num_cuts,cut_data *** cuts,char send_to_pool)1492 void generate_cgl_cuts(LPdata *lp_data, int *num_cuts, cut_data ***cuts,
1493 char send_to_pool)
1494 {
1495 return;
1496 }
1497 #endif
1498
1499 #endif /* __OSL__ */
1500
1501
1502 #ifdef __CPLEX__
1503
1504 /*****************************************************************************/
1505 /*****************************************************************************/
1506 /******* *******/
1507 /******* routines when CPLEX is used *******/
1508 /******* *******/
1509 /*****************************************************************************/
1510 /*****************************************************************************/
1511
1512 static int cpx_status;
1513
1514 #include <memory.h>
1515
1516 /*===========================================================================*/
1517
CPX_check_error(const char * erring_func)1518 void CPX_check_error(const char *erring_func)
1519 {
1520 if (cpx_status){
1521 printf("!!! Cplex status is nonzero !!! [%s, %i]\n",
1522 (char *)erring_func, cpx_status);
1523 }
1524 }
1525
1526 /*===========================================================================*/
1527
open_lp_solver(LPdata * lp_data)1528 void open_lp_solver(LPdata *lp_data)
1529 {
1530 int i;
1531
1532 i = CPX_OFF;
1533 lp_data->cpxenv = CPXopenCPLEX(&cpx_status);
1534 CPX_check_error("open_lp_solver - error opening environment");
1535 cpx_status = CPXsetintparam(lp_data->cpxenv, CPX_PARAM_SCRIND, i);
1536 CPX_check_error("open_lp_solver - CPXsetintparam, SCRIND");
1537 #if 0
1538 lp_data->lpetol = 1e-09;
1539 cpx_status = CPXsetdblparam(lp_data->cpxenv, CPX_PARAM_EPRHS,
1540 lp_data->lpetol);
1541 CPX_check_error("open_lp_solver - CPXsetdblparam");
1542 #else
1543 cpx_status = CPXgetdblparam(lp_data->cpxenv, CPX_PARAM_EPRHS,
1544 &lp_data->lpetol);
1545 CPX_check_error("open_lp_solver - CPXgetdblparam");
1546 #endif
1547 }
1548
1549 /*===========================================================================*/
1550
close_lp_solver(LPdata * lp_data)1551 void close_lp_solver(LPdata *lp_data)
1552 {
1553 if (lp_data->lp){
1554 cpx_status = CPXfreeprob(lp_data->cpxenv, &(lp_data->lp));
1555 CPX_check_error("close_lp_solver");
1556 lp_data->lp = NULL;
1557 }
1558 cpx_status = CPXcloseCPLEX(&lp_data->cpxenv);
1559 CPX_check_error("close_lp_solver");
1560 }
1561
1562 /*===========================================================================*/
1563
1564 /*===========================================================================*\
1565 * This function loads the data of an lp into the lp solver. This involves
1566 * transforming the data into CPLEX format and calling the CPLEX function
1567 * 'loadlp'.
1568 \*===========================================================================*/
1569
load_lp_prob(LPdata * lp_data,int scaling,int fastmip)1570 void load_lp_prob(LPdata *lp_data, int scaling, int fastmip)
1571 {
1572 int i, *matcnt, *matbeg;
1573
1574 /* realloc_lp_arrays(lp_data); */
1575
1576 matcnt = (int *) malloc (lp_data->n*ISIZE);
1577 matbeg = lp_data->mip->matbeg;
1578 for (i = lp_data->n - 1; i >= 0; i--)
1579 matcnt[i] = matbeg[i+1] - matbeg[i];
1580
1581 cpx_status = CPXsetintparam(lp_data->cpxenv, CPX_PARAM_SCAIND, -1);
1582 CPX_check_error("load_lp - CPXsetintparam - SCAIND");
1583
1584 cpx_status = CPXsetintparam(lp_data->cpxenv, CPX_PARAM_FASTMIP, fastmip);
1585 CPX_check_error("load_lp - CPXsetintparam - FASTMIP");
1586
1587 /* essentially disable basis snapshots */
1588 cpx_status =
1589 CPXsetintparam(lp_data->cpxenv, CPX_PARAM_BASINTERVAL, 2100000000);
1590 CPX_check_error("load_lp - CPXsetintparam - BASINTERVAL");
1591
1592 /* This is for the old memory model (user manages memory) */
1593 #if CPX_VERSION <= 600
1594 printf("\nSorry, CPLEX versions 6.0 and earlier are no longer supported");
1595 printf("due to incompatibilities in memory management.");
1596 printf("Please use SYMPHONY 3.0.1 or earlier.\n\n");
1597 FREE(matcnt);
1598 exit(-1);
1599 /* legacy code left for posterity */
1600 #else /* This is for the new memory model (CPLEX manages memory) */
1601 lp_data->lp = CPXcreateprob(lp_data->cpxenv,&cpx_status,(char *) "BB_prob");
1602 CPX_check_error("load_lp - CPXcreateprob");
1603 cpx_status = CPXcopylp(lp_data->cpxenv, lp_data->lp,
1604 lp_data->n, lp_data->m, 1, lp_data->mip->obj,
1605 lp_data->mip->rhs, lp_data->mip->sense,lp_data->mip->matbeg,
1606 matcnt, lp_data->mip->matind, lp_data->mip->matval,
1607 lp_data->mip->lb, lp_data->mip->ub, lp_data->mip->rngval);
1608 CPX_check_error("load_lp - CPXcopylp");
1609 FREE(matcnt);
1610 #endif
1611 }
1612
1613 /*===========================================================================*/
1614
unload_lp_prob(LPdata * lp_data)1615 void unload_lp_prob(LPdata *lp_data)
1616 {
1617 cpx_status = CPXfreeprob(lp_data->cpxenv, &lp_data->lp);
1618 CPX_check_error("unload_lp - CPXfreeprob");
1619 lp_data->lp = NULL;
1620
1621 lp_data->m = lp_data->n = lp_data->nz = 0;
1622 }
1623
1624 /*===========================================================================*/
1625
load_basis(LPdata * lp_data,int * cstat,int * rstat)1626 void load_basis(LPdata *lp_data, int *cstat, int *rstat)
1627 {
1628 cpx_status = CPXcopybase(lp_data->cpxenv, lp_data->lp, cstat, rstat);
1629 CPX_check_error("load_basis - CPXloadbase");
1630
1631 lp_data->lp_is_modified = LP_HAS_NOT_BEEN_MODIFIED;
1632 }
1633
1634 /*===========================================================================*/
1635
1636 /* There should be something nicer... */
refactorize(LPdata * lp_data)1637 void refactorize(LPdata *lp_data)
1638 {
1639 int itlim;
1640
1641 cpx_status = CPXgetintparam(lp_data->cpxenv, CPX_PARAM_ITLIM, &itlim);
1642 CPX_check_error("refactorize - CPXgetintparam");
1643 cpx_status = CPXsetintparam(lp_data->cpxenv, CPX_PARAM_ITLIM, 0);
1644 CPX_check_error("refactorize - CPXsetintparam");
1645 cpx_status = CPXprimopt(lp_data->cpxenv, lp_data->lp);
1646 CPX_check_error("refactorize - CPXoptimize");
1647 cpx_status = CPXsetintparam(lp_data->cpxenv, CPX_PARAM_ITLIM, itlim);
1648 CPX_check_error("refactorize - CPXsetintparam");
1649 }
1650
1651 /*===========================================================================*/
1652
add_rows(LPdata * lp_data,int rcnt,int nzcnt,double * rhs,char * sense,int * rmatbeg,int * rmatind,double * rmatval)1653 void add_rows(LPdata *lp_data, int rcnt, int nzcnt, double *rhs,
1654 char *sense, int *rmatbeg, int *rmatind, double *rmatval)
1655 {
1656 int i, j, indicator = FALSE;
1657
1658 if (indicator)
1659 for (i = 0; i < rcnt; i++){
1660 printf("\n");
1661 printf("%c %1f\n", sense[i], rhs[i]);
1662 for (j = rmatbeg[i]; j < rmatbeg[i+1]; j++){
1663 printf("%i ", rmatind[j]);
1664 }
1665 printf("\n");
1666 for (j = rmatbeg[i]; j < rmatbeg[i+1]; j++){
1667 printf("%1f ", rmatval[j]);
1668 }
1669 }
1670
1671 cpx_status = CPXaddrows(lp_data->cpxenv, lp_data->lp, 0, rcnt, nzcnt,
1672 rhs, sense, rmatbeg, rmatind, rmatval, NULL, NULL);
1673 CPX_check_error("add_rows");
1674 lp_data->m += rcnt;
1675 lp_data->nz += nzcnt;
1676 lp_data->lp_is_modified = LP_HAS_BEEN_MODIFIED;
1677 }
1678
1679 /*===========================================================================*/
1680
add_cols(LPdata * lp_data,int ccnt,int nzcnt,double * obj,int * cmatbeg,int * cmatind,double * cmatval,double * lb,double * ub,char * where_to_move)1681 void add_cols(LPdata *lp_data, int ccnt, int nzcnt, double *obj,
1682 int *cmatbeg, int *cmatind, double *cmatval,
1683 double *lb, double *ub, char *where_to_move)
1684 {
1685 cpx_status = CPXaddcols(lp_data->cpxenv, lp_data->lp, ccnt, nzcnt,
1686 obj, cmatbeg, cmatind, cmatval, lb, ub, NULL);
1687 CPX_check_error("add_cols");
1688 lp_data->n += ccnt;
1689 lp_data->nz += nzcnt;
1690 }
1691
1692 /*===========================================================================*/
1693
change_row(LPdata * lp_data,int row_ind,char sense,double rhs,double range)1694 void change_row(LPdata *lp_data, int row_ind,
1695 char sense, double rhs, double range)
1696 {
1697 cpx_status = CPXchgsense(lp_data->cpxenv, lp_data->lp, 1, &row_ind, &sense);
1698 CPX_check_error("change_row - CPXchgsense");
1699 cpx_status = CPXchgcoef(lp_data->cpxenv, lp_data->lp, row_ind, -1, rhs);
1700 CPX_check_error("change_row - CPXchgcoef");
1701 if (sense == 'R'){
1702 cpx_status = CPXchgcoef(lp_data->cpxenv, lp_data->lp, row_ind, -2,range);
1703 CPX_check_error("change_row - CPXchgcoef");
1704 }
1705 lp_data->lp_is_modified = LP_HAS_BEEN_MODIFIED;
1706 }
1707
1708 /*===========================================================================*/
1709
change_col(LPdata * lp_data,int col_ind,char sense,double lb,double ub)1710 void change_col(LPdata *lp_data, int col_ind,
1711 char sense, double lb, double ub)
1712 {
1713 switch (sense){
1714 case 'E': change_lbub(lp_data, col_ind, lb, ub); break;
1715 case 'R': change_lbub(lp_data, col_ind, lb, ub); break;
1716 case 'G': change_lb(lp_data, col_ind, lb); break;
1717 case 'L': change_ub(lp_data, col_ind, ub); break;
1718 }
1719 lp_data->lp_is_modified = LP_HAS_BEEN_MODIFIED;
1720 }
1721
1722 /*===========================================================================*/
1723
1724 /*===========================================================================*\
1725 * Solve the lp specified in lp_data->lp with dual simplex. The number of
1726 * iterations is returned in 'iterd'. The return value of the function is
1727 * the termination code of the dual simplex method.
1728 \*===========================================================================*/
1729
dual_simplex(LPdata * lp_data,int * iterd)1730 int dual_simplex(LPdata *lp_data, int *iterd)
1731 {
1732 int real_term, term, itlim, defit, minit, maxit;
1733 double objulim, objllim, defobj;
1734
1735 if (lp_data->lp_is_modified == LP_HAS_BEEN_ABANDONED){
1736 cpx_status = CPXsetintparam(lp_data->cpxenv, CPX_PARAM_ADVIND, CPX_OFF);
1737 CPX_check_error("dual_simplex - CPXsetintparam, ADVIND");
1738 }
1739
1740 term = CPXdualopt(lp_data->cpxenv, lp_data->lp);
1741 if (term == CPXERR_PRESLV_INForUNBD){
1742 cpx_status = CPXsetintparam(lp_data->cpxenv, CPX_PARAM_PREIND, CPX_OFF);
1743 CPX_check_error("dual_simplex - CPXsetintparam");
1744 term = CPXdualopt(lp_data->cpxenv, lp_data->lp);
1745 CPX_check_error("dual_simplex - CPXdualopt");
1746 cpx_status = CPXsetintparam(lp_data->cpxenv, CPX_PARAM_PREIND, CPX_ON);
1747 CPX_check_error("dual_simplex - CPXsetintparam");
1748 }
1749
1750 term = CPXgetstat(lp_data->cpxenv,lp_data->lp);
1751 #if CPX_VERSION >= 800
1752 if (term == CPX_STAT_UNBOUNDED){
1753 /* } to unconfuse vi */
1754 #else
1755 if (term == CPX_INFEASIBLE){
1756 #endif
1757 /* Dual infeas. This is impossible, so we must have had iteration
1758 * limit AND bound shifting AND dual feasibility not restored within
1759 * the given iteration limit. */
1760 cpx_status = CPXgetintparam(lp_data->cpxenv, CPX_PARAM_ITLIM, &itlim);
1761 CPX_check_error("dual_simplex - CPXgetintparam, ITLIM");
1762 cpx_status = CPXinfointparam(lp_data->cpxenv, CPX_PARAM_ITLIM,
1763 &defit, &minit, &maxit);
1764 CPX_check_error("dual_simplex - CPXinfointparam, ITLIM");
1765 cpx_status = CPXsetintparam(lp_data->cpxenv, CPX_PARAM_ITLIM, defit);
1766 CPX_check_error("dual_simplex - CPXsetintparam, ITLIM");
1767 cpx_status = CPXgetdblparam(lp_data->cpxenv, CPX_PARAM_OBJULIM,&objulim);
1768 CPX_check_error("dual_simplex - CPXgetdblparam, OBJULIM");
1769 cpx_status = CPXgetdblparam(lp_data->cpxenv, CPX_PARAM_OBJULIM,&objllim);
1770 CPX_check_error("dual_simplex - CPXgetdblparam, OBJULIM");
1771 defobj = 1e75;
1772 cpx_status = CPXsetdblparam(lp_data->cpxenv, CPX_PARAM_OBJULIM, defobj);
1773 CPX_check_error("dual_simplex - CPXsetdblparam, OBJULIM");
1774 defobj = -1e75;
1775 cpx_status = CPXsetdblparam(lp_data->cpxenv, CPX_PARAM_OBJLLIM, defobj);
1776 CPX_check_error("dual_simplex - CPXsetdblparam, OBJLLIM");
1777 term = CPXdualopt(lp_data->cpxenv, lp_data->lp);
1778 cpx_status = CPXsetdblparam(lp_data->cpxenv, CPX_PARAM_OBJULIM, objulim);
1779 CPX_check_error("dual_simplex - CPXsetdblparam, OBJULIM");
1780 cpx_status = CPXsetdblparam(lp_data->cpxenv, CPX_PARAM_OBJLLIM, objllim);
1781 CPX_check_error("dual_simplex - CPXsetdblparam, OBJLLIM");
1782 cpx_status = CPXsetintparam(lp_data->cpxenv, CPX_PARAM_ITLIM, itlim);
1783 CPX_check_error("dual_simplex - CPXsetintparam, ITLIM");
1784 }
1785
1786 #if CPX_VERSION >= 800
1787 switch (real_term = CPXgetstat(lp_data->cpxenv,lp_data->lp)){
1788 case CPX_STAT_OPTIMAL: term = LP_OPTIMAL; break;
1789 case CPX_STAT_INFEASIBLE: term = LP_D_UNBOUNDED; break;
1790 case CPX_STAT_UNBOUNDED: term = LP_D_INFEASIBLE; break;
1791 case CPX_STAT_ABORT_OBJ_LIM: term = LP_D_OBJLIM; break;
1792 case CPX_STAT_ABORT_IT_LIM: term = LP_D_ITLIM; break;
1793 default: term = LP_ABANDONED; break;
1794 }
1795 #else
1796 switch (real_term = CPXgetstat(lp_data->cpxenv,lp_data->lp)){
1797 case CPX_OPTIMAL: term = LP_OPTIMAL; break;
1798 case CPX_INFEASIBLE: term = LP_D_INFEASIBLE; break;
1799 case CPX_UNBOUNDED: term = LP_D_UNBOUNDED; break;
1800 case CPX_OBJ_LIM: term = LP_D_OBJLIM; break;
1801 case CPX_IT_LIM_FEAS: case CPX_IT_LIM_INFEAS: term = LP_D_ITLIM; break;
1802 default: term = LP_ABANDONED; break;
1803 }
1804 #endif
1805
1806 lp_data->termcode = term;
1807
1808 if (term != LP_ABANDONED){
1809 *iterd = CPXgetitcnt(lp_data->cpxenv, lp_data->lp);
1810 cpx_status = CPXgetobjval(lp_data->cpxenv,lp_data->lp, &lp_data->objval);
1811 CPX_check_error("dual_simplex - CPXgetobjval");
1812 cpx_status = CPXsetintparam(lp_data->cpxenv, CPX_PARAM_ADVIND, CPX_ON);
1813 CPX_check_error("dual_simplex - CPXsetintparam, ADVIND");
1814 lp_data->lp_is_modified = LP_HAS_NOT_BEEN_MODIFIED;
1815 }else{
1816 lp_data->lp_is_modified = LP_HAS_BEEN_ABANDONED;
1817 printf("CPLEX Abandoned calculation: Code %i \n\n", real_term);
1818 }
1819 return(term);
1820 }
1821
1822 /*===========================================================================*/
1823 int solve_hotstart(LPdata *lp_data, int *iterd)
1824 {
1825 return(dual_simplex(lp_data,iterd));
1826 }
1827 /*===========================================================================*/
1828 int unmark_hotstart(LPdata *lp_data)
1829 {
1830 /* only when using osi */
1831 return (0);
1832 }
1833
1834 /*===========================================================================*/
1835 int mark_hotstart(LPdata *lp_data)
1836 {
1837 /* only when using osi */
1838 return (0);
1839 }
1840
1841
1842
1843 /*===========================================================================*/
1844
1845 void btran(LPdata *lp_data, double *col)
1846 {
1847 cpx_status = CPXbtran(lp_data->cpxenv, lp_data->lp, col);
1848 CPX_check_error("btran");
1849 }
1850
1851 /*===========================================================================*/
1852
1853 void get_binvcol(LPdata *lp_data, int j, double *col)
1854 {
1855 cpx_status = CPXbinvcol(lp_data->cpxenv, lp_data->lp, j, col);
1856 CPX_check_error("get_binvcol");
1857 }
1858
1859 /*===========================================================================*/
1860
1861 void get_binvrow(LPdata *lp_data, int i, double *row)
1862 {
1863 cpx_status = CPXbinvrow(lp_data->cpxenv, lp_data->lp, i, row);
1864 CPX_check_error("get_binvrow");
1865 }
1866
1867 /*===========================================================================*/
1868
1869 void get_basis(LPdata *lp_data, int *cstat, int *rstat)
1870 {
1871 cpx_status = CPXgetbase(lp_data->cpxenv, lp_data->lp, cstat, rstat);
1872 CPX_check_error("get_basis");
1873 }
1874
1875 /*===========================================================================*/
1876
1877 /*===========================================================================*\
1878 * Set an upper limit on the objective function value. Call the 'setobjulim'
1879 * CPLEX function.
1880 \*===========================================================================*/
1881
1882 void set_obj_upper_lim(LPdata *lp_data, double lim)
1883 {
1884 cpx_status = CPXsetdblparam(lp_data->cpxenv, CPX_PARAM_OBJULIM, lim);
1885 CPX_check_error("set_obj_upper_lim");
1886 }
1887
1888 /*===========================================================================*/
1889
1890 /*===========================================================================*\
1891 * Set an upper limit on the number of iterations. If itlim < 0 then set
1892 * it to the maximum.
1893 \*===========================================================================*/
1894
1895 void set_itlim(LPdata *lp_data, int itlim)
1896 {
1897 if (itlim < 0)
1898 cpx_status = CPXinfointparam(lp_data->cpxenv,
1899 CPX_PARAM_ITLIM, &itlim, NULL, NULL);
1900 CPX_check_error("set_itlim - CPXinfointparam");
1901 cpx_status = CPXsetintparam(lp_data->cpxenv, CPX_PARAM_ITLIM, itlim);
1902 CPX_check_error("set_itlim - CPXsetintparam");
1903 }
1904 /*===========================================================================*/
1905 void set_itlim_hotstart(LPdata *lp_data, int itlim)
1906 {
1907 /* read being and nothingness -- Jean Paul Sartre */
1908 }
1909
1910
1911 /*===========================================================================*/
1912
1913 void get_column(LPdata *lp_data, int j,
1914 double *colval, int *colind, int *collen, double *cj)
1915 {
1916 int matbeg, surplus;
1917 /* If there was no scaling, then we could probably copy the data out
1918 * directly. Try sometime... */
1919 cpx_status = CPXgetcols(lp_data->cpxenv, lp_data->lp, collen, &matbeg,
1920 colind, colval, lp_data->m, &surplus, j, j);
1921 CPX_check_error("get_column - CPXgetcols");
1922 cpx_status = CPXgetobj(lp_data->cpxenv, lp_data->lp, cj, j, j);
1923 CPX_check_error("get_column - CPXgetobj");
1924 }
1925
1926 /*===========================================================================*/
1927
1928 void get_row(LPdata *lp_data, int i,
1929 double *rowval, int *rowind, int *rowlen,
1930 double *rowub, double *rowlb)
1931 {
1932 int rmatbeg, surplus;
1933 /* If there was no scaling, then we could probably copy the data out
1934 * directly. Try sometime... */
1935 cpx_status = CPXgetrows(lp_data->cpxenv, lp_data->lp, rowlen, &rmatbeg,
1936 rowind, rowval, lp_data->n, &surplus, i, i);
1937 CPX_check_error("get_row - CPXgetrows");
1938 }
1939
1940 /*===========================================================================*/
1941
1942 /* This routine returns the index of a row which proves the lp to be primal
1943 * infeasible. There must be one, or this function wouldn't be called. */
1944 /* There MUST be something better than this...
1945 * A function call perhaps... Ask CPLEX... */
1946 int get_proof_of_infeas(LPdata *lp_data, int *infind)
1947 {
1948 int idiv, jdiv;
1949 double bd;
1950
1951 #if 0
1952 /*something like this should work...*/
1953 CPXdualfarkas(lp_data->cpxenv, lp_data->lp, ...);
1954 CPX_check_error("get_proof_of_infeas - CPXdualfarkas");
1955 #endif
1956
1957 CPXgetijdiv(lp_data->cpxenv, lp_data->lp, &idiv, &jdiv);
1958 CPX_check_error("get_proof_of_infeas - CPXgetijdiv");
1959 cpx_status = CPXgetijrow(lp_data->cpxenv, lp_data->lp, idiv, jdiv, infind);
1960 CPX_check_error("get_proof_of_infeas - CPXgetijrow");
1961 if (cpx_status)
1962 return(0);
1963 if (jdiv < 0){ /* the diverging variable is a slack/range */
1964 if (lp_data->slacks)
1965 return(lp_data->slacks[idiv] < 0 ? LOWER_THAN_LB : HIGHER_THAN_UB);
1966 }else{ /* the diverging variable is structural */
1967 cpx_status = CPXgetlb(lp_data->cpxenv, lp_data->lp, &bd, jdiv, jdiv);
1968 CPX_check_error("get_proof_of_infeas - CPXgetlb");
1969 if (lp_data->x){
1970 return(bd < lp_data->x[jdiv] ? LOWER_THAN_LB : HIGHER_THAN_UB);
1971 }
1972 }
1973 return(0); /* fake return */
1974 }
1975
1976 /*===========================================================================*/
1977
1978 /*===========================================================================*\
1979 * Get the solution (values of the structural variables in an optimal
1980 * solution) to the lp (specified by lp_data->lp) into the vector
1981 * lp_data->x. This can be done by calling the 'getx' CPLEX function.
1982 \*===========================================================================*/
1983
1984 void get_x(LPdata *lp_data)
1985 {
1986 cpx_status = CPXgetx(lp_data->cpxenv, lp_data->lp, lp_data->x, 0,
1987 lp_data->n-1);
1988 CPX_check_error("get_x");
1989 }
1990
1991 /*===========================================================================*/
1992
1993 void get_dj_pi(LPdata *lp_data)
1994 {
1995 cpx_status = CPXgetpi(lp_data->cpxenv, lp_data->lp, lp_data->dualsol, 0,
1996 lp_data->m-1);
1997 CPX_check_error("get_dj_pi - CPXgetpi");
1998 cpx_status = CPXgetdj(lp_data->cpxenv, lp_data->lp, lp_data->dj, 0,
1999 lp_data->n-1);
2000 CPX_check_error("get_dj_pi - CPXgetdj");
2001 }
2002
2003 /*===========================================================================*/
2004
2005 void get_slacks(LPdata *lp_data)
2006 {
2007 row_data *rows = lp_data->rows;
2008 double *slacks = lp_data->slacks;
2009 int i, m = lp_data->m;
2010 cpx_status = CPXgetslack(lp_data->cpxenv, lp_data->lp, lp_data->slacks, 0,
2011 lp_data->m-1);
2012 CPX_check_error("get_slacks");
2013 /* Compute the real slacks for the free rows */
2014 for (i=m-1; i>=0; i--){
2015 if (rows[i].free){
2016 switch (rows[i].cut->sense){
2017 case 'E': slacks[i] += rows[i].cut->rhs - SYM_INFINITY; break;
2018 case 'L': slacks[i] += rows[i].cut->rhs - SYM_INFINITY; break;
2019 case 'G': slacks[i] += rows[i].cut->rhs + SYM_INFINITY; break;
2020 case 'R': slacks[i] += -rows[i].cut->rhs - SYM_INFINITY; break;
2021 }
2022 }
2023 }
2024
2025 }
2026
2027 /*===========================================================================*/
2028
2029 void change_range(LPdata *lp_data, int rowind, double value)
2030 {
2031 cpx_status = CPXchgcoef(lp_data->cpxenv, lp_data->lp, rowind, -2, value);
2032 CPX_check_error("change_range");
2033 lp_data->lp_is_modified = LP_HAS_BEEN_MODIFIED;
2034 }
2035
2036 /*===========================================================================*/
2037
2038 void change_rhs(LPdata *lp_data, int rownum, int *rhsind, double *rhsval)
2039 {
2040 cpx_status = CPXchgrhs(lp_data->cpxenv, lp_data->lp, rownum, rhsind,rhsval);
2041 CPX_check_error("change_rhs");
2042 lp_data->lp_is_modified = LP_HAS_BEEN_MODIFIED;
2043 }
2044
2045 /*===========================================================================*/
2046
2047 void change_sense(LPdata *lp_data, int cnt, int *index, char *sense)
2048 {
2049 cpx_status = CPXchgsense(lp_data->cpxenv, lp_data->lp, cnt, index, sense);
2050 CPX_check_error("change_sense");
2051 lp_data->lp_is_modified = LP_HAS_BEEN_MODIFIED;
2052 }
2053
2054 /*===========================================================================*/
2055
2056 void change_bounds(LPdata *lp_data, int cnt, int *index, char *lu, double *bd)
2057 {
2058 cpx_status = CPXchgbds(lp_data->cpxenv, lp_data->lp, cnt, index, lu, bd);
2059 CPX_check_error("change_bounds");
2060 lp_data->lp_is_modified = LP_HAS_BEEN_MODIFIED;
2061 }
2062
2063 /*===========================================================================*/
2064
2065 void change_lbub(LPdata *lp_data, int j, double lb, double ub)
2066 {
2067 int ind[2];
2068 double bd[2];
2069 ind[0] = ind[1] = j;
2070 bd[0] = lb; bd[1] = ub;
2071 cpx_status =
2072 CPXchgbds(lp_data->cpxenv, lp_data->lp, 2, ind, (char *)"LU", bd);
2073 CPX_check_error("change_lbub");
2074 lp_data->lp_is_modified = LP_HAS_BEEN_MODIFIED;
2075 }
2076
2077 /*===========================================================================*/
2078
2079 void change_ub(LPdata *lp_data, int j, double ub)
2080 {
2081 cpx_status = CPXchgbds(lp_data->cpxenv, lp_data->lp, 1, &j, (char *)"U",
2082 &ub);
2083 CPX_check_error("change_ub");
2084 lp_data->lp_is_modified = LP_HAS_BEEN_MODIFIED;
2085 }
2086
2087 /*===========================================================================*/
2088
2089 void change_lb(LPdata *lp_data, int j, double lb)
2090 {
2091 cpx_status = CPXchgbds(lp_data->cpxenv, lp_data->lp, 1, &j, (char *)"L",
2092 &lb);
2093 CPX_check_error("change_lb");
2094 lp_data->lp_is_modified = LP_HAS_BEEN_MODIFIED;
2095 }
2096
2097 /*===========================================================================*/
2098
2099 void get_ub(LPdata *lp_data, int j, double *ub)
2100 {
2101 cpx_status = CPXgetub(lp_data->cpxenv, lp_data->lp, ub, j, j);
2102 CPX_check_error("get_ub");
2103 }
2104
2105 /*===========================================================================*/
2106
2107 void get_lb(LPdata *lp_data, int j, double *lb)
2108 {
2109 cpx_status = CPXgetlb(lp_data->cpxenv, lp_data->lp, lb, j, j);
2110 CPX_check_error("get_lb");
2111 }
2112
2113 /*===========================================================================*/
2114
2115 void get_bounds(LPdata *lp_data)
2116 {
2117 if (!lp_data->lb){
2118 }
2119 cpx_status = CPXgetlb(lp_data->cpxenv, lp_data->lp, lp_data->lb, 0,
2120 lp_data->n-1);
2121 CPX_check_error("get_lb");
2122 cpx_status = CPXgetub(lp_data->cpxenv, lp_data->lp, lp_data->ub, 0,
2123 lp_data->n-1);
2124 CPX_check_error("get_ub");
2125
2126 }
2127
2128 /*===========================================================================*/
2129
2130 void get_objcoef(LPdata *lp_data, int j, double *objcoef)
2131 {
2132 cpx_status = CPXgetobj(lp_data->cpxenv, lp_data->lp, objcoef, j, j);
2133 CPX_check_error("get_objcoef");
2134 }
2135
2136 /*===========================================================================*/
2137
2138 void delete_rows(LPdata *lp_data, int deletable, int *free_rows)
2139 {
2140 cpx_status = CPXdelsetrows(lp_data->cpxenv, lp_data->lp, free_rows);
2141 CPX_check_error("delete_rows");
2142 lp_data->m -= deletable;
2143 lp_data->lp_is_modified = LP_HAS_BEEN_MODIFIED;
2144 }
2145
2146 /*===========================================================================*/
2147
2148 int delete_cols(LPdata *lp_data, int delnum, int *delstat)
2149 {
2150 double *dj = lp_data->dj;
2151 double *x = lp_data->x;
2152 char *status = lp_data->status;
2153 int i, num_to_keep;
2154
2155 cpx_status = CPXdelsetcols(lp_data->cpxenv, lp_data->lp, delstat);
2156 CPX_check_error("delete_cols - CPXdelsetcols");
2157 lp_data->nz = CPXgetnumnz(lp_data->cpxenv, lp_data->lp);
2158 CPX_check_error("delete_cols - CPXgetnumnz");
2159
2160 for (i = 0, num_to_keep = 0; i < lp_data->n; i++){
2161 if (delstat[i] != -1){
2162 dj[delstat[i]] = dj[i];
2163 x[delstat[i]] = x[i];
2164 status[delstat[i]] = status[i];
2165 num_to_keep++;
2166 }
2167 }
2168
2169 lp_data->n = num_to_keep;
2170
2171 return(delnum);
2172 }
2173
2174 /*===========================================================================*/
2175
2176 void release_var(LPdata *lp_data, int j, int where_to_move)
2177 {
2178 #if 0
2179 switch (where_to_move){
2180 case MOVE_TO_UB:
2181 lp_data->lpbas.cstat[j] = 2; break; /* non-basic at its upper bound */
2182 case MOVE_TO_LB:
2183 lp_data->lpbas.cstat[j] = 0; break; /* non-basic at its lower bound */
2184 }
2185 lp_data->lp_is_modified = LP_HAS_BEEN_MODIFIED;
2186 #endif
2187 }
2188
2189 /*===========================================================================*/
2190
2191 void free_row_set(LPdata *lp_data, int length, int *index)
2192 {
2193 int i, j;
2194 row_data *rows = lp_data->rows;
2195 double *rhsval = lp_data->tmp.d; /* m */
2196 int *ind_e = lp_data->tmp.i1 + 2 * lp_data->m; /* m (now) */
2197 /* See comment in check_row_effectiveness why the shift! */
2198 char *sen_e = lp_data->tmp.c; /* m (now) */
2199
2200 for (j=0, i=length-1; i>=0; i--){
2201 switch (rows[index[i]].cut->sense){
2202 case 'E': rhsval[i] = SYM_INFINITY; ind_e[j++] = index[i]; break;
2203 case 'L': rhsval[i] = SYM_INFINITY; break;
2204 case 'R':
2205 cpx_status = CPXchgcoef(lp_data->cpxenv, lp_data->lp, index[i], -2,
2206 2*SYM_INFINITY);
2207 CPX_check_error("free_row_set - CPXchgcoef");
2208 case 'G': rhsval[i] = -SYM_INFINITY; break;
2209 }
2210 }
2211 cpx_status = CPXchgrhs(lp_data->cpxenv, lp_data->lp, length, index, rhsval);
2212 CPX_check_error("free_row_set - CPXchgrhs");
2213 if (j > 0){
2214 memset(sen_e, 'L', j);
2215 cpx_status = CPXchgsense(lp_data->cpxenv, lp_data->lp, j, ind_e, sen_e);
2216 CPX_check_error("free_row_set - CPXchgsense");
2217 }
2218 lp_data->lp_is_modified = LP_HAS_BEEN_MODIFIED;
2219 }
2220
2221 /*===========================================================================*/
2222
2223 void constrain_row_set(LPdata *lp_data, int length, int *index)
2224 {
2225 int i;
2226 row_data *rows = lp_data->rows;
2227 cut_data *cut;
2228 double *rhsval = lp_data->tmp.d; /* m (now) */
2229 char *sense = lp_data->tmp.c + lp_data->m; /* m (now) */
2230 char range_constraint = FALSE;
2231
2232 for (i = length-1; i >= 0; i--){
2233 cut = rows[index[i]].cut;
2234 rhsval[i] = cut->rhs;
2235 if ((sense[i] = cut->sense) == 'R'){
2236 range_constraint = TRUE;
2237 }
2238 }
2239 cpx_status = CPXchgrhs(lp_data->cpxenv, lp_data->lp, length, index, rhsval);
2240 CPX_check_error("constrain_row_set - CPXchgrhs");
2241 cpx_status=CPXchgsense(lp_data->cpxenv, lp_data->lp, length, index, sense);
2242 CPX_check_error("constrain_row_set - CPXchgsense");
2243 if (range_constraint){
2244 for (i = length-1; i >= 0; i--){
2245 if (sense[i] == 'R'){
2246 cpx_status = CPXchgcoef(lp_data->cpxenv,lp_data->lp, index[i], -2,
2247 rows[index[i]].cut->range);
2248 CPX_check_error("constrain_row_set - CPXchgcoef");
2249 }
2250 }
2251 }
2252 lp_data->lp_is_modified = LP_HAS_BEEN_MODIFIED;
2253 }
2254
2255 /*===========================================================================*/
2256
2257 int read_mps(MIPdesc *desc, char *infile, char *probname, int verbosity)
2258 {
2259 printf("\nMps-format file can be read only through OSI interface.\n");
2260
2261 return(1);
2262 }
2263
2264 /*===========================================================================*/
2265
2266 int read_lp(MIPdesc *desc, char *infile, char *probname, int verbosity)
2267 {
2268 printf("\nLP-format file can be read only through OSI interface.\n");
2269
2270 return(1);
2271 }
2272
2273 /*===========================================================================*/
2274
2275 void write_mip_desc_mps(MIPdesc *mip, char *fname)
2276 {
2277 fprintf(stderr, "Function not implemented yet.");
2278 exit(-1);
2279 }
2280
2281 /*===========================================================================*/
2282
2283 void write_mip_desc_lp(MIPdesc *mip, char *fname)
2284 {
2285 fprintf(stderr, "Function not implemented yet.");
2286 exit(-1);
2287 }
2288
2289 /*===========================================================================*/
2290
2291 void write_mps(LPdata *lp_data, char *fname)
2292 {
2293 cpx_status = CPXmpswrite(lp_data->cpxenv, lp_data->lp, fname);
2294 CPX_check_error("write_mps");
2295 }
2296
2297 /*===========================================================================*/
2298
2299 void write_sav(LPdata *lp_data, char *fname)
2300 {
2301 cpx_status = CPXsavwrite(lp_data->cpxenv, lp_data->lp, fname);
2302 CPX_check_error("write_sav");
2303 }
2304
2305 /*===========================================================================*/
2306
2307 #ifdef USE_CGL_CUTS
2308 void generate_cgl_cuts(LPdata *lp_data, int *num_cuts, cut_data ***cuts,
2309 char send_to_pool)
2310 {
2311 return;
2312 }
2313 #endif
2314
2315 #endif /* __CPLEX__ */
2316
2317 #if defined(__OSI_CPLEX__) || defined(__OSI_OSL__) || defined(__OSI_CLP__) \
2318 || defined(__OSI_XPRESS__) || defined(__OSI_SOPLEX__) || defined(__OSI_VOL__) \
2319 || defined(__OSI_DYLP__) || defined (__OSI_GLPK__)
2320
2321 static bool retval = false;
2322
2323 void open_lp_solver(LPdata *lp_data)
2324 {
2325 lp_data->si = new OsiXSolverInterface();
2326
2327 /* Turn off the OSL messages (There are LOTS of them) */
2328 lp_data->si->setHintParam(OsiDoReducePrint);
2329 lp_data->si->messageHandler()->setLogLevel(0);
2330 #ifdef __OSI_CLP__
2331 lp_data->si->setupForRepeatedUse();
2332 //lp_data->si->setupForRepeatedUse(3,0);
2333 //lp_data->si->getModelPtr()->setFactorizationFrequency(200);
2334 //lp_data->si->getModelPtr()->setSparseFactorization(true);
2335 //lp_data->si->getModelPtr()->setSpecialOptions(524288);
2336 //lp_data->si->getModelPtr()->setSpecialOptions(4);
2337 lp_data->si->getModelPtr()->setPerturbation(50);
2338 //set cleanup param if unscaled primal is infeasible
2339 lp_data->si->setCleanupScaling(1);
2340 #endif
2341 #ifdef __OSI_GLPK__
2342 lp_data->lpetol = 1e-07; /* glpk doesn't return the value of this param */
2343 #else
2344 lp_data->si->getDblParam(OsiPrimalTolerance, lp_data->lpetol);
2345 #endif
2346 }
2347
2348 /*===========================================================================*/
2349
2350 void close_lp_solver(LPdata *lp_data)
2351 {
2352 delete lp_data->si;
2353 }
2354
2355 /*===========================================================================*/
2356
2357 /*===========================================================================*\
2358 * This function loads the data of an lp into the lp solver. This involves
2359 * transforming the data into CPLEX format and calling the CPLEX function
2360 * 'loadlp'.
2361 \*===========================================================================*/
2362
2363 void load_lp_prob(LPdata *lp_data, int scaling, int fastmip)
2364 {
2365
2366 /* Turn off scaling for CLP */
2367 //lp_data->si->setHintParam(OsiDoScale,false,OsiHintDo);
2368 MIPdesc *mip = lp_data->mip;
2369
2370 lp_data->si->loadProblem(lp_data->n, lp_data->m,
2371 mip->matbeg, mip->matind,
2372 mip->matval, mip->lb,
2373 mip->ub, mip->obj,
2374 mip->sense, mip->rhs,
2375 mip->rngval);
2376 }
2377
2378 /*===========================================================================*/
2379 int reset_lp_prob(LPdata *lp_data, int scaling, int fastmip)
2380 {
2381 lp_data->si->restoreBaseModel(lp_data->m);
2382 return 0;
2383 }
2384
2385 /*===========================================================================*/
2386 int save_lp(LPdata *lp_data)
2387 {
2388 lp_data->si->saveBaseModel();
2389 return 0;
2390 }
2391 /*===========================================================================*/
2392
2393 void unload_lp_prob(LPdata *lp_data)
2394 {
2395
2396 //lp_data->si->reset();
2397
2398 /* Set parameters as in open_lp_solver() (do these persist?) */
2399 lp_data->si->setHintParam(OsiDoReducePrint);
2400 lp_data->si->messageHandler()->setLogLevel(0);
2401 lp_data->m = lp_data->n = lp_data->nz = 0;
2402 }
2403
2404 /*===========================================================================*/
2405
2406 void load_basis(LPdata *lp_data, int *cstat, int *rstat)
2407 {
2408
2409 CoinWarmStartBasis *warmstart = new CoinWarmStartBasis;
2410
2411 int numcols = lp_data->n;
2412 int numrows = lp_data->m;
2413 int i;
2414
2415 warmstart->setSize(numcols, numrows);
2416
2417 for (i = 0; i < numrows; i++){
2418 switch (rstat[i]){
2419 case SLACK_AT_LB:
2420 warmstart->setArtifStatus(i,CoinWarmStartBasis::atLowerBound);
2421 break;
2422 case SLACK_BASIC:
2423 warmstart->setArtifStatus(i,CoinWarmStartBasis::basic);
2424 break;
2425 case SLACK_AT_UB:
2426 warmstart->setArtifStatus(i,CoinWarmStartBasis::atUpperBound);
2427 break;
2428 case SLACK_FREE:
2429 warmstart->setArtifStatus(i,CoinWarmStartBasis::isFree);
2430 break;
2431 default:
2432 break;
2433 }
2434 }
2435
2436 for (i = 0; i < numcols; i++){
2437 switch (cstat[i]){
2438 case VAR_AT_LB:
2439 warmstart->setStructStatus(i,CoinWarmStartBasis::atLowerBound);
2440 break;
2441 case VAR_BASIC:
2442 warmstart->setStructStatus(i,CoinWarmStartBasis::basic);
2443 break;
2444 case VAR_AT_UB:
2445 warmstart->setStructStatus(i,CoinWarmStartBasis::atUpperBound);
2446 break;
2447 case VAR_FREE:
2448 warmstart->setStructStatus(i,CoinWarmStartBasis::isFree);
2449 break;
2450 default:
2451 break;
2452 }
2453 }
2454
2455 retval = lp_data->si->setWarmStart(warmstart);
2456
2457 delete warmstart;
2458 }
2459
2460 /*===========================================================================*/
2461
2462 void add_rows(LPdata *lp_data, int rcnt, int nzcnt, double *rhs,
2463 char *sense, int *rmatbeg, int *rmatind, double *rmatval)
2464 {
2465 int i;// start, size;
2466 OsiXSolverInterface *si = lp_data->si;
2467 double *rlb = lp_data->tmp.d + rcnt;
2468 double *rub = lp_data->tmp.d + 2*rcnt;
2469 const double infinity = si->getInfinity();
2470
2471 /*
2472 for (i = 0; i < rcnt; i++){
2473 CoinPackedVector new_row;
2474 for (j = rmatbeg[i]; j < rmatbeg[i+1]; j++){
2475 new_row.insert(rmatind[j], rmatval[j]);
2476 }
2477 si->addRow(new_row, sense[i], rhs[i], 0);
2478 }
2479 */
2480 #if 0
2481 for (i = 0; i < rcnt; i++){
2482 start = rmatbeg[i];
2483 size = rmatbeg[i+1] - start;
2484 CoinPackedVector new_row(size, &rmatind[start], &rmatval[start], FALSE);
2485 si->addRow(new_row, sense[i], rhs[i], 0);
2486 }
2487 #else
2488 /* convert sense to bound */
2489 for(i = 0; i < rcnt; i++){
2490 switch (sense[i]){
2491 printf("%c \n", sense[i]);
2492 case 'E':
2493 rlb[i] = rub[i] = rhs[i];
2494 break;
2495 case 'L':
2496 rlb[i] = -infinity;
2497 rub[i] = rhs[i];
2498 break;
2499 case 'G':
2500 rlb[i] = rhs[i];
2501 rub[i] = infinity;
2502 break;
2503 case 'R': // we should not have this case
2504 rlb[i] = -infinity;
2505 rub[i] = rhs[i];
2506 break;
2507 case 'N':
2508 rlb[i] = -infinity;
2509 rub[i] = infinity;
2510 break;
2511 }
2512 }
2513
2514 si->addRows(rcnt, rmatbeg, rmatind, rmatval, rlb, rub);
2515 #endif
2516 lp_data->m += rcnt;
2517 lp_data->nz += nzcnt;
2518 lp_data->lp_is_modified=LP_HAS_BEEN_MODIFIED;
2519 }
2520
2521 /*===========================================================================*/
2522
2523 void add_cols(LPdata *lp_data, int ccnt, int nzcnt, double *obj,
2524 int *cmatbeg, int *cmatind, double *cmatval,
2525 double *lb, double *ub, char *where_to_move)
2526 {
2527 // TODO: eliminate the inner loop. its inefficient.
2528 int i, j;
2529 OsiXSolverInterface *si = lp_data->si;
2530 for (i = 0; i < ccnt; i++){
2531 CoinPackedVector col;
2532 for (j = cmatbeg[i]; j < cmatbeg[i+1]; j++)
2533 col.insert(cmatind[j], cmatval[j]);
2534 si->addCol(col, lb[i], ub[i], obj[i]);
2535 }
2536
2537 lp_data->n += ccnt;
2538 lp_data->nz += nzcnt;
2539 }
2540
2541 /*===========================================================================*/
2542
2543 void change_row(LPdata *lp_data, int row_ind,
2544 char sense, double rhs, double range)
2545 {
2546 lp_data->si->setRowType(row_ind, sense, rhs, range);
2547 }
2548
2549 /*===========================================================================*/
2550
2551 void change_col(LPdata *lp_data, int col_ind,
2552 char sense, double lb, double ub)
2553 {
2554 switch (sense){
2555 case 'E': change_lbub(lp_data, col_ind, lb, ub); break;
2556 case 'R': change_lbub(lp_data, col_ind, lb, ub); break;
2557 case 'G': change_lb(lp_data, col_ind, lb); break;
2558 case 'L': change_ub(lp_data, col_ind, ub); break;
2559 }
2560 }
2561
2562 /*===========================================================================*/
2563
2564 /*===========================================================================*\
2565 * Solve the lp specified in lp_data->lp with dual simplex. The number of
2566 * iterations is returned in 'iterd'. The return value of the function is
2567 * the termination code of the dual simplex method.
2568 \*===========================================================================*/
2569
2570 int initial_lp_solve (LPdata *lp_data, int *iterd)
2571 {
2572
2573 //int term = LP_ABANDONED;
2574 int term = 0;
2575 OsiXSolverInterface *si = lp_data->si;
2576
2577 //si->setHintParam(OsiDoPresolveInInitial, false, OsiHintDo);
2578
2579 si->initialSolve();
2580
2581 if (si->isProvenDualInfeasible()){
2582 term = LP_D_INFEASIBLE;
2583 }else if (si->isProvenPrimalInfeasible()){
2584 term = LP_D_UNBOUNDED;
2585 }else if (si->isDualObjectiveLimitReached()){
2586 term = LP_D_OBJLIM;
2587 }else if (si->isProvenOptimal()){
2588 term = LP_OPTIMAL;
2589 }else if (si->isIterationLimitReached()){
2590 term = LP_D_ITLIM;
2591 #ifdef __OSI_CLP__
2592 /* If max iterations and had switched to primal, bound is no good */
2593 if (si->getModelPtr()->secondaryStatus() == 10){
2594 term = LP_ABANDONED;
2595 }
2596 #endif
2597 }else if (si->isAbandoned()){
2598 term = LP_ABANDONED;
2599 }else{
2600 // Osi doesn't have a check for time limit or a way of setting it
2601 // This is the only posibility left.
2602 term = LP_TIME_LIMIT;
2603 }
2604
2605 lp_data->termcode = term;
2606
2607 if (term != LP_ABANDONED && term != LP_D_INFEASIBLE){
2608
2609 *iterd = si->getIterationCount();
2610
2611 lp_data->objval = si->getObjValue();
2612
2613 /* Get relevant data */
2614 if (lp_data->dualsol && lp_data->dj) {
2615 get_dj_pi(lp_data);
2616 }
2617 if (lp_data->slacks && term == LP_OPTIMAL) {
2618 get_slacks(lp_data);
2619 }
2620
2621 get_x(lp_data);
2622
2623 #ifdef CHECK_DUAL_SOLUTION
2624 if (term == LP_D_INFEASIBLE || term == LP_OPTIMAL) {
2625 //This code checks the dual solution values
2626 int t;
2627 double intercept = 0;
2628 double lb = 0;
2629
2630 for (t=0; t < lp_data->n; t++){
2631 intercept += lp_data->x[t]* lp_data->dj[t];
2632 }
2633 for (int i = 0; i <lp_data->m; i++){
2634 if (si->getRowUpper()[i] < 1000000){
2635 lb += si->getRowUpper()[i]*lp_data->dualsol[i];
2636 }else{
2637 lb += si->getRowLower()[i]*lp_data->dualsol[i];
2638 }
2639 }
2640
2641 if (fabs(intercept + lb - lp_data->objval) > 0.1){
2642 write_mps(lp_data, "lp.assert");
2643 }
2644
2645 assert(fabs(intercept + lb - lp_data->objval) <= 0.1);
2646 }
2647 #endif
2648
2649 lp_data->lp_is_modified = LP_HAS_NOT_BEEN_MODIFIED;
2650 }
2651 else{
2652 lp_data->lp_is_modified = LP_HAS_BEEN_ABANDONED;
2653 #ifdef __OSI_CLP__
2654 if (si->getModelPtr()->secondaryStatus() != 10)
2655 #endif
2656 printf("OSI Abandoned calculation: Code %i \n\n", term);
2657 }
2658
2659 /*
2660 si->getModelPtr()->tightenPrimalBounds(0.0,0,true);
2661 */
2662 return(term);
2663 }
2664
2665 /*===========================================================================*/
2666
2667 /*===========================================================================*\
2668 * Solve the lp specified in lp_data->lp with dual simplex. The number of
2669 * iterations is returned in 'iterd'. The return value of the function is
2670 * the termination code of the dual simplex method.
2671 \*===========================================================================*/
2672
2673
2674 int dual_simplex(LPdata *lp_data, int *iterd)
2675 {
2676
2677 int term = 0;
2678 OsiXSolverInterface *si = lp_data->si;
2679 #ifdef __OSI_CLP__
2680 int sp = si->specialOptions();
2681 if((sp&2) != 0) sp ^=2;
2682 si->setSpecialOptions(sp);
2683 //si->setSpecialOptions(0x80000000);
2684 si->getModelPtr()->setPerturbation(50);
2685 //si->getModelPtr()->setFactorizationFrequency(150);
2686 //si->getModelPtr()->setSubstitution(3);
2687 #endif
2688 si->resolve();
2689
2690 if (si->isProvenDualInfeasible()){
2691 term = LP_D_INFEASIBLE;
2692 }else if (si->isProvenPrimalInfeasible()){
2693 term = LP_D_UNBOUNDED;
2694 }else if (si->isDualObjectiveLimitReached()){
2695 term = LP_D_OBJLIM;
2696 }else if (si->isProvenOptimal()){
2697 term = LP_OPTIMAL;
2698 }else if (si->isIterationLimitReached()){
2699 term = LP_D_ITLIM;
2700 #ifdef __OSI_CLP__
2701 /* If max iterations and had switched to primal, bound is no good */
2702 if (si->getModelPtr()->secondaryStatus() == 10){
2703 term = LP_ABANDONED;
2704 }
2705 #endif
2706 }else if (si->isAbandoned()){
2707 term = LP_ABANDONED;
2708 }else{
2709 // Osi doesn't have a check for time limit or a way of setting it
2710 // This is the only posibility left.
2711 term = LP_TIME_LIMIT;
2712 }
2713
2714 lp_data->termcode = term;
2715
2716 if (term != LP_ABANDONED){
2717
2718 *iterd = si->getIterationCount();
2719
2720 lp_data->objval = si->getObjValue();
2721
2722 /* Get relevant data */
2723 if (lp_data->dualsol && lp_data->dj) {
2724 get_dj_pi(lp_data);
2725 }
2726 if (lp_data->slacks && term == LP_OPTIMAL) {
2727 get_slacks(lp_data);
2728 }
2729
2730 get_x(lp_data);
2731
2732 #ifdef CHECK_DUAL_SOLUTION
2733 if (term == LP_D_INFEASIBLE || term == LP_OPTIMAL) {
2734 //This code checks the dual solution values
2735 int t;
2736 double intercept = 0;
2737 double lb = 0;
2738
2739 for (t=0; t < lp_data->n; t++){
2740 intercept += lp_data->x[t]* lp_data->dj[t];
2741 }
2742 for (int i = 0; i <lp_data->m; i++){
2743 if (si->getRowUpper()[i] < 1000000){
2744 lb += si->getRowUpper()[i]*lp_data->dualsol[i];
2745 }else{
2746 lb += si->getRowLower()[i]*lp_data->dualsol[i];
2747 }
2748 }
2749
2750 if (fabs(intercept + lb - lp_data->objval) > 0.1){
2751 write_mps(lp_data, "lp.assert");
2752 }
2753
2754 assert(fabs(intercept + lb - lp_data->objval) <= 0.1);
2755 }
2756 #endif
2757
2758 lp_data->lp_is_modified = LP_HAS_NOT_BEEN_MODIFIED;
2759 }
2760 else{
2761 lp_data->lp_is_modified = LP_HAS_BEEN_ABANDONED;
2762 #ifdef __OSI_CLP__
2763 if (si->getModelPtr()->secondaryStatus() != 10)
2764 #endif
2765 printf("OSI Abandoned calculation: Code %i \n\n", term);
2766 }
2767
2768 /*
2769 si->getModelPtr()->tightenPrimalBounds(0.0,0,true);
2770 */
2771 return(term);
2772 }
2773
2774
2775 /*===========================================================================*/
2776 /*
2777 * Following hot-start functions make it faster for the lp solver to do strong
2778 * branching
2779 */
2780 /*===========================================================================*/
2781 int solve_hotstart(LPdata *lp_data, int *iterd)
2782 {
2783
2784 //int term = LP_ABANDONED;
2785 int term = 0;
2786 OsiXSolverInterface *si = lp_data->si;
2787
2788 si->solveFromHotStart();
2789
2790 if (si->isProvenDualInfeasible())
2791 term = LP_D_INFEASIBLE;
2792 else if (si->isProvenPrimalInfeasible())
2793 term = LP_D_UNBOUNDED;
2794 else if (si->isDualObjectiveLimitReached())
2795 term = LP_D_OBJLIM;
2796 else if (si->isProvenOptimal())
2797 term = LP_OPTIMAL;
2798 else if (si->isIterationLimitReached())
2799 term = LP_D_ITLIM;
2800 else if (si->isAbandoned())
2801 term = LP_ABANDONED;
2802
2803 /* if(term == D_UNBOUNDED){
2804 retval=si->getIntParam(OsiMaxNumIteration, itlim);
2805 CAN NOT GET DEFAULT, MIN VALUES in OSI of CPXinfointparam() */
2806 /* } to unconfuse vi */
2807
2808 lp_data->termcode = term;
2809
2810 if (term != LP_ABANDONED){
2811
2812 *iterd = si->getIterationCount();
2813
2814 lp_data->objval = si->getObjValue();
2815
2816 /* Get relevant data */
2817 if (lp_data->dualsol && lp_data->dj) {
2818 get_dj_pi(lp_data);
2819 }
2820 if (lp_data->slacks && term == LP_OPTIMAL) {
2821 get_slacks(lp_data);
2822 }
2823
2824 get_x(lp_data);
2825
2826 #ifdef CHECK_DUAL_SOLUTION
2827 if (term == LP_D_INFEASIBLE || term == LP_OPTIMAL) {
2828 //This code checks the dual solution values
2829 int t;
2830 double intercept = 0;
2831 double lb = 0;
2832
2833 for (t=0; t < lp_data->n; t++){
2834 intercept += lp_data->x[t]* lp_data->dj[t];
2835 }
2836 for (int i = 0; i <lp_data->m; i++){
2837 if (si->getRowUpper()[i] < 1000000){
2838 lb += si->getRowUpper()[i]*lp_data->dualsol[i];
2839 }else{
2840 lb += si->getRowLower()[i]*lp_data->dualsol[i];
2841 }
2842 }
2843
2844 if (fabs(intercept + lb - lp_data->objval) > 0.1){
2845 write_mps(lp_data, "lp.assert");
2846 }
2847
2848 assert(fabs(intercept + lb - lp_data->objval) <= 0.1);
2849 }
2850 #endif
2851
2852 lp_data->lp_is_modified = LP_HAS_NOT_BEEN_MODIFIED;
2853 }
2854 else{
2855 lp_data->lp_is_modified = LP_HAS_BEEN_ABANDONED;
2856 printf("OSI Abandoned calculation: Code %i \n\n", term);
2857 }
2858
2859 return(term);
2860 }
2861
2862 /*===========================================================================*/
2863 int mark_hotstart(LPdata *lp_data)
2864 {
2865 lp_data->si->markHotStart();
2866 return (0);
2867 }
2868
2869 /*===========================================================================*/
2870 int unmark_hotstart(LPdata *lp_data)
2871 {
2872 lp_data->si->unmarkHotStart();
2873 return (0);
2874 }
2875
2876 /*===========================================================================*/
2877 /*===========================================================================*/
2878 /* This function is used only together with get_proof_of_infeasibility... */
2879
2880 void get_binvrow(LPdata *lp_data, int i, double *row)
2881 {
2882 fprintf(stderr, "Function not implemented yet.");
2883 exit(-1);
2884 }
2885
2886 /*===========================================================================*/
2887
2888 void get_basis(LPdata *lp_data, int *cstat, int *rstat)
2889 {
2890 CoinWarmStart * warmstart = lp_data->si->getWarmStart();
2891
2892 CoinWarmStartBasis * ws = dynamic_cast<CoinWarmStartBasis*>(warmstart);
2893
2894 int numcols = ws->getNumStructural(); /* has to be <= lp_data->n */
2895 int numrows = ws->getNumArtificial(); /* has to be <= lp_data->m */
2896 int i; /* hence an assert? */
2897
2898 if (rstat){
2899 for (i = 0; i < numrows; i++){
2900 switch (ws->getArtifStatus(i)){
2901 case CoinWarmStartBasis::basic:
2902 rstat[i] = SLACK_BASIC;
2903 break;
2904 case CoinWarmStartBasis::atLowerBound:
2905 rstat[i] = SLACK_AT_LB;
2906 break;
2907 case CoinWarmStartBasis::atUpperBound:
2908 rstat[i] = SLACK_AT_UB;
2909 break;
2910 case CoinWarmStartBasis::isFree: //can it happen?
2911 rstat[i] = SLACK_FREE;
2912 break;
2913 default:
2914 break; //can it happen?
2915 }
2916 }
2917 }
2918
2919 if (cstat){
2920 for (i = 0; i < numcols; i++){
2921 switch (ws->getStructStatus(i)){
2922 case CoinWarmStartBasis::basic:
2923 cstat[i] = VAR_BASIC;
2924 break;
2925 case CoinWarmStartBasis::atLowerBound:
2926 cstat[i] = VAR_AT_LB;
2927 break;
2928 case CoinWarmStartBasis::atUpperBound:
2929 cstat[i] = VAR_AT_UB;
2930 break;
2931 case CoinWarmStartBasis::isFree:
2932 cstat[i] = VAR_FREE;
2933 break;
2934 default:
2935 break; //can it happen?
2936 }
2937 }
2938 }
2939
2940 delete ws;
2941 }
2942
2943 /*===========================================================================*/
2944
2945 /*===========================================================================*\
2946 * Set an upper limit on the objective function value. Call the 'setobjulim'
2947 * CPLEX function.
2948 \*===========================================================================*/
2949
2950 void set_obj_upper_lim(LPdata *lp_data, double lim)
2951 {
2952
2953 #ifndef __OSI_CPLEX__
2954 OsiDblParam key = OsiDualObjectiveLimit;
2955
2956 retval = lp_data->si->setDblParam(key, lim);
2957 #else
2958
2959 CPXsetdblparam(lp_data->si->getEnvironmentPtr(), CPX_PARAM_OBJULIM, lim);
2960
2961 #endif
2962 }
2963
2964 /*===========================================================================*/
2965
2966 void set_timelim(LPdata *lp_data, double timelim)
2967 {
2968
2969 #ifdef __OSI_CLP__
2970
2971 ClpDblParam key = ClpMaxWallSeconds;
2972 lp_data->si->getModelPtr()->setDblParam(key, timelim);
2973
2974 #endif
2975
2976 }
2977
2978 /*===========================================================================*\
2979 * Set an upper limit on the number of iterations. If itlim < 0 then set
2980 * it to the maximum.
2981 \*===========================================================================*/
2982
2983 void set_itlim(LPdata *lp_data, int itlim)
2984 {
2985 if (itlim < 0) itlim = LP_MAX_ITER;
2986
2987 OsiIntParam key = OsiMaxNumIteration;
2988
2989 retval = lp_data->si->setIntParam(key, itlim);
2990 }
2991
2992 /*===========================================================================*/
2993
2994 void set_itlim_hotstart(LPdata *lp_data, int itlim)
2995 {
2996 if (itlim < 0) itlim = LP_MAX_ITER;
2997
2998 OsiIntParam key = OsiMaxNumIterationHotStart;
2999
3000 retval = lp_data->si->setIntParam(key, itlim);
3001 }
3002
3003 /*===========================================================================*/
3004
3005 void get_column(LPdata *lp_data, int j,
3006 double *colval, int *colind, int *collen, double *cj)
3007 {
3008 const CoinPackedMatrix *matrixByCol = lp_data->si->getMatrixByCol();
3009
3010 int i;
3011
3012 const double *matval = matrixByCol->getElements();
3013 const int *matind = matrixByCol->getIndices();
3014 const int *matbeg = matrixByCol->getVectorStarts();
3015 const int matbeg_j = matbeg[j];
3016
3017 *collen = matrixByCol->getVectorSize(j);
3018
3019 for (i = 0; i < (*collen); i++){
3020 colval[i] = matval[matbeg_j + i];
3021 colind[i] = matind[matbeg_j + i];
3022 }
3023
3024 const double * objval = lp_data->si->getObjCoefficients();
3025
3026 *cj = objval[j];
3027 }
3028
3029 /*===========================================================================*/
3030
3031 void get_row(LPdata *lp_data, int i,
3032 double *rowval, int *rowind, int *rowlen,
3033 double *rowub, double *rowlb)
3034 {
3035 const CoinPackedMatrix * matrixByRow = lp_data->si->getMatrixByRow();
3036
3037 int j;
3038
3039 const double *matval = matrixByRow->getElements();
3040 const int *matind = matrixByRow->getIndices();
3041 const int *matbeg = matrixByRow->getVectorStarts();
3042 const int matbeg_i = matbeg[i];
3043
3044 *rowlen = matrixByRow->getVectorSize(i);
3045 *rowub = lp_data->si->getRowUpper()[i];
3046 *rowlb = lp_data->si->getRowLower()[i];
3047
3048 for (j = 0; j < (*rowlen); j++){
3049 rowval[j] = matval[matbeg_i + j];
3050 rowind[j] = matind[matbeg_i + j];
3051 }
3052 }
3053
3054 /*===========================================================================*/
3055
3056 /* This routine returns the index of a row which proves the lp to be primal
3057 * infeasible. There must be one, or this function wouldn't be called. */
3058
3059 int get_proof_of_infeas(LPdata *lp_data, int *infind)
3060 {
3061 fprintf(stderr, "Function not implemented yet.");
3062 return(0);
3063 }
3064
3065 /*===========================================================================*/
3066
3067 /*===========================================================================*\
3068 * Get the solution (values of the structural variables in an optimal
3069 * solution) to the lp (specified by lp_data->lp) into the vector
3070 * lp_data->x. This can be done by calling the 'getx' CPLEX function.
3071 \*===========================================================================*/
3072
3073 void get_x(LPdata *lp_data)
3074 {
3075 memcpy(lp_data->x, lp_data->si->getColSolution(), lp_data->n * DSIZE);
3076 }
3077
3078 /*===========================================================================*/
3079
3080 void get_dj_pi(LPdata *lp_data)
3081 {
3082 const double * pi;
3083 const CoinPackedMatrix * matrix = lp_data->si->getMatrixByCol();
3084 const int * row = matrix->getIndices();
3085 const int * columnLength = matrix->getVectorLengths();
3086 const CoinBigIndex * columnStart = matrix->getVectorStarts();
3087 const double * elementByColumn = matrix->getElements();
3088 const double * objective = lp_data->si->getObjCoefficients();
3089 const double * lower = lp_data->si->getColLower();
3090 const double * upper = lp_data->si->getColUpper();
3091 double * dj = lp_data->dj;
3092 int numberColumns = lp_data->n;
3093 int t;
3094 memcpy(lp_data->dualsol, lp_data->si->getRowPrice(), lp_data->m * DSIZE);
3095 pi=lp_data->dualsol;
3096 memcpy(dj, lp_data->si->getReducedCost(), lp_data->n * DSIZE);
3097 /* djs may not be correct on fixed variables */
3098 /* fix assumes minimization */
3099 for (t=0; t < numberColumns; t++) {
3100 if (lower[t] == upper[t]) {
3101 int k;
3102 double value=objective[t];
3103 for (k=columnStart[t];k<columnStart[t]+columnLength[t];k++) {
3104 int iRow=row[k];
3105 value -= elementByColumn[k]*pi[iRow];
3106 }
3107 dj[t] = value;
3108 }
3109 }
3110 }
3111
3112 /*===========================================================================*/
3113 /*===========================================================================*/
3114
3115 void get_slacks(LPdata *lp_data)
3116 {
3117 int m = lp_data->m, i = 0;
3118 double * slacks = lp_data->slacks;
3119 row_data *rows = lp_data->rows;
3120 cut_data *cut;
3121
3122 #ifndef __OSI_CPLEX__
3123
3124 const double * rowActivity = lp_data->si->getRowActivity();
3125
3126 for (i = m - 1; i >= 0; i--) {
3127 cut = rows[i].cut;
3128 if ((cut->sense == 'R') && (cut->range < 0) ) {
3129 slacks[i] = - cut->rhs + rowActivity[i];
3130 } else {
3131 slacks[i] = cut->rhs - rowActivity[i];
3132 }
3133 }
3134
3135 #else
3136
3137 CPXgetslack(lp_data->si->getEnvironmentPtr(), lp_data->si->getLpPtr(),
3138 lp_data->slacks, 0, lp_data->m-1);
3139 /* Compute the real slacks for the free rows */
3140 for (i = m - 1; i >= 0; i--){
3141 if (rows[i].free){
3142 switch (rows[i].cut->sense){
3143 case 'E': slacks[i] += rows[i].cut->rhs - SYM_INFINITY; break;
3144 case 'L': slacks[i] += rows[i].cut->rhs - SYM_INFINITY; break;
3145 case 'G': slacks[i] += rows[i].cut->rhs + SYM_INFINITY; break;
3146 case 'R': slacks[i] += -rows[i].cut->rhs - SYM_INFINITY; break;
3147 }
3148 }
3149 }
3150
3151 #endif
3152 }
3153
3154 /*===========================================================================*/
3155
3156 void change_range(LPdata *lp_data, int rowind, double value)
3157 {
3158
3159 double rhs = lp_data->si->getRightHandSide()[rowind];
3160
3161 lp_data->si->setRowType(rowind,'R', rhs, value);
3162 }
3163
3164 /*===========================================================================*/
3165
3166 void change_rhs(LPdata *lp_data, int rownum, int *rhsind, double *rhsval)
3167 {
3168 char *sense = lp_data->tmp.c;
3169 double *range = lp_data->tmp.d;
3170 OsiXSolverInterface *si = lp_data->si;
3171 int i;
3172 const char *si_sense = si->getRowSense();
3173 const double *si_range = si->getRowRange();
3174
3175 for (i = 0; i < rownum; i++){
3176 sense[i] = si_sense[rhsind[i]];
3177 if (sense[i] == 'R'){
3178 range[i] = si_range[rhsind[i]];
3179 }
3180 }
3181
3182 si->setRowSetTypes(rhsind, rhsind + rownum, sense, rhsval, range);
3183 }
3184
3185 /*===========================================================================*/
3186
3187 void change_sense(LPdata *lp_data, int cnt, int *index, char *sense)
3188 {
3189 double *rhs = lp_data->tmp.d;
3190 double *range = (double *) calloc(cnt, DSIZE);
3191 OsiXSolverInterface *si = lp_data->si;
3192 const double *si_rhs = si->getRightHandSide();
3193 const double *si_range = si->getRowRange();
3194 int i;
3195
3196 for (i = 0; i < cnt; i++){
3197 rhs[i] = si_rhs[index[i]];
3198 if (sense[i] == 'R')
3199 range[i] = si_range[index[i]];
3200 }
3201
3202 si->setRowSetTypes(index, index + cnt, sense, rhs, range);
3203
3204 FREE(range);
3205 }
3206
3207 /*===========================================================================*/
3208
3209 void change_bounds(LPdata *lp_data, int cnt, int *index, char *lu, double *bd)
3210 {
3211 int i;
3212 OsiXSolverInterface *si = lp_data->si;
3213
3214 for (i = 0; i < cnt; i++){
3215 switch (lu[i]){
3216 case 'L':
3217 si->setColLower(index[i], bd[i]);
3218 break;
3219 case 'U':
3220 si->setColUpper(index[i], bd[i]);
3221 break;
3222 default:
3223 /* default: can't happen */
3224 break;
3225 }
3226 }
3227
3228 lp_data->lp_is_modified = LP_HAS_BEEN_MODIFIED;
3229 }
3230
3231 /*===========================================================================*/
3232
3233 void change_lbub(LPdata *lp_data, int j, double lb, double ub)
3234 {
3235 lp_data->si->setColBounds(j,lb,ub);
3236 lp_data->lp_is_modified = LP_HAS_BEEN_MODIFIED;
3237 }
3238
3239 /*===========================================================================*/
3240
3241 void change_ub(LPdata *lp_data, int j, double ub)
3242 {
3243 lp_data->si->setColUpper(j,ub);
3244 lp_data->lp_is_modified = LP_HAS_BEEN_MODIFIED;
3245 }
3246
3247 /*===========================================================================*/
3248
3249 void change_lb(LPdata *lp_data, int j, double lb)
3250 {
3251 lp_data->si->setColLower(j,lb);
3252 lp_data->lp_is_modified = LP_HAS_BEEN_MODIFIED;
3253 }
3254
3255 /*===========================================================================*/
3256
3257 void get_ub(LPdata *lp_data, int j, double *ub)
3258 {
3259 *ub=lp_data->si->getColUpper()[j];
3260 }
3261
3262 /*===========================================================================*/
3263
3264 void get_lb(LPdata *lp_data, int j, double *lb)
3265 {
3266 *lb=lp_data->si->getColLower()[j];
3267 }
3268
3269 /*===========================================================================*/
3270
3271 void get_bounds(LPdata *lp_data)
3272 {
3273 lp_data->lb = const_cast<double *>(lp_data->si->getColLower());
3274 lp_data->ub = const_cast<double *>(lp_data->si->getColUpper());
3275 }
3276
3277 /*===========================================================================*/
3278
3279 void get_objcoef(LPdata *lp_data, int j, double *objcoef)
3280 {
3281 *objcoef = lp_data->si->getObjCoefficients()[j];
3282 }
3283
3284 /*===========================================================================*/
3285 void get_objcoeffs(LPdata *lp_data)
3286 {
3287 const double *si_objcoeffs = lp_data->si->getObjCoefficients();
3288 memcpy (lp_data->mip->obj,si_objcoeffs,lp_data->n*DSIZE);
3289 }
3290
3291 /*===========================================================================*/
3292 void change_objcoeff(LPdata *lp_data, const int* indexFirst,
3293 const int* indexLast, double *coeffs)
3294 {
3295 lp_data->si->setObjCoeffSet(indexFirst, indexLast, coeffs);
3296 }
3297
3298 /*===========================================================================*/
3299 void get_rhs_rng_sense(LPdata *lp_data)
3300 {
3301 const double *rowub = lp_data->si->getRowUpper();
3302 const double *rowlb = lp_data->si->getRowLower();
3303 double *mip_rhs = lp_data->mip->rhs;
3304 double *mip_rngval = lp_data->mip->rngval;
3305 char *mip_sense = lp_data->mip->sense;
3306
3307 for (int i=0;i<lp_data->m;i++) {
3308 if (rowub[i]>=SYM_INFINITY) {
3309 mip_sense[i] = 'G';
3310 mip_rhs[i] = rowlb[i];
3311 }
3312 else if (rowlb[i]<=-SYM_INFINITY) {
3313 mip_sense[i] = 'L';
3314 mip_rhs[i] = rowub[i];
3315 }
3316 else {
3317 mip_sense[i] = 'R';
3318 mip_rhs[i] = rowub[i];
3319 mip_rngval[i] = rowub[i]-rowlb[i];
3320 }
3321 }
3322 }
3323
3324 /*===========================================================================*/
3325 /*
3326 * copy everything from lp_data into new_data. it is assumed that new_data is
3327 * already mallocced.
3328 */
3329 int copy_lp_data(LPdata *lp_data, LPdata *new_data)
3330 {
3331 int termcode = FUNCTION_TERMINATED_NORMALLY;
3332 int n = lp_data->n;
3333 int m = lp_data->m;
3334 //double *lb, *ub;
3335 OsiXSolverInterface *si = lp_data->si;
3336
3337 if (!new_data) {
3338 return FUNCTION_TERMINATED_ABNORMALLY;
3339 }
3340
3341 new_data->lpetol = lp_data->lpetol;
3342 new_data->n = n;
3343 new_data->m = m;
3344 new_data->nz = lp_data->nz;
3345 new_data->maxn = lp_data->maxn;
3346 new_data->maxm = lp_data->maxm;
3347 new_data->maxnz = lp_data->maxnz;
3348
3349 //lb = (double *)malloc(n*DSIZE);
3350 //ub = (double *)malloc(n*DSIZE);
3351
3352 open_lp_solver(new_data);
3353 /* Turn off the OSI messages (There are LOTS of them) */
3354 new_data->si->setHintParam(OsiDoReducePrint);
3355 new_data->si->messageHandler()->setLogLevel(0);
3356
3357 new_data->si->loadProblem(*(si->getMatrixByRow()),
3358 si->getColLower(),
3359 si->getColUpper(),
3360 si->getObjCoefficients(),
3361 si->getRowLower(),
3362 si->getRowUpper()
3363 );
3364 /* get_bounds just returns a const pointer to si->ub, si->lb. we need to
3365 * memcpy because these pointers get changed when addCols is used */
3366 /* menal - I don't see where these pointers get changed, so disabling for now */
3367 get_bounds(new_data);
3368 //memcpy(lb,new_data->lb,DSIZE*n);
3369 //memcpy(ub,new_data->ub,DSIZE*n);
3370
3371 //new_data->lb = lb;
3372 //new_data->ub = ub;
3373
3374
3375 return termcode;
3376 }
3377
3378 /*===========================================================================*/
3379
3380 void delete_rows(LPdata *lp_data, int deletable, int *free_rows)
3381 {
3382
3383 int i, m = lp_data->m;
3384 int *which = lp_data->tmp.i1 + lp_data->m;
3385 int delnum = 0;
3386
3387 CoinFillN(which, deletable, 0);
3388
3389 for (i = 0; i < m; i++){
3390 if (free_rows[i]){
3391 which[delnum++] = i;
3392 }
3393 }
3394
3395 lp_data->si->deleteRows(delnum, which);
3396 lp_data->nz = lp_data->si->getNumElements();
3397 lp_data->m -= delnum;
3398 }
3399
3400 /*===========================================================================*/
3401
3402 void delete_rows_with_ind(LPdata *lp_data, int deletable, int *rowind)
3403 {
3404
3405 lp_data->si->deleteRows(deletable, rowind);
3406 lp_data->nz = lp_data->si->getNumElements();
3407 lp_data->m -= deletable;
3408 }
3409
3410 /*===========================================================================*/
3411
3412 int delete_cols(LPdata *lp_data, int delnum, int *delstat)
3413 {
3414 int i, n = lp_data->n;
3415 int *which = (int *) calloc(delnum, ISIZE);
3416 int num_to_delete = 0, num_to_keep = 0;
3417 double *dj = lp_data->dj;
3418 double *x = lp_data->x;
3419 char *status = lp_data->status;
3420
3421 for (i = n - 1; i >= 0; i--){
3422 if (delstat[i]){
3423 which[num_to_delete++]=i;
3424 }
3425 }
3426
3427 lp_data->si->deleteCols(num_to_delete, which);
3428 lp_data->nz = lp_data->si->getNumElements();
3429 FREE(which);
3430
3431 /* make result as CPLEX does */
3432 for (i = 0, num_to_keep = 0; i < lp_data->n; i++){
3433 if (delstat[i]){
3434 delstat[i] = -1;
3435 }else{
3436 delstat[i] = num_to_keep++;
3437 dj[delstat[i]] = dj[i];
3438 x[delstat[i]] = x[i];
3439 status[delstat[i]] = status[i];
3440 }
3441 }
3442
3443 lp_data->n = num_to_keep;
3444
3445 return(num_to_delete);
3446 }
3447
3448 /*===========================================================================*/
3449
3450 void release_var(LPdata *lp_data, int j, int where_to_move)
3451 {
3452 fprintf(stderr, "Function not implemented yet.");
3453 exit(-1);
3454 }
3455
3456 /*===========================================================================*/
3457
3458 void free_row_set(LPdata *lp_data, int length, int *index)
3459 {
3460 char *sense = lp_data->tmp.c; /* m (now) */
3461 double *rhs = lp_data->tmp.d; /* m */
3462 double *range = (double *) calloc(length, DSIZE);
3463 int i;
3464 OsiXSolverInterface *si = lp_data->si;
3465 const double infinity = si->getInfinity();
3466 const double *si_rhs = si->getRightHandSide();
3467 const double *si_rowrange = si->getRowRange();
3468 const char *si_rowsense = si->getRowSense();
3469
3470 for (i = 0; i < length; i++){
3471 rhs[i] = si_rhs[index[i]];
3472 sense[i] = si_rowsense[index[i]];
3473 if (sense[i] =='R'){
3474 range[i] = si_rowrange[index[i]];
3475 }
3476 }
3477
3478 for (i = 0; i < length; i++) {
3479 // range[i]=0;
3480 switch (sense[i]){
3481 case 'E':
3482 rhs[i] = infinity;
3483 sense[i] = 'L';
3484 break;
3485 case 'L':
3486 rhs[i] = infinity;
3487 break;
3488 case 'R':
3489 range[i] = 2*infinity;
3490 break;
3491 case 'G':
3492 rhs[i] = -infinity;
3493 }
3494 }
3495
3496 si->setRowSetTypes(index, index + length, sense, rhs, range);
3497
3498 FREE(range);
3499 }
3500
3501 /*===========================================================================*/
3502
3503 void constrain_row_set(LPdata *lp_data, int length, int *index)
3504 {
3505 char *sense = lp_data->tmp.c; /* m (now) */
3506 double *rhs = lp_data->tmp.d; /* m */
3507 double *range = (double *) calloc(length, DSIZE);
3508 row_data *rows = lp_data->rows;
3509 cut_data *cut;
3510 int i;
3511
3512 for (i = length - 1; i >= 0; i--){
3513 cut = rows[index[i]].cut;
3514 rhs[i] = cut->rhs;
3515 if ((sense[i] = cut->sense) == 'R'){
3516 range[i] = cut->range;
3517 }
3518 }
3519
3520 lp_data->si->setRowSetTypes(index, index + length, sense, rhs, range);
3521
3522 FREE(range);
3523
3524 }
3525
3526 /*===========================================================================*/
3527
3528 int read_mps(MIPdesc *mip, char *infile, char *probname, int versbotiy)
3529 {
3530 int j, errors;
3531 CoinMpsIO mps;
3532
3533 mps.messageHandler()->setLogLevel(0);
3534
3535 #if 0
3536
3537 int j, last_dot = 0, last_dir = 0;
3538 char fname[80] = "", ext[10] = "";
3539
3540 size_t size = 1000;
3541 char* buf = 0;
3542
3543 while (true) {
3544 buf = (char*)malloc(CSIZE*size);
3545 if (getcwd(buf, size))
3546 break;
3547 FREE(buf);
3548 buf = 0;
3549 size = 2*size;
3550 }
3551 char slash = buf[0] == '/' ? '/' : '\\';
3552 FREE(buf);
3553
3554 for (j = 0;; j++){
3555 if (infile[j] == '\0')
3556 break;
3557 if (infile[j] == '.') {
3558 last_dot = j;
3559 }
3560 if(infile[j] == slash){
3561 last_dir = j;
3562 }
3563 }
3564
3565 if(last_dir < last_dot){
3566 memcpy(fname, infile, CSIZE*last_dot);
3567 memcpy(ext, infile + last_dot + 1, CSIZE*(j - last_dot - 1));
3568 }
3569 else{
3570 memcpy(fname, infile, CSIZE*j);
3571 }
3572 #endif
3573
3574 mps.setInfinity(mps.getInfinity());
3575
3576 if ((errors = mps.readMps(infile,""))){
3577 return(errors);
3578 }
3579
3580 strncpy(probname, mps.getProblemName(), 80);
3581
3582 mip->m = mps.getNumRows();
3583 mip->n = mps.getNumCols();
3584 mip->nz = mps.getNumElements();
3585
3586 const CoinPackedMatrix * matrixByCol= mps.getMatrixByCol();
3587
3588 if (mip->n){
3589 mip->obj = (double *) malloc(DSIZE * mip->n);
3590 mip->obj1 = NULL;
3591 mip->obj2 = NULL;
3592 mip->ub = (double *) malloc(DSIZE * mip->n);
3593 mip->lb = (double *) malloc(DSIZE * mip->n);
3594 mip->is_int = (char *) calloc(CSIZE, mip->n);
3595 memcpy(mip->obj, mps.getObjCoefficients(), DSIZE * mip->n);
3596 memcpy(mip->ub, mps.getColUpper(), DSIZE * mip->n);
3597 memcpy(mip->lb, mps.getColLower(), DSIZE * mip->n);
3598
3599 mip->matbeg = (int *) malloc(ISIZE * (mip->n + 1));
3600 memcpy(mip->matbeg, matrixByCol->getVectorStarts(), ISIZE * (mip->n + 1));
3601
3602 mip->colname = (char **) malloc(sizeof(char *) * mip->n);
3603 }
3604 if (mip->m){
3605 mip->rhs = (double *) malloc(DSIZE * mip->m);
3606 mip->sense = (char *) malloc(CSIZE * mip->m);
3607 mip->rngval = (double *) malloc(DSIZE * mip->m);
3608 memcpy(mip->rhs, mps.getRightHandSide(), DSIZE * mip->m);
3609 memcpy(mip->sense, mps.getRowSense(), CSIZE * mip->m);
3610 memcpy(mip->rngval, mps.getRowRange(), DSIZE * mip->m);
3611 }
3612
3613 //user defined matind, matval, matbeg--fill as column ordered
3614
3615 if (mip->nz){
3616 mip->matval = (double *) malloc(DSIZE*mip->matbeg[mip->n]);
3617 mip->matind = (int *) malloc(ISIZE*mip->matbeg[mip->n]);
3618
3619 memcpy(mip->matval, matrixByCol->getElements(), DSIZE * mip->matbeg[mip->n]);
3620 memcpy(mip->matind, matrixByCol->getIndices(), ISIZE * mip->matbeg[mip->n]);
3621 }
3622
3623 for (j = 0; j < mip->n; j++){
3624 mip->is_int[j] = mps.isInteger(j);
3625 mip->colname[j] = (char *) malloc(CSIZE * MAX_NAME_SIZE);
3626 strncpy(mip->colname[j], mps.columnName(j), MAX_NAME_SIZE);
3627 mip->colname[j][MAX_NAME_SIZE-1] = 0;
3628 }
3629
3630 if (mip->obj_sense == SYM_MAXIMIZE){
3631 for (j = 0; j < mip->n; j++){
3632 mip->obj[j] *= -1.0;
3633 }
3634 }
3635
3636 mip->obj_offset = -mps.objectiveOffset();
3637
3638 return(errors);
3639 }
3640
3641 /*===========================================================================*/
3642
3643 int read_lp(MIPdesc *mip, char *infile, char *probname, int verbosity)
3644 {
3645
3646 int j;
3647 CoinLpIO lp;
3648
3649 lp.readLp(infile);
3650
3651 strncpy(probname, lp.getProblemName(), 80);
3652
3653 mip->m = lp.getNumRows();
3654 mip->n = lp.getNumCols();
3655 mip->nz = lp.getNumElements();
3656
3657 mip->obj = (double *) malloc(DSIZE * mip->n);
3658 mip->obj1 = NULL;
3659 mip->obj2 = NULL;
3660 mip->rhs = (double *) malloc(DSIZE * mip->m);
3661 mip->sense = (char *) malloc(CSIZE * mip->m);
3662 mip->rngval = (double *) malloc(DSIZE * mip->m);
3663 mip->ub = (double *) malloc(DSIZE * mip->n);
3664 mip->lb = (double *) malloc(DSIZE * mip->n);
3665 mip->is_int = (char *) calloc(CSIZE, mip->n);
3666
3667 if (lp.getNumObjectives() >= 2){
3668 mip->obj1 = (double *) calloc(mip->n, DSIZE);
3669 mip->obj2 = (double *) calloc(mip->n, DSIZE);
3670 memcpy(mip->obj, lp.getObjCoefficients(0), DSIZE * mip->n);
3671 memcpy(mip->obj1, lp.getObjCoefficients(0), DSIZE * mip->n);
3672 memcpy(mip->obj2, lp.getObjCoefficients(1), DSIZE * mip->n);
3673 if (lp.getNumObjectives() > 2){
3674 PRINT(verbosity, 2, ("Ignoring extra objectives...\n\n"));
3675 }
3676 }
3677 else{
3678 memcpy(mip->obj, lp.getObjCoefficients(), DSIZE * mip->n);
3679 }
3680
3681 memcpy(mip->rhs, lp.getRightHandSide(), DSIZE * mip->m);
3682 memcpy(mip->sense, lp.getRowSense(), CSIZE * mip->m);
3683 memcpy(mip->rngval, lp.getRowRange(), DSIZE * mip->m);
3684 memcpy(mip->ub, lp.getColUpper(), DSIZE * mip->n);
3685 memcpy(mip->lb, lp.getColLower(), DSIZE * mip->n);
3686
3687 //user defined matind, matval, matbeg--fill as column ordered
3688
3689 const CoinPackedMatrix * matrixByCol= lp.getMatrixByCol();
3690
3691 mip->matbeg = (int *) malloc(ISIZE * (mip->n + 1));
3692 memcpy(mip->matbeg, matrixByCol->getVectorStarts(), ISIZE * (mip->n + 1));
3693
3694 mip->matval = (double *) malloc(DSIZE*mip->matbeg[mip->n]);
3695 mip->matind = (int *) malloc(ISIZE*mip->matbeg[mip->n]);
3696
3697 memcpy(mip->matval, matrixByCol->getElements(), DSIZE * mip->matbeg[mip->n]);
3698 memcpy(mip->matind, matrixByCol->getIndices(), ISIZE * mip->matbeg[mip->n]);
3699
3700 mip->colname = (char **) malloc(sizeof(char *) * mip->n);
3701
3702 for (j = 0; j < mip->n; j++){
3703 mip->is_int[j] = lp.isInteger(j);
3704 mip->colname[j] = (char *) malloc(CSIZE * MAX_NAME_SIZE);
3705 strncpy(mip->colname[j], lp.columnName(j), MAX_NAME_SIZE);
3706 mip->colname[j][MAX_NAME_SIZE-1] = 0;
3707 }
3708
3709 if (mip->obj_sense == SYM_MAXIMIZE){
3710 for (j = 0; j < mip->n; j++){
3711 mip->obj[j] *= -1.0;
3712 }
3713 }
3714
3715 mip->obj_offset = -lp.objectiveOffset();
3716
3717 return 0;
3718 }
3719
3720 /*===========================================================================*/
3721
3722 void write_mps(LPdata *lp_data, char *fname)
3723 {
3724 const char * extension = "MPS";
3725 OsiXSolverInterface *si = lp_data->si;
3726 double ObjSense = si->getObjSense();
3727 int i;
3728
3729 for (i = 0; i < lp_data->n; i++) {
3730 si->setContinuous(i);
3731 }
3732
3733 si->writeMps(fname, extension, ObjSense);
3734 }
3735
3736 /*===========================================================================*/
3737
3738 void write_lp(LPdata *lp_data, char *fname)
3739 {
3740 const char * extension = "LP";
3741 OsiXSolverInterface *si = lp_data->si;
3742 double ObjSense = si->getObjSense();
3743 int i;
3744
3745 for (i = 0; i < lp_data->n; i++) {
3746 si->setContinuous(i);
3747 }
3748
3749 si->writeLp(fname, extension);
3750 }
3751
3752 /*===========================================================================*/
3753
3754 void write_mip_desc_mps(MIPdesc *mip, char *fname)
3755 {
3756 int i;
3757 double * obj;
3758 char filename[80] = "";
3759 CoinMpsIO mps;
3760 CoinPackedMatrix mip_matrix(true, mip->m, mip->n, mip->nz, mip->matval,
3761 mip->matind, mip->matbeg, 0);
3762
3763 obj = (double *) malloc(DSIZE*mip->n);
3764 memcpy(obj, mip->obj, DSIZE*mip->n);
3765 if (mip->obj_sense == SYM_MAXIMIZE){
3766 for (i = 0; i < mip->n; i++){
3767 obj[i] *= -1.0;
3768 }
3769 }
3770
3771 mps.setMpsData(mip_matrix, mps.getInfinity(), mip->lb, mip->ub, obj,
3772 mip->is_int, mip->sense, mip->rhs, mip->rngval,
3773 mip->colname, NULL);
3774 mps.setObjectiveOffset(mip->obj_offset);
3775
3776 sprintf(filename, "%s%s%s", fname, ".","MPS");
3777 mps.writeMps(filename);
3778 FREE(obj);
3779 }
3780
3781 /*===========================================================================*/
3782
3783 void write_mip_desc_lp(MIPdesc *mip, char *fname)
3784 {
3785 int i;
3786 double * obj, * rlb, * rub, infinity;
3787 char filename[80] = "";
3788 CoinLpIO lp;
3789 CoinPackedMatrix mip_matrix(true, mip->m, mip->n, mip->nz, mip->matval,
3790 mip->matind, mip->matbeg, 0);
3791
3792 obj = (double *) malloc(DSIZE*mip->n);
3793 memcpy(obj, mip->obj, DSIZE*mip->n);
3794 if (mip->obj_sense == SYM_MAXIMIZE){
3795 for (i = 0; i < mip->n; i++){
3796 obj[i] *= -1.0;
3797 }
3798 }
3799
3800 rlb = (double *) malloc(DSIZE*mip->m);
3801 rub = (double *) malloc(DSIZE*mip->m);
3802 infinity = lp.getInfinity();
3803
3804 /* convert sense to bound */
3805 for(i = 0; i < mip->m; i++){
3806 switch (mip->sense[i]){
3807 case 'E':
3808 rlb[i] = rub[i] = mip->rhs[i];
3809 break;
3810 case 'L':
3811 rlb[i] = -infinity;
3812 rub[i] = mip->rhs[i];
3813 break;
3814 case 'G':
3815 rlb[i] = mip->rhs[i];
3816 rub[i] = infinity;
3817 break;
3818 case 'R':
3819 rlb[i] = mip->rhs[i] - mip->rngval[i];
3820 rub[i] = mip->rhs[i];
3821 break;
3822 case 'N':
3823 rlb[i] = -infinity;
3824 rub[i] = infinity;
3825 break;
3826 }
3827 }
3828
3829 lp.setLpDataWithoutRowAndColNames(mip_matrix, mip->lb, mip->ub, obj,
3830 mip->is_int, rlb, rub);
3831 lp.setObjectiveOffset(mip->obj_offset);
3832 lp.setLpDataRowAndColNames(NULL, mip->colname);
3833 sprintf(filename, "%s%s%s", fname, ".","LPT");
3834 lp.writeLp(filename);
3835
3836 FREE(obj);
3837 FREE(rlb);
3838 FREE(rub);
3839 }
3840
3841
3842 /*===========================================================================*/
3843
3844 void write_sav(LPdata *lp_data, char *fname)
3845 {
3846 fprintf(stderr, "Function not implemented yet.");
3847 exit(-1);
3848 }
3849
3850 /*===========================================================================*/
3851
3852 #ifdef USE_CGL_CUTS
3853
3854 #include "sym_qsort.h"
3855
3856 void generate_cgl_cuts(LPdata *lp_data, int *num_cuts, cut_data ***cuts,
3857 char send_to_pool, int bc_index, int bc_level,
3858 int node_iter_num, int max_cuts_before_resolve,
3859 double ub, int *bnd_changes,
3860 lp_stat_desc *lp_stat, node_times *comp_times,
3861 int verbosity)
3862 {
3863 #if 0
3864 OsiCuts cutlist;
3865 OsiRowCut cut, cut2;
3866 int n = lp_data->n;
3867 int i = 0, j = 0, k = 0, l = 0;
3868 int *matind;
3869 double *matval, total_time=0, cut_time=0;
3870 cgl_params *par = &(lp_data->cgl);
3871 int termcode, iterd, cut_num = 0;
3872 int new_cut_num = 0;
3873 int is_top_iter = (lp_data->lp_count == 1) ? TRUE : FALSE;
3874 OsiXSolverInterface *si = lp_data->si;
3875 var_desc **vars = lp_data->vars;
3876 int is_rootnode = (bc_index>0) ? FALSE : TRUE;
3877 //double *newLower = lp_data->tmp.d;
3878 //double *newUpper = lp_data->tmp.d+n;
3879 int sizeColCuts, should_generate;
3880 int num_duplicate_cuts = 0;
3881 const double lpetol = lp_data->lpetol;
3882 const double etol1000 = lpetol*1000;
3883
3884 #ifndef COMPILE_IN_LP
3885 par->probing_generated_in_root = TRUE;
3886 par->gomory_generated_in_root = TRUE;
3887 par->redsplit_generated_in_root = FALSE;
3888 par->oddhole_generated_in_root = TRUE;
3889 par->mir_generated_in_root = TRUE;
3890 par->twomir_generated_in_root = FALSE;
3891 par->clique_generated_in_root = FALSE;
3892 par->flowcover_generated_in_root = TRUE;
3893 par->rounding_generated_in_root = FALSE;
3894 par->lift_and_project_generated_in_root = FALSE;
3895 par->landp_generated_in_root = FALSE;
3896 #endif
3897
3898 /* Set proper variables to be integer */
3899 /*
3900 * TODO: take this loop outside, should not be called in every call of
3901 * generate_cgl_cuts
3902 */
3903 if (node_iter_num < 2) {
3904 for (i = 0; i < n; i++) {
3905 if (vars[i]->is_int) { // integer or binary
3906 si->setInteger(i);
3907 }
3908 }
3909 }
3910 /* TODO: move these to vars[i]->... */
3911 //get_bounds(lp_data);
3912 //memcpy(newLower,lp_data->lb,DSIZE*n);
3913 //memcpy(newUpper,lp_data->ub,DSIZE*n);
3914
3915
3916 /* twice is necessary */
3917 cut_time = used_time(&total_time);
3918 cut_time = used_time(&total_time);
3919
3920 /* create CGL probing cuts */
3921 should_generate_this_cgl_cut(cut_num, max_cuts_before_resolve,
3922 par->generate_cgl_probing_cuts,
3923 par->generate_cgl_probing_cuts_freq, bc_level, bc_index,
3924 lp_stat->probing_cuts_root, &should_generate);
3925 if (should_generate==TRUE) {
3926 CglProbing *probe = new CglProbing;
3927 probe->setRowCuts(0);
3928 probe->setMode(2);
3929 if (ub < SYM_INFINITY/10) {
3930 probe->setUsingObjective(1);
3931 } else {
3932 probe->setUsingObjective(0);
3933 }
3934 if ((bc_level<6 && comp_times->probing_cuts>comp_times->lp) ||
3935 (bc_level>6 && comp_times->probing_cuts>comp_times->lp/10)) {
3936 /* since we are not using cgltreeinfo,
3937 * all nodes are like root nodes
3938 */
3939 probe->setMaxLookRoot(5);
3940 } else if (bc_level < 1) {
3941 probe->setMaxPass(10); /* default is 3 */
3942 probe->setMaxPassRoot(10); /* default is 3 */
3943 probe->setMaxElements(10000); /* default is 1000, 10000 for root */
3944 probe->setMaxElementsRoot(10000); /* default is 1000, 10000 for root */
3945 probe->setMaxLook(100); /* default is 50 */
3946 probe->setMaxLookRoot(100); /* default is 50 */
3947 probe->setMaxProbe(200); /* default is 100 */
3948 probe->setMaxProbeRoot(200); /* default is 100 */
3949 }
3950 probe->generateCuts(*(si), cutlist);
3951 if ((new_cut_num = cutlist.sizeCuts() - cut_num) > 0) {
3952 if (is_top_iter){
3953 par->probing_generated_in_root = TRUE;
3954 }
3955 PRINT(verbosity, 4,
3956 ("%i probing cuts added\n", new_cut_num));
3957 lp_stat->cuts_generated += new_cut_num;
3958 lp_stat->probing_cuts += new_cut_num;
3959 if (is_rootnode) {
3960 lp_stat->cuts_root += new_cut_num;
3961 lp_stat->probing_cuts_root += new_cut_num;
3962 }
3963 }
3964 cut_num = cutlist.sizeCuts();
3965 delete probe;
3966 cut_time = used_time(&total_time);
3967 comp_times->cuts += cut_time;
3968 comp_times->probing_cuts += cut_time;
3969 //cutlist.printCuts();
3970 }
3971 check_cuts(cutlist, lp_data, bc_level, num_cuts, cuts, send_to_pool,
3972 bnd_changes, lp_stat, comp_times, verbosity);
3973
3974 /* create CGL knapsack cuts */
3975 should_generate_this_cgl_cut(cut_num, max_cuts_before_resolve,
3976 par->generate_cgl_knapsack_cuts,
3977 par->generate_cgl_knapsack_cuts_freq, bc_level, bc_index,
3978 lp_stat->knapsack_cuts_root, &should_generate);
3979 if (should_generate==TRUE) {
3980 CglKnapsackCover *knapsack = new CglKnapsackCover;
3981 if (bc_level<6) {
3982 knapsack->setMaxInKnapsack(1000); // default is 50
3983 knapsack->switchOffExpensive(); // seems to get into infinite loop if
3984 // turned on
3985 }
3986 knapsack->generateCuts(*si, cutlist);
3987 if ((new_cut_num = cutlist.sizeCuts() - cut_num) > 0) {
3988 if (is_top_iter){
3989 par->knapsack_generated_in_root = TRUE;
3990 }
3991 PRINT(verbosity, 4,
3992 ("%i knapsack cuts added\n", new_cut_num));
3993 lp_stat->cuts_generated += new_cut_num;
3994 lp_stat->knapsack_cuts += new_cut_num;
3995 if (is_rootnode) {
3996 lp_stat->cuts_root += new_cut_num;
3997 lp_stat->knapsack_cuts_root += new_cut_num;
3998 }
3999 }
4000 cut_num = cutlist.sizeCuts();
4001 delete knapsack;
4002 cut_time = used_time(&total_time);
4003 comp_times->cuts += cut_time;
4004 comp_times->knapsack_cuts += cut_time;
4005 }
4006
4007 /* create CGL clique cuts */
4008 should_generate_this_cgl_cut(cut_num, max_cuts_before_resolve,
4009 par->generate_cgl_clique_cuts, par->generate_cgl_clique_cuts_freq,
4010 bc_level, bc_index, lp_stat->clique_cuts_root, &should_generate);
4011 if (should_generate==TRUE) {
4012 CglClique *clique = new CglClique;
4013 if (bc_level<6) {
4014 clique->setStarCliqueCandidateLengthThreshold(100); // default 12
4015 clique->setRowCliqueCandidateLengthThreshold(100); // default 12
4016 }
4017
4018 clique->setStarCliqueReport(FALSE);
4019 clique->setRowCliqueReport(FALSE);
4020 clique->generateCuts(*(lp_data->si), cutlist);
4021 if ((new_cut_num = cutlist.sizeCuts() - cut_num) > 0) {
4022 if (is_top_iter){
4023 par->clique_generated_in_root = TRUE;
4024 }
4025 PRINT(verbosity, 4,
4026 ("%i clique cuts added\n", new_cut_num));
4027 lp_stat->cuts_generated += new_cut_num;
4028 lp_stat->clique_cuts += new_cut_num;
4029 if (is_rootnode) {
4030 lp_stat->cuts_root += new_cut_num;
4031 lp_stat->clique_cuts_root += new_cut_num;
4032 }
4033 }
4034 cut_num = cutlist.sizeCuts();
4035 delete clique;
4036 cut_time = used_time(&total_time);
4037 comp_times->cuts += cut_time;
4038 comp_times->clique_cuts += cut_time;
4039 }
4040
4041 /* create CGL gomory cuts */
4042 should_generate_this_cgl_cut(cut_num, max_cuts_before_resolve,
4043 par->generate_cgl_gomory_cuts,
4044 par->generate_cgl_gomory_cuts_freq, bc_level, bc_index,
4045 lp_stat->gomory_cuts_root, &should_generate);
4046 if (should_generate==TRUE) {
4047 if ((bc_level<6 && comp_times->gomory_cuts<comp_times->lp) ||
4048 (bc_level>6 && comp_times->gomory_cuts<comp_times->lp/10)) {
4049 CglGomory *gomory = new CglGomory;
4050 // TODO: change this to something based on number of cols, sparsity
4051 // etc.
4052 if (bc_level<1) {
4053 gomory->setLimitAtRoot(100);
4054 gomory->setLimit(100);
4055 } else {
4056 gomory->setLimitAtRoot(100);
4057 gomory->setLimit(100);
4058 }
4059 gomory->generateCuts(*(lp_data->si), cutlist);
4060 if ((new_cut_num = cutlist.sizeCuts() - cut_num) > 0) {
4061 if (is_top_iter){
4062 par->gomory_generated_in_root = TRUE;
4063 }
4064 PRINT(verbosity, 4,
4065 ("%i Gomory cuts added\n", new_cut_num));
4066 lp_stat->cuts_generated += new_cut_num;
4067 lp_stat->gomory_cuts += new_cut_num;
4068 if (is_rootnode) {
4069 lp_stat->cuts_root += new_cut_num;
4070 lp_stat->gomory_cuts_root += new_cut_num;
4071 }
4072 }
4073 cut_num = cutlist.sizeCuts();
4074 delete gomory;
4075 cut_time = used_time(&total_time);
4076 comp_times->cuts += cut_time;
4077 comp_times->gomory_cuts += cut_time;
4078 }
4079 }
4080 //printf("gomory time = %f, lp time = %f\n",comp_times->gomory_cuts, comp_times->lp);
4081
4082 /* create CGL flow cover cuts */
4083 should_generate_this_cgl_cut(cut_num, max_cuts_before_resolve,
4084 par->generate_cgl_flowcover_cuts,
4085 par->generate_cgl_flowcover_cuts_freq, bc_level, bc_index,
4086 lp_stat->flowcover_cuts_root, &should_generate);
4087 if (should_generate==TRUE) {
4088 CglFlowCover *flow = new CglFlowCover;
4089 /* numFlowCuts_ is a static variable! needs to be reset */
4090 flow->setNumFlowCuts(0);
4091 flow->generateCuts(*(lp_data->si), cutlist);
4092 if ((new_cut_num = cutlist.sizeCuts() - cut_num) > 0) {
4093 if (is_top_iter){
4094 par->flowcover_generated_in_root = TRUE;
4095 }
4096 PRINT(verbosity, 4,
4097 ("%i flow cover cuts added\n", new_cut_num));
4098 lp_stat->cuts_generated += new_cut_num;
4099 lp_stat->flowcover_cuts += new_cut_num;
4100 if (is_rootnode) {
4101 lp_stat->cuts_root += new_cut_num;
4102 lp_stat->flowcover_cuts_root += new_cut_num;
4103 }
4104 }
4105 cut_num = cutlist.sizeCuts();
4106 delete flow;
4107 cut_time = used_time(&total_time);
4108 comp_times->cuts += cut_time;
4109 comp_times->flowcover_cuts += cut_time;
4110 }
4111
4112 /* create CGL twomir cuts */
4113 should_generate_this_cgl_cut(cut_num, max_cuts_before_resolve,
4114 par->generate_cgl_twomir_cuts,
4115 par->generate_cgl_twomir_cuts_freq,
4116 bc_level, bc_index,
4117 lp_stat->twomir_cuts_root, &should_generate);
4118 if (should_generate==TRUE) {
4119 CglTwomir *twomir = new CglTwomir;
4120 twomir->setMaxElements (100);
4121 twomir->setCutTypes (TRUE, TRUE, TRUE, TRUE);
4122 twomir->generateCuts(*(lp_data->si), cutlist);
4123 if ((new_cut_num = cutlist.sizeCuts() - cut_num) > 0) {
4124 if (is_top_iter){
4125 par->twomir_generated_in_root = TRUE;
4126 }
4127 PRINT(verbosity, 4,
4128 ("%i 2-MIR cuts added\n", new_cut_num));
4129 lp_stat->cuts_generated += new_cut_num;
4130 lp_stat->twomir_cuts += new_cut_num;
4131 if (is_rootnode) {
4132 lp_stat->cuts_root += new_cut_num;
4133 lp_stat->twomir_cuts_root += new_cut_num;
4134 }
4135 }
4136 cut_num = cutlist.sizeCuts();
4137 delete twomir;
4138 cut_time = used_time(&total_time);
4139 comp_times->cuts += cut_time;
4140 comp_times->twomir_cuts += cut_time;
4141 }
4142 }
4143
4144 /* create CGL redsplit cuts */
4145 if(par->generate_cgl_redsplit_cuts > -1 &&
4146 par->generate_cgl_redsplit_cuts_freq > 0){
4147 if(par->generate_cgl_redsplit_cuts == GENERATE_ALWAYS ||
4148 (par->generate_cgl_redsplit_cuts == GENERATE_ONLY_IN_ROOT &&
4149 is_rootnode && par->redsplit_generated_in_root) ||
4150 ((par->generate_cgl_redsplit_cuts == GENERATE_DEFAULT ||
4151 par->generate_cgl_redsplit_cuts == GENERATE_IF_IN_ROOT) &&
4152 par->redsplit_generated_in_root) ||
4153 (par->generate_cgl_redsplit_cuts == GENERATE_PERIODICALLY &&
4154 (lp_data->lp_count % par->generate_cgl_redsplit_cuts_freq == 0)) ||
4155 is_top_iter){
4156
4157 /* make basis ready first */
4158 termcode = dual_simplex(lp_data, &iterd);
4159 CglRedSplit *redsplit = new CglRedSplit;
4160 redsplit->generateCuts(*si, cutlist);
4161 if ((new_cut_num = cutlist.sizeCuts() - cut_num) > 0) {
4162 if (is_top_iter){
4163 par->redsplit_generated_in_root = TRUE;
4164 }
4165 PRINT(verbosity, 4,
4166 ("%i reduce and split cuts added\n", new_cut_num));
4167 lp_stat->cuts_generated += new_cut_num;
4168 lp_stat->redsplit_cuts += new_cut_num;
4169 if (is_rootnode) {
4170 lp_stat->cuts_root += new_cut_num;
4171 lp_stat->redsplit_cuts_root += new_cut_num;
4172 }
4173 }
4174 cut_num = cutlist.sizeCuts();
4175 delete redsplit;
4176 cut_time = used_time(&total_time);
4177 comp_times->cuts += cut_time;
4178 comp_times->redsplit_cuts += cut_time;
4179 }
4180 }
4181
4182 /* create CGL odd hole cuts */
4183 should_generate_this_cgl_cut(cut_num, max_cuts_before_resolve,
4184 par->generate_cgl_oddhole_cuts,
4185 par->generate_cgl_oddhole_cuts_freq,
4186 bc_level, bc_index,
4187 lp_stat->oddhole_cuts_root,
4188 &should_generate);
4189 if (should_generate==TRUE) {
4190 CglOddHole *oddhole = new CglOddHole;
4191 //#if 0
4192 oddhole->setMinimumViolation(0.005);
4193 oddhole->setMinimumViolationPer(0.00002);
4194 oddhole->setMaximumEntries(200);
4195 //#endif
4196 oddhole->generateCuts(*si, cutlist);
4197 if ((new_cut_num = cutlist.sizeCuts() - cut_num) > 0) {
4198 if (is_top_iter){
4199 par->oddhole_generated_in_root = TRUE;
4200 }
4201 PRINT(verbosity, 4,
4202 ("%i odd hole cuts added\n", new_cut_num));
4203 lp_stat->cuts_generated += new_cut_num;
4204 lp_stat->oddhole_cuts += new_cut_num;
4205 if (is_rootnode) {
4206 lp_stat->cuts_root += new_cut_num;
4207 lp_stat->oddhole_cuts_root += new_cut_num;
4208 }
4209 }
4210 cut_num = cutlist.sizeCuts();
4211 delete oddhole;
4212 cut_time = used_time(&total_time);
4213 comp_times->cuts += cut_time;
4214 comp_times->oddhole_cuts += cut_time;
4215 }
4216 }
4217
4218 /* create CGL mir cuts */
4219 if(par->generate_cgl_mir_cuts > -1 &&
4220 par->generate_cgl_mir_cuts_freq > 0){
4221 if(par->generate_cgl_mir_cuts == GENERATE_ALWAYS ||
4222 (par->generate_cgl_mir_cuts == GENERATE_ONLY_IN_ROOT &&
4223 is_rootnode && par->mir_generated_in_root) ||
4224 ((par->generate_cgl_mir_cuts == GENERATE_DEFAULT ||
4225 par->generate_cgl_mir_cuts == GENERATE_IF_IN_ROOT) &&
4226 par->mir_generated_in_root) ||
4227 (par->generate_cgl_mir_cuts == GENERATE_PERIODICALLY &&
4228 (lp_data->lp_count % par->generate_cgl_mir_cuts_freq == 0)) ||
4229 is_top_iter){
4230
4231 CglMixedIntegerRounding *mir = new CglMixedIntegerRounding;
4232 if (bc_level<6) {
4233 mir->setMAXAGGR_(5); // default __seems__ 1
4234 }
4235 mir->generateCuts(*si, cutlist);
4236 if ((new_cut_num = cutlist.sizeCuts() - cut_num) > 0) {
4237 if (is_top_iter){
4238 par->mir_generated_in_root = TRUE;
4239 }
4240 PRINT(verbosity, 4,
4241 ("%i MIR cuts added\n", new_cut_num));
4242 lp_stat->cuts_generated += new_cut_num;
4243 lp_stat->mir_cuts += new_cut_num;
4244 if (is_rootnode) {
4245 lp_stat->cuts_root += new_cut_num;
4246 lp_stat->mir_cuts_root += new_cut_num;
4247 }
4248 }
4249 cut_num = cutlist.sizeCuts();
4250 delete mir;
4251 cut_time = used_time(&total_time);
4252 comp_times->cuts += cut_time;
4253 comp_times->mir_cuts += cut_time;
4254 }
4255 }
4256
4257
4258 /* create CGL simple rounding cuts */
4259 if(par->generate_cgl_rounding_cuts > -1 &&
4260 par->generate_cgl_rounding_cuts_freq > 0){
4261 if(par->generate_cgl_rounding_cuts == GENERATE_ALWAYS ||
4262 (par->generate_cgl_rounding_cuts == GENERATE_ONLY_IN_ROOT &&
4263 is_rootnode && par->rounding_generated_in_root) ||
4264 ((par->generate_cgl_rounding_cuts == GENERATE_DEFAULT ||
4265 par->generate_cgl_rounding_cuts == GENERATE_IF_IN_ROOT) &&
4266 par->rounding_generated_in_root) ||
4267 (par->generate_cgl_rounding_cuts == GENERATE_PERIODICALLY &&
4268 (lp_data->lp_count % par->generate_cgl_rounding_cuts_freq == 0)) ||
4269 is_top_iter){
4270
4271 CglSimpleRounding * rounding = new CglSimpleRounding;
4272 rounding->generateCuts(*(lp_data->si), cutlist);
4273 if ((new_cut_num = cutlist.sizeCuts() - cut_num) > 0) {
4274 if (is_top_iter){
4275 par->rounding_generated_in_root = TRUE;
4276 }
4277 PRINT(verbosity, 4,
4278 ("%i rounding cuts added\n", new_cut_num));
4279 lp_stat->cuts_generated += new_cut_num;
4280 lp_stat->rounding_cuts += new_cut_num;
4281 if (is_rootnode) {
4282 lp_stat->cuts_root += new_cut_num;
4283 lp_stat->rounding_cuts_root += new_cut_num;
4284 }
4285 }
4286 cut_num = cutlist.sizeCuts();
4287 delete rounding;
4288 //printf("%i\n", cutlist.sizeCuts());
4289 cut_time = used_time(&total_time);
4290 comp_times->cuts += cut_time;
4291 comp_times->rounding_cuts += cut_time;
4292 }
4293 }
4294
4295 /* create CGL liftandproject cuts (currently buggy) */
4296 if(par->generate_cgl_lift_and_project_cuts > -1 &&
4297 par->generate_cgl_lift_and_project_cuts_freq > 0){
4298 if(par->generate_cgl_lift_and_project_cuts == GENERATE_ALWAYS ||
4299 (par->generate_cgl_lift_and_project_cuts == GENERATE_ONLY_IN_ROOT &&
4300 is_rootnode && par->lift_and_project_generated_in_root) ||
4301 ((par->generate_cgl_lift_and_project_cuts == GENERATE_DEFAULT ||
4302 par->generate_cgl_lift_and_project_cuts == GENERATE_IF_IN_ROOT) &&
4303 par->lift_and_project_generated_in_root) ||
4304 (par->generate_cgl_lift_and_project_cuts == GENERATE_PERIODICALLY &&
4305 (lp_data->lp_count % par->generate_cgl_lift_and_project_cuts_freq == 0)) ||
4306 is_top_iter){
4307
4308 CglLiftAndProject *liftandproject = new CglLiftAndProject;
4309 liftandproject->generateCuts(*(lp_data->si), cutlist);
4310 if ((new_cut_num = cutlist.sizeCuts() - cut_num) > 0) {
4311 if (is_top_iter){
4312 par->lift_and_project_generated_in_root = TRUE;
4313 }
4314 PRINT(verbosity, 4,
4315 ("%i lift and project cuts added\n", new_cut_num));
4316 lp_stat->cuts_generated += new_cut_num;
4317 lp_stat->lift_and_project_cuts += new_cut_num;
4318 if (is_rootnode) {
4319 lp_stat->cuts_root += new_cut_num;
4320 lp_stat->lift_and_project_cuts_root += new_cut_num;
4321 }
4322 }
4323 cut_num = cutlist.sizeCuts();
4324 delete liftandproject;
4325 cut_time = used_time(&total_time);
4326 comp_times->cuts += cut_time;
4327 comp_times->lift_and_project_cuts += cut_time;
4328 }
4329 }
4330
4331 /* create CGL LandP cuts */
4332 #ifndef __OSI_CLP__
4333 // PRINTF(verbosity, -1,
4334 // ("LandP cuts can be generated only with Clp...Skipping LandP cut generation..."));
4335 // /* { */
4336 // }
4337 par->generate_cgl_landp_cuts = DO_NOT_GENERATE;
4338 #else
4339 if(par->generate_cgl_landp_cuts > -1 &&
4340 par->generate_cgl_landp_cuts_freq > 0){
4341 if(par->generate_cgl_landp_cuts == GENERATE_ALWAYS ||
4342 (par->generate_cgl_landp_cuts == GENERATE_ONLY_IN_ROOT &&
4343 is_rootnode && par->landp_generated_in_root) ||
4344 ((par->generate_cgl_landp_cuts == GENERATE_DEFAULT ||
4345 par->generate_cgl_landp_cuts == GENERATE_IF_IN_ROOT) &&
4346 par->landp_generated_in_root) ||
4347 (par->generate_cgl_landp_cuts == GENERATE_PERIODICALLY &&
4348 (lp_data->lp_count % par->generate_cgl_landp_cuts_freq == 0)) ||
4349 is_top_iter){
4350 /* make basis ready first */
4351 termcode = dual_simplex(lp_data, &iterd);
4352 /* if(termcode != 0){
4353 write_mps(lp_data, "lanp.mps");
4354 si->initialSolve();
4355 lp_data->objval = si->getObjValue();
4356 dual_simplex(lp_data, &iterd);
4357 } */
4358 CglLandP *landp = new CglLandP;
4359 //landp->parameter().pivotLimit = 10;
4360 landp->parameter().maxCutPerRound = 30;
4361 landp->generateCuts(*si, cutlist);
4362 if ((new_cut_num = cutlist.sizeCuts() - cut_num) > 0) {
4363 if (is_top_iter){
4364 par->landp_generated_in_root = TRUE;
4365 }
4366 PRINT(verbosity, 4,
4367 ("%i landp cuts added\n", new_cut_num));
4368 lp_stat->cuts_generated += new_cut_num;
4369 lp_stat->landp_cuts += new_cut_num;
4370 if (is_rootnode) {
4371 lp_stat->cuts_root += new_cut_num;
4372 lp_stat->landp_cuts_root += new_cut_num;
4373 }
4374 }
4375 cut_num = cutlist.sizeCuts();
4376 delete landp;
4377 cut_time = used_time(&total_time);
4378 comp_times->cuts += cut_time;
4379 comp_times->landp_cuts += cut_time;
4380 }
4381 }
4382 #endif
4383
4384 #endif
4385
4386 return;
4387 }
4388
4389
4390 /*===========================================================================*/
should_generate_this_cgl_cut(int cut_num,int max_cuts_before_resolve,int generation_flag,int freq,int bc_level,int bc_index,int cuts_in_root,int * should_generate)4391 int should_generate_this_cgl_cut(int cut_num, int max_cuts_before_resolve,
4392 int generation_flag, int freq, int bc_level,
4393 int bc_index, int cuts_in_root, int *should_generate)
4394 {
4395 if (cut_num > max_cuts_before_resolve) {
4396 *should_generate = FALSE;
4397 return 0;
4398 }
4399 switch (generation_flag) {
4400 case (GENERATE_DEFAULT):
4401 if (freq>0 && (bc_level<6 || bc_index % freq == 0)) {
4402 *should_generate = TRUE;
4403 } else {
4404 *should_generate = FALSE;
4405 }
4406 break;
4407 case (GENERATE_ALWAYS):
4408 *should_generate = TRUE;
4409 break;
4410 case (GENERATE_ONLY_IN_ROOT):
4411 if (bc_level<1) {
4412 *should_generate = TRUE;
4413 } else {
4414 *should_generate = FALSE;
4415 }
4416 break;
4417 case (GENERATE_IF_IN_ROOT):
4418 if (bc_level<1) {
4419 *should_generate = TRUE;
4420 } else if (cuts_in_root>0 && bc_index % freq == 0) {
4421 *should_generate = TRUE;
4422 } else {
4423 *should_generate = FALSE;
4424 }
4425 break;
4426 case (GENERATE_PERIODICALLY):
4427 if (bc_index % freq == 0) {
4428 *should_generate = TRUE;
4429 } else {
4430 *should_generate = FALSE;
4431 }
4432 break;
4433 default:
4434 *should_generate = FALSE;
4435 }
4436 return 0;
4437 }
4438
4439
4440 /*===========================================================================*/
check_cuts(OsiCuts & cutlist,LPdata * lp_data,int bc_level,int * num_cuts,cut_data *** cuts,char send_to_pool,int * bnd_changes,lp_stat_desc * lp_stat,node_times * comp_times,int verbosity)4441 int check_cuts(OsiCuts &cutlist, LPdata *lp_data, int bc_level, int
4442 *num_cuts, cut_data ***cuts, char send_to_pool, int *bnd_changes,
4443 lp_stat_desc *lp_stat, node_times *comp_times, int verbosity)
4444 {
4445 #if 0
4446 //OsiCuts cutlist = *cutlist_p;
4447 OsiRowCut cut;
4448 int i, j, k, sizeColCuts;
4449 var_desc **vars = lp_data->vars;
4450 double cut_time, total_time;
4451 /* twice is necessary */
4452 cut_time = used_time(&total_time);
4453 cut_time = used_time(&total_time);
4454 if (cutlist.sizeRowCuts() > 0){
4455 int num_discarded_cuts = 0;
4456 int num_unviolated_cuts = 0;
4457 int num_duplicate_cuts = 0;
4458 int *tmp_matind = lp_data->tmp.i1;
4459 int *is_deleted = (int *) calloc(cutlist.sizeRowCuts(), ISIZE);
4460 double *hashes = (double *) malloc(cutlist.sizeRowCuts()* DSIZE);
4461 int num_elements, num_elements2;
4462 const int *indices, *indices2;
4463 const double *elements, *elements2;
4464 double min_coeff, max_coeff;
4465 int discard_cut, is_duplicate;
4466 double rhs, rhs2;
4467 const double max_elements = (bc_level < 1) ? 100 : 100;
4468 double hash_value, violation;
4469 double *random_hash = lp_data->random_hash;
4470 const double *x = lp_data->x;
4471 const double lpetol = lp_data->lpetol;
4472 const double etol1000 = lpetol*1000;
4473 double *matval;
4474 int *matind;
4475
4476 if (*cuts){
4477 *cuts = (cut_data **)realloc(*cuts, (*num_cuts+cutlist.sizeRowCuts())
4478 * sizeof(cut_data *));
4479 }else{
4480 *cuts = (cut_data **)malloc(cutlist.sizeRowCuts()*sizeof(cut_data *));
4481 }
4482
4483 for (i = 0, j = *num_cuts; i < cutlist.sizeRowCuts(); i++){
4484 cut = cutlist.rowCut(i);
4485 num_elements = cut.row().getNumElements();
4486 indices = cut.row().getIndices();
4487 elements = cut.row().getElements();
4488 rhs = cut.rhs();
4489 discard_cut = FALSE;
4490 max_coeff = 0;
4491 min_coeff = DBL_MAX;
4492
4493 if (num_elements > max_elements) {
4494 is_deleted[i] = TRUE;
4495 PRINT(verbosity,5,("Threw out cut because its length %d is too "
4496 "high.\n\n\n", num_elements));
4497 num_discarded_cuts++;
4498 continue;
4499 }
4500
4501 /*
4502 * Find the largest and the smallest non-zero coeffs to test the
4503 * numerical stability of the cut
4504 * also calculate the hash value
4505 */
4506
4507 hash_value = 0;
4508 violation = 0;
4509 for (int el_num=0; el_num<num_elements; el_num++) {
4510 if (fabs(elements[el_num])>max_coeff) {
4511 max_coeff = fabs(elements[el_num]);
4512 }
4513 if (fabs(elements[el_num]) < min_coeff) {
4514 min_coeff = fabs(elements[el_num]);
4515 }
4516 tmp_matind[el_num] = vars[indices[el_num]]->userind;
4517 hash_value += elements[el_num]*random_hash[tmp_matind[el_num]];
4518 violation += elements[el_num]*x[tmp_matind[el_num]];
4519 }
4520 /* see rhs as well */
4521 if (fabs(rhs) > lpetol) {
4522 if (fabs(rhs) < min_coeff) {
4523 min_coeff = fabs(rhs);
4524 }
4525 if (fabs(rhs) > max_coeff) {
4526 max_coeff = fabs(rhs);
4527 }
4528 }
4529 switch (cut.sense()) {
4530 case 'L':
4531 violation -= rhs;
4532 break;
4533 case 'G':
4534 violation = rhs - violation;
4535 break;
4536 case 'E':
4537 violation = fabs(rhs - violation);
4538 break;
4539 }
4540
4541 /*
4542 * display the cut
4543 */
4544 if (verbosity>11) {
4545 PRINT(12, 11, ("Cut #%i: rhs = %f sense = %c\n", i, rhs,
4546 cut.sense()));
4547 for (int el_num=0; el_num<num_elements; el_num++) {
4548 PRINT(12,11,("%d\t%f\n",indices[el_num],elements[el_num]));
4549 }
4550 }
4551 PRINT(verbosity,5,("generate_cgl_cuts: Number of Coefficients = "
4552 "%d\tMax = %f, Min = %f\n",num_elements,max_coeff,
4553 min_coeff));
4554
4555 if (violation < lpetol) {
4556 num_unviolated_cuts++;
4557 discard_cut = TRUE;
4558 PRINT(verbosity,5,("violation = %f. Threw out cut.\n\n",
4559 violation));
4560 }
4561
4562 if (discard_cut != TRUE && num_elements>0) {
4563 if ( (max_coeff > 0 && min_coeff/max_coeff < etol1000)||
4564 (min_coeff > 0 && min_coeff < etol1000) ) {
4565 num_discarded_cuts++;
4566 discard_cut = TRUE;
4567 PRINT(verbosity,5,("Threw out cut because of bad coeffs.\n\n"));
4568 }
4569 }
4570 if (discard_cut==TRUE) {
4571 is_deleted[i] = TRUE;
4572 continue;
4573 }
4574
4575 /* check for duplicates */
4576 if (num_elements>0) {
4577 is_duplicate = FALSE;
4578 /* check against last 50 cuts only. otherwise, takes a lot of time
4579 */
4580 /* for (k = j-1; k > MAX(-1,j-51); k--) */
4581 for (k=j-1; k>-1; k--) {
4582 num_elements2 = ((int *) ((*cuts)[k]->coef))[0];
4583 rhs2 = (*cuts)[k]->rhs;
4584 if (num_elements2 != num_elements ||
4585 fabs(rhs2 - rhs) > lpetol ||
4586 fabs(hashes[k]-hash_value) > lpetol) {
4587 continue;
4588 } else {
4589 break;
4590 }
4591 }
4592 /* if (k>MAX(-1,i-51)) */
4593 if (k>-1) {
4594 is_deleted[i] = TRUE;
4595 PRINT(verbosity,5,("cut #%d is same as accepted cut #%d\n",i,k));
4596 num_duplicate_cuts++;
4597 continue;
4598 }
4599 }
4600
4601 (*cuts)[j] = (cut_data *) calloc(1, sizeof(cut_data));
4602 (*cuts)[j]->type = EXPLICIT_ROW;
4603 if (((*cuts)[j]->sense = cut.sense()) == 'R'){
4604 FREE((*cuts)[j]);
4605 continue; /* This must be a bug. */
4606 }
4607 hashes[j] = hash_value;
4608 PRINT(verbosity, 12, ("Cut #%i: accepted as cut number %i\n", i, j));
4609
4610 (*cuts)[j]->rhs = rhs;
4611 (*cuts)[j]->range = cut.range();
4612 (*cuts)[j]->size = (num_elements + 1) * (ISIZE + DSIZE);
4613 (*cuts)[j]->coef = (char *) malloc ((*cuts)[j]->size);
4614 ((int *) ((*cuts)[j]->coef))[0] = num_elements;
4615 //Here, we have to pad the initial int to avoid misalignment, so we
4616 //add DSIZE bytes to get to a double boundary
4617 matval = (double *) ((*cuts)[j]->coef + DSIZE);
4618 matind = (int *) ((*cuts)[j]->coef + (num_elements + 1)*DSIZE);
4619 memcpy((char *)matval, (char *)elements, num_elements * DSIZE);
4620 memcpy((char*)matind, (char *)tmp_matind, num_elements * ISIZE);
4621 qsort_id(matind, matval, num_elements);
4622
4623 (*cuts)[j]->branch = DO_NOT_BRANCH_ON_THIS_ROW;
4624
4625 (*cuts)[j]->deletable = TRUE;
4626 if (send_to_pool){
4627 (*cuts)[j++]->name = CUT__SEND_TO_CP;
4628 }else{
4629 (*cuts)[j++]->name = CUT__DO_NOT_SEND_TO_CP;
4630 }
4631 }
4632 *num_cuts = j;
4633 if (num_discarded_cuts>0) {
4634 PRINT(verbosity,3,("generate_cgl_cuts: Number of discarded cuts = %d\n",num_discarded_cuts));
4635 }
4636 lp_stat->num_poor_cuts += num_discarded_cuts;
4637 lp_stat->num_unviolated_cuts += num_unviolated_cuts;
4638 if (num_duplicate_cuts>0) {
4639 PRINT(verbosity,3,("generate_cgl_cuts: Number of duplicate cuts = %d\n",num_duplicate_cuts));
4640 }
4641 lp_stat->num_duplicate_cuts += num_duplicate_cuts;
4642 FREE(is_deleted);
4643 FREE(hashes);
4644 }
4645
4646 sizeColCuts = cutlist.sizeColCuts();
4647 if (sizeColCuts > 0){
4648 PRINT(verbosity,3,("cgl_generate_cuts: %d colCuts found\n",sizeColCuts));
4649 OsiColCut colCut;
4650 const int *indices;
4651 const double *elements;
4652 for (i=0;i<sizeColCuts;i++) {
4653 colCut = cutlist.colCut(i);
4654 if (verbosity>10) {
4655 colCut.print();
4656 }
4657 indices = colCut.lbs().getIndices();
4658 elements = colCut.lbs().getElements();
4659 for (j=0;j<colCut.lbs().getNumElements();j++) {
4660 if (vars[indices[j]]->new_lb < elements[j]) {
4661 vars[indices[j]]->new_lb = elements[j];
4662 change_lbub(lp_data, indices[j], elements[j],
4663 vars[indices[j]]->new_ub);
4664 (*bnd_changes)++;
4665 }
4666 }
4667 indices = colCut.ubs().getIndices();
4668 elements = colCut.ubs().getElements();
4669 for (j=0;j<colCut.ubs().getNumElements();j++) {
4670 if (vars[indices[j]]->new_ub > elements[j]) {
4671 vars[indices[j]]->new_ub = elements[j];
4672 change_lbub(lp_data, indices[j], vars[indices[j]]->new_lb,
4673 elements[j]);
4674 (*bnd_changes)++;
4675 }
4676 }
4677 }
4678 //exit(0);
4679 }
4680 cut_time = used_time(&total_time);
4681 comp_times->dupes_and_bad_coeffs_in_cuts += cut_time;
4682 #endif
4683 return 0;
4684 }
4685 #endif
4686 #endif /* __OSI_xxx__ */
4687 /*===========================================================================*/
4688 /*===========================================================================*/
4689