1/* $Id: CoinPresolveSubst.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#include "CoinPresolveMatrix.hpp"
10#include "CoinPresolveEmpty.hpp" // for DROP_COL/DROP_ROW
11#include "CoinPresolvePsdebug.hpp"
12#include "CoinPresolveFixed.hpp"
13#include "CoinPresolveZeros.hpp"
14#include "CoinPresolveSubst.hpp"
15#include "CoinMessage.hpp"
16#include "CoinHelperFunctions.hpp"
17#include "CoinSort.hpp"
18#include "CoinError.hpp"
19#include "CoinFinite.hpp"
20
21#if PRESOLVE_DEBUG || PRESOLVE_CONSISTENCY
22#include "CoinPresolvePsdebug.hpp"
23#endif
24
25namespace { // begin unnamed file-local namespace
26
27inline void prepend_elem(int jcol, double coeff, int irow,
28 CoinBigIndex *mcstrt,
29 double *colels,
30 int *hrow,
31 int *link, CoinBigIndex *free_listp)
32{
33 CoinBigIndex kk = *free_listp;
34 assert(kk >= 0) ;
35 *free_listp = link[*free_listp];
36 link[kk] = mcstrt[jcol];
37 mcstrt[jcol] = kk;
38 colels[kk] = coeff;
39 hrow[kk] = irow;
40}
41
42// add coeff_factor * rowy to rowx
43static bool add_row(CoinBigIndex *mrstrt,
44 double *rlo, double * acts, double *rup,
45 double *rowels,
46 int *hcol,
47 int *hinrow,
48 presolvehlink *rlink, int nrows,
49 double coeff_factor,
50 int irowx, int irowy,
51 int *x_to_y)
52{
53 CoinBigIndex krs = mrstrt[irowy];
54 CoinBigIndex kre = krs + hinrow[irowy];
55 CoinBigIndex krsx = mrstrt[irowx];
56 CoinBigIndex krex = krsx + hinrow[irowx];
57 // const int maxk = mrstrt[nrows]; // (22)
58
59 // if irowx is very long, the searching gets very slow,
60 // so we always sort.
61 // whatever sorts rows should handle almost-sorted data efficiently
62 // (quicksort may not)
63 CoinSort_2(hcol+krsx,hcol+krsx+hinrow[irowx],rowels+krsx);
64 CoinSort_2(hcol+krs,hcol+krs+hinrow[irowy],rowels+krs);
65 //ekk_sort2(hcol+krsx, rowels+krsx, hinrow[irowx]);
66 //ekk_sort2(hcol+krs, rowels+krs, hinrow[irowy]);
67
68 //printf("%s x=%d y=%d cf=%g nx=%d ny=%d\n",
69 // "ADD_ROW:",
70 // irowx, irowy, coeff_factor, hinrow[irowx], hinrow[irowy]);
71
72# if PRESOLVE_DEBUG
73 printf("%s x=%d y=%d cf=%g nx=%d ycols=(",
74 "ADD_ROW:",
75 irowx, irowy, coeff_factor, hinrow[irowx]);
76# endif
77
78 // adjust row bounds of rowx;
79 // analogous to adjusting bounds info of colx in doubleton,
80 // or perhaps adjustment to rlo/rup in elim_doubleton
81 //
82 // I believe that since we choose a column that is implied free,
83 // no other column bounds need to be updated.
84 // This is what would happen in doubleton if y's bounds were implied free;
85 // in that case,
86 // lo1 would never improve clo[icolx] and
87 // up1 would never improve cup[icolx].
88 {
89 double rhsy = rlo[irowy];
90
91 // (1)
92 if (-PRESOLVE_INF < rlo[irowx]) {
93# if PRESOLVE_DEBUG
94 if (rhsy * coeff_factor)
95 printf("ELIM_ROW RLO: %g -> %g\n",
96 rlo[irowx],
97 rlo[irowx] + rhsy * coeff_factor);
98# endif
99 rlo[irowx] += rhsy * coeff_factor;
100 }
101 // (2)
102 if (rup[irowx] < PRESOLVE_INF) {
103# if PRESOLVE_DEBUG
104 if (rhsy * coeff_factor)
105 printf("ELIM_ROW RUP: %g -> %g\n",
106 rup[irowx],
107 rup[irowx] + rhsy * coeff_factor);
108# endif
109 rup[irowx] += rhsy * coeff_factor;
110 }
111 if (acts)
112 { acts[irowx] += rhsy * coeff_factor ; }
113 }
114
115 CoinBigIndex kcolx = krsx;
116 CoinBigIndex krex0 = krex;
117 int x_to_y_i = 0;
118
119 for (CoinBigIndex krowy=krs; krowy<kre; krowy++) {
120 int jcol = hcol[krowy];
121
122 // even though these values are updated, they remain consistent
123 PRESOLVEASSERT(krex == krsx + hinrow[irowx]);
124
125 // see if row appears in colx
126 // do NOT look beyond the original elements of rowx
127 //CoinBigIndex kcolx = presolve_find_col1(jcol, krsx, krex, hcol);
128 while (kcolx < krex0 && hcol[kcolx] < jcol)
129 kcolx++;
130
131# if PRESOLVE_DEBUG
132 printf("%d%s ", jcol, (kcolx < krex0 && hcol[kcolx] == jcol) ? "+" : "");
133# endif
134
135 if (kcolx < krex0 && hcol[kcolx] == jcol) {
136 // before: both x and y are in the jcol
137 // after: only x is in the jcol
138 // so: number of elems in col x unchanged, and num elems in jcol is one less
139
140 // update row rep - just modify coefficent
141 // column y is deleted as a whole at the end of the loop
142# if PRESOLVE_DEBUG
143 printf("CHANGING %g + %g -> %g\n",
144 rowels[kcolx],
145 rowels[krowy],
146 rowels[kcolx] + rowels[krowy] * coeff_factor);
147# endif
148 rowels[kcolx] += rowels[krowy] * coeff_factor;
149
150 // this is where this element in rowy ended up
151 x_to_y[x_to_y_i++] = kcolx - krsx;
152 kcolx++;
153 } else {
154 // before: only y is in the jcol
155 // after: only x is in the jcol
156 // so: number of elems in col x is one greater, but num elems in jcol remains same
157 {
158 bool outOfSpace=presolve_expand_row(mrstrt,rowels,hcol,hinrow,rlink,nrows,irowx) ;
159 if (outOfSpace)
160 return true;
161 // this may force a compaction
162 // this will be called excessively if the rows are packed too tightly
163
164 // have to adjust various induction variables
165 krowy = mrstrt[irowy] + (krowy - krs);
166 krs = mrstrt[irowy]; // do this for ease of debugging
167 kre = mrstrt[irowy] + hinrow[irowy];
168
169 kcolx = mrstrt[irowx] + (kcolx - krsx); // don't really need to do this
170 krex0 = mrstrt[irowx] + (krex0 - krsx);
171 krsx = mrstrt[irowx];
172 krex = mrstrt[irowx] + hinrow[irowx];
173 }
174 // this is where this element in rowy ended up
175 x_to_y[x_to_y_i++] = krex - krsx;
176
177 // there is now an unused entry in the memory after the column - use it
178 // mrstrt[nrows] == penultimate index of arrays hcol/rowels
179 hcol[krex] = jcol;
180 rowels[krex] = rowels[krowy] * coeff_factor;
181 hinrow[irowx]++, krex++; // expand the col
182
183 // do NOT increment kcolx
184 }
185 }
186
187# if PRESOLVE_DEBUG
188 printf(")\n");
189# endif
190 return false;
191}
192
193
194// It is common in osl to copy from one representation to another
195// (say from a col rep to a row rep).
196// One such routine is ekkclcp.
197// This is similar, except that it does not assume that the
198// representation is packed, and it adds some slack space
199// in the target rep.
200// It assumes both hincol/hinrow are correct.
201// Note that such routines automatically sort the target rep by index,
202// because they sweep the rows in ascending order.
203void copyrep(const int * mrstrt, const int *hcol, const double *rowels,
204 const int *hinrow, int nrows,
205 int *mcstrt, int *hrow, double *colels,
206 int *hincol, int ncols)
207{
208 int pos = 0;
209 for (int j = 0; j < ncols; ++j) {
210 mcstrt[j] = pos;
211 pos += hincol[j];
212 pos += CoinMin(hincol[j], 10); // slack
213 hincol[j] = 0;
214 }
215
216 for (int i = 0; i < nrows; ++i) {
217 CoinBigIndex krs = mrstrt[i];
218 CoinBigIndex kre = krs + hinrow[i];
219 for (CoinBigIndex kr = krs; kr < kre; ++kr) {
220 int icol = hcol[kr];
221 int iput = hincol[icol];
222 hincol[icol] = iput + 1;
223 iput += mcstrt[icol];
224
225 hrow[iput] = i;
226 colels[iput] = rowels[kr];
227 }
228 }
229}
230
231} // end unnamed file-local namespace
232
233
234const char *subst_constraint_action::name() const
235{
236 return ("subst_constraint_action");
237}
238
239// add -x/y times row y to row x, thus cancelling out one column of rowx;
240// afterwards, that col will be singleton for rowy, so we drop the row.
241//
242// This no longer maintains the col rep as it goes along.
243// Instead, it reconstructs it from scratch afterward.
244//
245// This implements the functionality of ekkrdc3.
246
247/*
248 This routine is called only from implied_free_action. There are several
249 oddities and redundancies in the relationship. The two routines need a good
250 grooming.
251
252 try_fill_level limits the allowable number of coefficients in a column
253 under consideration for substitution. There's some sort of hack going on
254 that has the following effect: if try_fill_level comes in as 2, and that
255 seems overly limiting (number of substitutions < 30), try increasing it to
256 3. To trigger a wider examination of columns, this is actually passed back
257 as -3. The next entry of implied_free_action (and then this routine) will
258 override ColsToDo and examine all columns.
259
260 Hence the initial loop triggered when try_fill_level < 0. Other positive
261 values of fill_level will have no effect. A value of -3 will be converted
262 (and passed back out) as +3. Arbitrary negative values of try_fill_level
263 will also trigger the expansion of search and be converted to positive
264 values.
265
266 I would have thought that the columns considered by implied_free_action
267 should also be limited by fill_level, but that's not currently the case.
268 It's hard-wired to consider columns with 1 to 3 coefficients.
269
270 There must be a better way. -- lh, 040818 --
271*/
272
273const CoinPresolveAction *
274subst_constraint_action::presolve(CoinPresolveMatrix *prob,
275 const int *implied_free,
276 const int * whichFree,
277 int numberFree,
278 const CoinPresolveAction *next,
279 int &try_fill_level)
280{
281 double *colels = prob->colels_;
282 int *hrow = prob->hrow_;
283 CoinBigIndex *mcstrt = prob->mcstrt_;
284 int *hincol = prob->hincol_;
285 const int ncols = prob->ncols_;
286
287 double *rowels = prob->rowels_;
288 int *hcol = prob->hcol_;
289 CoinBigIndex *mrstrt = prob->mrstrt_;
290 int *hinrow = prob->hinrow_;
291 const int nrows = prob->nrows_;
292
293 double *rlo = prob->rlo_;
294 double *rup = prob->rup_;
295 double *acts = prob->acts_;
296
297 double *dcost = prob->cost_;
298
299 presolvehlink *clink = prob->clink_;
300 presolvehlink *rlink = prob->rlink_;
301
302 const double tol = prob->feasibilityTolerance_;
303
304 action *actions = new action [ncols];
305# ifdef ZEROFAULT
306 CoinZeroN(reinterpret_cast<char *>(actions),ncols*sizeof(action)) ;
307# endif
308 int nactions = 0;
309
310 int *zerocols = new int[ncols];
311 int nzerocols = 0;
312
313 int *x_to_y = new int[ncols];
314
315#if 0
316 // follmer.mps presents a challenge, since it has some very
317 // long rows. I started experimenting with how to deal with it,
318 // but haven't yet finished.
319 // The idea was to space out the rows to add some padding between them.
320 // Ideally, we shouldn't have to do this just here, but could try to
321 // do it a little everywhere.
322
323 // sort the row rep by reconstructing from col rep
324 copyrep(mcstrt, hrow, colels, hincol, ncols,
325 mrstrt, hcol, rowels, hinrow, nrows);
326 presolve_make_memlists(/*mrstrt,*/ hinrow, rlink, nrows);
327 // NEED SOME ASSERTION ABOUT NELEMS
328
329 copyrep(mrstrt, hcol, rowels, hinrow, nrows,
330 mcstrt, hrow, colels, hincol, ncols);
331 presolve_make_memlists(/*mcstrt,*/ hincol, clink, ncols);
332#endif
333
334 // in the original presolve, I don't think the two representations were
335 // kept in sync. It may be useful not to do that here, either,
336 // but rather just keep the columns with nfill_level rows accurate
337 // and resync at the end of the function.
338
339 // DEBUGGING
340#if PRESOLVE_DEBUGx
341 int maxsubst = atoi(getenv("MAXSUBST"));
342#else
343 const int maxsubst = 1000000;
344#endif
345
346 int nsubst = 0;
347
348 // This loop does very nearly the same thing as
349 // the first loop in implied_free_action::presolve.
350 // We have to do it again in case constraints change while we process
351 // them (???).
352/*
353 No --- given the hack with -3 coming in to implied_free_action and overriding
354 ColsToDo, we could have columns in implied_free that aren't in ColsToDo.
355 -- lh, 040818 --
356*/
357 int numberLook = prob->numberColsToDo_;
358 int iLook;
359 int * look = prob->colsToDo_;
360 int fill_level = try_fill_level;
361 int * look2 = NULL;
362 // if gone from 2 to 3 look at all
363 if (fill_level<0) {
364 //abort();
365 fill_level=-fill_level;
366 try_fill_level=fill_level;
367 look2 = new int[ncols];
368 look=look2;
369 if (!prob->anyProhibited()) {
370 for (iLook=0;iLook<ncols;iLook++)
371 look[iLook]=iLook;
372 numberLook=ncols;
373 } else {
374 // some prohibited
375 numberLook=0;
376 for (iLook=0;iLook<ncols;iLook++)
377 if (!prob->colProhibited(iLook))
378 look[numberLook++]=iLook;
379 }
380 }
381
382
383 int * rowsUsed = prob->usefulRowInt_+prob->nrows_;
384 int nRowsUsed=0;
385 for (iLook=0;iLook<numberFree;iLook++) {
386 int jcoly=whichFree[iLook];
387 int whichRow = implied_free[iLook];
388 if (hincol[jcoly] <2 || hincol[jcoly] > fill_level)
389 continue;
390 CoinBigIndex kcs = mcstrt[jcoly];
391 CoinBigIndex kce = kcs + hincol[jcoly];
392
393 int bestrowy_size = 0;
394 int bestrowy_row=-1;
395 int bestrowy_k=-1;
396 double bestrowy_coeff=0.0;
397 CoinBigIndex k;
398 for (k=kcs; k<kce; ++k) {
399 int row = hrow[k];
400 double coeffj = colels[k];
401
402 // we don't clean up zeros in the middle of the routine.
403 // if there is one, skip this candidate.
404 if (fabs(coeffj) <= ZTOLDP2 || prob->rowUsed(row)) {
405 bestrowy_size = 0;
406 break;
407 }
408
409 if (row==whichRow) {
410 // if its row is an equality constraint...
411 if (hinrow[row] > 1 && // don't bother with singleton rows
412
413 fabs(rlo[row] - rup[row]) < tol &&
414 !prob->rowUsed(row)) {
415 // both column bounds implied by the constraint bounds
416
417 // we want coeffy to be smaller than x, BACKWARDS from in doubleton
418 bestrowy_size = hinrow[row];
419 bestrowy_row = row;
420 bestrowy_coeff = coeffj;
421 bestrowy_k = k;
422 }
423 }
424 }
425
426 if (bestrowy_size == 0)
427 continue;
428
429 bool all_ok = true;
430 for (k=kcs; k<kce; ++k) {
431 double coeff_factor = fabs(colels[k] / bestrowy_coeff);
432 if (fabs(coeff_factor) > 10.0)
433 all_ok = false;
434 }
435#if 0 // block A
436 // check fill-in
437 if (all_ok && hincol[jcoly] == 3) {
438 // compute fill-in
439 int row1 = -1;
440 int row2=-1;
441 CoinBigIndex kk;
442 for (kk=kcs; kk<kce; ++kk)
443 if (kk != bestrowy_k) {
444 if (row1 == -1)
445 row1 = hrow[kk];
446 else
447 row2 = hrow[kk];
448 }
449
450
451 CoinBigIndex krs = mrstrt[bestrowy_row];
452 CoinBigIndex kre = krs + hinrow[bestrowy_row];
453 CoinBigIndex krs1 = mrstrt[row1];
454 CoinBigIndex kre1 = krs + hinrow[row1];
455 CoinBigIndex krs2 = mrstrt[row2];
456 CoinBigIndex kre2 = krs + hinrow[row2];
457
458 CoinSort_2(hcol+krs,hcol+krs+hinrow[bestrowy_row],rowels+krs);
459 CoinSort_2(hcol+krs1,hcol+krs1+hinrow[row1],rowels+krs1);
460 CoinSort_2(hcol+krs2,hcol+krs2+hinrow[row2],rowels+krs2);
461 //ekk_sort2(hcol+krs, rowels+krs, hinrow[bestrowy_row]);
462 //ekk_sort2(hcol+krs1, rowels+krs1, hinrow[row1]);
463 //ekk_sort2(hcol+krs2, rowels+krs2, hinrow[row2]);
464
465 int nfill = -hinrow[bestrowy_row];
466 CoinBigIndex kcol1 = krs1;
467 for (kk=krs; kk<kre; ++kk) {
468 int jcol = hcol[kk];
469
470 while (kcol1 < kre1 && hcol[kcol1] < jcol)
471 kcol1++;
472 if (! (kcol1 < kre1 && hcol[kcol1] == jcol))
473 nfill++;
474 }
475 CoinBigIndex kcol2 = krs2;
476 for (kk=krs; kk<kre; ++kk) {
477 int jcol = hcol[kk];
478
479 while (kcol2 < kre2 && hcol[kcol2] < jcol)
480 kcol2++;
481 if (! (kcol2 < kre2 && hcol[kcol2] == jcol))
482 nfill++;
483 }
484#if PRESOLVE_DEBUG
485 printf("FILL: %d\n", nfill);
486#endif
487
488#if 0
489 static int maxfill = atoi(getenv("MAXFILL"));
490
491 if (nfill > maxfill)
492 all_ok = false;
493#endif
494
495 // not too much
496 if (nfill <= 0)
497 ngood++;
498
499#if 0
500 static int nts = 0;
501 if (++nts > atoi(getenv("NTS")))
502 all_ok = false;
503 else
504 nt++;
505#endif
506 }
507#endif // end block A
508 // probably never happens
509 if (all_ok && nzerocols + hinrow[bestrowy_row] >= ncols)
510 all_ok = false;
511
512 if (nsubst >= maxsubst) {
513 all_ok = false;
514 }
515
516 if (all_ok) {
517 nsubst++;
518#if 0
519 // debug
520 if (numberLook<ncols&&iLook==numberLook-1) {
521 printf("found last one?? %d\n", jcoly);
522 }
523#endif
524
525 CoinBigIndex kcs = mcstrt[jcoly];
526 int rowy = bestrowy_row;
527 double coeffy = bestrowy_coeff;
528
529 PRESOLVEASSERT(fabs(colels[kcs]) > ZTOLDP);
530 PRESOLVEASSERT(fabs(colels[kcs+1]) > ZTOLDP);
531
532 PRESOLVEASSERT(hinrow[rowy] > 1);
533
534 const bool nonzero_cost = (fabs(dcost[jcoly]) > tol);
535
536 double *costsx = nonzero_cost ? new double[hinrow[rowy]] : 0;
537
538 int ntotels = 0;
539 for (k=kcs; k<kce; ++k) {
540 int irow = hrow[k];
541 ntotels += hinrow[irow];
542 // mark row as contaminated
543 assert (!prob->rowUsed(irow));
544 prob->setRowUsed(irow);
545 rowsUsed[nRowsUsed++]=irow;
546 }
547
548 {
549 action *ap = &actions[nactions++];
550 int nincol = hincol[jcoly];
551
552 ap->col = jcoly;
553 ap->rowy = rowy;
554 PRESOLVE_DETAIL_PRINT(printf("pre_subst %dC %dR E\n",jcoly,rowy));
555
556 ap->nincol = nincol;
557 ap->rows = new int[nincol];
558 ap->rlos = new double[nincol];
559 ap->rups = new double[nincol];
560
561 // coefficients in deleted col
562 ap->coeffxs = new double[nincol];
563
564 ap->ninrowxs = new int[nincol];
565 ap->rowcolsxs = new int[ntotels];
566 ap->rowelsxs = new double[ntotels];
567
568 ap->costsx = costsx;
569
570 // copy all the rows for restoring later - wasteful
571 {
572 int nel = 0;
573 for (CoinBigIndex k=kcs; k<kce; ++k) {
574 int irow = hrow[k];
575 CoinBigIndex krs = mrstrt[irow];
576 //#define COIN_SAFE_SUBST
577#ifdef COIN_SAFE_SUBST
578 CoinBigIndex kre = krs + hinrow[irow];
579 for (CoinBigIndex k1=krs; k1<kre; ++k1) {
580 int jcol = hcol[k1];
581 if (jcol != jcoly) {
582 CoinBigIndex kcs = mcstrt[jcol];
583 CoinBigIndex kce = kcs + hincol[jcol];
584 for (CoinBigIndex k2=kcs; k2<kce; ++k2) {
585 int irow = hrow[k2];
586 if (!prob->rowUsed(irow)) {
587 // mark row as contaminated
588 prob->setRowUsed(irow);
589 rowsUsed[nRowsUsed++]=irow;
590 }
591 }
592 }
593 }
594#endif
595
596 prob->addRow(irow);
597 ap->rows[k-kcs] = irow;
598 ap->ninrowxs[k-kcs] = hinrow[irow];
599 ap->rlos[k-kcs] = rlo[irow];
600 ap->rups[k-kcs] = rup[irow];
601
602 ap->coeffxs[k-kcs] = colels[k];
603
604 CoinMemcpyN( &hcol[krs],hinrow[irow], &ap->rowcolsxs[nel]);
605 CoinMemcpyN( &rowels[krs],hinrow[irow], &ap->rowelsxs[nel]);
606 nel += hinrow[irow];
607 }
608 }
609 }
610
611 // rowy is supposed to be an equality row
612 PRESOLVEASSERT(fabs(rup[rowy] - rlo[rowy]) < ZTOLDP);
613
614 // now adjust for the implied free row - COPIED
615 if (nonzero_cost) {
616#if 0&&PRESOLVE_DEBUG
617 printf("NONZERO SUBST COST: %d %g\n", jcoly, dcost[jcoly]);
618#endif
619 double *cost = dcost;
620 double *save_costs = costsx;
621 double coeffj = coeffy;
622 CoinBigIndex krs = mrstrt[rowy];
623 CoinBigIndex kre = krs + hinrow[rowy];
624
625 double rhs = rlo[rowy];
626 double costj = cost[jcoly];
627
628 for (CoinBigIndex k=krs; k<kre; k++) {
629 int jcol = hcol[k];
630 prob->addCol(jcol);
631 save_costs[k-krs] = cost[jcol];
632
633 if (jcol != jcoly) {
634 double coeff = rowels[k];
635
636 /*
637 * Similar to eliminating doubleton:
638 * cost1 x = cost1 (c - b y) / a = (c cost1)/a - (b cost1)/a
639 * cost[icoly] += cost[icolx] * (-coeff2 / coeff1);
640 */
641 cost[jcol] += costj * (-coeff / coeffj);
642 }
643 }
644
645 // I'm not sure about this
646 prob->change_bias(costj * rhs / coeffj);
647
648 // ??
649 cost[jcoly] = 0.0;
650 }
651
652#if 0&&PRESOLVE_DEBUG
653 if (hincol[jcoly] == 3) {
654 CoinBigIndex krs = mrstrt[rowy];
655 CoinBigIndex kre = krs + hinrow[rowy];
656 printf("HROW0 (%d): ", rowy);
657 for (CoinBigIndex k=krs; k<kre; ++k) {
658 int jcol = hcol[k];
659 double coeff = rowels[k];
660 printf("%d:%g (%d) ", jcol, coeff, hincol[jcol]);
661 }
662 printf("\n");
663 }
664#endif
665
666 if (hincol[jcoly] != 2) {
667 CoinBigIndex krs = mrstrt[rowy];
668 // CoinBigIndex kre = krs + hinrow[rowy];
669 CoinSort_2(hcol+krs,hcol+krs+hinrow[rowy],rowels+krs);
670 //ekk_sort2(hcol+krs, rowels+krs, hinrow[rowy]);
671 }
672
673 // substitute away jcoly in the other rows
674 // Use ap as mcstrt etc may move if compacted
675 kce = hincol[jcoly];
676 action *ap = &actions[nactions-1];
677 for (k=0; k<kce; ++k) {
678 int rowx = ap->rows[k];
679 //assert(rowx==hrow[k+kcs]);
680 //assert(ap->coeffxs[k]==colels[k+kcs]);
681 if (rowx != rowy) {
682 double coeffx = ap->coeffxs[k];
683 double coeff_factor = -coeffx / coeffy; // backwards from doubleton
684
685#if 0&&PRESOLVE_DEBUG
686 {
687 CoinBigIndex krs = mrstrt[rowx];
688 CoinBigIndex kre = krs + hinrow[rowx];
689 printf("HROW (%d %d %d): ", rowx, hinrow[rowx], jcoly);
690 for (CoinBigIndex k=krs; k<kre; ++k) {
691 int jcol = hcol[k];
692 double coeff = rowels[k];
693 printf("%d ", jcol);
694 }
695 printf("\n");
696#if 0
697 for (CoinBigIndex k=krs; k<kre; ++k) {
698 int jcol = hcol[k];
699 prob->addCol(jcol);
700 double coeff = rowels[k];
701 printf("%g ", coeff);
702 }
703 printf("\n");
704#endif
705 }
706#endif
707 {
708 CoinBigIndex krsx = mrstrt[rowx];
709 CoinBigIndex krex = krsx + hinrow[rowx];
710 int i;
711 for (i=krsx;i<krex;i++)
712 prob->addCol(hcol[i]);
713 if (hincol[jcoly] != 2)
714 CoinSort_2(hcol+krsx,hcol+krsx+hinrow[rowx],rowels+krsx);
715 //ekk_sort2(hcol+krsx, rowels+krsx, hinrow[rowx]);
716 }
717
718 // add (coeff_factor * <rowy>) to rowx
719 // does not affect rowy
720 // may introduce (or cancel) elements in rowx
721 bool outOfSpace = add_row(mrstrt,
722 rlo, acts, rup,
723 rowels, hcol,
724 hinrow,
725 rlink, nrows,
726 coeff_factor,
727 rowx, rowy,
728 x_to_y);
729 if (outOfSpace)
730 throwCoinError("out of memory",
731 "CoinImpliedFree::presolve");
732
733 // update col rep of rowx from row rep:
734 // for every col in rowy, copy the elem for that col in rowx
735 // from the row rep to the col rep
736 {
737 CoinBigIndex krs = mrstrt[rowy];
738 // CoinBigIndex kre = krs + hinrow[rowy];
739 int niny = hinrow[rowy];
740
741 CoinBigIndex krsx = mrstrt[rowx];
742 // CoinBigIndex krex = krsx + hinrow[rowx];
743 for (CoinBigIndex ki=0; ki<niny; ++ki) {
744 CoinBigIndex k = krs + ki;
745 int jcol = hcol[k];
746 prob->addCol(jcol);
747 CoinBigIndex kcs = mcstrt[jcol];
748 CoinBigIndex kce = kcs + hincol[jcol];
749
750 //double coeff = rowels[presolve_find_col(jcol, krsx, krex, hcol)];
751 if (hcol[krsx + x_to_y[ki]] != jcol)
752 abort();
753 double coeff = rowels[krsx + x_to_y[ki]];
754
755 // see if rowx appeared in jcol in the col rep
756 CoinBigIndex k2 = presolve_find_row1(rowx, kcs, kce, hrow);
757
758 //PRESOLVEASSERT(fabs(coeff) > ZTOLDP);
759
760 if (k2 < kce) {
761 // yes - just update the entry
762 colels[k2] = coeff;
763 } else {
764 // no - make room, then append
765 bool outOfSpace=presolve_expand_row(mcstrt,colels,hrow,hincol,
766 clink,ncols,jcol) ;
767 if (outOfSpace)
768 throwCoinError("out of memory",
769 "CoinImpliedFree::presolve");
770 krsx = mrstrt[rowx];
771 krs = mrstrt[rowy];
772 kcs = mcstrt[jcol];
773 kce = kcs + hincol[jcol];
774
775 hrow[kce] = rowx;
776 colels[kce] = coeff;
777 hincol[jcol]++;
778 }
779 }
780 }
781 // now colels[k] == 0.0
782
783#if 1
784 // now remove jcoly from rowx in the row rep
785 // better if this were first
786 presolve_delete_from_row(rowx, jcoly, mrstrt, hinrow, hcol, rowels);
787#endif
788#if 0&&PRESOLVE_DEBUG
789 {
790 CoinBigIndex krs = mrstrt[rowx];
791 CoinBigIndex kre = krs + hinrow[rowx];
792 printf("HROW (%d %d %d): ", rowx, hinrow[rowx], jcoly);
793 for (CoinBigIndex k=krs; k<kre; ++k) {
794 int jcol = hcol[k];
795 double coeff = rowels[k];
796 printf("%d ", jcol);
797 }
798 printf("\n");
799#if 0
800 for (CoinBigIndex k=krs; k<kre; ++k) {
801 int jcol = hcol[k];
802 double coeff = rowels[k];
803 printf("%g ", coeff);
804 }
805 printf("\n");
806#endif
807 }
808#endif
809
810 // don't have to update col rep, since entire col deleted later
811 }
812 }
813
814#if 0&&PRESOLVE_DEBUG
815 printf("\n");
816#endif
817
818 // the addition of rows may have created zero coefficients
819 CoinMemcpyN( &hcol[mrstrt[rowy]],hinrow[rowy], &zerocols[nzerocols]);
820 nzerocols += hinrow[rowy];
821
822 // delete rowy in col rep
823 {
824 CoinBigIndex krs = mrstrt[rowy];
825 CoinBigIndex kre = krs + hinrow[rowy];
826 for (CoinBigIndex k=krs; k<kre; ++k) {
827 int jcol = hcol[k];
828
829 // delete rowy from the jcol
830 presolve_delete_from_col(rowy,jcol,mcstrt,hincol,hrow,colels) ;
831 if (hincol[jcol] == 0)
832 { PRESOLVE_REMOVE_LINK(clink,jcol) ; }
833 }
834 }
835 // delete rowy in row rep
836 hinrow[rowy] = 0;
837
838 // This last is entirely dual to doubleton, but for the cost adjustment
839
840 // eliminate col entirely from the col rep
841 PRESOLVE_REMOVE_LINK(clink, jcoly);
842 hincol[jcoly] = 0;
843
844 // eliminate rowy entirely from the row rep
845 PRESOLVE_REMOVE_LINK(rlink, rowy);
846 //cost[irowy] = 0.0;
847
848 rlo[rowy] = 0.0;
849 rup[rowy] = 0.0;
850
851#if 0 && PRESOLVE_DEBUG
852 printf("ROWY COLS: ");
853 for (CoinBigIndex k=0; k<save_ninrowy; ++k)
854 if (rowycols[k] != col) {
855 printf("%d ", rowycols[k]);
856 (void)presolve_find_col(rowycols[k], mrstrt[rowx], mrstrt[rowx]+hinrow[rowx],
857 hcol);
858 }
859 printf("\n");
860#endif
861# if PRESOLVE_CONSISTENCY
862 presolve_links_ok(prob) ;
863 presolve_consistent(prob) ;
864# endif
865 }
866
867 }
868
869 // Clear row used flags
870 for (int i=0;i<nRowsUsed;i++)
871 prob->unsetRowUsed(rowsUsed[i]);
872 // general idea - only do doubletons until there are almost none left
873 if (nactions < 30&&fill_level<prob->maxSubstLevel_)
874 try_fill_level = -fill_level-1;
875 if (nactions) {
876# if PRESOLVE_SUMMARY
877 printf("NSUBSTS: %d\n", nactions);
878 //printf("NT: %d NGOOD: %d FILL_LEVEL: %d\n", nt, ngood, fill_level);
879# endif
880 next = new subst_constraint_action(nactions, CoinCopyOfArray(actions,nactions), next);
881
882 next = drop_zero_coefficients_action::presolve(prob, zerocols, nzerocols, next);
883 }
884 delete [] look2;
885 deleteAction(actions,action*);
886
887 delete[]x_to_y;
888 delete[]zerocols;
889
890 return (next);
891}
892
893void subst_constraint_action::postsolve(CoinPostsolveMatrix *prob) const
894{
895 const action *const actions = actions_;
896 const int nactions = nactions_;
897
898 double *colels = prob->colels_;
899 int *hrow = prob->hrow_;
900 CoinBigIndex *mcstrt = prob->mcstrt_;
901 int *hincol = prob->hincol_;
902 int *link = prob->link_;
903 // int ncols = prob->ncols_;
904
905 double *rlo = prob->rlo_;
906 double *rup = prob->rup_;
907
908 double *dcost = prob->cost_;
909
910 double *sol = prob->sol_;
911 double *rcosts = prob->rcosts_;
912
913 double *acts = prob->acts_;
914 double *rowduals = prob->rowduals_;
915
916# if PRESOLVE_DEBUG || PRESOLVE_CONSISTENCY
917 char *cdone = prob->cdone_;
918 char *rdone = prob->rdone_;
919# endif
920
921 CoinBigIndex &free_list = prob->free_list_;
922
923 // const double ztoldj = prob->ztoldj_;
924 const double maxmin = prob->maxmin_;
925 int k;
926
927 for (const action *f = &actions[nactions-1]; actions<=f; f--) {
928 int icol = f->col;
929
930 int nincoly = f->nincol;
931 double *rlos = f->rlos;
932 double *rups = f->rups;
933 int *rows = f->rows;
934
935 double *coeffxs = f->coeffxs;
936
937 int jrowy = f->rowy;
938
939 int *ninrowxs = f->ninrowxs;
940 const int *rowcolsxs = f->rowcolsxs;
941 const double *rowelsxs = f->rowelsxs;
942
943 /* the row was in the reduced problem */
944 for (int i=0; i<nincoly; ++i) {
945 if (rows[i] != jrowy)
946 PRESOLVEASSERT(rdone[rows[i]]);
947 }
948 PRESOLVEASSERT(cdone[icol]==DROP_COL);
949 PRESOLVEASSERT(rdone[jrowy]==DROP_ROW);
950
951 // DEBUG CHECK
952#if 1 && PRESOLVE_DEBUG
953 {
954 double actx = 0.0;
955 const double ztolzb = prob->ztolzb_;
956 for (int j=0; j<prob->ncols_; ++j)
957 if (hincol[j] > 0 && cdone[j]) {
958 CoinBigIndex krow = presolve_find_row1(jrowy, mcstrt[j], mcstrt[j] + hincol[j], hrow);
959 if (krow < mcstrt[j] + hincol[j])
960 actx += colels[krow] * sol[j];
961 }
962 if (! (fabs(acts[jrowy] - actx) < 100*ztolzb))
963 printf("BAD ACTSX: acts[%d]==%g != %g\n",
964 jrowy, acts[jrowy], actx);
965 if (! (rlo[jrowy] - 100*ztolzb <= actx && actx <= rup[jrowy] + 100*ztolzb))
966 printf("ACTSX NOT IN RANGE: %d %g %g %g\n",
967 jrowy, rlo[jrowy], actx, rup[jrowy]);
968 }
969#endif
970
971 int ninrowy=-1;
972 const int *rowcolsy=NULL;
973 const double *rowelsy=NULL;
974 double coeffy=0.0;
975
976 double rloy=1.0e50;
977 {
978 int nel = 0;
979 for (int i=0; i<nincoly; ++i) {
980 int row = rows[i];
981 rlo[row] = rlos[i];
982 rup[row] = rups[i];
983 if (row == jrowy) {
984 ninrowy = ninrowxs[i];
985 rowcolsy = &rowcolsxs[nel];
986 rowelsy = &rowelsxs[nel];
987
988 coeffy = coeffxs[i];
989 rloy = rlo[row];
990
991 }
992 nel += ninrowxs[i];
993 }
994 }
995 double rhsy = rloy;
996
997 // restore costs
998 {
999 const double *costs = f->costsx;
1000 if (costs)
1001 for (int i = 0; i<ninrowy; ++i) {
1002 dcost[rowcolsy[i]] = costs[i];
1003 }
1004 }
1005
1006 // solve for the equality to find the solution for the eliminated col
1007 // this is why we want coeffx < coeffy (55)
1008 {
1009 double sol0 = rloy;
1010 sol[icol] = 0.0; // to avoid condition in loop
1011 for (k = 0; k<ninrowy; ++k) {
1012 int jcolx = rowcolsy[k];
1013 double coeffx = rowelsy[k];
1014 sol0 -= coeffx * sol[jcolx];
1015 }
1016 sol[icol] = sol0 / coeffy;
1017
1018# if PRESOLVE_DEBUG
1019 const double ztolzb = prob->ztolzb_;
1020 double *clo = prob->clo_;
1021 double *cup = prob->cup_;
1022
1023 if (! (sol[icol] > clo[icol] - ztolzb &&
1024 cup[icol] + ztolzb > sol[icol]))
1025 printf("NEW SOL OUT-OF-TOL: %g %g %g\n", clo[icol],
1026 sol[icol], cup[icol]);
1027# endif
1028 }
1029
1030 // since this row is fixed
1031 acts[jrowy] = rloy;
1032
1033 // acts[irow] always ok, since slack is fixed
1034 prob->setRowStatus(jrowy,CoinPrePostsolveMatrix::atLowerBound);
1035
1036 // remove old rowx from col rep
1037 // we don't explicitly store what the current rowx is;
1038 // however, after the presolve, rowx contains a col for every
1039 // col in either the original rowx or the original rowy.
1040 // If there were cancellations, those were handled in subsequent
1041 // presolves.
1042 {
1043 // erase those cols in the other rows that occur in rowy
1044 // (with the exception of icol, which was deleted);
1045 // the other rows *must* contain these cols
1046 for (k = 0; k<ninrowy; ++k) {
1047 int col = rowcolsy[k];
1048
1049 // remove jrowx from col in the col rep
1050 // icol itself was deleted, so won't be there
1051 if (col != icol)
1052 for (int i = 0; i<nincoly; ++i) {
1053 if (rows[i] != jrowy)
1054 presolve_delete_from_col2(rows[i],col,mcstrt,hincol,hrow,/*colels,*/
1055 link,&free_list) ;
1056 }
1057 }
1058# if PRESOLVE_CONSISTENCY
1059 presolve_check_free_list(prob) ;
1060# endif
1061
1062 // initialize this for loops below
1063 hincol[icol] = 0;
1064
1065 // now restore the original rows (other than rowy).
1066 // those cols that were also in rowy were just removed;
1067 // otherwise, they *must* already be there.
1068 // This loop and the next automatically create the rep for the new col.
1069 {
1070 const int *rowcolsx = rowcolsxs;
1071 const double *rowelsx = rowelsxs;
1072
1073 for (int i = 0; i<nincoly; ++i) {
1074 int ninrowx = ninrowxs[i];
1075 int jrowx = rows[i];
1076
1077 if (jrowx != jrowy)
1078 for (k = 0; k<ninrowx; ++k) {
1079 int col = rowcolsx[k];
1080 CoinBigIndex kcolx = presolve_find_row3(jrowx, mcstrt[col], hincol[col], hrow, link);
1081
1082 if (kcolx != -1) {
1083 PRESOLVEASSERT(presolve_find_col1(col, 0, ninrowy, rowcolsy) == ninrowy);
1084 // overwrite the existing entry
1085 colels[kcolx] = rowelsx[k];
1086 } else {
1087 PRESOLVEASSERT(presolve_find_col1(col, 0, ninrowy, rowcolsy) < ninrowy);
1088
1089 {
1090 CoinBigIndex kk = free_list;
1091 assert(kk >= 0 && kk < prob->bulk0_) ;
1092 free_list = link[free_list];
1093
1094 link[kk] = mcstrt[col];
1095 mcstrt[col] = kk;
1096 colels[kk] = rowelsx[k];
1097 hrow[kk] = jrowx;
1098 }
1099 ++hincol[col];
1100 }
1101 }
1102 rowcolsx += ninrowx;
1103 rowelsx += ninrowx;
1104 }
1105# if PRESOLVE_CONSISTENCY
1106 presolve_check_free_list(prob) ;
1107# endif
1108 }
1109
1110 // finally, add original rowy elements
1111 for (k = 0; k<ninrowy; ++k) {
1112 int col = rowcolsy[k];
1113
1114 {
1115 prepend_elem(col, rowelsy[k], jrowy, mcstrt, colels, hrow, link, &free_list);
1116 ++hincol[col];
1117 }
1118 }
1119# if PRESOLVE_CONSISTENCY
1120 presolve_check_free_list(prob) ;
1121# endif
1122 }
1123
1124 // my guess is that the CLAIM in doubleton generalizes to
1125 // equations with more than one x-style variable.
1126 // Since I can't see how to distinguish among them,
1127 // I assume that any of them will do.
1128
1129 {
1130 // CoinBigIndex k;
1131 double dj = maxmin*dcost[icol];
1132 double bounds_factor = rhsy/coeffy;
1133 for (int i=0; i<nincoly; ++i)
1134 if (rows[i] != jrowy) {
1135 int row = rows[i];
1136 double coeff = coeffxs[i];
1137
1138 // PROBABLY DOESN'T MAKE SENSE
1139 acts[row] += coeff * bounds_factor;
1140
1141 dj -= rowduals[row] * coeff;
1142 }
1143
1144 // DEBUG CHECK
1145 double acty = 0.0;
1146 for (k = 0; k<ninrowy; ++k) {
1147 int col = rowcolsy[k];
1148 acty += rowelsy[k] * sol[col];
1149 }
1150
1151 PRESOLVEASSERT(fabs(acty - acts[jrowy]) < 100*ZTOLDP);
1152
1153 // RECOMPUTING
1154 {
1155 const int *rowcolsx = rowcolsxs;
1156 const double *rowelsx = rowelsxs;
1157
1158 for (int i=0; i<nincoly; ++i) {
1159 int ninrowx = ninrowxs[i];
1160
1161 if (rows[i] != jrowy) {
1162 int jrowx = rows[i];
1163
1164 double actx = 0.0;
1165 for (k = 0; k<ninrowx; ++k) {
1166 int col = rowcolsx[k];
1167 actx += rowelsx[k] * sol[col];
1168 }
1169 PRESOLVEASSERT(rlo[jrowx] - prob->ztolzb_ <= actx
1170 && actx <= rup[jrowx] + prob->ztolzb_);
1171 acts[jrowx] = actx;
1172 if (prob->getRowStatus(jrowx)!=CoinPrePostsolveMatrix::basic) {
1173 if (actx-rlo[jrowx]<rup[jrowx]-actx)
1174 prob->setRowStatus(jrowx,CoinPrePostsolveMatrix::atLowerBound);
1175 else
1176 prob->setRowStatus(jrowx,CoinPrePostsolveMatrix::atUpperBound);
1177 }
1178 }
1179 rowcolsx += ninrowx;
1180 rowelsx += ninrowx;
1181 }
1182 }
1183
1184 // this is the coefficient we need to force col y's reduced cost to 0.0;
1185 // for example, this is obviously true if y is a singleton column
1186 rowduals[jrowy] = dj / coeffy;
1187 rcosts[icol] = 0.0;
1188
1189 // furthermore, this should leave rcosts[colx] for all colx
1190 // in jrowx unchanged (????).
1191 }
1192
1193 // Unlike doubleton, there should never be a problem with keeping
1194 // the reduced costs the way they were, because the other
1195 // variable's bounds are never changed, since col was implied free.
1196 //rowstat[jrowy] = 0;
1197 prob->setColumnStatus(icol,CoinPrePostsolveMatrix::basic);
1198
1199# if PRESOLVE_DEBUG || PRESOLVE_CONSISTENCY
1200 cdone[icol] = SUBST_ROW;
1201 rdone[jrowy] = SUBST_ROW;
1202# endif
1203 }
1204
1205# if PRESOLVE_CONSISTENCY
1206 presolve_check_threads(prob) ;
1207# endif
1208
1209 return ;
1210}
1211
1212
1213
1214subst_constraint_action::~subst_constraint_action()
1215{
1216 const action *actions = actions_;
1217
1218 for (int i=0; i<nactions_; ++i) {
1219 delete[]actions[i].rows;
1220 delete[]actions[i].rlos;
1221 delete[]actions[i].rups;
1222 delete[]actions[i].coeffxs;
1223 delete[]actions[i].ninrowxs;
1224 delete[]actions[i].rowcolsxs;
1225 delete[]actions[i].rowelsxs;
1226
1227
1228 //delete [](double*)actions[i].costsx;
1229 deleteAction(actions[i].costsx,double*);
1230 }
1231
1232 // Must add cast to placate MS compiler
1233 //delete [] (subst_constraint_action::action*)actions_;
1234 deleteAction(actions_,subst_constraint_action::action*);
1235}
1236