1/* $Id: CoinPresolveDupcol.cpp 1448 2011-06-19 15:34:41Z stefan $ */
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
26void 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
40namespace { // 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
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{ for (int cndx = 0 ; cndx < nlook ; ++cndx)
54 { int i = majcands[cndx] ;
55 PRESOLVEASSERT(majlens[i] > 0) ;
56
57 CoinBigIndex kcs = majstrts[i] ;
58 CoinBigIndex kce = kcs + majlens[i] ;
59
60 double value = 0.0 ;
61
62 for (CoinBigIndex k = kcs ; k < kce ; k++)
63 { int irow = minndxs[k] ;
64 value += minmuls[irow]*elems[k] ; }
65
66 majsums[cndx] = value ; }
67
68 return ; }
69
70
71void create_col (int col, int n, double *els,
72 CoinBigIndex *mcstrt, double *colels, int *hrow, int *link,
73 CoinBigIndex *free_listp)
74{
75 int *rows = reinterpret_cast<int *>(els+n) ;
76 CoinBigIndex free_list = *free_listp;
77 int xstart = NO_LINK;
78 for (int i=0; i<n; ++i) {
79 CoinBigIndex k = free_list;
80 assert(k >= 0) ;
81 free_list = link[free_list];
82 hrow[k] = rows[i];
83 colels[k] = els[i];
84 link[k] = xstart;
85 xstart = k;
86 }
87 mcstrt[col] = xstart;
88 *free_listp = free_list;
89}
90
91
92} // end unnamed file-local namespace
93
94
95
96const char *dupcol_action::name () const
97{
98 return ("dupcol_action");
99}
100
101
102/*
103 Original comment: This is just ekkredc5, adapted into the new framework.
104 The datasets scorpion.mps and allgrade.mps have duplicate columns.
105
106 In case you don't have your OSL manual handy, a somewhat more informative
107 explanation: We're looking for an easy-to-detect special case of linearly
108 dependent columns, where the coefficients of the duplicate columns are
109 exactly equal. The idea for locating such columns is to generate pseudo-
110 random weights for each row and then calculate the weighted sum of
111 coefficients of each column. Columns with equal sums are checked more
112 thoroughly.
113
114 Analysis of the situation says there are two major cases:
115 * If the columns have equal objective coefficients, we can combine
116 them.
117 * If the columns have unequal objective coefficients, we may be able to
118 fix one at bound. If the required bound doesn't exist, we have dual
119 infeasibility (hence one of primal infeasibility or unboundedness).
120
121 In the comments below are a few fragments of code from the original
122 routine. I don't think they make sense, but I've left them for the nonce in
123 case someone else recognises the purpose. -- lh, 040909 --
124*/
125
126const CoinPresolveAction
127 *dupcol_action::presolve (CoinPresolveMatrix *prob,
128 const CoinPresolveAction *next)
129{
130 double startTime = 0.0;
131 int startEmptyRows=0;
132 int startEmptyColumns = 0;
133 if (prob->tuning_) {
134 startTime = CoinCpuTime();
135 startEmptyRows = prob->countEmptyRows();
136 startEmptyColumns = prob->countEmptyCols();
137 }
138
139 double maxmin = prob->maxmin_ ;
140
141 double *colels = prob->colels_ ;
142 int *hrow = prob->hrow_ ;
143 CoinBigIndex *mcstrt = prob->mcstrt_ ;
144 int *hincol = prob->hincol_ ;
145 int ncols = prob->ncols_ ;
146 int nrows = prob->nrows_ ;
147
148 double *clo = prob->clo_ ;
149 double *cup = prob->cup_ ;
150 double *sol = prob->sol_ ;
151 double *rlo = prob->rlo_ ;
152 double *rup = prob->rup_ ;
153
154 // If all coefficients positive do more simply
155 bool allPositive=true;
156 double * rhs = prob->usefulRowDouble_; //new double[nrows];
157 CoinMemcpyN(rup,nrows,rhs);
158/*
159 Scan the columns for candidates, and write the indices into sort. We're not
160 interested in columns that are empty, prohibited, or integral.
161
162 Question: Should we exclude singletons, which are useful in other transforms?
163 Question: Why are we excluding integral columns?
164*/
165 // allow integral columns if asked for
166 bool allowIntegers = ( prob->presolveOptions_&1)!=0;
167 int *sort = prob->usefulColumnInt_; //new int[ncols] ;
168 int nlook = 0 ;
169 for (int j = 0 ; j < ncols ; j++) {
170 if (hincol[j] == 0) continue ;
171 // sort
172 CoinSort_2(hrow+mcstrt[j],hrow+mcstrt[j]+hincol[j],
173 colels+mcstrt[j]);
174 // check all positive and adjust rhs
175 if (allPositive) {
176 double lower = clo[j];
177 if (lower<cup[j]) {
178 for (int k=mcstrt[j];k<mcstrt[j]+hincol[j];k++) {
179 double value=colels[k];
180 if (value<0.0)
181 allPositive=false;
182 else
183 rhs[hrow[k]] -= lower*value;
184 }
185 } else {
186 for (int k=mcstrt[j];k<mcstrt[j]+hincol[j];k++) {
187 double value=colels[k];
188 rhs[hrow[k]] -= lower*value;
189 }
190 }
191 }
192 if (prob->colProhibited2(j)) continue ;
193 //#define PRESOLVE_INTEGER_DUPCOL
194#ifndef PRESOLVE_INTEGER_DUPCOL
195 if (prob->isInteger(j)&&!allowIntegers) continue ;
196#endif
197 sort[nlook++] = j ; }
198 if (nlook == 0)
199 { //delete[] sort ;
200 //delete [] rhs;
201 return (next) ; }
202/*
203 Prep: add the coefficients of each candidate column. To reduce false
204 positives, multiply each row by a `random' multiplier when forming the
205 sums. On return from compute_sums, sort and colsum are loaded with the
206 indices and column sums, respectively, of candidate columns. The pair of
207 arrays are then sorted by sum so that equal sums are adjacent.
208*/
209 double *colsum = prob->usefulColumnDouble_; //new double[ncols] ;
210 double *rowmul;
211 if (!prob->randomNumber_) {
212 rowmul = new double[nrows] ;
213 coin_init_random_vec(rowmul,nrows) ;
214 } else {
215 rowmul = prob->randomNumber_;
216 }
217 compute_sums(ncols,hincol,mcstrt,hrow,colels,rowmul,sort,colsum,nlook) ;
218 CoinSort_2(colsum,colsum+nlook,sort) ;
219/*
220 General prep --- unpack the various vectors we'll need, and allocate arrays
221 to record the results.
222*/
223 presolvehlink *clink = prob->clink_ ;
224
225 double *rowels = prob->rowels_ ;
226 int *hcol = prob->hcol_ ;
227 const CoinBigIndex *mrstrt = prob->mrstrt_ ;
228 int *hinrow = prob->hinrow_ ;
229
230 double *dcost = prob->cost_ ;
231
232 action *actions = new action [nlook] ;
233 int nactions = 0 ;
234# ifdef ZEROFAULT
235 memset(actions,0,nlook*sizeof(action)) ;
236# endif
237
238 int *fixed_down = new int[nlook] ;
239 int nfixed_down = 0 ;
240 int *fixed_up = new int[nlook] ;
241 int nfixed_up = 0 ;
242
243#if 0
244 // Excluded in the original routine. I'm guessing it's excluded because
245 // it's just not cost effective to worry about this. -- lh, 040908 --
246
247 // It may be the case that several columns are duplicate.
248 // If not all have the same cost, then we have to make sure
249 // that we set the most expensive one to its minimum
250 // now sort in each class by cost
251 {
252 double dval = colsum[0] ;
253 int first = 0 ;
254 for (int jj = 1; jj < nlook; jj++) {
255 while (colsum[jj]==dval)
256 jj++ ;
257
258 if (first + 1 < jj) {
259 double buf[jj - first] ;
260 for (int i=first; i<jj; ++i)
261 buf[i-first] = dcost[sort[i]]*maxmin ;
262
263 CoinSort_2(buf,buf+jj-first,sort+first) ;
264 //ekk_sortonDouble(buf,&sort[first],jj-first) ;
265 }
266 }
267 }
268#endif
269 // We will get all min/max but only if needed
270 bool gotStuff=false;
271/*
272 Original comment: It appears to be the case that this loop is finished,
273 there may still be duplicate cols left. I haven't done anything
274 about that yet.
275
276 Open the main loop to compare column pairs. We'll compare sort[jj] to
277 sort[tgt]. This allows us to accumulate multiple columns into one. But
278 we don't manage all-pairs comparison when we can't combine columns.
279
280 We can quickly dismiss pairs which have unequal sums or lengths.
281*/
282 int isorted = -1 ;
283 int tgt = 0 ;
284 for (int jj = 1 ; jj < nlook ; jj++)
285 { if (colsum[jj] != colsum[jj-1]) {
286 tgt = jj; // Must update before continuing
287 continue ;
288 }
289
290 int j2 = sort[jj] ;
291 int j1 = sort[tgt] ;
292 int len2 = hincol[j2] ;
293 int len1 = hincol[j1] ;
294
295 if (len2 != len1)
296 { tgt = jj ;
297 continue ; }
298/*
299 The final test: sort the columns by row index and compare each index and
300 coefficient.
301*/
302 CoinBigIndex kcs = mcstrt[j2] ;
303 CoinBigIndex kce = kcs+hincol[j2] ;
304 int ishift = mcstrt[j1]-kcs ;
305
306 if (len1 > 1 && isorted < j1)
307 { CoinSort_2(hrow+mcstrt[j1],hrow+mcstrt[j1]+len1,
308 colels+mcstrt[j1]) ;
309 isorted = j1 ; }
310 if (len2 > 1 && isorted < j2)
311 { CoinSort_2(hrow+kcs,hrow+kcs+len2,colels+kcs) ;
312 isorted = j2 ; }
313
314 CoinBigIndex k ;
315 for (k = kcs ; k < kce ; k++)
316 { if (hrow[k] != hrow[k+ishift] || colels[k] != colels[k+ishift])
317 { break ; } }
318 if (k != kce)
319 { tgt = jj ;
320 continue ; }
321/*
322 These really are duplicate columns. Grab values for convenient reference.
323 Convert the objective coefficients for minimization.
324*/
325 double clo1 = clo[j1] ;
326 double cup1 = cup[j1] ;
327 double clo2 = clo[j2] ;
328 double cup2 = cup[j2] ;
329 double c1 = dcost[j1]*maxmin ;
330 double c2 = dcost[j2]*maxmin ;
331 PRESOLVEASSERT(!(clo1 == cup1 || clo2 == cup2)) ;
332 // Get reasonable bounds on sum of two variables
333 double lowerBound=-COIN_DBL_MAX;
334 double upperBound=COIN_DBL_MAX;
335 // For now only if lower bounds are zero
336 if (!clo1&&!clo2) {
337 // Only need bounds if c1 != c2
338 if (c1!=c2) {
339 if (!allPositive) {
340#if 0
341
342 for (k=kcs;k<kce;k++) {
343 int iRow = hrow[k];
344 bool posinf = false;
345 bool neginf = false;
346 double maxup = 0.0;
347 double maxdown = 0.0;
348
349 // compute sum of all bounds except for j1,j2
350 CoinBigIndex kk;
351 CoinBigIndex kre = mrstrt[iRow]+hinrow[iRow];
352 double value1=0.0;
353 for (kk=mrstrt[iRow]; kk<kre; kk++) {
354 int col = hcol[kk];
355 if (col == j1||col==j2) {
356 value1=rowels[kk];
357 continue;
358 }
359 double coeff = rowels[kk];
360 double lb = clo[col];
361 double ub = cup[col];
362
363 if (coeff > 0.0) {
364 if (PRESOLVE_INF <= ub) {
365 posinf = true;
366 if (neginf)
367 break; // pointless
368 } else {
369 maxup += ub * coeff;
370 }
371 if (lb <= -PRESOLVE_INF) {
372 neginf = true;
373 if (posinf)
374 break; // pointless
375 } else {
376 maxdown += lb * coeff;
377 }
378 } else {
379 if (PRESOLVE_INF <= ub) {
380 neginf = true;
381 if (posinf)
382 break; // pointless
383 } else {
384 maxdown += ub * coeff;
385 }
386 if (lb <= -PRESOLVE_INF) {
387 posinf = true;
388 if (neginf)
389 break; // pointless
390 } else {
391 maxup += lb * coeff;
392 }
393 }
394 }
395
396 if (kk==kre) {
397 assert (value1);
398 if (value1>1.0e-5) {
399 if (!neginf&&rup[iRow]<1.0e10)
400 if (upperBound*value1>rup[iRow]-maxdown)
401 upperBound = (rup[iRow]-maxdown)/value1;
402 if (!posinf&&rlo[iRow]>-1.0e10)
403 if (lowerBound*value1<rlo[iRow]-maxup)
404 lowerBound = (rlo[iRow]-maxup)/value1;
405 } else if (value1<-1.0e-5) {
406 if (!neginf&&rup[iRow]<1.0e10)
407 if (lowerBound*value1>rup[iRow]-maxdown) {
408#ifndef NDEBUG
409 double x=lowerBound;
410#endif
411 lowerBound = (rup[iRow]-maxdown)/value1;
412 assert (lowerBound == CoinMax(x,(rup[iRow]-maxdown)/value1));
413 }
414 if (!posinf&&rlo[iRow]>-1.0e10)
415 if (upperBound*value1<rlo[iRow]-maxup) {
416#ifndef NDEBUG
417 double x=upperBound;
418#endif
419 upperBound = (rlo[iRow]-maxup)/value1;
420 assert(upperBound == CoinMin(x,(rlo[iRow]-maxup)/value1));
421 }
422 }
423 }
424 }
425 double l=lowerBound;
426 double u=upperBound;
427#endif
428 if (!gotStuff) {
429 prob->recomputeSums(-1); // get min max
430 gotStuff=true;
431 }
432 int positiveInf=0;
433 int negativeInf=0;
434 double lo=0;
435 double up=0.0;
436 if (clo1<-PRESOLVE_INF)
437 negativeInf++;
438 else
439 lo+=clo1;
440 if (clo2<-PRESOLVE_INF)
441 negativeInf++;
442 else
443 lo+=clo2;
444 if (cup1>PRESOLVE_INF)
445 positiveInf++;
446 else
447 up+=cup1;
448 if (cup2>PRESOLVE_INF)
449 positiveInf++;
450 else
451 up+=cup2;
452 for (k=kcs;k<kce;k++) {
453 int iRow = hrow[k];
454 double value = colels[k];
455 int pInf = (value>0.0) ? positiveInf : negativeInf;
456 int nInf = (value>0.0) ? negativeInf : positiveInf;
457 int posinf = prob->infiniteUp_[iRow]-pInf;
458 int neginf = prob->infiniteDown_[iRow]-nInf;
459 if (posinf>0&&neginf>0)
460 continue; // this row can't bound
461 double maxup = prob->sumUp_[iRow];
462 double maxdown = prob->sumDown_[iRow];
463
464 if (value>0.0) {
465 maxdown -= value*lo;
466 maxup -= value*up;
467 } else {
468 maxdown -= value*up;
469 maxup -= value*lo;
470 }
471 if (value>1.0e-5) {
472 if (!neginf&&rup[iRow]<1.0e10)
473 if (upperBound*value>rup[iRow]-maxdown)
474 upperBound = (rup[iRow]-maxdown)/value;
475 if (!posinf&&rlo[iRow]>-1.0e10)
476 if (lowerBound*value<rlo[iRow]-maxup)
477 lowerBound = (rlo[iRow]-maxup)/value;
478 } else if (value<-1.0e-5) {
479 if (!neginf&&rup[iRow]<1.0e10)
480 if (lowerBound*value>rup[iRow]-maxdown) {
481 lowerBound = (rup[iRow]-maxdown)/value;
482 }
483 if (!posinf&&rlo[iRow]>-1.0e10)
484 if (upperBound*value<rlo[iRow]-maxup) {
485 upperBound = (rlo[iRow]-maxup)/value;
486 }
487 }
488 }
489 //assert (fabs(l-lowerBound)<1.0e-5&&fabs(u-upperBound)<1.0e-5);
490 } else {
491 // can do faster
492 lowerBound=0.0;
493 for (k=kcs;k<kce;k++) {
494 int iRow = hrow[k];
495 double value=colels[k];
496 if (upperBound*value>rhs[iRow])
497 upperBound = rhs[iRow]/value;
498 }
499 }
500 }
501 // relax a bit
502 upperBound -= 1.0e-9;
503 } else {
504 // Not sure what to do so give up
505 continue;
506 }
507/*
508 There are two main cases: The objective coefficients are equal or unequal.
509
510 For equal objective coefficients c1 == c2, we can combine the columns by
511 making the substitution x<j1> = x'<j1> - x<j2>. This will eliminate column
512 sort[jj] = j2 and leave the new variable x' in column sort[tgt] = j1. tgt
513 doesn't move.
514*/
515 if (c1 == c2)
516 {
517#ifdef PRESOLVE_INTEGER_DUPCOL
518 if (!allowIntegers) {
519 if (prob->isInteger(j1)) {
520 if (!prob->isInteger(j2)) {
521 if (cup2 < upperBound) //if (!prob->colInfinite(j2))
522 continue;
523 else
524 cup2 = COIN_DBL_MAX;
525 }
526 } else if (prob->isInteger(j2)) {
527 if (cup1 < upperBound) //if (!prob->colInfinite(j1))
528 continue;
529 else
530 cup1 = COIN_DBL_MAX;
531 }
532 //printf("TakingINTeq\n");
533 }
534#endif
535/*
536 As far as the presolved lp, there's no difference between columns. But we
537 need this relation to hold in order to guarantee that we can split the
538 value of the combined column during postsolve without damaging the basis.
539 (The relevant case is when the combined column is basic --- we need to be
540 able to retain column j1 in the basis and make column j2 nonbasic.)
541*/
542 if (!(clo2+cup1 <= clo1+cup2))
543 { CoinSwap(j1,j2) ;
544 CoinSwap(clo1,clo2) ;
545 CoinSwap(cup1,cup2) ;
546 tgt = jj ; }
547/*
548 Create the postsolve action before we start to modify the columns.
549*/
550 PRESOLVE_STMT(printf("DUPCOL: (%d,%d) %d += %d\n",j1,j2,j1,j2)) ;
551 PRESOLVE_DETAIL_PRINT(printf("pre_dupcol %dC %dC E\n",j2,j1));
552
553 action *s = &actions[nactions++] ;
554 s->thislo = clo[j2] ;
555 s->thisup = cup[j2] ;
556 s->lastlo = clo[j1] ;
557 s->lastup = cup[j1] ;
558 s->ithis = j2 ;
559 s->ilast = j1 ;
560 s->nincol = hincol[j2] ;
561 s->colels = presolve_dupmajor(colels,hrow,hincol[j2],mcstrt[j2]) ;
562/*
563 Combine the columns into column j1. Upper and lower bounds and solution
564 simply add, and the coefficients are unchanged.
565
566 I'm skeptical of pushing a bound to infinity like this, but leave it for now.
567 -- lh, 040908 --
568*/
569 clo1 += clo2 ;
570 if (clo1 < -1.0e20)
571 { clo1 = -PRESOLVE_INF ; }
572 clo[j1] = clo1 ;
573 cup1 += cup2 ;
574 if (cup1 > 1.0e20)
575 { cup1 = PRESOLVE_INF ; }
576 cup[j1] = cup1 ;
577 if (sol)
578 { sol[j1] += sol[j2] ; }
579 if (prob->colstat_)
580 { if (prob->getColumnStatus(j1) == CoinPrePostsolveMatrix::basic ||
581 prob->getColumnStatus(j2) == CoinPrePostsolveMatrix::basic)
582 { prob->setColumnStatus(j1,CoinPrePostsolveMatrix::basic); } }
583/*
584 Empty column j2.
585*/
586 dcost[j2] = 0.0 ;
587 if (sol)
588 { sol[j2] = clo2 ; }
589 CoinBigIndex k2cs = mcstrt[j2] ;
590 CoinBigIndex k2ce = k2cs + hincol[j2] ;
591 for (CoinBigIndex k = k2cs ; k < k2ce ; ++k)
592 { presolve_delete_from_row(hrow[k],j2,mrstrt,hinrow,hcol,rowels) ; }
593 hincol[j2] = 0 ;
594 PRESOLVE_REMOVE_LINK(clink,j2) ;
595 continue ; }
596/*
597 Unequal reduced costs. In this case, we may be able to fix one of the columns
598 or prove dual infeasibility. Given column a_k, duals y, objective
599 coefficient c_k, the reduced cost cbar_k = c_k - dot(y,a_k). Given
600 a_1 = a_2, substitute for dot(y,a_1) in the relation for cbar_2 to get
601 cbar_2 = (c_2 - c_1) + cbar_1
602 Independent elements here are variable bounds l_k, u_k, and difference
603 (c_2 - c_1). Infinite bounds for l_k, u_k will constrain the sign of cbar_k.
604 Assume minimization. If you do the case analysis, you find these cases of
605 interest:
606
607 l_1 u_1 l_2 u_2 cbar_1 c_2-c_1 cbar_2 result
608
609 A any finite -inf any <= 0 > 0 <= 0 x_1 -> NBUB
610 B -inf any any finite <= 0 < 0 < 0 x_2 -> NBUB
611
612 C finite any any +inf >= 0 < 0 >= 0 x_1 -> NBLB
613 D any +inf finite any >= 0 > 0 >= 0 x_2 -> NBLB
614
615 E -inf any any +inf <= 0 < 0 >= 0 dual infeas
616 F any inf -inf any >= 0 > 0 <= 0 dual infeas
617
618 G any finite finite any > 0 no inference
619 H finite any any finite < 0 no inference
620
621 The cases labelled dual infeasible are primal unbounded.
622
623 To keep the code compact, we'll always aim to take x_2 to bound. In the cases
624 where x_1 should go to bound, we'll swap. The implementation is boolean
625 algebra. Define bits for infinite bounds and (c_2 > c_1), then look for the
626 correct patterns.
627*/
628 else
629 { int minterm = 0 ;
630#ifdef PRESOLVE_INTEGER_DUPCOL
631 if (!allowIntegers) {
632 if (c2 > c1) {
633 if (cup1 < upperBound/*!prob->colInfinite(j1)*/ && (prob->isInteger(j1)||prob->isInteger(j2)))
634 continue ;
635 } else {
636 if (cup2 < upperBound/*!prob->colInfinite(j2)*/ && (prob->isInteger(j1)||prob->isInteger(j2)))
637 continue ;
638 }
639 //printf("TakingINTne\n");
640 }
641#endif
642 bool swapped = false ;
643#if PRESOLVE_DEBUG
644 printf("bounds %g %g\n",lowerBound,upperBound);
645#endif
646 if (c2 > c1) minterm |= 1<<0 ;
647 if (cup2 >= PRESOLVE_INF/*prob->colInfinite(j2)*/) minterm |= 1<<1 ;
648 if (clo2 <= -PRESOLVE_INF) minterm |= 1<<2 ;
649 if (cup1 >= PRESOLVE_INF/*prob->colInfinite(j1)*/) minterm |= 1<<3 ;
650 if (clo1 <= -PRESOLVE_INF) minterm |= 1<<4 ;
651 // for now be careful - just one special case
652 if (!clo1&&!clo2) {
653 if (c2 > c1 && cup1 >= upperBound)
654 minterm |= 1<<3;
655 else if (c2 < c1 && cup2 >= upperBound)
656 minterm |= 1<<1;
657 }
658/*
659 The most common case in a well-formed system should be no inference. We're
660 looking for x00x1 (case G) and 0xx00 (case H). This is where we have the
661 potential to miss inferences: If there are three or more columns with the
662 same sum, sort[tgt] == j1 will only be compared to the second in the
663 group.
664*/
665 if ((minterm&0x0d) == 0x1 || (minterm&0x13) == 0)
666 { tgt = jj ;
667 continue ; }
668/*
669 Next remove the unbounded cases, 1xx10 and x11x1.
670*/
671 if ((minterm&0x13) == 0x12 || (minterm&0x0d) == 0x0d)
672 { prob->setStatus(2) ;
673 PRESOLVE_STMT(printf("DUPCOL: (%d,%d) Unbounded\n",j1,j2)) ;
674 break ; }
675/*
676 With the no inference and unbounded cases removed, all that's left are the
677 cases where we can push a variable to bound. Swap if necessary (x01x1 or
678 0xx10) so that we're always fixing index j2. This means that column
679 sort[tgt] = j1 will be fixed. Unswapped, we fix column sort[jj] = j2.
680*/
681 if ((minterm&0x0d) == 0x05 || (minterm&0x13) == 0x02)
682 { CoinSwap(j1, j2) ;
683 CoinSwap(clo1, clo2) ;
684 CoinSwap(cup1, cup2) ;
685 CoinSwap(c1, c2) ;
686 int tmp1 = minterm&0x18 ;
687 int tmp2 = minterm&0x06 ;
688 int tmp3 = minterm&0x01 ;
689 minterm = (tmp1>>2)|(tmp2<<2)|(tmp3^0x01) ;
690 swapped = true ; }
691/*
692 Force x_2 to upper bound? (Case B, boolean 1X100, where X == don't care.)
693*/
694 if ((minterm&0x13) == 0x10)
695 { fixed_up[nfixed_up++] = j2 ;
696 PRESOLVE_STMT(printf("DUPCOL: (%d,%d) %d -> NBUB\n",j1,j2,j2)) ;
697 if (prob->colstat_)
698 { if (prob->getColumnStatus(j1) == CoinPrePostsolveMatrix::basic ||
699 prob->getColumnStatus(j2) == CoinPrePostsolveMatrix::basic)
700 { prob->setColumnStatus(j1,CoinPrePostsolveMatrix::basic) ; }
701 prob->setColumnStatus(j2,CoinPrePostsolveMatrix::atUpperBound) ; }
702 if (sol)
703 { double delta2 = cup2-sol[j2] ;
704 sol[j2] = cup2 ;
705 sol[j1] -= delta2 ; }
706 if (swapped)
707 { tgt = jj ; }
708 continue ; }
709/*
710 Force x_2 to lower bound? (Case C, boolean X1011.)
711*/
712 if ((minterm&0x0d) == 0x09)
713 { fixed_down[nfixed_down++] = j2 ;
714 PRESOLVE_STMT(printf("DUPCOL: (%d,%d) %d -> NBLB\n",j1,j2,j2)) ;
715 if (prob->colstat_)
716 { if (prob->getColumnStatus(j1) == CoinPrePostsolveMatrix::basic ||
717 prob->getColumnStatus(j2) == CoinPrePostsolveMatrix::basic)
718 { prob->setColumnStatus(j1,CoinPrePostsolveMatrix::basic) ; }
719 prob->setColumnStatus(j2,CoinPrePostsolveMatrix::atLowerBound) ; }
720 if (sol)
721 { double delta2 = clo2-sol[j2] ;
722 sol[j2] = clo2 ;
723 sol[j1] -= delta2 ; }
724 if (swapped)
725 { tgt = jj ; }
726 continue ; } }
727/*
728 We should never reach this point in the loop --- all cases force a new
729 iteration or loop termination. If we get here, something happened that we
730 didn't anticipate.
731*/
732 PRESOLVE_STMT(printf("DUPCOL: (%d,%d) UNEXPECTED!\n",j1,j2)) ; }
733/*
734 What's left? Deallocate vectors, and call make_fixed_action to handle any
735 variables that were fixed to bound.
736*/
737 if (rowmul != prob->randomNumber_)
738 delete[] rowmul ;
739 //delete[] colsum ;
740 //delete[] sort ;
741 //delete [] rhs;
742
743# if PRESOLVE_SUMMARY || PRESOLVE_DEBUG
744 if (nactions+nfixed_down+nfixed_up > 0)
745 { printf("DUPLICATE COLS: %d combined, %d lb, %d ub\n",
746 nactions,nfixed_down,nfixed_up) ; }
747# endif
748 if (nactions)
749 { next = new dupcol_action(nactions,CoinCopyOfArray(actions,nactions),next) ;
750 // we can't go round again in integer
751 prob->presolveOptions_ |= 0x80000000;
752}
753 deleteAction(actions,action*) ;
754
755 if (nfixed_down)
756 { next =
757 make_fixed_action::presolve(prob,fixed_down,nfixed_down,true,next) ; }
758 delete[]fixed_down ;
759
760 if (nfixed_up)
761 { next =
762 make_fixed_action::presolve(prob,fixed_up,nfixed_up,false,next) ; }
763 delete[]fixed_up ;
764
765 if (prob->tuning_) {
766 double thisTime=CoinCpuTime();
767 int droppedRows = prob->countEmptyRows() - startEmptyRows ;
768 int droppedColumns = prob->countEmptyCols() - startEmptyColumns;
769 printf("CoinPresolveDupcol(128) - %d rows, %d columns dropped in time %g, total %g\n",
770 droppedRows,droppedColumns,thisTime-startTime,thisTime-prob->startTime_);
771 }
772 return (next) ; }
773
774
775void dupcol_action::postsolve(CoinPostsolveMatrix *prob) const
776{
777 const action *const actions = actions_;
778 const int nactions = nactions_;
779
780 double *clo = prob->clo_;
781 double *cup = prob->cup_;
782
783 double *sol = prob->sol_;
784 double *dcost = prob->cost_;
785
786 double *colels = prob->colels_;
787 int *hrow = prob->hrow_;
788 CoinBigIndex *mcstrt = prob->mcstrt_;
789 int *hincol = prob->hincol_;
790 int *link = prob->link_;
791
792 double *rcosts = prob->rcosts_;
793 double tolerance = prob->ztolzb_;
794
795 for (const action *f = &actions[nactions-1]; actions<=f; f--) {
796 int icol = f->ithis; // was fixed
797 int icol2 = f->ilast; // was kept
798
799 dcost[icol] = dcost[icol2];
800 clo[icol] = f->thislo;
801 cup[icol] = f->thisup;
802 clo[icol2] = f->lastlo;
803 cup[icol2] = f->lastup;
804
805 create_col(icol,f->nincol,f->colels,mcstrt,colels,hrow,link,
806 &prob->free_list_) ;
807# if PRESOLVE_CONSISTENCY
808 presolve_check_free_list(prob) ;
809# endif
810 // hincol[icol] = hincol[icol2]; // right? - no - has to match number in create_col
811 hincol[icol] = f->nincol;
812
813 double l_j = f->thislo;
814 double u_j = f->thisup;
815 double l_k = f->lastlo;
816 double u_k = f->lastup;
817 double x_k_sol = sol[icol2];
818 PRESOLVE_DETAIL_PRINT(printf("post icol %d %g %g %g icol2 %d %g %g %g\n",
819 icol,clo[icol],sol[icol],cup[icol],
820 icol2,clo[icol2],sol[icol2],cup[icol2]));
821 if (l_j>-PRESOLVE_INF&& x_k_sol-l_j>=l_k-tolerance&&x_k_sol-l_j<=u_k+tolerance) {
822 // j at lb, leave k
823 prob->setColumnStatus(icol,CoinPrePostsolveMatrix::atLowerBound);
824 sol[icol] = l_j;
825 sol[icol2] = x_k_sol - sol[icol];
826 } else if (u_j<PRESOLVE_INF&& x_k_sol-u_j>=l_k-tolerance&&x_k_sol-u_j<=u_k+tolerance) {
827 // j at ub, leave k
828 prob->setColumnStatus(icol,CoinPrePostsolveMatrix::atUpperBound);
829 sol[icol] = u_j;
830 sol[icol2] = x_k_sol - sol[icol];
831 } else if (l_k>-PRESOLVE_INF&& x_k_sol-l_k>=l_j-tolerance&&x_k_sol-l_k<=u_j+tolerance) {
832 // k at lb make j basic
833 prob->setColumnStatus(icol,prob->getColumnStatus(icol2));
834 sol[icol2] = l_k;
835 sol[icol] = x_k_sol - l_k;
836 prob->setColumnStatus(icol2,CoinPrePostsolveMatrix::atLowerBound);
837 } else if (u_k<PRESOLVE_INF&& x_k_sol-u_k>=l_j-tolerance&&x_k_sol-u_k<=u_j+tolerance) {
838 // k at ub make j basic
839 prob->setColumnStatus(icol,prob->getColumnStatus(icol2));
840 sol[icol2] = u_k;
841 sol[icol] = x_k_sol - u_k;
842 prob->setColumnStatus(icol2,CoinPrePostsolveMatrix::atUpperBound);
843 } else {
844 // both free! superbasic time
845 sol[icol] = 0.0; // doesn't matter
846 prob->setColumnStatus(icol,CoinPrePostsolveMatrix::isFree);
847 }
848 PRESOLVE_DETAIL_PRINT(printf("post2 icol %d %g icol2 %d %g\n",
849 icol,sol[icol],
850 icol2,sol[icol2]));
851 // row activity doesn't change
852 // dj of both variables is the same
853 rcosts[icol] = rcosts[icol2];
854 // leave until destructor
855 // deleteAction(f->colels,double *);
856
857# if PRESOLVE_DEBUG
858 const double ztolzb = prob->ztolzb_;
859 if (! (clo[icol] - ztolzb <= sol[icol] && sol[icol] <= cup[icol] + ztolzb))
860 printf("BAD DUPCOL BOUNDS: %g %g %g\n", clo[icol], sol[icol], cup[icol]);
861 if (! (clo[icol2] - ztolzb <= sol[icol2] && sol[icol2] <= cup[icol2] + ztolzb))
862 printf("BAD DUPCOL BOUNDS: %g %g %g\n", clo[icol2], sol[icol2], cup[icol2]);
863# endif
864 }
865 // leave until desctructor
866 // deleteAction(actions_,action *);
867}
868
869dupcol_action::~dupcol_action()
870{
871 for (int i = nactions_-1; i >= 0; --i) {
872 deleteAction(actions_[i].colels, double *);
873 }
874 deleteAction(actions_, action*);
875}
876
877
878
879/*
880 Routines for duplicate rows. This is definitely unfinished --- there's no
881 postsolve action.
882*/
883
884const char *duprow_action::name () const
885{
886 return ("duprow_action");
887}
888
889// This is just ekkredc4, adapted into the new framework.
890/*
891 I've made minimal changes for compatibility with dupcol: An initial scan to
892 accumulate rows of interest in sort.
893 -- lh, 040909 --
894*/
895const CoinPresolveAction
896 *duprow_action::presolve (CoinPresolveMatrix *prob,
897 const CoinPresolveAction *next)
898{
899 double startTime = 0.0;
900 int startEmptyRows=0;
901 int startEmptyColumns = 0;
902 if (prob->tuning_) {
903 startTime = CoinCpuTime();
904 startEmptyRows = prob->countEmptyRows();
905 startEmptyColumns = prob->countEmptyCols();
906 }
907 double *rowels = prob->rowels_;
908 int *hcol = prob->hcol_;
909 CoinBigIndex *mrstrt = prob->mrstrt_;
910 int *hinrow = prob->hinrow_;
911 int ncols = prob->ncols_;
912 int nrows = prob->nrows_;
913
914/*
915 Scan the rows for candidates, and write the indices into sort. We're not
916 interested in rows that are empty or prohibited.
917
918 Question: Should we exclude singletons, which are useful in other transforms?
919 Question: Why are we excluding integral columns?
920*/
921 int *sort = new int[nrows] ;
922 int nlook = 0 ;
923 for (int i = 0 ; i < nrows ; i++)
924 { if (hinrow[i] == 0) continue ;
925 if (prob->rowProhibited2(i)) continue ;
926 // sort
927 CoinSort_2(hcol+mrstrt[i],hcol+mrstrt[i]+hinrow[i],
928 rowels+mrstrt[i]);
929 sort[nlook++] = i ; }
930 if (nlook == 0)
931 { delete[] sort ;
932 return (next) ; }
933
934 double * workrow = new double[nrows+1];
935
936 double * workcol;
937 if (!prob->randomNumber_) {
938 workcol = new double[ncols+1];
939 coin_init_random_vec(workcol, ncols);
940 } else {
941 workcol = prob->randomNumber_;
942 }
943 compute_sums(nrows,hinrow,mrstrt,hcol,rowels,workcol,sort,workrow,nlook);
944 CoinSort_2(workrow,workrow+nlook,sort);
945
946 double *rlo = prob->rlo_;
947 double *rup = prob->rup_;
948
949 int nuseless_rows = 0;
950 bool fixInfeasibility = (prob->presolveOptions_&16384)!=0;
951 bool allowIntersection = ( prob->presolveOptions_&16)!=0;
952 double tolerance = prob->feasibilityTolerance_;
953
954 double dval = workrow[0];
955 for (int jj = 1; jj < nlook; jj++) {
956 if (workrow[jj]==dval) {
957 int ithis=sort[jj];
958 int ilast=sort[jj-1];
959 CoinBigIndex krs = mrstrt[ithis];
960 CoinBigIndex kre = krs + hinrow[ithis];
961 if (hinrow[ithis] == hinrow[ilast]) {
962 int ishift = mrstrt[ilast] - krs;
963 CoinBigIndex k;
964 for (k=krs;k<kre;k++) {
965 if (hcol[k] != hcol[k+ishift] ||
966 rowels[k] != rowels[k+ishift]) {
967 break;
968 }
969 }
970 if (k == kre) {
971 /* now check rhs to see what is what */
972 double rlo1=rlo[ilast];
973 double rup1=rup[ilast];
974 double rlo2=rlo[ithis];
975 double rup2=rup[ithis];
976
977 int idelete=-1;
978 if (rlo1<=rlo2) {
979 if (rup2<=rup1) {
980 /* this is strictly tighter than last */
981 idelete=ilast;
982 PRESOLVE_DETAIL_PRINT(printf("pre_duprow %dR %dR E\n",ilast,ithis));
983 } else if (fabs(rlo1-rlo2)<1.0e-12) {
984 /* last is strictly tighter than this */
985 idelete=ithis;
986 PRESOLVE_DETAIL_PRINT(printf("pre_duprow %dR %dR E\n",ithis,ilast));
987 // swap so can carry on deleting
988 sort[jj-1]=ithis;
989 sort[jj]=ilast;
990 } else {
991 if (rup1<rlo2-tolerance&&!fixInfeasibility) {
992 // infeasible
993 prob->status_|= 1;
994 // wrong message - correct if works
995 prob->messageHandler()->message(COIN_PRESOLVE_ROWINFEAS,
996 prob->messages())
997 <<ithis
998 <<rlo[ithis]
999 <<rup[ithis]
1000 <<CoinMessageEol;
1001 break;
1002 } else if (allowIntersection/*||fabs(rup1-rlo2)<tolerance*/) {
1003 /* overlapping - could merge */
1004#ifdef CLP_INVESTIGATE7
1005 printf("overlapping duplicate row %g %g, %g %g\n",
1006 rlo1,rup1,rlo2,rup2);
1007# endif
1008 // pretend this is stricter than last
1009 idelete=ilast;
1010 PRESOLVE_DETAIL_PRINT(printf("pre_duprow %dR %dR E\n",ilast,ithis));
1011 rup[ithis]=rup1;
1012 }
1013 }
1014 } else {
1015 // rlo1>rlo2
1016 if (rup1<=rup2) {
1017 /* last is strictly tighter than this */
1018 idelete=ithis;
1019 PRESOLVE_DETAIL_PRINT(printf("pre_duprow %dR %dR E\n",ithis,ilast));
1020 // swap so can carry on deleting
1021 sort[jj-1]=ithis;
1022 sort[jj]=ilast;
1023 } else {
1024 /* overlapping - could merge */
1025 // rlo1>rlo2
1026 // rup1>rup2
1027 if (rup2<rlo1-tolerance&&!fixInfeasibility) {
1028 // infeasible
1029 prob->status_|= 1;
1030 // wrong message - correct if works
1031 prob->messageHandler()->message(COIN_PRESOLVE_ROWINFEAS,
1032 prob->messages())
1033 <<ithis
1034 <<rlo[ithis]
1035 <<rup[ithis]
1036 <<CoinMessageEol;
1037 break;
1038 } else if (allowIntersection/*||fabs(rup2-rlo1)<tolerance*/) {
1039#ifdef CLP_INVESTIGATE7
1040 printf("overlapping duplicate row %g %g, %g %g\n",
1041 rlo1,rup1,rlo2,rup2);
1042# endif
1043 // pretend this is stricter than last
1044 idelete=ilast;
1045 PRESOLVE_DETAIL_PRINT(printf("pre_duprow %dR %dR E\n",ilast,ithis));
1046 rlo[ithis]=rlo1;
1047 }
1048 }
1049 }
1050 if (idelete>=0)
1051 sort[nuseless_rows++]=idelete;
1052 }
1053 }
1054 }
1055 dval=workrow[jj];
1056 }
1057
1058 delete[]workrow;
1059 if(workcol != prob->randomNumber_)
1060 delete[]workcol;
1061
1062
1063 if (nuseless_rows) {
1064# if PRESOLVE_SUMMARY
1065 printf("DUPLICATE ROWS: %d\n", nuseless_rows);
1066# endif
1067 next = useless_constraint_action::presolve(prob,
1068 sort, nuseless_rows,
1069 next);
1070 }
1071 delete[]sort;
1072
1073 if (prob->tuning_) {
1074 double thisTime=CoinCpuTime();
1075 int droppedRows = prob->countEmptyRows() - startEmptyRows ;
1076 int droppedColumns = prob->countEmptyCols() - startEmptyColumns;
1077 printf("CoinPresolveDuprow(256) - %d rows, %d columns dropped in time %g, total %g\n",
1078 droppedRows,droppedColumns,thisTime-startTime,thisTime-prob->startTime_);
1079 }
1080 return (next);
1081}
1082
1083void duprow_action::postsolve(CoinPostsolveMatrix *) const
1084{
1085 printf("STILL NO POSTSOLVE FOR DUPROW!\n");
1086 abort();
1087}
1088
1089
1090
1091/*
1092 Routines for gub rows. This is definitely unfinished --- there's no
1093 postsolve action.
1094*/
1095
1096const char *gubrow_action::name () const
1097{
1098 return ("gubrow_action");
1099}
1100
1101
1102const CoinPresolveAction
1103 *gubrow_action::presolve (CoinPresolveMatrix *prob,
1104 const CoinPresolveAction *next)
1105{
1106 double startTime = 0.0;
1107 int droppedElements=0;
1108 int affectedRows=0;
1109 if (prob->tuning_) {
1110 startTime = CoinCpuTime();
1111 }
1112 double *rowels = prob->rowels_;
1113 int *hcol = prob->hcol_;
1114 CoinBigIndex *mrstrt = prob->mrstrt_;
1115 int *hinrow = prob->hinrow_;
1116 double *colels = prob->colels_ ;
1117 int *hrow = prob->hrow_ ;
1118 CoinBigIndex *mcstrt = prob->mcstrt_ ;
1119 int *hincol = prob->hincol_ ;
1120 int ncols = prob->ncols_;
1121 int nrows = prob->nrows_;
1122 double *rlo = prob->rlo_;
1123 double *rup = prob->rup_;
1124
1125/*
1126 Scan the rows. We're not
1127 interested in rows that are empty or prohibited.
1128
1129*/
1130 int *which = prob->usefulRowInt_;
1131 int * number = which + nrows;
1132 double * els = prob->usefulRowDouble_;
1133 char * markCol = reinterpret_cast<char *> (prob->usefulColumnInt_);
1134 memset(markCol,0,ncols);
1135 CoinZeroN(els,nrows);
1136 for (int i = 0 ; i < nrows ; i++) {
1137 int nInRow = hinrow[i];
1138 if (nInRow>1 &&!prob->rowProhibited2(i)&&rlo[i]==rup[i]) {
1139 CoinBigIndex rStart = mrstrt[i];
1140 CoinBigIndex k = rStart;
1141 CoinBigIndex rEnd = rStart+nInRow;
1142 double value1=rowels[k];
1143 k++;
1144 for (;k<rEnd;k++) {
1145 if (rowels[k]!=value1)
1146 break;
1147 }
1148 if (k==rEnd) {
1149 // Gub row
1150 int nLook = 0 ;
1151 for (k=rStart;k<rEnd;k++) {
1152 int iColumn = hcol[k];
1153 markCol[iColumn]=1;
1154 CoinBigIndex kk = mcstrt[iColumn];
1155 CoinBigIndex cEnd = kk+hincol[iColumn];
1156 for (;kk<cEnd;kk++) {
1157 int iRow = hrow[kk];
1158 double value = colels[kk];
1159 if (iRow!=i) {
1160 double value2 = els[iRow];
1161 if (value2) {
1162 if (value==value2)
1163 number[iRow]++;
1164 } else {
1165 // first
1166 els[iRow]=value;
1167 number[iRow]=1;
1168 which[nLook++]=iRow;
1169 }
1170 }
1171 }
1172 }
1173 // Now see if any promising
1174 for (int j=0;j<nLook;j++) {
1175 int iRow = which[j];
1176 if (number[iRow]==nInRow) {
1177 // can delete elements and adjust rhs
1178 affectedRows++;
1179 droppedElements += nInRow;
1180 for (CoinBigIndex kk=rStart; kk<rEnd; kk++)
1181 presolve_delete_from_col(iRow,hcol[kk],mcstrt,hincol,hrow,colels) ;
1182 int nInRow2 = hinrow[iRow];
1183 CoinBigIndex rStart2 = mrstrt[iRow];
1184 CoinBigIndex rEnd2 = rStart2+nInRow2;
1185 for (CoinBigIndex kk=rStart2; kk<rEnd2; kk++) {
1186 int iColumn = hcol[kk];
1187 if (markCol[iColumn]==0) {
1188 hcol[rStart2]=iColumn;
1189 rowels[rStart2++]=rowels[kk];
1190 }
1191 }
1192 hinrow[iRow] = nInRow2-nInRow;
1193 if (!hinrow[iRow])
1194 PRESOLVE_REMOVE_LINK(prob->rlink_,iRow) ;
1195 double value =(rlo[i]/value1)*els[iRow];
1196 // correct rhs
1197 if (rlo[iRow]>-1.0e20)
1198 rlo[iRow] -= value;
1199 if (rup[iRow]<1.0e20)
1200 rup[iRow] -= value;
1201 }
1202 els[iRow]=0.0;
1203 }
1204 for (k=rStart;k<rEnd;k++) {
1205 int iColumn = hcol[k];
1206 markCol[iColumn]=0;
1207 }
1208 }
1209 }
1210 }
1211 if (prob->tuning_) {
1212 double thisTime=CoinCpuTime();
1213 printf("CoinPresolveGubrow(1024) - %d elements dropped (%d rows) in time %g, total %g\n",
1214 droppedElements,affectedRows,thisTime-startTime,thisTime-prob->startTime_);
1215 } else if (droppedElements) {
1216#ifdef CLP_INVESTIGATE
1217 printf("CoinPresolveGubrow(1024) - %d elements dropped (%d rows)\n",
1218 droppedElements,affectedRows);
1219#endif
1220 }
1221 return (next);
1222}
1223
1224void gubrow_action::postsolve(CoinPostsolveMatrix *) const
1225{
1226 printf("STILL NO POSTSOLVE FOR GUBROW!\n");
1227 abort();
1228}
1229