1 /* $Id: CoinPresolveDupcol.cpp 2083 2019-01-06 19:38:09Z unxusr $ */
2 // Copyright (C) 2002, International Business Machines
3 // Corporation and others. All Rights Reserved.
4 // This code is licensed under the terms of the Eclipse Public License (EPL).
5
6 #include <stdio.h>
7 #include <math.h>
8
9 //#define PRESOLVE_DEBUG 1
10 // Debugging macros/functions
11 //#define PRESOLVE_DETAIL 1
12 #include "CoinPresolveMatrix.hpp"
13 #include "CoinPresolveFixed.hpp"
14 #include "CoinPresolveDupcol.hpp"
15 #include "CoinSort.hpp"
16 #include "CoinFinite.hpp"
17 #include "CoinHelperFunctions.hpp"
18 #include "CoinPresolveUseless.hpp"
19 #include "CoinMessage.hpp"
20 #if PRESOLVE_DEBUG || PRESOLVE_CONSISTENCY
21 #include "CoinPresolvePsdebug.hpp"
22 #endif
23
24 #define DSEED2 2147483647.0
25 // Can be used from anywhere
coin_init_random_vec(double * work,int n)26 void coin_init_random_vec(double *work, int n)
27 {
28 double deseed = 12345678.0;
29
30 for (int i = 0; i < n; ++i) {
31 deseed *= 16807.;
32 int jseed = static_cast< int >(deseed / DSEED2);
33 deseed -= static_cast< double >(jseed) * DSEED2;
34 double random = deseed / DSEED2;
35
36 work[i] = random;
37 }
38 }
39
40 namespace { // begin unnamed file-local namespace
41
42 /*
43 For each candidate major-dimension vector in majcands, calculate the sum
44 over the vector, with each minor dimension weighted by a random amount.
45 (E.g., calculate column sums with each row weighted by a random amount.)
46 The sums are returned in the corresponding entries of majsums.
47 */
48
compute_sums(int,const int * majlens,const CoinBigIndex * majstrts,int * minndxs,double * elems,const double * minmuls,int * majcands,double * majsums,int nlook)49 void compute_sums(int /*n*/, const int *majlens, const CoinBigIndex *majstrts,
50 int *minndxs, double *elems, const double *minmuls,
51 int *majcands, double *majsums, int nlook)
52
53 {
54 for (int cndx = 0; cndx < nlook; ++cndx) {
55 int i = majcands[cndx];
56 PRESOLVEASSERT(majlens[i] > 0);
57
58 CoinBigIndex kcs = majstrts[i];
59 CoinBigIndex kce = kcs + majlens[i];
60
61 double value = 0.0;
62
63 for (CoinBigIndex k = kcs; k < kce; k++) {
64 int irow = minndxs[k];
65 value += minmuls[irow] * elems[k];
66 }
67
68 majsums[cndx] = value;
69 }
70
71 return;
72 }
73
create_col(int col,int n,double * els,CoinBigIndex * mcstrt,double * colels,int * hrow,CoinBigIndex * link,CoinBigIndex * free_listp)74 void create_col(int col, int n, double *els,
75 CoinBigIndex *mcstrt, double *colels, int *hrow, CoinBigIndex *link,
76 CoinBigIndex *free_listp)
77 {
78 int *rows = reinterpret_cast< int * >(els + n);
79 CoinBigIndex free_list = *free_listp;
80 CoinBigIndex xstart = NO_LINK;
81 for (int i = 0; i < n; ++i) {
82 CoinBigIndex k = free_list;
83 assert(k >= 0);
84 free_list = link[free_list];
85 hrow[k] = rows[i];
86 colels[k] = els[i];
87 link[k] = xstart;
88 xstart = k;
89 }
90 mcstrt[col] = xstart;
91 *free_listp = free_list;
92 }
93
94 } // end unnamed file-local namespace
95
name() const96 const char *dupcol_action::name() const
97 {
98 return ("dupcol_action");
99 }
100
101 /*
102 Original comment: This is just ekkredc5, adapted into the new framework.
103 The datasets scorpion.mps and allgrade.mps have duplicate columns.
104
105 In case you don't have your OSL manual handy, a somewhat more informative
106 explanation: We're looking for an easy-to-detect special case of linearly
107 dependent columns, where the coefficients of the duplicate columns are
108 exactly equal. The idea for locating such columns is to generate pseudo-
109 random weights for each row and then calculate the weighted sum of
110 coefficients of each column. Columns with equal sums are checked more
111 thoroughly.
112
113 Analysis of the situation says there are two major cases:
114 * If the columns have equal objective coefficients, we can combine
115 them.
116 * If the columns have unequal objective coefficients, we may be able to
117 fix one at bound. If the required bound doesn't exist, we have dual
118 infeasibility (hence one of primal infeasibility or unboundedness).
119
120 In the comments below are a few fragments of code from the original
121 routine. I don't think they make sense, but I've left them for the nonce in
122 case someone else recognises the purpose. -- lh, 040909 --
123 */
124
125 const CoinPresolveAction
126 *
presolve(CoinPresolveMatrix * prob,const CoinPresolveAction * next)127 dupcol_action::presolve(CoinPresolveMatrix *prob,
128 const CoinPresolveAction *next)
129 {
130 #if PRESOLVE_DEBUG > 0 || PRESOLVE_CONSISTENCY > 0
131 #if PRESOLVE_DEBUG > 0
132 std::cout
133 << "Entering dupcol_action::presolve." << std::endl;
134 #endif
135 presolve_consistent(prob);
136 presolve_links_ok(prob);
137 presolve_check_sol(prob);
138 presolve_check_nbasic(prob);
139 #endif
140
141 #if PRESOLVE_DEBUG > 0 || COIN_PRESOLVE_TUNING > 0
142 int startEmptyRows = 0;
143 int startEmptyColumns = 0;
144 startEmptyRows = prob->countEmptyRows();
145 startEmptyColumns = prob->countEmptyCols();
146 #if COIN_PRESOLVE_TUNING > 0
147 double startTime = 0.0;
148 if (prob->tuning_)
149 startTime = CoinCpuTime();
150 #endif
151 #endif
152
153 double maxmin = prob->maxmin_;
154
155 double *colels = prob->colels_;
156 int *hrow = prob->hrow_;
157 CoinBigIndex *mcstrt = prob->mcstrt_;
158 int *hincol = prob->hincol_;
159 int ncols = prob->ncols_;
160 int nrows = prob->nrows_;
161
162 double *clo = prob->clo_;
163 double *cup = prob->cup_;
164 double *sol = prob->sol_;
165 double *rlo = prob->rlo_;
166 double *rup = prob->rup_;
167
168 // If all coefficients positive do more simply
169 bool allPositive = true;
170 double *rhs = prob->usefulRowDouble_; //new double[nrows];
171 CoinMemcpyN(rup, nrows, rhs);
172 /*
173 Scan the columns for candidates, and write the indices into sort. We're not
174 interested in columns that are empty, prohibited, or integral.
175
176 Question: Should we exclude singletons, which are useful in other transforms?
177 Question: Why are we excluding integral columns?
178 */
179 // allow integral columns if asked for
180 bool allowIntegers = ((prob->presolveOptions_ & 0x01) != 0);
181 int *sort = prob->usefulColumnInt_; //new int[ncols] ;
182 int nlook = 0;
183 for (int j = 0; j < ncols; j++) {
184 if (hincol[j] == 0)
185 continue;
186 // sort
187 CoinSort_2(hrow + mcstrt[j], hrow + mcstrt[j] + hincol[j],
188 colels + mcstrt[j]);
189 // check all positive and adjust rhs
190 if (allPositive) {
191 double lower = clo[j];
192 if (lower < cup[j]) {
193 for (CoinBigIndex k = mcstrt[j]; k < mcstrt[j] + hincol[j]; k++) {
194 double value = colels[k];
195 if (value < 0.0)
196 allPositive = false;
197 else
198 rhs[hrow[k]] -= lower * value;
199 }
200 } else {
201 for (CoinBigIndex k = mcstrt[j]; k < mcstrt[j] + hincol[j]; k++) {
202 double value = colels[k];
203 rhs[hrow[k]] -= lower * value;
204 }
205 }
206 }
207 if (prob->colProhibited2(j))
208 continue;
209 //#define PRESOLVE_INTEGER_DUPCOL
210 #ifndef PRESOLVE_INTEGER_DUPCOL
211 if (prob->isInteger(j) && !allowIntegers)
212 continue;
213 #endif
214 sort[nlook++] = j;
215 }
216 if (nlook == 0) { //delete[] sort ;
217 //delete [] rhs;
218 return (next);
219 }
220 /*
221 Prep: add the coefficients of each candidate column. To reduce false
222 positives, multiply each row by a `random' multiplier when forming the
223 sums. On return from compute_sums, sort and colsum are loaded with the
224 indices and column sums, respectively, of candidate columns. The pair of
225 arrays are then sorted by sum so that equal sums are adjacent.
226 */
227 double *colsum = prob->usefulColumnDouble_; //new double[ncols] ;
228 double *rowmul;
229 if (!prob->randomNumber_) {
230 rowmul = new double[nrows];
231 coin_init_random_vec(rowmul, nrows);
232 } else {
233 rowmul = prob->randomNumber_;
234 }
235 compute_sums(ncols, hincol, mcstrt, hrow, colels, rowmul, sort, colsum, nlook);
236 #define USE_LBS 0
237 #define SWAP_SIGNS 0
238 #if SWAP_SIGNS
239 int nPiece = 0;
240 // array to chain piecewise linear
241 int *piece = new int[ncols];
242 for (int i = 0; i < ncols; i++)
243 piece[i] = -1;
244 if (!allPositive) {
245 // swap signs
246 memcpy(sort + nlook, sort, nlook * sizeof(int));
247 for (int i = 0; i < nrows; i++)
248 rowmul[i] = -rowmul[i];
249 compute_sums(ncols, hincol, mcstrt, hrow, colels, rowmul, sort + nlook, colsum + nlook, nlook);
250 nlook += nlook;
251 }
252 #endif
253 CoinSort_2(colsum, colsum + nlook, sort);
254 /*
255 General prep --- unpack the various vectors we'll need, and allocate arrays
256 to record the results.
257 */
258 presolvehlink *clink = prob->clink_;
259
260 double *rowels = prob->rowels_;
261 int *hcol = prob->hcol_;
262 const CoinBigIndex *mrstrt = prob->mrstrt_;
263 int *hinrow = prob->hinrow_;
264
265 double *dcost = prob->cost_;
266
267 action *actions = new action[nlook];
268 int nactions = 0;
269 #ifdef ZEROFAULT
270 memset(actions, 0, nlook * sizeof(action));
271 #endif
272
273 int *fixed_down = new int[nlook];
274 int nfixed_down = 0;
275 int *fixed_up = new int[nlook];
276 int nfixed_up = 0;
277
278 #if 0
279 // Excluded in the original routine. I'm guessing it's excluded because
280 // it's just not cost effective to worry about this. -- lh, 040908 --
281
282 // It may be the case that several columns are duplicate.
283 // If not all have the same cost, then we have to make sure
284 // that we set the most expensive one to its minimum
285 // now sort in each class by cost
286 {
287 double dval = colsum[0] ;
288 int first = 0 ;
289 for (int jj = 1; jj < nlook; jj++) {
290 while (colsum[jj]==dval)
291 jj++ ;
292
293 if (first + 1 < jj) {
294 double buf[jj - first] ;
295 for (int i=first; i<jj; ++i)
296 buf[i-first] = dcost[sort[i]]*maxmin ;
297
298 CoinSort_2(buf,buf+jj-first,sort+first) ;
299 //ekk_sortonDouble(buf,&sort[first],jj-first) ;
300 }
301 }
302 }
303 #endif
304 // We will get all min/max but only if needed
305 bool gotStuff = false;
306 /*
307 Original comment: It appears to be the case that this loop is finished,
308 there may still be duplicate cols left. I haven't done anything
309 about that yet.
310
311 Open the main loop to compare column pairs. We'll compare sort[jj] to
312 sort[tgt]. This allows us to accumulate multiple columns into one. But
313 we don't manage all-pairs comparison when we can't combine columns.
314
315 We can quickly dismiss pairs which have unequal sums or lengths.
316 */
317 int isorted = -1;
318 int tgt = 0;
319 for (int jj = 1; jj < nlook; jj++) {
320 if (colsum[jj] != colsum[jj - 1]) {
321 tgt = jj; // Must update before continuing
322 continue;
323 }
324
325 int j2 = sort[jj];
326 int j1 = sort[tgt];
327 int len2 = hincol[j2];
328 int len1 = hincol[j1];
329
330 if (len2 != len1) {
331 tgt = jj;
332 continue;
333 }
334 /*
335 The final test: sort the columns by row index and compare each index and
336 coefficient.
337 */
338 CoinBigIndex kcs = mcstrt[j2];
339 CoinBigIndex kce = kcs + hincol[j2];
340 CoinBigIndex ishift = mcstrt[j1] - kcs;
341
342 if (len1 > 1 && isorted < j1) {
343 CoinSort_2(hrow + mcstrt[j1], hrow + mcstrt[j1] + len1,
344 colels + mcstrt[j1]);
345 isorted = j1;
346 }
347 if (len2 > 1 && isorted < j2) {
348 CoinSort_2(hrow + kcs, hrow + kcs + len2, colels + kcs);
349 isorted = j2;
350 }
351
352 CoinBigIndex k;
353 for (k = kcs; k < kce; k++) {
354 if (hrow[k] != hrow[k + ishift] || colels[k] != colels[k + ishift]) {
355 break;
356 }
357 }
358 #if SWAP_SIGNS
359 if (k != kce) {
360 if (allPositive) {
361 tgt = jj;
362 continue;
363 } else {
364 // try negative
365 for (k = kcs; k < kce; k++) {
366 if (hrow[k] != hrow[k + ishift] || colels[k] != -colels[k + ishift]) {
367 break;
368 }
369 }
370 if (k == kce) {
371 nPiece++;
372 if (piece[j1] < 0 && piece[j2] < 0) {
373 piece[j1] = j2;
374 piece[j2] = j1;
375 } else if (piece[j1] < 0) {
376 int j3 = piece[j2];
377 piece[j2] = j1;
378 piece[j1] = j3;
379 }
380 }
381 tgt = jj;
382 continue;
383 }
384 }
385 #else
386 if (k != kce) {
387 tgt = jj;
388 continue;
389 }
390 #endif
391 // check both continuous or both integer
392 if (prob->isInteger(j1) != prob->isInteger(j2))
393 continue;
394 /*
395 These really are duplicate columns. Grab values for convenient reference.
396 Convert the objective coefficients for minimization.
397 */
398 double clo1 = clo[j1];
399 double cup1 = cup[j1];
400 double clo2 = clo[j2];
401 double cup2 = cup[j2];
402 double c1 = dcost[j1] * maxmin;
403 double c2 = dcost[j2] * maxmin;
404 PRESOLVEASSERT(!(clo1 == cup1 || clo2 == cup2));
405 // Get reasonable bounds on sum of two variables
406 double lowerBound = -COIN_DBL_MAX;
407 double upperBound = COIN_DBL_MAX;
408 #if USE_LBS == 0
409 // For now only if lower bounds are zero
410 bool takeColumn = (!clo1 && !clo2);
411 #else
412 bool takeColumn = true;
413 if (clo1 || clo2)
414 printf("DUPCOL %d and %d have nonzero lbs %g %g\n",
415 j1, j2, clo1, clo2);
416 #endif
417 if (takeColumn) {
418 // Only need bounds if c1 != c2
419 if (c1 != c2) {
420 if (!allPositive) {
421 #if 0
422
423 for (k=kcs;k<kce;k++) {
424 int iRow = hrow[k];
425 bool posinf = false;
426 bool neginf = false;
427 double maxup = 0.0;
428 double maxdown = 0.0;
429
430 // compute sum of all bounds except for j1,j2
431 CoinBigIndex kk;
432 CoinBigIndex kre = mrstrt[iRow]+hinrow[iRow];
433 double value1=0.0;
434 for (kk=mrstrt[iRow]; kk<kre; kk++) {
435 int col = hcol[kk];
436 if (col == j1||col==j2) {
437 value1=rowels[kk];
438 continue;
439 }
440 double coeff = rowels[kk];
441 double lb = clo[col];
442 double ub = cup[col];
443
444 if (coeff > 0.0) {
445 if (PRESOLVE_INF <= ub) {
446 posinf = true;
447 if (neginf)
448 break; // pointless
449 } else {
450 maxup += ub * coeff;
451 }
452 if (lb <= -PRESOLVE_INF) {
453 neginf = true;
454 if (posinf)
455 break; // pointless
456 } else {
457 maxdown += lb * coeff;
458 }
459 } else {
460 if (PRESOLVE_INF <= ub) {
461 neginf = true;
462 if (posinf)
463 break; // pointless
464 } else {
465 maxdown += ub * coeff;
466 }
467 if (lb <= -PRESOLVE_INF) {
468 posinf = true;
469 if (neginf)
470 break; // pointless
471 } else {
472 maxup += lb * coeff;
473 }
474 }
475 }
476
477 if (kk==kre) {
478 assert (value1);
479 if (value1>1.0e-5) {
480 if (!neginf&&rup[iRow]<1.0e10)
481 if (upperBound*value1>rup[iRow]-maxdown)
482 upperBound = (rup[iRow]-maxdown)/value1;
483 if (!posinf&&rlo[iRow]>-1.0e10)
484 if (lowerBound*value1<rlo[iRow]-maxup)
485 lowerBound = (rlo[iRow]-maxup)/value1;
486 } else if (value1<-1.0e-5) {
487 if (!neginf&&rup[iRow]<1.0e10)
488 if (lowerBound*value1>rup[iRow]-maxdown) {
489 #ifndef NDEBUG
490 double x=lowerBound;
491 #endif
492 lowerBound = (rup[iRow]-maxdown)/value1;
493 assert (lowerBound == CoinMax(x,(rup[iRow]-maxdown)/value1));
494 }
495 if (!posinf&&rlo[iRow]>-1.0e10)
496 if (upperBound*value1<rlo[iRow]-maxup) {
497 #ifndef NDEBUG
498 double x=upperBound;
499 #endif
500 upperBound = (rlo[iRow]-maxup)/value1;
501 assert(upperBound == CoinMin(x,(rlo[iRow]-maxup)/value1));
502 }
503 }
504 }
505 }
506 double l=lowerBound;
507 double u=upperBound;
508 #endif
509 if (!gotStuff) {
510 prob->recomputeSums(-1); // get min max
511 gotStuff = true;
512 }
513 int positiveInf = 0;
514 int negativeInf = 0;
515 double lo = 0;
516 double up = 0.0;
517 if (clo1 < -PRESOLVE_INF)
518 negativeInf++;
519 else
520 lo += clo1;
521 if (clo2 < -PRESOLVE_INF)
522 negativeInf++;
523 else
524 lo += clo2;
525 if (cup1 > PRESOLVE_INF)
526 positiveInf++;
527 else
528 up += cup1;
529 if (cup2 > PRESOLVE_INF)
530 positiveInf++;
531 else
532 up += cup2;
533 for (k = kcs; k < kce; k++) {
534 int iRow = hrow[k];
535 double value = colels[k];
536 int pInf = (value > 0.0) ? positiveInf : negativeInf;
537 int nInf = (value > 0.0) ? negativeInf : positiveInf;
538 int posinf = prob->infiniteUp_[iRow] - pInf;
539 int neginf = prob->infiniteDown_[iRow] - nInf;
540 if (posinf > 0 && neginf > 0)
541 continue; // this row can't bound
542 double maxup = prob->sumUp_[iRow];
543 double maxdown = prob->sumDown_[iRow];
544
545 if (value > 0.0) {
546 maxdown -= value * lo;
547 maxup -= value * up;
548 } else {
549 maxdown -= value * up;
550 maxup -= value * lo;
551 }
552 if (value > 1.0e-5) {
553 if (!neginf && rup[iRow] < 1.0e10)
554 if (upperBound * value > rup[iRow] - maxdown)
555 upperBound = (rup[iRow] - maxdown) / value;
556 if (!posinf && rlo[iRow] > -1.0e10)
557 if (lowerBound * value < rlo[iRow] - maxup)
558 lowerBound = (rlo[iRow] - maxup) / value;
559 } else if (value < -1.0e-5) {
560 if (!neginf && rup[iRow] < 1.0e10)
561 if (lowerBound * value > rup[iRow] - maxdown) {
562 lowerBound = (rup[iRow] - maxdown) / value;
563 }
564 if (!posinf && rlo[iRow] > -1.0e10)
565 if (upperBound * value < rlo[iRow] - maxup) {
566 upperBound = (rlo[iRow] - maxup) / value;
567 }
568 }
569 }
570 //assert (fabs(l-lowerBound)<1.0e-5&&fabs(u-upperBound)<1.0e-5);
571 } else {
572 // can do faster
573 lowerBound = 0.0;
574 for (k = kcs; k < kce; k++) {
575 int iRow = hrow[k];
576 double value = colels[k];
577 if (upperBound * value > rhs[iRow])
578 upperBound = rhs[iRow] / value;
579 }
580 }
581 }
582 // relax a bit
583 upperBound -= 1.0e-9;
584 } else {
585 // Not sure what to do so give up
586 continue;
587 }
588 /*
589 There are two main cases: The objective coefficients are equal or unequal.
590
591 For equal objective coefficients c1 == c2, we can combine the columns by
592 making the substitution x<j1> = x'<j1> - x<j2>. This will eliminate column
593 sort[jj] = j2 and leave the new variable x' in column sort[tgt] = j1. tgt
594 doesn't move.
595 */
596 if (c1 == c2) {
597 #ifdef PRESOLVE_INTEGER_DUPCOL
598 if (!allowIntegers) {
599 if (prob->isInteger(j1)) {
600 if (!prob->isInteger(j2)) {
601 if (cup2 < upperBound) //if (!prob->colInfinite(j2))
602 continue;
603 else
604 cup2 = COIN_DBL_MAX;
605 }
606 } else if (prob->isInteger(j2)) {
607 if (cup1 < upperBound) //if (!prob->colInfinite(j1))
608 continue;
609 else
610 cup1 = COIN_DBL_MAX;
611 }
612 //printf("TakingINTeq\n");
613 }
614 #endif
615 /*
616 As far as the presolved lp, there's no difference between columns. But we
617 need this relation to hold in order to guarantee that we can split the
618 value of the combined column during postsolve without damaging the basis.
619 (The relevant case is when the combined column is basic --- we need to be
620 able to retain column j1 in the basis and make column j2 nonbasic.)
621 */
622 if (!(clo2 + cup1 <= clo1 + cup2)) {
623 CoinSwap(j1, j2);
624 CoinSwap(clo1, clo2);
625 CoinSwap(cup1, cup2);
626 tgt = jj;
627 }
628 /*
629 Create the postsolve action before we start to modify the columns.
630 */
631 #if PRESOLVE_DEBUG > 1
632 PRESOLVE_STMT(printf("DUPCOL: (%d,%d) %d += %d\n", j1, j2, j1, j2));
633 PRESOLVE_DETAIL_PRINT(printf("pre_dupcol %dC %dC E\n", j2, j1));
634 #endif
635
636 action *s = &actions[nactions++];
637 s->thislo = clo[j2];
638 s->thisup = cup[j2];
639 s->lastlo = clo[j1];
640 s->lastup = cup[j1];
641 s->ithis = j2;
642 s->ilast = j1;
643 s->nincol = hincol[j2];
644 s->colels = presolve_dupmajor(colels, hrow, hincol[j2], mcstrt[j2]);
645 /*
646 Combine the columns into column j1. Upper and lower bounds and solution
647 simply add, and the coefficients are unchanged.
648
649 I'm skeptical of pushing a bound to infinity like this, but leave it for now.
650 -- lh, 040908 --
651 */
652 clo1 += clo2;
653 if (clo1 < -1.0e20) {
654 clo1 = -PRESOLVE_INF;
655 }
656 clo[j1] = clo1;
657 cup1 += cup2;
658 if (cup1 > 1.0e20) {
659 cup1 = PRESOLVE_INF;
660 }
661 cup[j1] = cup1;
662 if (sol) {
663 sol[j1] += sol[j2];
664 }
665 if (prob->colstat_) {
666 if (prob->getColumnStatus(j1) == CoinPrePostsolveMatrix::basic || prob->getColumnStatus(j2) == CoinPrePostsolveMatrix::basic) {
667 prob->setColumnStatus(j1, CoinPrePostsolveMatrix::basic);
668 }
669 }
670 /*
671 Empty column j2.
672 */
673 dcost[j2] = 0.0;
674 if (sol) {
675 sol[j2] = clo2;
676 }
677 CoinBigIndex k2cs = mcstrt[j2];
678 CoinBigIndex k2ce = k2cs + hincol[j2];
679 for (CoinBigIndex k = k2cs; k < k2ce; ++k) {
680 presolve_delete_from_row(hrow[k], j2, mrstrt, hinrow, hcol, rowels);
681 }
682 hincol[j2] = 0;
683 PRESOLVE_REMOVE_LINK(clink, j2);
684 continue;
685 }
686 /*
687 Unequal reduced costs. In this case, we may be able to fix one of the columns
688 or prove dual infeasibility. Given column a_k, duals y, objective
689 coefficient c_k, the reduced cost cbar_k = c_k - dot(y,a_k). Given
690 a_1 = a_2, substitute for dot(y,a_1) in the relation for cbar_2 to get
691 cbar_2 = (c_2 - c_1) + cbar_1
692 Independent elements here are variable bounds l_k, u_k, and difference
693 (c_2 - c_1). Infinite bounds for l_k, u_k will constrain the sign of cbar_k.
694 Assume minimization. If you do the case analysis, you find these cases of
695 interest:
696
697 l_1 u_1 l_2 u_2 cbar_1 c_2-c_1 cbar_2 result
698
699 A any finite -inf any <= 0 > 0 <= 0 x_1 -> NBUB
700 B -inf any any finite <= 0 < 0 < 0 x_2 -> NBUB
701
702 C finite any any +inf >= 0 < 0 >= 0 x_1 -> NBLB
703 D any +inf finite any >= 0 > 0 >= 0 x_2 -> NBLB
704
705 E -inf any any +inf <= 0 < 0 >= 0 dual infeas
706 F any inf -inf any >= 0 > 0 <= 0 dual infeas
707
708 G any finite finite any > 0 no inference
709 H finite any any finite < 0 no inference
710
711 The cases labelled dual infeasible are primal unbounded.
712
713 To keep the code compact, we'll always aim to take x_2 to bound. In the cases
714 where x_1 should go to bound, we'll swap. The implementation is boolean
715 algebra. Define bits for infinite bounds and (c_2 > c_1), then look for the
716 correct patterns.
717 */
718 else {
719 int minterm = 0;
720 #ifdef PRESOLVE_INTEGER_DUPCOL
721 if (!allowIntegers) {
722 if (c2 > c1) {
723 if (cup1 < upperBound /*!prob->colInfinite(j1)*/ && (prob->isInteger(j1) || prob->isInteger(j2)))
724 continue;
725 } else {
726 if (cup2 < upperBound /*!prob->colInfinite(j2)*/ && (prob->isInteger(j1) || prob->isInteger(j2)))
727 continue;
728 }
729 //printf("TakingINTne\n");
730 }
731 #endif
732 bool swapped = false;
733 #if PRESOLVE_DEBUG > 1
734 printf("bounds %g %g\n", lowerBound, upperBound);
735 #endif
736 if (c2 > c1)
737 minterm |= 1 << 0;
738 if (cup2 >= PRESOLVE_INF /*prob->colInfinite(j2)*/)
739 minterm |= 1 << 1;
740 if (clo2 <= -PRESOLVE_INF)
741 minterm |= 1 << 2;
742 if (cup1 >= PRESOLVE_INF /*prob->colInfinite(j1)*/)
743 minterm |= 1 << 3;
744 if (clo1 <= -PRESOLVE_INF)
745 minterm |= 1 << 4;
746 // for now be careful - just one special case
747 if (!clo1 && !clo2) {
748 if (c2 > c1 && cup1 >= upperBound)
749 minterm |= 1 << 3;
750 else if (c2 < c1 && cup2 >= upperBound)
751 minterm |= 1 << 1;
752 }
753 /*
754 The most common case in a well-formed system should be no inference. We're
755 looking for x00x1 (case G) and 0xx00 (case H). This is where we have the
756 potential to miss inferences: If there are three or more columns with the
757 same sum, sort[tgt] == j1 will only be compared to the second in the
758 group.
759 */
760 if ((minterm & 0x0d) == 0x1 || (minterm & 0x13) == 0) {
761 tgt = jj;
762 continue;
763 }
764 /*
765 Next remove the unbounded cases, 1xx10 and x11x1.
766 */
767 if ((minterm & 0x13) == 0x12 || (minterm & 0x0d) == 0x0d) {
768 prob->setStatus(2);
769 #if PRESOLVE_DEBUG > 1
770 PRESOLVE_STMT(printf("DUPCOL: (%d,%d) Unbounded\n", j1, j2));
771 #endif
772 break;
773 }
774 /*
775 With the no inference and unbounded cases removed, all that's left are the
776 cases where we can push a variable to bound. Swap if necessary (x01x1 or
777 0xx10) so that we're always fixing index j2. This means that column
778 sort[tgt] = j1 will be fixed. Unswapped, we fix column sort[jj] = j2.
779 */
780 if ((minterm & 0x0d) == 0x05 || (minterm & 0x13) == 0x02) {
781 CoinSwap(j1, j2);
782 CoinSwap(clo1, clo2);
783 CoinSwap(cup1, cup2);
784 CoinSwap(c1, c2);
785 int tmp1 = minterm & 0x18;
786 int tmp2 = minterm & 0x06;
787 int tmp3 = minterm & 0x01;
788 minterm = (tmp1 >> 2) | (tmp2 << 2) | (tmp3 ^ 0x01);
789 swapped = true;
790 }
791 /*
792 Force x_2 to upper bound? (Case B, boolean 1X100, where X == don't care.)
793 */
794 if ((minterm & 0x13) == 0x10) {
795 fixed_up[nfixed_up++] = j2;
796 #if PRESOLVE_DEBUG > 1
797 PRESOLVE_STMT(printf("DUPCOL: (%d,%d) %d -> NBUB\n", j1, j2, j2));
798 #endif
799 if (prob->colstat_) {
800 if (prob->getColumnStatus(j1) == CoinPrePostsolveMatrix::basic || prob->getColumnStatus(j2) == CoinPrePostsolveMatrix::basic) {
801 prob->setColumnStatus(j1, CoinPrePostsolveMatrix::basic);
802 }
803 prob->setColumnStatus(j2, CoinPrePostsolveMatrix::atUpperBound);
804 }
805 if (sol) {
806 double delta2 = cup2 - sol[j2];
807 sol[j2] = cup2;
808 sol[j1] -= delta2;
809 }
810 if (swapped) {
811 tgt = jj;
812 }
813 continue;
814 }
815 /*
816 Force x_2 to lower bound? (Case C, boolean X1011.)
817 */
818 if ((minterm & 0x0d) == 0x09) {
819 fixed_down[nfixed_down++] = j2;
820 #if PRESOLVE_DEBUG > 1
821 PRESOLVE_STMT(printf("DUPCOL: (%d,%d) %d -> NBLB\n", j1, j2, j2));
822 #endif
823 if (prob->colstat_) {
824 if (prob->getColumnStatus(j1) == CoinPrePostsolveMatrix::basic || prob->getColumnStatus(j2) == CoinPrePostsolveMatrix::basic) {
825 prob->setColumnStatus(j1, CoinPrePostsolveMatrix::basic);
826 }
827 prob->setColumnStatus(j2, CoinPrePostsolveMatrix::atLowerBound);
828 }
829 if (sol) {
830 double delta2 = clo2 - sol[j2];
831 sol[j2] = clo2;
832 sol[j1] -= delta2;
833 }
834 if (swapped) {
835 tgt = jj;
836 }
837 continue;
838 }
839 }
840 /*
841 We should never reach this point in the loop --- all cases force a new
842 iteration or loop termination. If we get here, something happened that we
843 didn't anticipate.
844 */
845 PRESOLVE_STMT(printf("DUPCOL: (%d,%d) UNEXPECTED!\n", j1, j2));
846 }
847 /*
848 What's left? Deallocate vectors, and call make_fixed_action to handle any
849 variables that were fixed to bound.
850 */
851 if (rowmul != prob->randomNumber_)
852 delete[] rowmul;
853 #if SWAP_SIGNS
854 if (nPiece) {
855 nPiece = 0;
856 int nTotal = 0;
857 for (int i = 0; i < ncols; i++) {
858 if ((piece[i] & 0x80000000) == 0) {
859 int j2 = i;
860 int number = 1;
861 while (piece[j2] != i) {
862 int jNext = piece[j2];
863 assert(jNext != -1);
864 if (jNext < i) {
865 // already done
866 number = 0;
867 break;
868 }
869 if ((jNext & 0x80000000) != 0 || clo[jNext] == cup[jNext]) {
870 // no good
871 piece[j2] |= 0x80000000;
872 jNext &= 0x7fffffff;
873 while (jNext != j2) {
874 piece[jNext] |= 0x80000000;
875 jNext = piece[jNext] & 0x7fffffff;
876 }
877 piece[j2] |= 0x80000000;
878 number = -1;
879 break;
880 } else {
881 j2 = jNext;
882 number++;
883 }
884 }
885 if (number > 0) {
886 nPiece++;
887 nTotal += number;
888 }
889 }
890 }
891 printf("%d odd duplicates in %d groups\n", nTotal, nPiece);
892 // now build real structures
893 // to test fix and adjust
894 nPiece = 0;
895 nTotal = 0;
896 for (int i = 0; i < ncols; i++) {
897 if ((piece[i] & 0x80000000) == 0) {
898 int number = 1;
899 double lo = CoinMax(clo[i], -1.0e100);
900 double up = CoinMin(cup[i], 1.0e100);
901 // get first
902 double value = colels[mcstrt[i]];
903 int iNext = piece[i];
904 piece[i] |= 0x80000000;
905 while (iNext != i) {
906 number++;
907 if (value == colels[mcstrt[iNext]]) {
908 lo += CoinMax(clo[iNext], -1.0e100);
909 up += CoinMin(cup[iNext], 1.0e100);
910 } else {
911 lo -= CoinMin(cup[iNext], 1.0e100);
912 up -= CoinMax(clo[iNext], -1.0e100);
913 }
914 clo[iNext] = 0.0;
915 cup[iNext] = 0.0;
916 fixed_down[nfixed_down++] = iNext;
917 piece[iNext] |= 0x80000000;
918 iNext = piece[iNext] & 0x7fffffff;
919 }
920 if (lo < -1.0e50)
921 lo = -COIN_DBL_MAX;
922 if (up > 1.0e50)
923 up = COIN_DBL_MAX;
924 if (lo == -COIN_DBL_MAX && up == COIN_DBL_MAX) {
925 printf("Help - free variable %d\n", i);
926 }
927 nPiece++;
928 nTotal += number;
929 }
930 }
931 printf("Pass two %d odd duplicates in %d groups\n", nTotal, nPiece);
932 }
933 delete[] piece;
934 #endif
935 //delete[] colsum ;
936 //delete[] sort ;
937 //delete [] rhs;
938
939 #if PRESOLVE_SUMMARY || PRESOLVE_DEBUG
940 if (nactions + nfixed_down + nfixed_up > 0) {
941 printf("DUPLICATE COLS: %d combined, %d lb, %d ub\n",
942 nactions, nfixed_down, nfixed_up);
943 }
944 #endif
945 if (nactions) {
946 next = new dupcol_action(nactions, CoinCopyOfArray(actions, nactions), next);
947 // we can't go round again in integer
948 prob->presolveOptions_ |= 0x80000000;
949 }
950 deleteAction(actions, action *);
951
952 if (nfixed_down) {
953 next = make_fixed_action::presolve(prob, fixed_down, nfixed_down, true, next);
954 }
955 delete[] fixed_down;
956
957 if (nfixed_up) {
958 next = make_fixed_action::presolve(prob, fixed_up, nfixed_up, false, next);
959 }
960 delete[] fixed_up;
961
962 #if COIN_PRESOLVE_TUNING > 0
963 double thisTime = 0.0;
964 if (prob->tuning_)
965 thisTime = CoinCpuTime();
966 #endif
967 #if PRESOLVE_CONSISTENCY > 0 || PRESOLVE_DEBUG > 0
968 presolve_check_sol(prob);
969 #endif
970 #if PRESOLVE_DEBUG > 0 || COIN_PRESOLVE_TUNING > 0
971 int droppedRows = prob->countEmptyRows() - startEmptyRows;
972 int droppedColumns = prob->countEmptyCols() - startEmptyColumns;
973 std::cout
974 << "Leaving dupcol_action::presolve, "
975 << droppedRows << " rows, " << droppedColumns << " columns dropped";
976 #if COIN_PRESOLVE_TUNING > 0
977 std::cout << " in " << thisTime - startTime << "s";
978 #endif
979 std::cout << "." << std::endl;
980 #endif
981
982 return (next);
983 }
984
postsolve(CoinPostsolveMatrix * prob) const985 void dupcol_action::postsolve(CoinPostsolveMatrix *prob) const
986 {
987 const action *const actions = actions_;
988 const int nactions = nactions_;
989
990 double *clo = prob->clo_;
991 double *cup = prob->cup_;
992
993 double *sol = prob->sol_;
994 double *dcost = prob->cost_;
995
996 double *colels = prob->colels_;
997 int *hrow = prob->hrow_;
998 CoinBigIndex *mcstrt = prob->mcstrt_;
999 int *hincol = prob->hincol_;
1000 CoinBigIndex *link = prob->link_;
1001
1002 double *rcosts = prob->rcosts_;
1003 double tolerance = prob->ztolzb_;
1004
1005 for (const action *f = &actions[nactions - 1]; actions <= f; f--) {
1006 int icol = f->ithis; // was fixed
1007 int icol2 = f->ilast; // was kept
1008
1009 dcost[icol] = dcost[icol2];
1010 clo[icol] = f->thislo;
1011 cup[icol] = f->thisup;
1012 clo[icol2] = f->lastlo;
1013 cup[icol2] = f->lastup;
1014
1015 create_col(icol, f->nincol, f->colels, mcstrt, colels, hrow, link,
1016 &prob->free_list_);
1017 #if PRESOLVE_CONSISTENCY
1018 presolve_check_free_list(prob);
1019 #endif
1020 // hincol[icol] = hincol[icol2]; // right? - no - has to match number in create_col
1021 hincol[icol] = f->nincol;
1022
1023 double l_j = f->thislo;
1024 double u_j = f->thisup;
1025 double l_k = f->lastlo;
1026 double u_k = f->lastup;
1027 double x_k_sol = sol[icol2];
1028 PRESOLVE_DETAIL_PRINT(printf("post icol %d %g %g %g icol2 %d %g %g %g\n",
1029 icol, clo[icol], sol[icol], cup[icol],
1030 icol2, clo[icol2], sol[icol2], cup[icol2]));
1031 if (l_j > -PRESOLVE_INF && x_k_sol - l_j >= l_k - tolerance && x_k_sol - l_j <= u_k + tolerance) {
1032 // j at lb, leave k
1033 prob->setColumnStatus(icol, CoinPrePostsolveMatrix::atLowerBound);
1034 sol[icol] = l_j;
1035 sol[icol2] = x_k_sol - sol[icol];
1036 } else if (u_j < PRESOLVE_INF && x_k_sol - u_j >= l_k - tolerance && x_k_sol - u_j <= u_k + tolerance) {
1037 // j at ub, leave k
1038 prob->setColumnStatus(icol, CoinPrePostsolveMatrix::atUpperBound);
1039 sol[icol] = u_j;
1040 sol[icol2] = x_k_sol - sol[icol];
1041 } else if (l_k > -PRESOLVE_INF && x_k_sol - l_k >= l_j - tolerance && x_k_sol - l_k <= u_j + tolerance) {
1042 // k at lb make j basic
1043 prob->setColumnStatus(icol, prob->getColumnStatus(icol2));
1044 sol[icol2] = l_k;
1045 sol[icol] = x_k_sol - l_k;
1046 prob->setColumnStatus(icol2, CoinPrePostsolveMatrix::atLowerBound);
1047 } else if (u_k < PRESOLVE_INF && x_k_sol - u_k >= l_j - tolerance && x_k_sol - u_k <= u_j + tolerance) {
1048 // k at ub make j basic
1049 prob->setColumnStatus(icol, prob->getColumnStatus(icol2));
1050 sol[icol2] = u_k;
1051 sol[icol] = x_k_sol - u_k;
1052 prob->setColumnStatus(icol2, CoinPrePostsolveMatrix::atUpperBound);
1053 } else {
1054 // both free! superbasic time
1055 sol[icol] = 0.0; // doesn't matter
1056 prob->setColumnStatus(icol, CoinPrePostsolveMatrix::isFree);
1057 }
1058 PRESOLVE_DETAIL_PRINT(printf("post2 icol %d %g icol2 %d %g\n",
1059 icol, sol[icol],
1060 icol2, sol[icol2]));
1061 // row activity doesn't change
1062 // dj of both variables is the same
1063 rcosts[icol] = rcosts[icol2];
1064 // leave until destructor
1065 // deleteAction(f->colels,double *);
1066
1067 #if PRESOLVE_DEBUG > 0
1068 const double ztolzb = prob->ztolzb_;
1069 if (!(clo[icol] - ztolzb <= sol[icol] && sol[icol] <= cup[icol] + ztolzb))
1070 printf("BAD DUPCOL BOUNDS: %g %g %g\n", clo[icol], sol[icol], cup[icol]);
1071 if (!(clo[icol2] - ztolzb <= sol[icol2] && sol[icol2] <= cup[icol2] + ztolzb))
1072 printf("BAD DUPCOL BOUNDS: %g %g %g\n", clo[icol2], sol[icol2], cup[icol2]);
1073 #endif
1074 }
1075 // leave until destructor
1076 // deleteAction(actions_,action *);
1077 }
1078
~dupcol_action()1079 dupcol_action::~dupcol_action()
1080 {
1081 for (int i = nactions_ - 1; i >= 0; --i) {
1082 deleteAction(actions_[i].colels, double *);
1083 }
1084 deleteAction(actions_, action *);
1085 }
1086
1087 /*
1088 Routines for duplicate rows. This is definitely unfinished --- there's no
1089 postsolve action.
1090 */
1091
name() const1092 const char *duprow_action::name() const
1093 {
1094 return ("duprow_action");
1095 }
1096
1097 // This is just ekkredc4, adapted into the new framework.
1098 /*
1099 I've made minimal changes for compatibility with dupcol: An initial scan to
1100 accumulate rows of interest in sort.
1101 -- lh, 040909 --
1102 */
1103 const CoinPresolveAction
1104 *
presolve(CoinPresolveMatrix * prob,const CoinPresolveAction * next)1105 duprow_action::presolve(CoinPresolveMatrix *prob,
1106 const CoinPresolveAction *next)
1107 {
1108 double startTime = 0.0;
1109 int startEmptyRows = 0;
1110 int startEmptyColumns = 0;
1111 if (prob->tuning_) {
1112 startTime = CoinCpuTime();
1113 startEmptyRows = prob->countEmptyRows();
1114 startEmptyColumns = prob->countEmptyCols();
1115 }
1116 double *rowels = prob->rowels_;
1117 int *hcol = prob->hcol_;
1118 CoinBigIndex *mrstrt = prob->mrstrt_;
1119 int *hinrow = prob->hinrow_;
1120 int ncols = prob->ncols_;
1121 int nrows = prob->nrows_;
1122
1123 /*
1124 Scan the rows for candidates, and write the indices into sort. We're not
1125 interested in rows that are empty or prohibited.
1126
1127 Question: Should we exclude singletons, which are useful in other transforms?
1128 Question: Why are we excluding integral columns?
1129 */
1130 int *sort = new int[nrows];
1131 int nlook = 0;
1132 for (int i = 0; i < nrows; i++) {
1133 if (hinrow[i] == 0)
1134 continue;
1135 if (prob->rowProhibited2(i))
1136 continue;
1137 // sort
1138 CoinSort_2(hcol + mrstrt[i], hcol + mrstrt[i] + hinrow[i],
1139 rowels + mrstrt[i]);
1140 sort[nlook++] = i;
1141 }
1142 if (nlook == 0) {
1143 delete[] sort;
1144 return (next);
1145 }
1146
1147 double *workrow = new double[nrows + 1];
1148
1149 double *workcol;
1150 if (!prob->randomNumber_) {
1151 workcol = new double[ncols + 1];
1152 coin_init_random_vec(workcol, ncols);
1153 } else {
1154 workcol = prob->randomNumber_;
1155 }
1156 compute_sums(nrows, hinrow, mrstrt, hcol, rowels, workcol, sort, workrow, nlook);
1157 CoinSort_2(workrow, workrow + nlook, sort);
1158
1159 double *rlo = prob->rlo_;
1160 double *rup = prob->rup_;
1161
1162 int nuseless_rows = 0;
1163 bool fixInfeasibility = ((prob->presolveOptions_ & 0x4000) != 0);
1164 bool allowIntersection = ((prob->presolveOptions_ & 0x10) != 0);
1165 double tolerance = prob->feasibilityTolerance_;
1166
1167 if (0) {
1168 static int xxxxxx = 0;
1169 int n = 0;
1170 double dval = workrow[0];
1171 int ilastX = 0;
1172 xxxxxx++;
1173 for (int jj = 1; jj < nlook; jj++) {
1174 if (workrow[jj] != dval) {
1175 n += jj - ilastX - 1;
1176 //if (jj>ilastX+2)
1177 //printf("%d matches\n",jj-ilastX);
1178 ilastX = jj;
1179 dval = workrow[jj];
1180 }
1181 }
1182 printf("may be able to delete %d rows - pass %d\n", n, xxxxxx);
1183 }
1184 double dval = workrow[0];
1185 for (int jj = 1; jj < nlook; jj++) {
1186 if (workrow[jj] == dval) {
1187 int ithis = sort[jj];
1188 int ilast = sort[jj - 1];
1189 CoinBigIndex krs = mrstrt[ithis];
1190 CoinBigIndex kre = krs + hinrow[ithis];
1191 if (hinrow[ithis] == hinrow[ilast]) {
1192 CoinBigIndex ishift = mrstrt[ilast] - krs;
1193 CoinBigIndex k;
1194 for (k = krs; k < kre; k++) {
1195 if (hcol[k] != hcol[k + ishift] || fabs(rowels[k] - rowels[k + ishift]) > 1.0e-14) {
1196 break;
1197 }
1198 }
1199 if (k != kre && false) {
1200 int rows[2];
1201 rows[0] = ithis;
1202 rows[1] = ilast;
1203 for (int k = 0; k < 2; k++) {
1204 int kRow = rows[k];
1205 CoinBigIndex krs = mrstrt[kRow];
1206 CoinBigIndex kre = krs + hinrow[kRow];
1207 printf("%g <=", rlo[kRow]);
1208 for (CoinBigIndex k1 = krs; k1 < kre; k1++)
1209 printf(" (%d,%g)", hcol[k1], rowels[k1]);
1210 printf(" <= %g\n", rup[kRow]);
1211 }
1212 }
1213 if (k == kre) {
1214 /* now check rhs to see what is what */
1215 double rlo1 = rlo[ilast];
1216 double rup1 = rup[ilast];
1217 double rlo2 = rlo[ithis];
1218 double rup2 = rup[ithis];
1219
1220 int idelete = -1;
1221 if (rlo1 <= rlo2) {
1222 if (rup2 <= rup1) {
1223 /* this is strictly tighter than last */
1224 idelete = ilast;
1225 PRESOLVE_DETAIL_PRINT(printf("pre_duprow %dR %dR E\n", ilast, ithis));
1226 } else if (fabs(rlo1 - rlo2) < 1.0e-12) {
1227 /* last is strictly tighter than this */
1228 idelete = ithis;
1229 PRESOLVE_DETAIL_PRINT(printf("pre_duprow %dR %dR E\n", ithis, ilast));
1230 // swap so can carry on deleting
1231 sort[jj - 1] = ithis;
1232 sort[jj] = ilast;
1233 } else {
1234 if (rup1 < rlo2 - tolerance && !fixInfeasibility) {
1235 // infeasible
1236 prob->status_ |= 1;
1237 // wrong message - correct if works
1238 prob->messageHandler()->message(COIN_PRESOLVE_ROWINFEAS,
1239 prob->messages())
1240 << ithis
1241 << rlo[ithis]
1242 << rup[ithis]
1243 << CoinMessageEol;
1244 break;
1245 } else if (allowIntersection /*||fabs(rup1-rlo2)<tolerance*/) {
1246 /* overlapping - could merge */
1247 #ifdef CLP_INVESTIGATE7
1248 printf("overlapping duplicate row %g %g, %g %g\n",
1249 rlo1, rup1, rlo2, rup2);
1250 #endif
1251 // pretend this is stricter than last
1252 idelete = ilast;
1253 PRESOLVE_DETAIL_PRINT(printf("pre_duprow %dR %dR E\n", ilast, ithis));
1254 rup[ithis] = rup1;
1255 }
1256 }
1257 } else {
1258 // rlo1>rlo2
1259 if (rup1 <= rup2) {
1260 /* last is strictly tighter than this */
1261 idelete = ithis;
1262 PRESOLVE_DETAIL_PRINT(printf("pre_duprow %dR %dR E\n", ithis, ilast));
1263 // swap so can carry on deleting
1264 sort[jj - 1] = ithis;
1265 sort[jj] = ilast;
1266 } else {
1267 /* overlapping - could merge */
1268 // rlo1>rlo2
1269 // rup1>rup2
1270 if (rup2 < rlo1 - tolerance && !fixInfeasibility) {
1271 // infeasible
1272 prob->status_ |= 1;
1273 // wrong message - correct if works
1274 prob->messageHandler()->message(COIN_PRESOLVE_ROWINFEAS,
1275 prob->messages())
1276 << ithis
1277 << rlo[ithis]
1278 << rup[ithis]
1279 << CoinMessageEol;
1280 break;
1281 } else if (allowIntersection /*||fabs(rup2-rlo1)<tolerance*/) {
1282 #ifdef CLP_INVESTIGATE7
1283 printf("overlapping duplicate row %g %g, %g %g\n",
1284 rlo1, rup1, rlo2, rup2);
1285 #endif
1286 // pretend this is stricter than last
1287 idelete = ilast;
1288 PRESOLVE_DETAIL_PRINT(printf("pre_duprow %dR %dR E\n", ilast, ithis));
1289 rlo[ithis] = rlo1;
1290 }
1291 }
1292 }
1293 if (idelete >= 0)
1294 sort[nuseless_rows++] = idelete;
1295 }
1296 }
1297 }
1298 dval = workrow[jj];
1299 }
1300
1301 delete[] workrow;
1302 if (workcol != prob->randomNumber_)
1303 delete[] workcol;
1304
1305 if (nuseless_rows) {
1306 #if PRESOLVE_SUMMARY
1307 printf("DUPLICATE ROWS: %d\n", nuseless_rows);
1308 #endif
1309 next = useless_constraint_action::presolve(prob,
1310 sort, nuseless_rows,
1311 next);
1312 }
1313 delete[] sort;
1314
1315 if (prob->tuning_) {
1316 double thisTime = CoinCpuTime();
1317 int droppedRows = prob->countEmptyRows() - startEmptyRows;
1318 int droppedColumns = prob->countEmptyCols() - startEmptyColumns;
1319 printf("CoinPresolveDuprow(256) - %d rows, %d columns dropped in time %g, total %g\n",
1320 droppedRows, droppedColumns, thisTime - startTime, thisTime - prob->startTime_);
1321 }
1322 return (next);
1323 }
1324
postsolve(CoinPostsolveMatrix *) const1325 void duprow_action::postsolve(CoinPostsolveMatrix *) const
1326 {
1327 printf("STILL NO POSTSOLVE FOR DUPROW!\n");
1328 abort();
1329 }
1330
1331 #include "CoinFactorization.hpp"
1332
name() const1333 const char *duprow3_action::name() const
1334 {
1335 return ("duprow3_action");
1336 }
1337
1338 // This is just ekkredc4, adapted into the new framework.
1339 /*
1340 I've made minimal changes for compatibility with dupcol: An initial scan to
1341 accumulate rows of interest in sort.
1342 -- lh, 040909 --
1343 */
1344 const CoinPresolveAction
1345 *
presolve(CoinPresolveMatrix * prob,const CoinPresolveAction * next)1346 duprow3_action::presolve(CoinPresolveMatrix *prob,
1347 const CoinPresolveAction *next)
1348 {
1349 double startTime = 0.0;
1350 if (prob->tuning_) {
1351 startTime = CoinCpuTime();
1352 }
1353 //double *rowels = prob->rowels_;
1354 //int *hcol = prob->hcol_;
1355 //CoinBigIndex *mrstrt = prob->mrstrt_;
1356 int *hinrow = prob->hinrow_;
1357 double *colels = prob->colels_;
1358 int *hrow = prob->hrow_;
1359 CoinBigIndex *mcstrt = prob->mcstrt_;
1360 int *hincol = prob->hincol_;
1361 double *rowLower = prob->rlo_;
1362 double *rowUpper = prob->rup_;
1363 double *columnLower = prob->clo_;
1364 double *columnUpper = prob->cup_;
1365 int ncols = prob->ncols_;
1366 int nrows = prob->nrows_;
1367 int numberDropped = 0;
1368 int numberRows = 0;
1369 int *columnMap = prob->usefulColumnInt_; //new int[ncols] ;
1370 int *columnBack = columnMap + ncols;
1371 int *rowMap = new int[2 * nrows];
1372 int *rowBack = rowMap + nrows;
1373 #define SKIP_RHS
1374 int numberRhs = 0;
1375 for (int i = 0; i < nrows; i++) {
1376 if (rowLower[i] == rowUpper[i] && hinrow[i] > 1) {
1377 if (rowLower[i]) {
1378 numberRhs++;
1379 #ifdef SKIP_RHS
1380 rowBack[i] = -1;
1381 continue;
1382 #endif
1383 }
1384 rowBack[i] = numberRows;
1385 rowMap[numberRows++] = i;
1386 } else {
1387 rowBack[i] = -1;
1388 }
1389 }
1390 if (numberRows < CoinMax(100, nrows / 4)) {
1391 //printf("equality rows %d - %d with rhs - unlikely to help much\n",numberRows,numberRhs);
1392 //numberRows=0;
1393 } else {
1394 //printf("equality rows %d - %d with rhs\n",numberRows,numberRhs);
1395 }
1396 #ifdef SKIP_RHS
1397 numberRhs = 0;
1398 #endif
1399 if (numberRows) {
1400 CoinIndexedVector one;
1401 one.reserve(numberRows);
1402 double *minValue = one.denseVector();
1403 CoinIndexedVector two;
1404 two.reserve(numberRows);
1405 double *maxValue = two.denseVector();
1406 for (int i = 0; i < numberRows; i++) {
1407 minValue[i] = COIN_DBL_MAX;
1408 maxValue[i] = 0.0;
1409 }
1410 CoinBigIndex nElements = 0;
1411 int numberColumns = 0;
1412 for (int i = 0; i < ncols; i++) {
1413 if (columnLower[i] < columnUpper[i]) {
1414 int n = 0;
1415 for (CoinBigIndex j = mcstrt[i]; j < mcstrt[i] + hincol[i]; j++) {
1416 int iRow = hrow[j];
1417 iRow = rowBack[iRow];
1418 if (iRow >= 0) {
1419 n++;
1420 double value = fabs(colels[j]);
1421 minValue[iRow] = CoinMin(minValue[iRow], value);
1422 maxValue[iRow] = CoinMax(maxValue[iRow], value);
1423 }
1424 }
1425 if (n) {
1426 nElements += n;
1427 columnBack[i] = numberColumns;
1428 columnMap[numberColumns++] = i;
1429 } else {
1430 columnBack[i] = -1;
1431 }
1432 }
1433 }
1434 //printf("%d columns %d elements\n",numberColumns,nElements);
1435 CoinFactorization factorization;
1436 factorization.setDenseThreshold(0);
1437 CoinPackedMatrix matrix(true, 0.0, 0.0);
1438 matrix.reserve(numberColumns, nElements);
1439 int maxDimension = CoinMax(numberRows, numberColumns);
1440 matrix.setDimensions(maxDimension, numberColumns);
1441 double *element = matrix.getMutableElements();
1442 int *row = matrix.getMutableIndices();
1443 CoinBigIndex *columnStart = matrix.getMutableVectorStarts();
1444 int *columnLength = matrix.getMutableVectorLengths();
1445 // scale rows
1446 for (int i = 0; i < numberRows; i++) {
1447 minValue[i] = 1.0 / sqrt(minValue[i] * maxValue[i]);
1448 }
1449 nElements = 0;
1450 columnStart[0] = 0;
1451 for (int i = 0; i < numberColumns; i++) {
1452 int iColumn = columnMap[i];
1453 for (CoinBigIndex j = mcstrt[iColumn]; j < mcstrt[iColumn] + hincol[iColumn]; j++) {
1454 int iRow = hrow[j];
1455 iRow = rowBack[iRow];
1456 if (iRow >= 0) {
1457 row[nElements] = iRow;
1458 element[nElements++] = minValue[iRow] * colels[j];
1459 }
1460 }
1461 columnLength[i] = static_cast< int >(nElements - columnStart[i]);
1462 columnStart[i + 1] = nElements;
1463 }
1464 matrix.setNumElements(nElements);
1465 int *rowIsBasic = new int[maxDimension];
1466 int *columnIsBasic = new int[maxDimension];
1467 //int numberBasic=0;
1468 for (int i = 0; i < maxDimension; i++) {
1469 rowIsBasic[i] = -1;
1470 }
1471 for (int i = 0; i < numberColumns; i++) {
1472 columnIsBasic[i] = 1;
1473 }
1474 //returns 0 -okay, -1 singular, -2 too many in basis, -99 memory */
1475 int status = factorization.factorize(matrix, rowIsBasic, columnIsBasic, 5.0);
1476 //printf("factorization status %d\n",status);
1477 if (status == -1) {
1478 const int *permute = factorization.permute();
1479 const int *pivotColumn = factorization.pivotColumn();
1480 int numberBasic = factorization.numberGoodColumns();
1481 //printf("%d good columns out of %d - %d rows\n",numberBasic,
1482 // numberColumns,numberRows);
1483 if (numberBasic < numberRows - CoinMax(20, nrows / 10)) {
1484 // get list of real rows which are dependent (maybe infeasible)
1485 int *badRow = new int[numberRows - numberBasic];
1486 int nBad = 0;
1487 numberBasic = 0;
1488 for (int i = 0; i < numberColumns; i++) {
1489 //printf("%d permute %d pivotColumn %d\n",
1490 // i,permute[i],pivotColumn[i]);
1491 if (pivotColumn[i] < 0)
1492 columnIsBasic[i] = -1;
1493 else
1494 numberBasic++;
1495 }
1496 for (int i = 0; i < numberRows; i++) {
1497 if (permute[i] < 0) {
1498 badRow[nBad++] = rowMap[i];
1499 rowIsBasic[i] = 1;
1500 }
1501 }
1502 assert(numberBasic + nBad == numberRows);
1503 // factorize again
1504 if (numberRows < maxDimension) {
1505 int nDelete = maxDimension - numberRows;
1506 int *del = new int[nDelete];
1507 for (int i = 0; i < nDelete; i++)
1508 del[i] = i + numberRows;
1509 matrix.deleteRows(nDelete, del);
1510 delete[] del;
1511 }
1512 if (numberRhs) {
1513 status = factorization.factorize(matrix, rowIsBasic, columnIsBasic, 5.0);
1514 //printf("second factorization status %d\n",status);
1515 } else {
1516 //printf("skipping second factorization as zero Rhs\n");
1517 status = 0;
1518 }
1519 if (status == 0) {
1520 // check feasible
1521 numberDropped = 0;
1522 if (numberRhs) {
1523 memset(minValue, 0, numberRows * sizeof(double));
1524 memset(maxValue, 0, numberRows * sizeof(double));
1525 int *index = one.getIndices();
1526 for (int i = 0; i < numberRows; i++) {
1527 int iRow = rowMap[i];
1528 double rhs = rowLower[iRow];
1529 minValue[i] = rhs;
1530 }
1531 for (int i = 0; i < nBad; i++) {
1532 int iRow = badRow[i];
1533 iRow = rowBack[iRow];
1534 minValue[iRow] = 0.0;
1535 }
1536 int n = 0;
1537 for (int i = 0; i < numberRows; i++) {
1538 if (minValue[i])
1539 index[n++] = i;
1540 }
1541 one.setNumElements(n);
1542 factorization.updateColumn(&two, &one);
1543 for (int i = 0; i < nBad; i++) {
1544 int iRealRow = badRow[i];
1545 int iRow = rowBack[iRealRow];
1546 if (fabs(rowLower[iRealRow] - minValue[iRow]) > 1.0e-7) {
1547 //printf("bad %d %d real row %d rhs %g - computed %g\n",
1548 // i,iRow,iRealRow,rowLower[iRealRow],minValue[iRow]);
1549 } else {
1550 badRow[numberDropped++] = iRealRow;
1551 }
1552 }
1553 } else {
1554 for (int i = 0; i < nBad; i++) {
1555 int iRealRow = badRow[i];
1556 badRow[numberDropped++] = iRealRow;
1557 }
1558 }
1559 if (numberDropped) {
1560 next = useless_constraint_action::presolve(prob,
1561 badRow, numberDropped,
1562 next);
1563 #ifdef CLP_USEFUL_PRINTOUT
1564 printf("Dependency checking dropped %d rows\n", numberDropped);
1565 #endif
1566 }
1567 }
1568 delete[] badRow;
1569 }
1570 }
1571 delete[] rowIsBasic;
1572 delete[] columnIsBasic;
1573 }
1574 delete[] rowMap;
1575 if (prob->tuning_) {
1576 double thisTime = CoinCpuTime();
1577 printf("CoinPresolveDuprow3 - %d rows dropped in time %g, total %g\n",
1578 numberDropped, thisTime - startTime, thisTime - prob->startTime_);
1579 }
1580 #ifdef CLP_INVESTIGATE
1581 printf("CoinPresolveDuprow3 - %d rows dropped\n", numberDropped);
1582 #endif
1583 return (next);
1584 }
1585
postsolve(CoinPostsolveMatrix *) const1586 void duprow3_action::postsolve(CoinPostsolveMatrix *) const
1587 {
1588 printf("STILL NO POSTSOLVE FOR DUPROW3!\n");
1589 abort();
1590 }
1591
1592 /*
1593 Routines for gub rows. This is definitely unfinished --- there's no
1594 postsolve action.
1595
1596 This is potentially called from ClpPresolve and OsiPresolve. Unclear that
1597 it can be backed out --- there's no postsolve.
1598 */
1599
name() const1600 const char *gubrow_action::name() const
1601 {
1602 return ("gubrow_action");
1603 }
1604
1605 const CoinPresolveAction
1606 *
presolve(CoinPresolveMatrix * prob,const CoinPresolveAction * next)1607 gubrow_action::presolve(CoinPresolveMatrix *prob,
1608 const CoinPresolveAction *next)
1609 {
1610 double startTime = 0.0;
1611 int droppedElements = 0;
1612 int affectedRows = 0;
1613 if (prob->tuning_) {
1614 startTime = CoinCpuTime();
1615 }
1616 double *rowels = prob->rowels_;
1617 int *hcol = prob->hcol_;
1618 CoinBigIndex *mrstrt = prob->mrstrt_;
1619 int *hinrow = prob->hinrow_;
1620 double *colels = prob->colels_;
1621 int *hrow = prob->hrow_;
1622 CoinBigIndex *mcstrt = prob->mcstrt_;
1623 int *hincol = prob->hincol_;
1624 int ncols = prob->ncols_;
1625 int nrows = prob->nrows_;
1626 double *rlo = prob->rlo_;
1627 double *rup = prob->rup_;
1628
1629 // maximum number of records
1630 action *actions = new action[nrows];
1631 int nactions = 0;
1632 /*
1633 Scan the rows. We're not
1634 interested in rows that are empty or prohibited.
1635
1636 */
1637 int *which = prob->usefulRowInt_;
1638 int *number = which + nrows;
1639 double *els = prob->usefulRowDouble_;
1640 char *markCol = reinterpret_cast< char * >(prob->usefulColumnInt_);
1641 memset(markCol, 0, ncols);
1642 CoinZeroN(els, nrows);
1643 for (int i = 0; i < nrows; i++) {
1644 int nInRow = hinrow[i];
1645 if (nInRow > 1 && !prob->rowProhibited2(i) && rlo[i] == rup[i]) {
1646 CoinBigIndex rStart = mrstrt[i];
1647 CoinBigIndex k = rStart;
1648 CoinBigIndex rEnd = rStart + nInRow;
1649 double value1 = rowels[k];
1650 k++;
1651 for (; k < rEnd; k++) {
1652 if (rowels[k] != value1)
1653 break;
1654 }
1655 if (k == rEnd) {
1656 // Gub row
1657 int nLook = 0;
1658 for (k = rStart; k < rEnd; k++) {
1659 int iColumn = hcol[k];
1660 markCol[iColumn] = 1;
1661 CoinBigIndex kk = mcstrt[iColumn];
1662 CoinBigIndex cEnd = kk + hincol[iColumn];
1663 for (; kk < cEnd; kk++) {
1664 int iRow = hrow[kk];
1665 double value = colels[kk];
1666 if (iRow != i) {
1667 double value2 = els[iRow];
1668 if (value2) {
1669 if (value == value2)
1670 number[iRow]++;
1671 } else {
1672 // first
1673 els[iRow] = value;
1674 number[iRow] = 1;
1675 which[nLook++] = iRow;
1676 }
1677 }
1678 }
1679 }
1680 // Now see if any promising
1681 int nDrop = 0;
1682 for (int j = 0; j < nLook; j++) {
1683 int iRow = which[j];
1684 if (number[iRow] == nInRow) {
1685 // can delete elements and adjust rhs
1686 for (CoinBigIndex kk = rStart; kk < rEnd; kk++)
1687 presolve_delete_from_col(iRow, hcol[kk], mcstrt, hincol, hrow, colels);
1688 int nInRow2 = hinrow[iRow];
1689 CoinBigIndex rStart2 = mrstrt[iRow];
1690 CoinBigIndex rEnd2 = rStart2 + nInRow2;
1691 for (CoinBigIndex kk = rStart2; kk < rEnd2; kk++) {
1692 int iColumn = hcol[kk];
1693 if (markCol[iColumn] == 0) {
1694 hcol[rStart2] = iColumn;
1695 rowels[rStart2++] = rowels[kk];
1696 }
1697 }
1698 hinrow[iRow] = nInRow2 - nInRow;
1699 nDrop++;
1700 if (!hinrow[iRow])
1701 PRESOLVE_REMOVE_LINK(prob->rlink_, iRow);
1702 double value = (rlo[i] / value1) * els[iRow];
1703 // correct rhs
1704 if (rlo[iRow] > -1.0e20)
1705 rlo[iRow] -= value;
1706 if (rup[iRow] < 1.0e20)
1707 rup[iRow] -= value;
1708 } else {
1709 number[iRow] = 0;
1710 }
1711 }
1712 if (nDrop) {
1713 affectedRows += nDrop;
1714 droppedElements += nDrop * nInRow;
1715 action &thisAction = actions[nactions];
1716 int *deletedRow = new int[nDrop + 1];
1717 int *indices = CoinCopyOfArray(hcol + rStart, nInRow);
1718 thisAction.indices = indices;
1719 double *rowels = new double[nDrop + 1];
1720 thisAction.rhs = rlo[i];
1721 deletedRow[nDrop] = i;
1722 rowels[nDrop] = value1;
1723 nDrop = 0;
1724 for (int j = 0; j < nLook; j++) {
1725 int iRow = which[j];
1726 if (number[iRow]) {
1727 deletedRow[nDrop] = iRow;
1728 rowels[nDrop++] = els[iRow];
1729 }
1730 }
1731 thisAction.nDrop = nDrop;
1732 thisAction.ninrow = nInRow;
1733 thisAction.deletedRow = deletedRow;
1734 thisAction.rowels = rowels;
1735 nactions++;
1736 }
1737 for (int j = 0; j < nLook; j++) {
1738 int iRow = which[j];
1739 els[iRow] = 0.0;
1740 }
1741 for (k = rStart; k < rEnd; k++) {
1742 int iColumn = hcol[k];
1743 markCol[iColumn] = 0;
1744 }
1745 }
1746 }
1747 }
1748 if (nactions) {
1749 next = new gubrow_action(nactions,
1750 CoinCopyOfArray(actions, nactions), next);
1751 }
1752 deleteAction(actions, action *);
1753 if (prob->tuning_) {
1754 double thisTime = CoinCpuTime();
1755 printf("CoinPresolveGubrow(1024) - %d elements dropped (%d rows) in time %g, total %g\n",
1756 droppedElements, affectedRows, thisTime - startTime, thisTime - prob->startTime_);
1757 } else if (droppedElements) {
1758 #ifdef CLP_INVESTIGATE
1759 printf("CoinPresolveGubrow(1024) - %d elements dropped (%d rows)\n",
1760 droppedElements, affectedRows);
1761 #endif
1762 }
1763 return (next);
1764 }
1765
postsolve(CoinPostsolveMatrix * prob) const1766 void gubrow_action::postsolve(CoinPostsolveMatrix *prob) const
1767 {
1768 #if PRESOLVE_DEBUG > 0 || PRESOLVE_CONSISTENCY > 0
1769 #if PRESOLVE_DEBUG > 0
1770 std::cout
1771 << "Entering gubrow_action::postsolve, "
1772 << nactions_ << " constraints to process." << std::endl;
1773 #endif
1774 //int ncols = prob->ncols_ ;
1775 //char *cdone = prob->cdone_ ;
1776 //char *rdone = prob->rdone_ ;
1777 //const double ztolzb = prob->ztolzb_ ;
1778
1779 presolve_check_threads(prob);
1780 presolve_check_free_list(prob);
1781 presolve_check_reduced_costs(prob);
1782 presolve_check_duals(prob);
1783 presolve_check_sol(prob, 2, 2, 2);
1784 presolve_check_nbasic(prob);
1785 #endif
1786
1787 /*
1788 Unpack the column-major representation.
1789 */
1790 CoinBigIndex *colStarts = prob->mcstrt_;
1791 int *colLengths = prob->hincol_;
1792 int *rowIndices = prob->hrow_;
1793 double *colCoeffs = prob->colels_;
1794 /*
1795 Rim vectors, solution, reduced costs, duals, row activity.
1796 */
1797 double *rlo = prob->rlo_;
1798 double *rup = prob->rup_;
1799 //double *cost = prob->cost_ ;
1800 //double *sol = prob->sol_ ;
1801 //double *rcosts = prob->rcosts_ ;
1802 double *acts = prob->acts_;
1803 double *rowduals = prob->rowduals_;
1804
1805 CoinBigIndex *link = prob->link_;
1806 CoinBigIndex &free_list = prob->free_list_;
1807
1808 //const double maxmin = prob->maxmin_ ;
1809
1810 const action *const actions = actions_;
1811 const int nactions = nactions_;
1812
1813 /*
1814 Open the main loop to step through the postsolve objects.
1815
1816 */
1817 for (const action *f = &actions[nactions - 1]; actions <= f; f--) {
1818 int *deletedRow = f->deletedRow;
1819 double *els = f->rowels;
1820 int *indices = f->indices;
1821 int nDrop = f->nDrop;
1822 int ninrow = f->ninrow;
1823 double rhs = f->rhs;
1824 double coeff = els[nDrop];
1825 int iRow = deletedRow[nDrop];
1826 for (int i = 0; i < nDrop; i++) {
1827 int tgtrow = deletedRow[i];
1828 double coeffy = els[i];
1829 double dualValue = rowduals[tgtrow];
1830 #if PRESOLVE_DEBUG > 2
1831 std::cout << " restoring row " << tgtrow << " coeff " << coeffy
1832 << " dual " << dualValue
1833 << " dual on " << iRow << " goes from "
1834 << rowduals[iRow] << " to "
1835 << rowduals[iRow] - (coeffy * dualValue) / coeff
1836 << std::endl;
1837 #endif
1838 // adjust dual on iRow
1839 rowduals[iRow] -= (coeffy * dualValue) / coeff;
1840 #if 0 //PRESOLVE_DEBUG > 2
1841 int * rowStart;
1842 int * column;
1843 double * element;
1844 postsolve_get_rowcopy(prob,rowStart,column,element);
1845 for (int k=rowStart[tgtrow];k<rowStart[tgtrow+1];k++) {
1846 printf("(%d, %g, %s - rcost %g) ",column[k],element[k],
1847 statusName(prob->getColumnStatus(column[k])),rcosts[column[k]]);
1848 if ((k-rowStart[tgtrow])%10==9)
1849 printf("\n");
1850 }
1851 printf("\n");
1852 #endif
1853
1854 for (int rndx = 0; rndx < ninrow; ++rndx) {
1855 int j = indices[rndx];
1856 #if PRESOLVE_DEBUG > 2
1857 std::cout << " a(" << j << " rcost " << rcosts[j]
1858 << " -> " << rcosts[j] - dualValue * coeffy << " " << statusName(prob->getColumnStatus(j)) << ")";
1859 #endif
1860 CoinBigIndex kk = free_list;
1861 assert(kk >= 0 && kk < prob->bulk0_);
1862 free_list = link[free_list];
1863 link[kk] = colStarts[j];
1864 colStarts[j] = kk;
1865 colCoeffs[kk] = coeffy;
1866 rowIndices[kk] = tgtrow;
1867 ++colLengths[j];
1868 }
1869 double value = (rhs / coeff) * coeffy;
1870 acts[tgtrow] += value;
1871 ;
1872 // correct rhs
1873 if (rlo[tgtrow] > -1.0e20)
1874 rlo[tgtrow] += value;
1875 if (rup[tgtrow] < 1.0e20)
1876 rup[tgtrow] += value;
1877 #if PRESOLVE_DEBUG > 2
1878 std::cout << std::endl;
1879 #endif
1880
1881 #if PRESOLVE_CONSISTENCY > 0
1882 presolve_check_threads(prob);
1883 presolve_check_free_list(prob);
1884 #endif
1885 }
1886 }
1887
1888 #if PRESOLVE_DEBUG > 0 || PRESOLVE_CONSISTENCY > 0
1889 presolve_check_threads(prob);
1890 presolve_check_free_list(prob);
1891 presolve_check_reduced_costs(prob);
1892 presolve_check_duals(prob);
1893 presolve_check_sol(prob, 2, 2, 2);
1894 presolve_check_nbasic(prob);
1895 #if PRESOLVE_DEBUG > 0
1896 std::cout << "Leaving gubrow_action::postsolve." << std::endl;
1897 #endif
1898 #endif
1899
1900 return;
1901 }
1902
~gubrow_action()1903 gubrow_action::~gubrow_action()
1904 {
1905 const action *actions = actions_;
1906
1907 for (int i = 0; i < nactions_; ++i) {
1908 delete[] actions[i].rowels;
1909 delete[] actions[i].deletedRow;
1910 delete[] actions[i].indices;
1911 }
1912
1913 // Must add cast to placate MS compiler
1914 deleteAction(actions_, gubrow_action::action *);
1915 }
1916
1917 /*
1918 Routines for two by two blocks. This is definitely unfinished --- there's no
1919 postsolve action.
1920 */
1921
name() const1922 const char *twoxtwo_action::name() const
1923 {
1924 return ("twoxtwo_action");
1925 }
1926
1927 const CoinPresolveAction
1928 *
presolve(CoinPresolveMatrix * prob,const CoinPresolveAction * next)1929 twoxtwo_action::presolve(CoinPresolveMatrix *prob,
1930 const CoinPresolveAction *next)
1931 {
1932 double startTime = 0.0;
1933 int startEmptyRows = 0;
1934 int startEmptyColumns = 0;
1935 if (prob->tuning_) {
1936 startTime = CoinCpuTime();
1937 startEmptyRows = prob->countEmptyRows();
1938 startEmptyColumns = prob->countEmptyCols();
1939 }
1940 // maximum number of records
1941 action *boundRecords = new action[(prob->nrows_ + 1) >> 1];
1942 int nactions = 0;
1943 // column-major representation
1944 const int ncols = prob->ncols_;
1945 const CoinBigIndex *const mcstrt = prob->mcstrt_;
1946 const int *const hincol = prob->hincol_;
1947 const int *const hrow = prob->hrow_;
1948 const double *colels = prob->colels_;
1949 double *cost = prob->cost_;
1950
1951 // column type, bounds, solution, and status
1952 const unsigned char *const integerType = prob->integerType_;
1953 double *clo = prob->clo_;
1954 double *cup = prob->cup_;
1955 // row-major representation
1956 //const int nrows = prob->nrows_ ;
1957 const CoinBigIndex *const mrstrt = prob->mrstrt_;
1958 const int *const hinrow = prob->hinrow_;
1959 const int *const hcol = prob->hcol_;
1960 const double *rowels = prob->rowels_;
1961
1962 // row bounds
1963 double *rlo = prob->rlo_;
1964 double *rup = prob->rup_;
1965
1966 // tolerances
1967 //const double ekkinf2 = PRESOLVE_SMALL_INF ;
1968 //const double ekkinf = ekkinf2*1.0e8 ;
1969 //const double ztolcbarj = prob->ztoldj_ ;
1970 //const CoinRelFltEq relEq(prob->ztolzb_) ;
1971 double bound[2];
1972 double alpha[2] = { 0.0, 0.0 };
1973 double offset = 0.0;
1974
1975 for (int icol = 0; icol < ncols; icol++) {
1976 if (hincol[icol] == 2) {
1977 CoinBigIndex start = mcstrt[icol];
1978 int row0 = hrow[start];
1979 if (hinrow[row0] != 2)
1980 continue;
1981 int row1 = hrow[start + 1];
1982 if (hinrow[row1] != 2)
1983 continue;
1984 double element0 = colels[start];
1985 double rowUpper0 = rup[row0];
1986 bool swapSigns0 = false;
1987 if (rlo[row0] > -1.0e30) {
1988 if (rup[row0] > 1.0e30) {
1989 swapSigns0 = true;
1990 rowUpper0 = -rlo[row0];
1991 element0 = -element0;
1992 } else {
1993 // range or equality
1994 continue;
1995 }
1996 } else if (rup[row0] > 1.0e30) {
1997 // free
1998 continue;
1999 }
2000 #if 0
2001 // skip here for speed
2002 // skip if no cost (should be able to get rid of)
2003 if (!cost[icol]) {
2004 PRESOLVE_DETAIL_PRINT(printf("should be able to get rid of %d with no cost\n",icol));
2005 continue;
2006 }
2007 // skip if negative cost for now
2008 if (cost[icol]<0.0) {
2009 PRESOLVE_DETAIL_PRINT(printf("code for negative cost\n"));
2010 continue;
2011 }
2012 #endif
2013 double element1 = colels[start + 1];
2014 double rowUpper1 = rup[row1];
2015 bool swapSigns1 = false;
2016 if (rlo[row1] > -1.0e30) {
2017 if (rup[row1] > 1.0e30) {
2018 swapSigns1 = true;
2019 rowUpper1 = -rlo[row1];
2020 element1 = -element1;
2021 } else {
2022 // range or equality
2023 continue;
2024 }
2025 } else if (rup[row1] > 1.0e30) {
2026 // free
2027 continue;
2028 }
2029 double lowerX = clo[icol];
2030 double upperX = cup[icol];
2031 int otherCol = -1;
2032 CoinBigIndex startRow = mrstrt[row0];
2033 for (CoinBigIndex j = startRow; j < startRow + 2; j++) {
2034 int jcol = hcol[j];
2035 if (jcol != icol) {
2036 alpha[0] = swapSigns0 ? -rowels[j] : rowels[j];
2037 otherCol = jcol;
2038 }
2039 }
2040 startRow = mrstrt[row1];
2041 bool possible = true;
2042 for (CoinBigIndex j = startRow; j < startRow + 2; j++) {
2043 int jcol = hcol[j];
2044 if (jcol != icol) {
2045 if (jcol == otherCol) {
2046 alpha[1] = swapSigns1 ? -rowels[j] : rowels[j];
2047 } else {
2048 possible = false;
2049 }
2050 }
2051 }
2052 if (possible) {
2053 // skip if no cost (should be able to get rid of)
2054 if (!cost[icol]) {
2055 PRESOLVE_DETAIL_PRINT(printf("should be able to get rid of %d with no cost\n", icol));
2056 continue;
2057 }
2058 // skip if negative cost for now
2059 if (cost[icol] < 0.0) {
2060 PRESOLVE_DETAIL_PRINT(printf("code for negative cost\n"));
2061 continue;
2062 }
2063 bound[0] = clo[otherCol];
2064 bound[1] = cup[otherCol];
2065 double lowestLowest = COIN_DBL_MAX;
2066 double highestLowest = -COIN_DBL_MAX;
2067 double lowestHighest = COIN_DBL_MAX;
2068 double highestHighest = -COIN_DBL_MAX;
2069 int binding0 = 0;
2070 int binding1 = 0;
2071 for (int k = 0; k < 2; k++) {
2072 bool infLow0 = false;
2073 bool infLow1 = false;
2074 double sum0 = 0.0;
2075 double sum1 = 0.0;
2076 double value = bound[k];
2077 if (fabs(value) < 1.0e30) {
2078 sum0 += alpha[0] * value;
2079 sum1 += alpha[1] * value;
2080 } else {
2081 if (alpha[0] > 0.0) {
2082 if (value < 0.0)
2083 infLow0 = true;
2084 } else if (alpha[0] < 0.0) {
2085 if (value > 0.0)
2086 infLow0 = true;
2087 }
2088 if (alpha[1] > 0.0) {
2089 if (value < 0.0)
2090 infLow1 = true;
2091 } else if (alpha[1] < 0.0) {
2092 if (value > 0.0)
2093 infLow1 = true;
2094 }
2095 }
2096 /* Got sums
2097 */
2098 double thisLowest0 = -COIN_DBL_MAX;
2099 double thisHighest0 = COIN_DBL_MAX;
2100 if (element0 > 0.0) {
2101 // upper bound unless inf&2 !=0
2102 if (!infLow0)
2103 thisHighest0 = (rowUpper0 - sum0) / element0;
2104 } else {
2105 // lower bound unless inf&2 !=0
2106 if (!infLow0)
2107 thisLowest0 = (rowUpper0 - sum0) / element0;
2108 }
2109 double thisLowest1 = -COIN_DBL_MAX;
2110 double thisHighest1 = COIN_DBL_MAX;
2111 if (element1 > 0.0) {
2112 // upper bound unless inf&2 !=0
2113 if (!infLow1)
2114 thisHighest1 = (rowUpper1 - sum1) / element1;
2115 } else {
2116 // lower bound unless inf&2 !=0
2117 if (!infLow1)
2118 thisLowest1 = (rowUpper1 - sum1) / element1;
2119 }
2120 if (thisLowest0 > thisLowest1 + 1.0e-12) {
2121 if (thisLowest0 > lowerX + 1.0e-12)
2122 binding0 |= 1 << k;
2123 } else if (thisLowest1 > thisLowest0 + 1.0e-12) {
2124 if (thisLowest1 > lowerX + 1.0e-12)
2125 binding1 |= 1 << k;
2126 thisLowest0 = thisLowest1;
2127 }
2128 if (thisHighest0 < thisHighest1 - 1.0e-12) {
2129 if (thisHighest0 < upperX - 1.0e-12)
2130 binding0 |= 1 << k;
2131 } else if (thisHighest1 < thisHighest0 - 1.0e-12) {
2132 if (thisHighest1 < upperX - 1.0e-12)
2133 binding1 |= 1 << k;
2134 thisHighest0 = thisHighest1;
2135 }
2136 lowestLowest = CoinMin(lowestLowest, thisLowest0);
2137 highestHighest = CoinMax(highestHighest, thisHighest0);
2138 lowestHighest = CoinMin(lowestHighest, thisHighest0);
2139 highestLowest = CoinMax(highestLowest, thisLowest0);
2140 }
2141 // see if any good
2142 //#define PRINT_VALUES
2143 if (!binding0 || !binding1) {
2144 PRESOLVE_DETAIL_PRINT(printf("Row redundant for column %d\n", icol));
2145 } else {
2146 PRESOLVE_DETAIL_PRINT(printf("Column %d bounds %g,%g lowest %g,%g highest %g,%g\n",
2147 icol, lowerX, upperX, lowestLowest, lowestHighest,
2148 highestLowest, highestHighest));
2149 // if integer adjust
2150 if (integerType[icol]) {
2151 lowestLowest = ceil(lowestLowest - 1.0e-5);
2152 highestLowest = ceil(highestLowest - 1.0e-5);
2153 lowestHighest = floor(lowestHighest + 1.0e-5);
2154 highestHighest = floor(highestHighest + 1.0e-5);
2155 }
2156 // if costed may be able to adjust
2157 if (cost[icol] >= 0.0) {
2158 if (highestLowest < upperX && highestLowest >= lowerX && highestHighest < 1.0e30) {
2159 highestHighest = CoinMin(highestHighest, highestLowest);
2160 }
2161 }
2162 if (cost[icol] <= 0.0) {
2163 if (lowestHighest > lowerX && lowestHighest <= upperX && lowestHighest > -1.0e30) {
2164 lowestLowest = CoinMax(lowestLowest, lowestHighest);
2165 }
2166 }
2167 #if 1
2168 if (lowestLowest > lowerX + 1.0e-8) {
2169 PRESOLVE_DETAIL_PRINT(printf("Can increase lower bound on %d from %g to %g\n",
2170 icol, lowerX, lowestLowest));
2171 lowerX = lowestLowest;
2172 }
2173 if (highestHighest < upperX - 1.0e-8) {
2174 PRESOLVE_DETAIL_PRINT(printf("Can decrease upper bound on %d from %g to %g\n",
2175 icol, upperX, highestHighest));
2176 upperX = highestHighest;
2177 }
2178 #endif
2179 // see if we can move costs
2180 double xValue;
2181 double yValue0;
2182 double yValue1;
2183 double newLower = COIN_DBL_MAX;
2184 double newUpper = -COIN_DBL_MAX;
2185 double costEqual;
2186 double slope[2];
2187 assert(binding0 + binding1 == 3);
2188 // get where equal
2189 xValue = (rowUpper0 * element1 - rowUpper1 * element0) / (alpha[0] * element1 - alpha[1] * element0);
2190 yValue0 = (rowUpper0 - xValue * alpha[0]) / element0;
2191 yValue1 = (rowUpper1 - xValue * alpha[1]) / element1;
2192 newLower = CoinMin(newLower, CoinMax(yValue0, yValue1));
2193 newUpper = CoinMax(newUpper, CoinMax(yValue0, yValue1));
2194 double xValueEqual = xValue;
2195 double yValueEqual = yValue0;
2196 costEqual = xValue * cost[otherCol] + yValueEqual * cost[icol];
2197 if (binding0 == 1) {
2198 // take x 1.0 down
2199 double x = xValue - 1.0;
2200 double y = (rowUpper0 - x * alpha[0]) / element0;
2201 double costTotal = x * cost[otherCol] + y * cost[icol];
2202 slope[0] = costEqual - costTotal;
2203 // take x 1.0 up
2204 x = xValue + 1.0;
2205 y = (rowUpper1 - x * alpha[1]) / element0;
2206 costTotal = x * cost[otherCol] + y * cost[icol];
2207 slope[1] = costTotal - costEqual;
2208 } else {
2209 // take x 1.0 down
2210 double x = xValue - 1.0;
2211 double y = (rowUpper1 - x * alpha[1]) / element0;
2212 double costTotal = x * cost[otherCol] + y * cost[icol];
2213 slope[1] = costEqual - costTotal;
2214 // take x 1.0 up
2215 x = xValue + 1.0;
2216 y = (rowUpper0 - x * alpha[0]) / element0;
2217 costTotal = x * cost[otherCol] + y * cost[icol];
2218 slope[0] = costTotal - costEqual;
2219 }
2220 PRESOLVE_DETAIL_PRINT(printf("equal value of %d is %g, value of %d is max(%g,%g) - %g\n",
2221 otherCol, xValue, icol, yValue0, yValue1, CoinMax(yValue0, yValue1)));
2222 PRESOLVE_DETAIL_PRINT(printf("Cost at equality %g for constraint 0 slope %g for constraint 1 slope %g\n",
2223 costEqual, slope[0], slope[1]));
2224 xValue = bound[0];
2225 yValue0 = (rowUpper0 - xValue * alpha[0]) / element0;
2226 yValue1 = (rowUpper1 - xValue * alpha[1]) / element1;
2227 PRESOLVE_DETAIL_PRINT(printf("value of %d is %g, value of %d is max(%g,%g) - %g\n",
2228 otherCol, xValue, icol, yValue0, yValue1, CoinMax(yValue0, yValue1)));
2229 newLower = CoinMin(newLower, CoinMax(yValue0, yValue1));
2230 // cost>0 so will be at lower
2231 //double yValueAtBound0=newLower;
2232 newUpper = CoinMax(newUpper, CoinMax(yValue0, yValue1));
2233 xValue = bound[1];
2234 yValue0 = (rowUpper0 - xValue * alpha[0]) / element0;
2235 yValue1 = (rowUpper1 - xValue * alpha[1]) / element1;
2236 PRESOLVE_DETAIL_PRINT(printf("value of %d is %g, value of %d is max(%g,%g) - %g\n",
2237 otherCol, xValue, icol, yValue0, yValue1, CoinMax(yValue0, yValue1)));
2238 newLower = CoinMin(newLower, CoinMax(yValue0, yValue1));
2239 // cost>0 so will be at lower
2240 //double yValueAtBound1=newLower;
2241 newUpper = CoinMax(newUpper, CoinMax(yValue0, yValue1));
2242 lowerX = CoinMax(lowerX, newLower - 1.0e-12 * fabs(newLower));
2243 upperX = CoinMin(upperX, newUpper + 1.0e-12 * fabs(newUpper));
2244 // Now make duplicate row
2245 // keep row 0 so need to adjust costs so same
2246 PRESOLVE_DETAIL_PRINT(printf("Costs for x %g,%g,%g are %g,%g,%g\n",
2247 xValueEqual - 1.0, xValueEqual, xValueEqual + 1.0,
2248 costEqual - slope[0], costEqual, costEqual + slope[1]));
2249 double costOther = cost[otherCol] + slope[1];
2250 double costThis = cost[icol] + slope[1] * (element0 / alpha[0]);
2251 xValue = xValueEqual;
2252 yValue0 = CoinMax((rowUpper0 - xValue * alpha[0]) / element0, lowerX);
2253 double thisOffset = costEqual - (costOther * xValue + costThis * yValue0);
2254 offset += thisOffset;
2255 PRESOLVE_DETAIL_PRINT(printf("new cost at equal %g\n", costOther * xValue + costThis * yValue0 + thisOffset));
2256 xValue = xValueEqual - 1.0;
2257 yValue0 = CoinMax((rowUpper0 - xValue * alpha[0]) / element0, lowerX);
2258 PRESOLVE_DETAIL_PRINT(printf("new cost at -1 %g\n", costOther * xValue + costThis * yValue0 + thisOffset));
2259 assert(fabs((costOther * xValue + costThis * yValue0 + thisOffset) - (costEqual - slope[0])) < 1.0e-5);
2260 xValue = xValueEqual + 1.0;
2261 yValue0 = CoinMax((rowUpper0 - xValue * alpha[0]) / element0, lowerX);
2262 PRESOLVE_DETAIL_PRINT(printf("new cost at +1 %g\n", costOther * xValue + costThis * yValue0 + thisOffset));
2263 assert(fabs((costOther * xValue + costThis * yValue0 + thisOffset) - (costEqual + slope[1])) < 1.0e-5);
2264 action &boundRecord = boundRecords[nactions++];
2265 boundRecord.row = row1;
2266 boundRecord.col = icol;
2267 boundRecord.othercol = otherCol;
2268 boundRecord.lbound_row = rlo[row1];
2269 boundRecord.ubound_row = rup[row1];
2270 boundRecord.lbound_col = clo[icol];
2271 boundRecord.ubound_col = cup[icol];
2272 boundRecord.cost_col = cost[icol];
2273 boundRecord.cost_othercol = cost[otherCol];
2274 cost[otherCol] = costOther;
2275 cost[icol] = costThis;
2276 clo[icol] = lowerX;
2277 cup[icol] = upperX;
2278 // make row useless
2279 rlo[row1] = -COIN_DBL_MAX;
2280 rup[row1] = COIN_DBL_MAX;
2281 }
2282 }
2283 }
2284 }
2285 if (nactions) {
2286 #if PRESOLVE_SUMMARY
2287 printf("Cost offset %g - from %d blocks\n", offset, nactions);
2288 printf("TWO by TWO blocks: %d - offset %g\n", nactions, offset);
2289 #endif
2290 action *actions = new action[nactions];
2291 memcpy(actions, boundRecords, nactions * sizeof(action));
2292 next = new twoxtwo_action(nactions, actions, next);
2293 int *sort = prob->usefulColumnInt_;
2294 for (int i = 0; i < nactions; i++)
2295 sort[i] = boundRecords[i].row;
2296 next = useless_constraint_action::presolve(prob,
2297 sort, nactions,
2298 next);
2299 // adjust offset
2300 prob->change_bias(offset);
2301 }
2302 delete[] boundRecords;
2303 if (prob->tuning_) {
2304 double thisTime = CoinCpuTime();
2305 int droppedRows = prob->countEmptyRows() - startEmptyRows;
2306 int droppedColumns = prob->countEmptyCols() - startEmptyColumns;
2307 printf("CoinPresolveTwoxtwo(2048) - %d rows, %d columns dropped in time %g, total %g\n",
2308 droppedRows, droppedColumns, thisTime - startTime, thisTime - prob->startTime_);
2309 }
2310 return (next);
2311 }
postsolve(CoinPostsolveMatrix * prob) const2312 void twoxtwo_action::postsolve(CoinPostsolveMatrix *prob) const
2313 {
2314 const CoinBigIndex *const mcstrt = prob->mcstrt_;
2315 const int *const hincol = prob->hincol_;
2316 const int *const hrow = prob->hrow_;
2317 const double *colels = prob->colels_;
2318 CoinBigIndex *link = prob->link_;
2319 double *cost = prob->cost_;
2320
2321 // column type, bounds, solution, and status
2322 //const unsigned char *const integerType = prob->integerType_ ;
2323 double *clo = prob->clo_;
2324 double *cup = prob->cup_;
2325
2326 // row bounds
2327 double *rlo = prob->rlo_;
2328 double *rup = prob->rup_;
2329 double *sol = prob->sol_;
2330 double *rcosts = prob->rcosts_;
2331 double *rowacts = prob->acts_;
2332 double *dual = prob->rowduals_;
2333 double tolerance = prob->ztolzb_;
2334 const double maxmin = prob->maxmin_;
2335 for (int iAction = 0; iAction < nactions_; iAction++) {
2336 const action &boundRecord = actions_[iAction];
2337 int row1 = boundRecord.row;
2338 int icol = boundRecord.col;
2339 int otherCol = boundRecord.othercol;
2340 CoinBigIndex start = mcstrt[icol];
2341 CoinBigIndex nextEl = link[start];
2342 int row0;
2343 // first is otherCol
2344 double els0[2] = { 0.0, 0.0 };
2345 double els1[2] = { 0.0, 0.0 };
2346 if (hrow[start] == row1) {
2347 row0 = hrow[nextEl];
2348 els0[1] = colels[nextEl];
2349 els1[1] = colels[start];
2350 } else {
2351 row0 = hrow[start];
2352 els0[1] = colels[start];
2353 els1[1] = colels[nextEl];
2354 }
2355 nextEl = mcstrt[otherCol];
2356 for (CoinBigIndex j = 0; j < hincol[otherCol]; j++) {
2357 if (hrow[nextEl] == row0)
2358 els0[0] = colels[nextEl];
2359 else if (hrow[nextEl] == row1)
2360 els1[0] = colels[nextEl];
2361 nextEl = link[nextEl];
2362 }
2363 prob->setRowStatus(row1, CoinPrePostsolveMatrix::basic);
2364 // put stuff back
2365 rlo[row1] = boundRecord.lbound_row;
2366 rup[row1] = boundRecord.ubound_row;
2367 clo[icol] = boundRecord.lbound_col;
2368 cup[icol] = boundRecord.ubound_col;
2369 double oldCost = cost[icol];
2370 //double oldOtherCost=cost[otherCol];
2371 cost[icol] = boundRecord.cost_col;
2372 cost[otherCol] = boundRecord.cost_othercol;
2373 double els0real[2];
2374 double els1real[2];
2375 els0real[0] = els0[0];
2376 els0real[1] = els0[1];
2377 els1real[0] = els1[0];
2378 els1real[1] = els1[1];
2379 // make <= rows
2380 double rowUpper0 = rup[row0];
2381 if (rlo[row0] > -1.0e30) {
2382 rowUpper0 = -rlo[row0];
2383 els0[0] = -els0[0];
2384 els0[1] = -els0[1];
2385 }
2386 double rowUpper1 = rup[row1];
2387 bool swapSigns1 = false;
2388 if (rlo[row1] > -1.0e30) {
2389 swapSigns1 = true;
2390 rowUpper1 = -rlo[row1];
2391 els1[0] = -els1[0];
2392 els1[1] = -els1[1];
2393 }
2394 // compute feasible value for icol
2395 double valueOther = sol[otherCol];
2396 double value;
2397 // first see if at bound is OK
2398 bool lowerBoundPossible = clo[icol] > -1.0e30;
2399 value = clo[icol];
2400 if (lowerBoundPossible) {
2401 double value0 = els0[0] * valueOther + els0[1] * value;
2402 if (value0 > rowUpper0 + tolerance)
2403 lowerBoundPossible = false;
2404 double value1 = els1[0] * valueOther + els1[1] * value;
2405 if (value1 > rowUpper1 + tolerance)
2406 lowerBoundPossible = false;
2407 }
2408 bool upperBoundPossible = cup[icol] < 1.0e30;
2409 value = cup[icol];
2410 if (upperBoundPossible) {
2411 double value0 = els0[0] * valueOther + els0[1] * value;
2412 if (value0 > rowUpper0 + tolerance)
2413 upperBoundPossible = false;
2414 double value1 = els1[0] * valueOther + els1[1] * value;
2415 if (value1 > rowUpper1 + tolerance)
2416 upperBoundPossible = false;
2417 }
2418 if (lowerBoundPossible && cost[icol] >= 0.0) {
2419 // set to lower bound
2420 prob->setColumnStatus(icol, CoinPrePostsolveMatrix::atLowerBound);
2421 sol[icol] = clo[icol];
2422 rcosts[icol] = maxmin * cost[icol] - dual[row0] * els0real[1];
2423 } else if (upperBoundPossible && cost[icol] <= 0.0) {
2424 // set to upper bound
2425 prob->setColumnStatus(icol, CoinPrePostsolveMatrix::atUpperBound);
2426 sol[icol] = cup[icol];
2427 rcosts[icol] = maxmin * cost[icol] - dual[row0] * els0real[1];
2428 } else {
2429 // need to make basic
2430 // we shouldn't get here (at present) if zero cost
2431 assert(cost[icol]);
2432 double value0 = (rowUpper0 - els0[0] * valueOther) / els0[1];
2433 double value1 = (rowUpper1 - els1[0] * valueOther) / els1[1];
2434 //bool binding0=true;
2435 double value;
2436 if (cost[icol] > 0) {
2437 if (value0 > value1) {
2438 value = value0;
2439 } else {
2440 value = value1;
2441 //binding0=false;
2442 }
2443 } else {
2444 if (value0 < value1) {
2445 value = value0;
2446 } else {
2447 value = value1;
2448 //binding0=false;
2449 }
2450 }
2451 sol[icol] = value;
2452 #if 0
2453 printf("row %d status %d, row %d status %d, col %d status %d, col %d status %d - binding0 %c\n",
2454 row0,prob->getRowStatus(row0),
2455 row1,prob->getRowStatus(row1),
2456 otherCol,prob->getColumnStatus(otherCol),
2457 icol,prob->getColumnStatus(icol),binding0 ? 'T' : 'F');
2458 #endif
2459 if (prob->getColumnStatus(icol) == CoinPrePostsolveMatrix::basic) {
2460 //printf("col %d above was basic\n",icol);
2461 if (prob->getRowStatus(row0) != CoinPrePostsolveMatrix::basic) {
2462 // adjust dual
2463 dual[row0] = maxmin * ((cost[icol] - oldCost) / els0real[1]);
2464 }
2465 continue;
2466 }
2467 //if (binding0)
2468 //printf("Says row0 %d binding?\n",row0);
2469 prob->setColumnStatus(icol, CoinPrePostsolveMatrix::basic);
2470 rcosts[icol] = 0.0;
2471 //printf("row1 %d taken out of basis\n",row1);
2472 if (!swapSigns1) {
2473 prob->setRowStatus(row1, CoinPrePostsolveMatrix::atUpperBound);
2474 rowacts[row1] = rup[row1];
2475 } else {
2476 prob->setRowStatus(row1, CoinPrePostsolveMatrix::atLowerBound);
2477 rowacts[row1] = rlo[row1];
2478 }
2479 dual[row1] = maxmin * ((cost[icol] - oldCost) / els1real[1]);
2480 if (iAction == -1)
2481 abort();
2482 }
2483 }
2484 }
2485
2486 /* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
2487 */
2488