1 /*===========================================================================*/
2 /* */
3 /* This file is part of the SYMPHONY Branch, Cut, and Price Library. */
4 /* */
5 /* SYMPHONY was jointly developed by Ted Ralphs (ted@lehigh.edu) and */
6 /* Laci Ladanyi (ladanyi@us.ibm.com). */
7 /* */
8 /* The author of this file is Menal Guzelsoy */
9 /* */
10 /* (c) Copyright 2006-2019 Lehigh University. All Rights Reserved. */
11 /* */
12 /* This software is licensed under the Eclipse Public License. Please see */
13 /* accompanying file for terms. */
14 /* */
15 /*===========================================================================*/
16 /* last modified: June 09, menal*/
17
18 #include <string.h>
19 #include <stdio.h>
20 #include <stdlib.h>
21
22 #include "sym_qsort.h"
23 #include "sym_macros.h"
24 #include "sym_constants.h"
25 #include "sym_prep.h"
26
27 /* Under Development */
28
29 /*===========================================================================*/
30 /*===========================================================================*/
sr_initialize(SRdesc ** sr,int n)31 void sr_initialize(SRdesc **sr, int n){
32
33 int do_clean = FALSE;
34
35 if(!(*sr)){
36 *sr = (SRdesc *)calloc(1, sizeof(SRdesc));
37 do_clean = TRUE;
38 }
39
40 if(!do_clean){
41 (*sr)->prob_type = 0;
42 (*sr)->max_n = (*sr)->min_n = 0;
43 (*sr)->ub = (*sr)->lb = 0.0;
44 (*sr)->ub_offset = (*sr)->lb_offset = 0.0;
45 (*sr)->sum_a_max = (*sr)->sum_a_min = 0.0;
46 (*sr)->sum_c_max = (*sr)->sum_c_min = 0.0;
47 (*sr)->ub_updated = (*sr)->lb_updated = FALSE;
48 (*sr)->rhs = (*sr)->rhs_max = (*sr)->rhs_min = 0.0;
49 (*sr)->sense = ' ';
50 if((*sr)->obj_max){
51 memset((*sr)->reversed_max, FALSE, CSIZE*n);
52 memset((*sr)->reversed_min, FALSE, CSIZE*n);
53 memset((*sr)->var_stat_max, SR_VAR_IN, ISIZE*n);
54 memset((*sr)->var_stat_min, SR_VAR_IN, ISIZE*n);
55 }
56 }
57 }
58
59
60
61 /*===========================================================================*/
62 /*===========================================================================*/
63
sr_allocate(SRdesc ** sr,int n)64 void sr_allocate(SRdesc **sr, int n){
65
66 int k;
67 (*sr)->obj_max = (double *)malloc(DSIZE*n);
68 (*sr)->matval_max = (double *)malloc(DSIZE*n);
69 (*sr)->matind_max = (int *)malloc(ISIZE*n);
70 (*sr)->ratio_max = (double *)malloc(DSIZE*n);
71 (*sr)->reversed_max = (char *)malloc(CSIZE*n);
72
73 (*sr)->obj_min = (double *)malloc(DSIZE*n);
74 (*sr)->matval_min = (double *)malloc(DSIZE*n);
75 (*sr)->matind_min = (int *)malloc(ISIZE*n);
76 (*sr)->ratio_min = (double *)malloc(DSIZE*n);
77 (*sr)->reversed_min = (char *)malloc(CSIZE*n);
78
79 /* for variable fixing, tightening etc... */
80
81 (*sr)->var_max_opt = (double *)malloc(n* DSIZE);
82 (*sr)->var_min_opt = (double *)malloc(n* DSIZE);
83 (*sr)->var_stat_max = (int *)malloc(ISIZE*n);
84 (*sr)->var_stat_min = (int *)malloc(n* ISIZE);
85 (*sr)->var_obj_max = (double *)malloc(n* DSIZE);
86 (*sr)->var_obj_min = (double *)malloc(n* DSIZE);
87 (*sr)->var_matval_max = (double *)malloc(n* DSIZE);
88 (*sr)->var_matval_min = (double *)malloc(n* DSIZE);
89
90 /* debug, get something smart instead of these */
91 (*sr)->tmp_ind = (int *)malloc(ISIZE*n);
92 (*sr)->fixed_ind = (int *)malloc(ISIZE*n);
93
94 for(k = 0; k < n; k++){
95 (*sr)->fixed_ind[k] = k;
96 }
97 }
98
99 /*===========================================================================*/
100 /*===========================================================================*/
prep_solve_sr_rlx(PREPdesc * P,int row_cnt,int * row_indices)101 int prep_solve_sr_rlx(PREPdesc *P, int row_cnt, int *row_indices)
102 {
103
104 int i, j, k, l;
105 int termcode = SR_NO_UPDATES;
106
107 MIPdesc * mip = P->mip;
108 prep_params params = P->params;
109 MIPinfo *mip_inf = mip->mip_inf;
110
111 COLinfo *cols = mip_inf->cols;
112 ROWinfo *rows = mip_inf->rows;
113
114 int n = mip->n, m = mip->m;
115 int *c_matbeg = mip->matbeg;
116 int *c_matind = mip->matind;
117
118 int * r_matbeg = mip->row_matbeg;
119 int * r_matind = mip->row_matind;
120 double * r_matval = mip->row_matval;
121 double *rhs = mip->rhs;
122 char *sense = mip->sense;
123
124 double *ub = mip->ub;
125 double *lb = mip->lb;
126
127 int max_sr_cnt, max_aggr_cnt, verbosity; //max_aggr_row_num, verbosity;
128 int p_level, do_sr_rlx, do_aggr_row_rlx;
129 int obj_ind, tot_sub_pr;
130 char can_iterate = TRUE;
131 double etol;
132
133 do_sr_rlx = params.do_single_row_rlx;
134 do_aggr_row_rlx = params.do_aggregate_row_rlx;
135 p_level = params.level;
136 etol = params.etol;
137 verbosity = params.verbosity;
138
139 /* get the max iteration limits */
140 max_sr_cnt = params.max_sr_cnt;
141 max_aggr_cnt = 0;
142
143 /* initialize arrays to be used for each subproblem*/
144
145
146 SRdesc * sr, *d_sr;
147
148 if(!(P->rows_checked)){
149 P->rows_checked = (char *)malloc(m* CSIZE);
150 }
151
152 char *rows_checked = P->rows_checked;
153 double old_bound;
154 //char no_upper, no_lower;
155 int row_ind;// const_row_ind;
156 //int updated_lb_cnt = 0, updated_ub_cnt = 0;
157 int last_col_loc, last_row_loc;
158
159 /*max_sr_cnt should be up to some ratio!!!! may be too many to handle*/
160 tot_sub_pr = max_sr_cnt + max_aggr_cnt;
161
162 /* do this only for all unbounded or all bounded rows for now...*/
163 /* fixme... extend this if results happen to be good...*/
164
165 for(i = 0; i < row_cnt; i++){
166
167 obj_ind = row_indices[i];
168
169 if(rows[obj_ind].bound_type == MIXED_BOUNDED_ROW ||
170 rows[obj_ind].is_redundant){
171 continue;
172 }
173
174 rows[obj_ind].orig_ub = rows[obj_ind].sr_ub = rows[obj_ind].ub;
175 rows[obj_ind].orig_lb = rows[obj_ind].sr_lb = rows[obj_ind].lb;
176
177 if(verbosity >=4){
178 printf("init bounds: row: %i", i);
179 printf("\told_lb:");
180 if(rows[obj_ind].sr_lb > -INF){
181 printf("%f", rows[obj_ind].sr_lb);
182 }else{
183 printf("-inf");
184 }
185 printf("\told_ub:");
186 if(rows[obj_ind].sr_ub < INF){
187 printf("%f", rows[obj_ind].sr_ub);
188 }else{
189 printf("inf");
190 }
191 printf("\n");
192 }
193
194
195 // srows[i] = (SRrlx *)calloc(tot_sub_pr, sizeof(SRrlx));
196 memset(rows_checked, FALSE, CSIZE*m);
197 last_col_loc = r_matbeg[obj_ind];
198 last_row_loc = c_matbeg[r_matind[last_col_loc]];
199
200 for(j = 0; j < tot_sub_pr; j++){
201 /* first do single row stuff if can*/
202 // is_open_prob = 1;
203 row_ind = -1;
204 can_iterate = TRUE;
205 if(j < max_sr_cnt){
206 //int max_shared_row_ind = -1;
207 //int max_shared_size = 0;
208 //int row_search_iter_cnt = 0;
209 /* get something smarter here */
210 /*find a row that has the most common shared vars with this
211 one */
212 /*find a row to be used as a constraint*/
213 for(k = last_col_loc; k < r_matbeg[obj_ind+1]; k++){
214 for(l = last_row_loc; l < c_matbeg[r_matind[k]+1];
215 l++){
216 if(!rows[c_matind[l]].is_redundant &&
217 !rows_checked[c_matind[l]]){
218 rows_checked[c_matind[l]] = TRUE;
219 if(rows[obj_ind].bound_type ==
220 rows[c_matind[l]].bound_type
221 && c_matind[l] != obj_ind) {
222 row_ind = c_matind[l];
223 break;
224 }
225 }
226 }
227
228 if(row_ind >= 0){
229 last_col_loc = k;
230 last_row_loc = l;
231 break;
232 }
233 }
234
235 if(row_ind >= 0){
236
237 sr_initialize(&(P->sr), n);
238 sr = P->sr;
239 sr->prob_type = rows[obj_ind].bound_type;
240 sr->rhs = rhs[row_ind];
241 sr->sense = sense[row_ind];
242
243 /* convert the problem to <= constraint*/
244 /* or solve it if it is unbounded */
245
246 switch(rows[obj_ind].bound_type){
247 case OPEN_ROW: /* easiest case */
248
249 sr->rhs_max = sr->rhs_min = sr->rhs;
250
251 sr_solve_open_prob(P, sr, obj_ind, row_ind, r_matbeg,
252 r_matind, r_matval, cols, ub, lb, etol);
253
254 break;
255 case ALL_BOUNDED_ROW:
256 if(rows[obj_ind].ub_inf_var_num +
257 rows[obj_ind].lb_inf_var_num +
258 rows[obj_ind].free_var_num > 0 ||
259 rows[row_ind].ub_inf_var_num +
260 rows[row_ind].lb_inf_var_num +
261 rows[row_ind].free_var_num > 0){
262 /* debug */
263 /* fixme, get rid of this bug*/
264 printf("something is wrong -case all_bounded_row-"
265 "prep_solve_sr_rlx(), exiting...\n");
266 return PREP_OTHER_ERROR;
267 }
268
269 /* always convert the problem into ax <= b ineq for max
270 ax >= b for min */
271
272 /* so,
273
274 _max arrays for solving [max cx st. ax <= b]
275 _min arrays for solving [min cx st. ax >= b]
276
277 if sense = E;
278
279 d_max arrays for solving [max cx st. -ax <= -b]
280 d_min arrays for solving [min cx st. -ax >= -b]
281
282 */
283
284 if(!sr->obj_max && rows[obj_ind].bound_type != OPEN_ROW){
285 sr_allocate(&sr, n);
286 }
287 switch(sr->sense){
288 case 'G':
289 sr->rhs_max = -sr->rhs;
290 sr->rhs_min = sr->rhs;
291 break;
292 case 'L':
293 sr->rhs_max = sr->rhs;
294 sr->rhs_min = -sr->rhs;
295 break;
296 case 'E':
297 sr->rhs_max = sr->rhs;
298 sr->rhs_min = -sr->rhs;
299
300 sr_initialize(&(P->d_sr), n);
301 d_sr = P->d_sr;
302 d_sr->prob_type = rows[obj_ind].bound_type;
303 d_sr->rhs = rhs[row_ind];
304 d_sr->sense = sense[row_ind];
305
306 d_sr->rhs_max = -d_sr->rhs;
307 d_sr->rhs_min = d_sr->rhs;
308
309 if(!d_sr->obj_max){
310 sr_allocate(&d_sr, n);
311 }
312
313 break;
314 }
315
316 sr_solve_bounded_prob(P, sr, d_sr, obj_ind, row_ind,
317 r_matbeg, r_matind, r_matval,
318 cols, ub, lb, etol);
319 if(!rows[obj_ind].is_redundant){
320 if(sr->sense == 'E'){
321 if(sr->ub > d_sr->ub){
322 sr->ub = d_sr->ub;
323 }
324 if(sr->lb < d_sr->lb){
325 sr->lb = d_sr->lb;
326 }
327 }
328
329 sr->lb_updated = sr->ub_updated = TRUE;
330 }else{
331 break;
332 }
333 }
334
335 /* check for any progress */
336 if(sr->lb_updated){
337 if(rows[obj_ind].sr_lb < sr->lb){
338 old_bound = rows[obj_ind].sr_lb;
339 rows[obj_ind].sr_lb = sr->lb;
340 /* debug */
341 if(termcode != SR_BOUNDS_UPDATED){
342 termcode = SR_BOUNDS_UPDATED;
343 }
344
345 if(verbosity >=5){
346 printf("lb improved, "
347 "row: %i \told_lb:%f \tnew_lb:%f\n",
348 obj_ind, old_bound <= -INF ? 1 : old_bound, sr->lb);
349 }
350 }else if (rows[obj_ind].orig_lb > sr->lb + etol){
351 /* debug */
352 printf("error-lb, row: %i \told_lb:%f \tnew_lb:%f\n",
353 obj_ind, rows[obj_ind].orig_lb, sr->lb);
354 }
355 }
356 if(sr->ub_updated){
357 if(rows[obj_ind].sr_ub > sr->ub){
358 old_bound = rows[obj_ind].sr_ub;
359 rows[obj_ind].sr_ub = sr->ub;
360 /* debug */
361 if(termcode != SR_BOUNDS_UPDATED){
362 termcode = SR_BOUNDS_UPDATED;
363 }
364 if(verbosity >=5){
365 printf("ub improved, "
366 "row: %i \told_ub:%f \tnew_ub:%f\n",
367 obj_ind, old_bound >= INF ? -1 : old_bound, sr->ub);
368 }
369 }else if(rows[obj_ind].orig_ub < sr->ub - etol){
370 /*debug*/
371 // if(verbosity >=5){
372 printf("error-ub, row: %i \told_ub:%f \tnew_ub:%f\n",
373 obj_ind, rows[obj_ind].orig_ub, sr->ub);
374 // }
375 }
376 if(sr->lb_updated){
377 if(sr->ub < sr->lb - etol){
378 /* debug */
379 printf("bounds err : "
380 "row: %i \tnew_ub:%f \tnew_lb:%f\n",
381 obj_ind, sr->ub, sr->lb);
382 termcode = SR_INFEAS;
383 break;
384 }
385 }
386 }
387
388 }
389 }
390 }
391
392 /* debug */
393 if(termcode == SR_INFEAS){
394 break;
395 }
396
397 if(verbosity >=4){
398 printf("finl bounds: row: %i", i);
399 printf("\tnew_lb:");
400 if(rows[obj_ind].sr_lb > -INF){
401 printf("%f", rows[obj_ind].sr_lb);
402 }else{
403 printf("-inf");
404 }
405 printf("\tnew_ub:");
406 if(rows[obj_ind].sr_ub < INF){
407 printf("%f", rows[obj_ind].sr_ub);
408 }else{
409 printf("inf");
410 }
411 printf("\n\n");
412
413 }
414 }
415
416 return termcode;
417 }
418
419 /*===========================================================================*/
420 /*===========================================================================*/
421
sr_solve_bounded_prob(PREPdesc * P,SRdesc * sr,SRdesc * d_sr,int obj_ind,int row_ind,int * r_matbeg,int * r_matind,double * r_matval,COLinfo * cols,double * ub,double * lb,double etol)422 int sr_solve_bounded_prob(PREPdesc *P, SRdesc *sr, SRdesc *d_sr,
423 int obj_ind, int row_ind,
424 int *r_matbeg, int *r_matind, double *r_matval,
425 COLinfo *cols, double *ub, double *lb, double etol)
426 {
427
428 int k, l, col_ind;
429 double c_val, a_val;
430
431 for( k = r_matbeg[obj_ind], l = r_matbeg[row_ind];;){
432 if(k < r_matbeg[obj_ind + 1] &&
433 (r_matind[k] < r_matind[l] ||
434 l >= r_matbeg[row_ind + 1])){
435 c_val = r_matval[k];
436 col_ind = r_matind[k];
437 sr_add_new_col(sr, d_sr, c_val, 0.0, col_ind,
438 cols[col_ind].var_type, ub[col_ind], lb[col_ind],
439 sr->sense, 1, 1);
440 k++;
441 }else if(l < r_matbeg[row_ind + 1] &&
442 (r_matind[k] > r_matind[l] ||
443 k >= r_matbeg[obj_ind+1])){
444 a_val = r_matval[l];
445 col_ind = r_matind[l];
446
447 sr_add_new_col(sr, d_sr, 0.0, a_val, col_ind,
448 cols[col_ind].var_type, ub[col_ind], lb[col_ind],
449 sr->sense, 0, 1);
450 l++;
451 }else{
452 /* now the indices are equal, fill in the arrays */
453 c_val = r_matval[k];
454 a_val = r_matval[l];
455 col_ind = r_matind[k];
456
457 if(c_val == 0.0 || a_val == 0.0){
458 printf("not nonzero???"
459 "numerical issues -case bounded row-"
460 "sr_solve_bounded_prob(), exiting...\n");
461 return PREP_OTHER_ERROR;
462 }
463
464 sr_add_new_col(sr, d_sr, c_val, a_val, col_ind,
465 cols[col_ind].var_type, ub[col_ind], lb[col_ind],
466 sr->sense, 2, 1);
467 k++;
468 l++;
469 }
470 if(k == r_matbeg[obj_ind + 1] && l == r_matbeg[row_ind + 1]){
471 break;
472 }
473 }
474
475 /* now solve the problem */
476 if(!P->mip->mip_inf->rows[obj_ind].is_redundant){
477 sr_find_opt_bounded(P, sr, obj_ind, ub, lb);
478 }
479
480 if(!P->mip->mip_inf->rows[obj_ind].is_redundant){
481 if(sr->sense == 'E'){
482 sr_find_opt_bounded(P, d_sr, obj_ind, ub, lb);
483 }
484 }
485
486
487 int termcode = 0;
488 ROWinfo *rows = P->mip->mip_inf->rows;
489 double min_ub = sr->ub;
490 double max_lb = sr->lb;
491
492 if(sr->sense == 'E'){
493 if(!P->mip->mip_inf->rows[obj_ind].is_redundant){
494 if(min_ub > d_sr->ub){
495 min_ub = d_sr->ub;
496 }
497
498 if(max_lb < d_sr->lb){
499 max_lb = d_sr->lb;
500 }
501 }
502 }
503 if(rows[obj_ind].ub > min_ub || rows[obj_ind].lb < max_lb){
504 termcode = prep_check_redundancy(P, obj_ind, TRUE, min_ub, max_lb,
505 FALSE, 0);
506 }
507
508 if(PREP_QUIT(termcode)){
509 return termcode;
510 }
511
512 return(0);
513
514 }
515
516 /*===========================================================================*/
517 /*===========================================================================*/
518
sr_find_opt_bounded(PREPdesc * P,SRdesc * sr,int obj_ind,double * ub,double * lb)519 int sr_find_opt_bounded(PREPdesc *P, SRdesc *sr, int obj_ind,
520 double *ub, double *lb)
521
522 {
523 int i, last_ind, col_loc, col_ind, *var_stat; //,j, var_ind;
524 char max_solved = FALSE, min_solved = FALSE;
525 double lhs, ax, var_frac_val;
526 /* get opt for each column (col fixed ub in min solved and
527 lb in max solved - check also a_vals)*/
528 double bound;
529
530 int * tmp_ind = sr->tmp_ind;
531 double etol = P->params.etol;
532
533 if(sr->sum_a_max < sr->rhs_max +etol || sr->max_n <= 0){
534 sr->ub += sr->sum_c_max + sr->ub_offset;
535 max_solved = TRUE;
536 }
537
538 if(sr->sum_a_min > sr->rhs_min - etol|| sr->min_n <= 0){
539 sr->lb += sr->sum_c_min + sr->lb_offset;
540 min_solved = TRUE;
541 }
542
543 if(max_solved && min_solved){
544 /* check redundancy */
545 /* redundant and useless row*/
546 return PREP_UNMODIFIED;
547 }
548
549 if(!max_solved){ /* otherwise, the row is redundant and useless */
550
551 var_stat = sr->var_stat_max;
552 memcpy(tmp_ind, sr->fixed_ind, ISIZE*sr->max_n);
553 qsort_di(sr->ratio_max, tmp_ind, sr->max_n);
554 //CoinSort_2(sr->ratio_max, sr->ratio_max + sr->max_n, tmp_ind);
555
556 /* now fill in knapsack */
557 lhs = 0;
558 for(i = sr->max_n - 1; i >=0; i--){
559 col_loc = tmp_ind[i];
560 col_ind = sr->matind_max[col_loc];
561 bound = ub[col_ind] - lb[col_ind];
562
563 if(lhs > sr->rhs_max - etol){
564 break;
565 }
566
567 ax = sr->matval_max[col_loc] * bound;
568
569 if(lhs + ax < sr->rhs_max - etol){
570 sr->ub += bound *
571 sr->obj_max[col_loc];
572 lhs += ax;
573 var_stat[col_ind] = SR_VAR_IN_FIXED_UB;
574 }else{
575 var_frac_val = sr->obj_max[col_loc] *
576 (sr->rhs_max - lhs)/sr->matval_max[col_loc];
577 sr->ub += var_frac_val; //sr->obj_max[col_loc] * var_frac_val;
578 var_stat[col_ind] = SR_VAR_IN_FRAC;
579 last_ind = i;
580 break;
581 }
582 }
583 sr->ub += sr->ub_offset;
584 }
585
586 if(!min_solved){ /* otherwise this row is redundant and useless */
587 memcpy(tmp_ind, sr->fixed_ind, ISIZE*sr->min_n);
588 qsort_di(sr->ratio_min, tmp_ind, sr->min_n);
589 //CoinSort_2(sr->ratio_min, sr->ratio_min + sr->min_n, tmp_ind);
590 /* now fill in knapsack */
591 lhs = 0;
592 var_stat = sr->var_stat_min;
593 for(i = 0; i < sr->min_n; i++){
594 col_loc = tmp_ind[i];
595 col_ind = sr->matind_min[col_loc];
596 bound = ub[col_ind] - lb[col_ind];
597 ax = sr->matval_min[col_loc] * bound;
598
599 if(lhs > sr->rhs_min - etol){
600 break;
601 }
602
603 if(lhs + ax < sr->rhs_min - etol){
604 sr->lb += bound *
605 sr->obj_min[col_loc];
606 lhs += ax;
607 var_stat[col_ind] = SR_VAR_IN_FIXED_UB;
608 }else{
609 // ax = (sr->rhs_max - lhs)/sr->matval_max[col_loc];
610 sr->lb += sr->obj_min[col_loc] *
611 (sr->rhs_min - lhs)/sr->matval_min[col_loc];
612 var_stat[col_ind] = SR_VAR_IN_FIXED_UB;
613 last_ind = i;
614 break;
615 }
616 }
617 sr->lb += sr->lb_offset;
618 }
619
620 return 0;
621
622 }
623
624 /*===========================================================================*/
625 /*===========================================================================*/
626
627 /* will add the column to problem if necessary */
628
sr_add_new_col(SRdesc * sr,SRdesc * d_sr,double c_val,double a_val,int col_ind,char var_type,double col_ub,double col_lb,char sense,int col_type,int col_bound_type)629 int sr_add_new_col(SRdesc *sr, SRdesc *d_sr, double c_val, double a_val,
630 int col_ind, char var_type, double col_ub,
631 double col_lb, char sense,
632 int col_type, int col_bound_type)
633 {
634 /* col_type =
635 0 => c_val = 0, a_val != 0
636 1 => c_val != 0, a_val = 0
637 2 => c_val != 0, a_val != 0
638 */
639
640 /* col_bound_type =
641 0 => open row
642 1 => all bounded row
643 2 => mixed bounded row
644 */
645
646 double rhs_ub_offset = a_val * col_ub;
647 double rhs_lb_offset = a_val * col_lb;
648
649 double obj_ub_offset = c_val * col_ub;
650 double obj_lb_offset = c_val * col_lb;
651
652 if(col_bound_type == 1){
653 if(col_type >= 0){
654 if(var_type != 'F'){
655 switch(sense){
656 case 'L':
657 sr_add_new_bounded_col(sr, c_val, a_val, col_ind,
658 rhs_ub_offset, rhs_lb_offset,
659 obj_ub_offset, obj_lb_offset,
660 col_ub, col_lb,
661 SR_MAX, var_type);
662 sr_add_new_bounded_col(sr, c_val, -a_val, col_ind,
663 -rhs_ub_offset, -rhs_lb_offset,
664 obj_ub_offset, obj_lb_offset,
665 col_ub, col_lb,
666 SR_MIN, var_type);
667 break;
668 case 'G':
669 sr_add_new_bounded_col(sr, c_val, -a_val, col_ind,
670 -rhs_ub_offset, -rhs_lb_offset,
671 obj_ub_offset, obj_lb_offset,
672 col_ub, col_lb,
673 SR_MAX, var_type);
674 sr_add_new_bounded_col(sr, c_val, a_val, col_ind,
675 rhs_ub_offset, rhs_lb_offset,
676 obj_ub_offset, obj_lb_offset,
677 col_ub, col_lb,
678 SR_MIN, var_type);
679 break;
680 case 'E':
681 sr_add_new_bounded_col(sr, c_val, a_val, col_ind,
682 rhs_ub_offset, rhs_lb_offset,
683 obj_ub_offset, obj_lb_offset,
684 col_ub, col_lb,
685 SR_MAX, var_type);
686 sr_add_new_bounded_col(sr, c_val, -a_val, col_ind,
687 -rhs_ub_offset, -rhs_lb_offset,
688 obj_ub_offset, obj_lb_offset,
689 col_ub, col_lb,
690 SR_MIN, var_type);
691 sr_add_new_bounded_col(d_sr, c_val, -a_val, col_ind,
692 -rhs_ub_offset, -rhs_lb_offset,
693 obj_ub_offset, obj_lb_offset,
694 col_ub, col_lb,
695 SR_MAX, var_type);
696 sr_add_new_bounded_col(d_sr, c_val, a_val, col_ind,
697 rhs_ub_offset, rhs_lb_offset,
698 obj_ub_offset, obj_lb_offset,
699 col_ub, col_lb,
700 SR_MIN, var_type);
701 break;
702 }
703 }else{
704
705 sr->ub_offset += obj_ub_offset;
706 sr->lb_offset += obj_ub_offset;
707 sr->rhs_max -= rhs_ub_offset;
708 sr->rhs_min -= rhs_ub_offset;
709
710 if(sense == 'E'){
711 d_sr->ub_offset += obj_ub_offset;
712 d_sr->lb_offset += obj_ub_offset;
713 d_sr->rhs_max -= rhs_ub_offset;
714 d_sr->rhs_min -= rhs_ub_offset;
715 }
716 }
717 }
718 }
719 return 0;
720 }
721
722 /*===========================================================================*/
723 /*===========================================================================*/
724
725 /* will add the column to problem if necessary, otherwise will update the
726 offset values.
727 will assume the sense is L for max and G for min, so
728 the a_val has to be sent after being updated.
729 For E, this function will be called twice each for max and min*/
730
sr_add_new_bounded_col(SRdesc * sr,double c_val,double a_val,int col_ind,double rhs_ub_offset,double rhs_lb_offset,double obj_ub_offset,double obj_lb_offset,double col_ub,double col_lb,int obj_sense,char var_type)731 int sr_add_new_bounded_col(SRdesc *sr, double c_val, double a_val,
732 int col_ind,
733 double rhs_ub_offset, double rhs_lb_offset,
734 double obj_ub_offset, double obj_lb_offset,
735 double col_ub, double col_lb, int obj_sense,
736 char var_type)
737 {
738 /*
739 ratio_type =
740 0 c_val >0, a_val>0
741 1 c_val >= 0, a_val <= 0
742 2 c_val <= 0, a_val >= 0
743 3 c_val < 0, a_val < 0
744 */
745
746 // int n;// = sr->max_n, min_n = sr->min_n;
747
748 /* we will convert the vars so that u-l >= x >= 0 */
749
750
751 int ratio_type = 0;
752
753 if(c_val > 0.0){
754 if(a_val <= 0.0){
755 ratio_type = 1;
756 }
757 }else if(c_val < 0.0){
758 if(a_val >= 0.0){
759 ratio_type = 2;
760 }else{
761 ratio_type = 3;
762 }
763 }else{
764 if(a_val <= 0.0){
765 ratio_type = 1;
766 }else{
767 ratio_type = 2;
768 }
769 }
770
771 int *n, *matind, *var_stat;
772 double *obj, *matval, *rhs, *obj_offset, *sum, *obj_sum, *ratios;
773 double *var_matval, *var_obj;
774 char *is_reversed;
775 if(obj_sense == SR_MAX){
776 n = &(sr->max_n);
777 obj_offset = &(sr->ub_offset);
778 sum = &(sr->sum_a_max);
779 obj_sum = &(sr->sum_c_max);
780 rhs = &(sr->rhs_max);
781 obj = sr->obj_max;
782 matind = sr->matind_max;
783 matval = sr->matval_max;
784 ratios = sr->ratio_max;
785 is_reversed = sr->reversed_max;
786 var_stat = sr->var_stat_max;
787 var_matval = sr->var_matval_max;
788 var_obj = sr->var_obj_max;
789 //var_lhs_offset = sr->opt_ub_var_offset;
790 }else{
791 n = &(sr->min_n);
792 obj_offset = &(sr->lb_offset);
793 sum = &(sr->sum_a_min);
794 obj_sum = &(sr->sum_c_min);
795 rhs = &(sr->rhs_min);
796 obj = sr->obj_min;
797 matind = sr->matind_min;
798 matval = sr->matval_min;
799 ratios = sr->ratio_min;
800 is_reversed = sr->reversed_min;
801 var_stat = sr->var_stat_min;
802 var_matval = sr->var_matval_min;
803 var_obj = sr->var_obj_min;
804 }
805
806 #if 0
807 sr->var_ub_lhs_offset[col_ind] = -rhs_ub_offset;
808 sr->var_lb_lhs_offset[col_ind] = -rhs_lb_offset;
809 sr->var_lb_obj_offset[col_ind] = obj_lb_offset;
810 sr->var_ub_obj_offset[col_ind] = obj_ub_offset;
811 #endif
812
813 if(ratio_type == 0){
814 obj[*n] = c_val;
815 matval[*n] = a_val;
816 matind[*n] = col_ind;
817 ratios[*n] = c_val/a_val;
818 if(obj_sense == SR_MAX){
819 *sum += (rhs_ub_offset - rhs_lb_offset);
820 *obj_sum += (obj_ub_offset - obj_ub_offset);
821 }else{
822 /* since bounds are converted to be
823 u - l >= x >= 0 */
824 *sum += 0.0;//rhs_lb_offset;
825 *obj_sum += 0.0;//obj_lb_offset;
826 }
827 (*n)++;
828 /* to solve bounds problem */
829 *rhs += -(rhs_lb_offset); /* conversion by x = y + l */
830 *obj_offset += obj_lb_offset;
831
832 }else if((ratio_type == 1 && obj_sense == SR_MAX) ||
833 (ratio_type == 2 && obj_sense == SR_MIN)){
834 *rhs += -rhs_ub_offset;
835 *obj_offset += obj_ub_offset;
836 var_stat[col_ind] = SR_VAR_FIXED_UB;
837 var_matval[col_ind] = a_val;
838 var_obj[col_ind] = c_val;
839 }else if((ratio_type == 1 && obj_sense == SR_MIN) ||
840 (ratio_type == 2 && obj_sense == SR_MAX)){
841 *rhs += -rhs_lb_offset;
842 *obj_offset += obj_lb_offset;
843 var_stat[col_ind] = SR_VAR_FIXED_LB;
844 var_matval[col_ind] = a_val;
845 var_obj[col_ind] = c_val;
846 }else{
847 obj[*n] = -c_val;
848 matval[*n] = -a_val;
849 matind[*n] = col_ind;
850 ratios[*n] = c_val/a_val;
851 is_reversed[*n] = TRUE;
852 if(obj_sense == SR_MAX){
853 *sum += -rhs_ub_offset + +rhs_lb_offset;
854 *obj_sum += -obj_ub_offset + rhs_lb_offset;
855 }else{
856 *sum += 0.0; //-rhs_lb_offset;
857 *obj_sum += 0.0; //-obj_lb_offset;
858 }
859 (*n)++;
860 /* to solve bounds problem */
861 *rhs += -(rhs_ub_offset); /* conversion by x = -y + u */
862 *obj_offset += obj_ub_offset;
863 }
864
865 return 0;
866 }
867 /*===========================================================================*/
868 /*===========================================================================*/
869
870 /* will modify the constraint to E and solve it */
871
sr_solve_open_prob(PREPdesc * P,SRdesc * sr,int obj_ind,int row_ind,int * r_matbeg,int * r_matind,double * r_matval,COLinfo * cols,double * ub,double * lb,double etol)872 int sr_solve_open_prob(PREPdesc *P, SRdesc *sr, int obj_ind,
873 int row_ind, int *r_matbeg,
874 int *r_matind, double *r_matval, COLinfo *cols,
875 double *ub, double *lb, double etol)
876 {
877
878 int l, k, col_ind;
879
880 double max_dual_ub = INF, min_dual_ub = INF;
881 double max_dual_lb = -INF, min_dual_lb = -INF;
882 double d_ratio, obj_val, a_val;
883
884 char no_upper = FALSE, no_lower = FALSE, is_free_column;
885 char is_fixed_column = FALSE;
886 char can_iterate = TRUE, prob_infeasible = FALSE, is_null_obj;
887
888 double *ub_offset = &(sr->ub_offset);
889 double *lb_offset = &(sr->lb_offset);
890 double rhs = sr->rhs;
891 char sense = sr->sense;
892
893 double obj_ub_offset;
894 double obj_lb_offset;
895
896
897 // sr->prob_type = OPEN_PROB;
898
899 for( k = r_matbeg[obj_ind], l = r_matbeg[row_ind];;){
900 if(k < r_matbeg[obj_ind + 1] &&
901 (r_matind[k] < r_matind[l] ||
902 l >= r_matbeg[row_ind + 1])){
903 if(r_matval[k] > 0.0){
904 if(!no_upper){
905 if(ub[r_matind[k]] >= INF){
906 no_upper = TRUE;
907 }else{
908 *ub_offset += ub[r_matind[k]] * r_matval[k];
909 }
910 }
911 if(!no_lower){
912 if(lb[r_matind[k]] <= -INF){
913 no_lower = TRUE;
914 }else{
915 *lb_offset += lb[r_matind[k]] * r_matval[k];
916 }
917 }
918 }else if (r_matval[k] < 0.0){
919 if(!no_lower){
920 if(ub[r_matind[k]] >= INF){
921 no_lower = TRUE;
922 }else{
923 *lb_offset += ub[r_matind[k]] * r_matval[k];
924 }
925 }
926 if(!no_upper){
927 if(lb[r_matind[k]] <= -INF){
928 no_upper = TRUE;
929 }else{
930 *ub_offset += lb[r_matind[k]] * r_matval[k];
931 }
932 }
933 }
934 k++;
935 }else{
936 if(l < r_matbeg[row_ind + 1] &&
937 (r_matind[k] > r_matind[l] ||
938 k >= r_matbeg[obj_ind+1])) {
939 is_null_obj = TRUE;
940 obj_val = 0.0;
941 }else{
942 is_null_obj = FALSE;
943 obj_val = r_matval[k];
944 }
945
946 /* first convert the column into 0 <= x <= inf */
947
948 a_val = r_matval[l];
949 col_ind = r_matind[l];
950 is_free_column = FALSE;
951 is_fixed_column = FALSE;
952 if(ub[col_ind] < INF && lb[col_ind] > -INF){
953 /* debug - get vars.type here */
954 if(ub[col_ind] > lb[col_ind] + etol){
955 /* debug */
956 printf("bounded column -case all open row-"
957 "sr_solve_open_prob(), exiting...\n");
958 return PREP_OTHER_ERROR;
959 }else{
960 /* fix column, take care of it here */
961 if(!is_null_obj){
962 obj_lb_offset = obj_val * lb[col_ind];
963 if(!no_upper){
964 *ub_offset += obj_lb_offset;
965 }
966 if(!no_lower){
967 *lb_offset += obj_lb_offset;
968 }
969 }
970 rhs += -(a_val *lb[col_ind]);
971 is_fixed_column = TRUE;
972 }
973 }else if(ub[col_ind] >= INF){
974 if(lb[col_ind] > -INF){
975 if(!is_null_obj){
976 obj_lb_offset = obj_val * lb[col_ind];
977 if(!no_upper){
978 *ub_offset += obj_lb_offset;
979 }
980 if(!no_lower){
981 *lb_offset += obj_lb_offset;
982 }
983 }
984 rhs += -(a_val * lb[col_ind]);
985 }else{
986 is_free_column = TRUE;
987 }
988 }else{
989 if(!is_null_obj){
990 obj_ub_offset = obj_val * ub[col_ind];
991 if(!no_upper){
992 *ub_offset += obj_ub_offset;
993 }
994 if(!no_lower){
995 *lb_offset += obj_ub_offset;
996 }
997 }
998 rhs += -(a_val * ub[col_ind]);
999 obj_val = -obj_val;
1000 a_val = -a_val;
1001 }
1002
1003 /* now get dual bounds */
1004 if(!is_fixed_column){
1005 if(a_val > 0.0 || a_val < 0.0){
1006 d_ratio = obj_val/a_val;
1007 if(a_val > 0.0){
1008 if(d_ratio < min_dual_ub){
1009 min_dual_ub = d_ratio;
1010 }
1011 if(-d_ratio < max_dual_ub){
1012 max_dual_ub = -d_ratio;
1013 }
1014 if(is_free_column){
1015 if(d_ratio > min_dual_lb){
1016 min_dual_lb = d_ratio;
1017 }
1018 if(-d_ratio > max_dual_lb){
1019 max_dual_lb = -d_ratio;
1020 }
1021 }
1022 }else{
1023 if(d_ratio > min_dual_lb){
1024 min_dual_lb = d_ratio;
1025 }
1026 if(-d_ratio > max_dual_lb){
1027 max_dual_lb = -d_ratio;
1028 }
1029
1030 if(is_free_column){
1031 if(d_ratio < min_dual_ub){
1032 min_dual_ub = d_ratio;
1033 }
1034 if(-d_ratio < max_dual_ub){
1035 max_dual_ub = -d_ratio;
1036 }
1037 }
1038 }
1039
1040 if(min_dual_lb > min_dual_ub){
1041 no_lower = TRUE;
1042 }
1043 if(max_dual_lb > max_dual_ub){
1044 no_upper = TRUE;
1045 /* debug */
1046 //printf("unbounded or infeasible problem?"
1047 // "-case all open row-"
1048 // "prep_solve_sr_rlx(), exiting...\n");
1049 //return PREP_OTHER_ERROR;
1050 }
1051 }else{
1052 /* debug */
1053 printf("not nonzero???"
1054 "numerical issues -case all open row-"
1055 "prep_solve_sr_rlx(), exiting...\n");
1056 return PREP_OTHER_ERROR;
1057 }
1058 }
1059 l++;
1060 if(!is_null_obj){
1061 k++;
1062 }
1063 }
1064
1065 if((no_upper && no_lower)){
1066 can_iterate = FALSE;
1067 break;
1068 }
1069 if(k == r_matbeg[obj_ind + 1] && l == r_matbeg[row_ind + 1]){
1070 break;
1071 }
1072 }
1073
1074 if(can_iterate){
1075 /* update the bounds for this row */
1076
1077 switch(sense){
1078 case 'L':
1079 if(max_dual_ub > 0.0){
1080 max_dual_ub = 0.0;
1081 }
1082 if(min_dual_ub > 0.0){
1083 min_dual_ub = 0.0;
1084 }
1085 break;
1086 case 'G':
1087 if(max_dual_lb < 0.0){
1088 max_dual_lb = 0.0;
1089 }
1090 if(min_dual_lb < 0.0){
1091 min_dual_lb = 0.0;
1092 }
1093 break;
1094 }
1095
1096 /* check again */
1097 // if(min_dual_lb > min_dual_ub ||/
1098 // max_dual_lb > max_dual_ub){/
1099 // printf("unbounded or infeasible problem?"
1100 // "-case all open row-"
1101 // "prep_solve_sr_rlx(), exiting...\n");
1102 // return PREP_OTHER_ERROR;
1103 // }
1104
1105 if(!no_lower){
1106 if(rhs >= 0){
1107 if(min_dual_ub < INF){
1108 sr->lb = min_dual_ub * rhs;
1109 }else{
1110 prob_infeasible = TRUE;
1111 }
1112 }else{
1113 if(min_dual_lb > -INF){
1114 sr->lb = min_dual_lb *rhs;
1115 }else{
1116 prob_infeasible = TRUE;
1117 }
1118 }
1119 if(!prob_infeasible){
1120 sr->lb += *lb_offset;
1121 sr->lb_updated = TRUE;
1122 }
1123
1124 }
1125
1126 if(!prob_infeasible){
1127 if(!no_upper){
1128 if(rhs >= 0){
1129 if(max_dual_ub < INF){
1130 sr->ub = -(max_dual_ub * rhs);
1131 }else{
1132 prob_infeasible = TRUE;
1133 }
1134 }else{
1135 if(max_dual_lb > -INF){
1136 sr->ub = -(max_dual_lb *rhs);
1137 }else{
1138 prob_infeasible = TRUE;
1139 }
1140 }
1141 if(!prob_infeasible){
1142 sr->ub += *ub_offset;
1143 sr->ub_updated = TRUE;
1144 }
1145
1146 }
1147 }
1148 }
1149
1150 // return(sr->lb_updated || sr->ub_updated);
1151 return(prob_infeasible);
1152 }
1153
1154 /*===========================================================================*/
1155 /*===========================================================================*/
free_sr_desc(SRdesc * sr)1156 void free_sr_desc(SRdesc *sr)
1157 {
1158 if(sr){
1159 FREE(sr->obj_max);
1160 FREE(sr->matval_max);
1161 FREE(sr->matind_max);
1162 FREE(sr->ratio_max);
1163
1164 FREE(sr->obj_min);
1165 FREE(sr->matval_min);
1166 FREE(sr->matind_min);
1167 FREE(sr->ratio_min);
1168
1169 FREE(sr->fixed_ind);
1170 FREE(sr->tmp_ind);
1171
1172 FREE(sr);
1173 }
1174 }
1175 /*===========================================================================*/
1176 /*===========================================================================*/
1177