1/* $Id: CoinPresolveTripleton.cpp 1448 2011-06-19 15:34:41Z stefan $ */
2// Copyright (C) 2003, 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 "CoinFinite.hpp"
10#include "CoinHelperFunctions.hpp"
11#include "CoinPresolveMatrix.hpp"
12
13#include "CoinPresolveEmpty.hpp" // for DROP_COL/DROP_ROW
14#include "CoinPresolveZeros.hpp"
15#include "CoinPresolveFixed.hpp"
16#include "CoinPresolveTripleton.hpp"
17
18#include "CoinPresolvePsdebug.hpp"
19#include "CoinMessage.hpp"
20
21#if PRESOLVE_DEBUG || PRESOLVE_CONSISTENCY
22#include "CoinPresolvePsdebug.hpp"
23#endif
24
25/*
26 * Substituting y away:
27 *
28 * y = (c - a x - d z) / b
29 *
30 * so adjust bounds by: c/b
31 * and x by: -a/b
32 * and z by: -d/b
33 *
34 * This affects both the row and col representations.
35 *
36 * mcstrt only modified if the column must be moved.
37 *
38 * for every row in icoly
39 * if icolx is also has an entry for row
40 * modify the icolx entry for row
41 * drop the icoly entry from row and modify the icolx entry
42 * else
43 * add a new entry to icolx column
44 * create a new icolx entry
45 * (this may require moving the column in memory)
46 * replace icoly entry from row and replace with icolx entry
47 *
48 * same for icolz
49 * The row and column reps are inconsistent during the routine,
50 * because icolx in the column rep is updated, and the entries corresponding
51 * to icolx in the row rep are updated, but nothing concerning icoly
52 * in the col rep is changed. icoly entries in the row rep are deleted,
53 * and icolx entries in both reps are consistent.
54 * At the end, we set the length of icoly to be zero, so the reps would
55 * be consistent if the row were deleted from the row rep.
56 * Both the row and icoly must be removed from both reps.
57 * In the col rep, icoly will be eliminated entirely, at the end of the routine;
58 * irow occurs in just two columns, one of which (icoly) is eliminated
59 * entirely, the other is icolx, which is not deleted here.
60 * In the row rep, irow will be eliminated entirely, but not here;
61 * icoly is removed from the rows it occurs in.
62 */
63static bool elim_tripleton(const char *
64#ifdef PRESOLVE_DEBUG
65msg
66#endif
67 ,
68 CoinBigIndex *mcstrt,
69 double *rlo, double * acts, double *rup,
70 double *colels,
71 int *hrow, int *hcol,
72 int *hinrow, int *hincol,
73 presolvehlink *clink, int ncols,
74 presolvehlink *rlink, int nrows,
75 CoinBigIndex *mrstrt, double *rowels,
76 //double a, double b, double c,
77 double coeff_factorx,double coeff_factorz,
78 double bounds_factor,
79 int row0, int icolx, int icoly, int icolz)
80{
81 CoinBigIndex kcs = mcstrt[icoly];
82 CoinBigIndex kce = kcs + hincol[icoly];
83 CoinBigIndex kcsx = mcstrt[icolx];
84 CoinBigIndex kcex = kcsx + hincol[icolx];
85 CoinBigIndex kcsz = mcstrt[icolz];
86 CoinBigIndex kcez = kcsz + hincol[icolz];
87
88# if PRESOLVE_DEBUG
89 printf("%s %d x=%d y=%d z=%d cfx=%g cfz=%g nx=%d yrows=(", msg,
90 row0,icolx,icoly,icolz,coeff_factorx,coeff_factorz,hincol[icolx]) ;
91# endif
92 for (CoinBigIndex kcoly=kcs; kcoly<kce; kcoly++) {
93 int row = hrow[kcoly];
94
95 // even though these values are updated, they remain consistent
96 PRESOLVEASSERT(kcex == kcsx + hincol[icolx]);
97 PRESOLVEASSERT(kcez == kcsz + hincol[icolz]);
98
99 // we don't need to update the row being eliminated
100 if (row != row0/* && hinrow[row] > 0*/) {
101 if (bounds_factor != 0.0) {
102 // (1)
103 if (-PRESOLVE_INF < rlo[row])
104 rlo[row] -= colels[kcoly] * bounds_factor;
105
106 // (2)
107 if (rup[row] < PRESOLVE_INF)
108 rup[row] -= colels[kcoly] * bounds_factor;
109
110 // and solution
111 if (acts)
112 { acts[row] -= colels[kcoly] * bounds_factor; }
113 }
114 // see if row appears in colx
115 CoinBigIndex kcolx = presolve_find_row1(row, kcsx, kcex, hrow);
116# if PRESOLVE_DEBUG
117 printf("%d%s ",row,(kcolx<kcex)?"x+":"") ;
118# endif
119 // see if row appears in colz
120 CoinBigIndex kcolz = presolve_find_row1(row, kcsz, kcez, hrow);
121# if PRESOLVE_DEBUG
122 printf("%d%s ",row,(kcolz<kcez)?"x+":"") ;
123# endif
124
125 if (kcolx>=kcex&&kcolz<kcez) {
126 // swap
127 int iTemp;
128 iTemp=kcolx;
129 kcolx=kcolz;
130 kcolz=iTemp;
131 iTemp=kcsx;
132 kcsx=kcsz;
133 kcsz=iTemp;
134 iTemp=kcex;
135 kcex=kcez;
136 kcez=iTemp;
137 iTemp=icolx;
138 icolx=icolz;
139 icolz=iTemp;
140 double dTemp=coeff_factorx;
141 coeff_factorx=coeff_factorz;
142 coeff_factorz=dTemp;
143 }
144 if (kcolx<kcex) {
145 // before: both x and y are in the row
146 // after: only x is in the row
147 // so: number of elems in col x unchanged, and num elems in row is one less
148
149 // update col rep - just modify coefficent
150 // column y is deleted as a whole at the end of the loop
151 colels[kcolx] += colels[kcoly] * coeff_factorx;
152 // update row rep
153 // first, copy new value for col x into proper place in rowels
154 CoinBigIndex k2 = presolve_find_col(icolx, mrstrt[row], mrstrt[row]+hinrow[row], hcol);
155 rowels[k2] = colels[kcolx];
156 if (kcolz<kcez) {
157 // before: both z and y are in the row
158 // after: only z is in the row
159 // so: number of elems in col z unchanged, and num elems in row is one less
160
161 // update col rep - just modify coefficent
162 // column y is deleted as a whole at the end of the loop
163 colels[kcolz] += colels[kcoly] * coeff_factorz;
164 // update row rep
165 // first, copy new value for col z into proper place in rowels
166 CoinBigIndex k2 = presolve_find_col(icolz, mrstrt[row], mrstrt[row]+hinrow[row], hcol);
167 rowels[k2] = colels[kcolz];
168 // now delete col y from the row; this changes hinrow[row]
169 presolve_delete_from_row(row, icoly, mrstrt, hinrow, hcol, rowels);
170 } else {
171 // before: only y is in the row
172 // after: only z is in the row
173 // so: number of elems in col z is one greater, but num elems in row remains same
174 // update entry corresponding to icolz in row rep
175 // by just overwriting the icoly entry
176 {
177 CoinBigIndex k2 = presolve_find_col(icoly, mrstrt[row], mrstrt[row]+hinrow[row], hcol);
178 hcol[k2] = icolz;
179 rowels[k2] = colels[kcoly] * coeff_factorz;
180 }
181
182 {
183 bool no_mem = presolve_expand_col(mcstrt,colels,hrow,hincol,
184 clink,ncols,icolz);
185 if (no_mem)
186 return (true);
187
188 // have to adjust various induction variables
189 kcolx = mcstrt[icolx] + (kcolx - kcsx);
190 kcsx = mcstrt[icolx];
191 kcex = mcstrt[icolx] + hincol[icolx];
192 kcoly = mcstrt[icoly] + (kcoly - kcs);
193 kcs = mcstrt[icoly]; // do this for ease of debugging
194 kce = mcstrt[icoly] + hincol[icoly];
195
196 kcolz = mcstrt[icolz] + (kcolz - kcs); // don't really need to do this
197 kcsz = mcstrt[icolz];
198 kcez = mcstrt[icolz] + hincol[icolz];
199 }
200
201 // there is now an unused entry in the memory after the column - use it
202 // mcstrt[ncols] == penultimate index of arrays hrow/colels
203 hrow[kcez] = row;
204 colels[kcez] = colels[kcoly] * coeff_factorz; // y factor is 0.0
205 hincol[icolz]++, kcez++; // expand the col
206 }
207 } else {
208 // before: only y is in the row
209 // after: only x and z are in the row
210 // update entry corresponding to icolx in row rep
211 // by just overwriting the icoly entry
212 {
213 CoinBigIndex k2 = presolve_find_col(icoly, mrstrt[row], mrstrt[row]+hinrow[row], hcol);
214 hcol[k2] = icolx;
215 rowels[k2] = colels[kcoly] * coeff_factorx;
216 }
217 presolve_expand_row(mrstrt,rowels,hcol,hinrow,rlink,nrows,row) ;
218 // there is now an unused entry in the memory after the column - use it
219 int krez = mrstrt[row]+hinrow[row];
220 hcol[krez] = icolz;
221 rowels[krez] = colels[kcoly] * coeff_factorz;
222 hinrow[row]++;
223
224 {
225 bool no_mem = presolve_expand_col(mcstrt,colels,hrow,hincol,
226 clink,ncols,icolx) ;
227 if (no_mem)
228 return (true);
229
230 // have to adjust various induction variables
231 kcoly = mcstrt[icoly] + (kcoly - kcs);
232 kcs = mcstrt[icoly]; // do this for ease of debugging
233 kce = mcstrt[icoly] + hincol[icoly];
234
235 kcolx = mcstrt[icolx] + (kcolx - kcs); // don't really need to do this
236 kcsx = mcstrt[icolx];
237 kcex = mcstrt[icolx] + hincol[icolx];
238 kcolz = mcstrt[icolz] + (kcolz - kcs); // don't really need to do this
239 kcsz = mcstrt[icolz];
240 kcez = mcstrt[icolz] + hincol[icolz];
241 }
242
243 // there is now an unused entry in the memory after the column - use it
244 hrow[kcex] = row;
245 colels[kcex] = colels[kcoly] * coeff_factorx; // y factor is 0.0
246 hincol[icolx]++, kcex++; // expand the col
247
248 {
249 bool no_mem = presolve_expand_col(mcstrt,colels,hrow,hincol,clink,
250 ncols,icolz);
251 if (no_mem)
252 return (true);
253
254 // have to adjust various induction variables
255 kcoly = mcstrt[icoly] + (kcoly - kcs);
256 kcs = mcstrt[icoly]; // do this for ease of debugging
257 kce = mcstrt[icoly] + hincol[icoly];
258
259 kcsx = mcstrt[icolx];
260 kcex = mcstrt[icolx] + hincol[icolx];
261 kcsz = mcstrt[icolz];
262 kcez = mcstrt[icolz] + hincol[icolz];
263 }
264
265 // there is now an unused entry in the memory after the column - use it
266 hrow[kcez] = row;
267 colels[kcez] = colels[kcoly] * coeff_factorz; // y factor is 0.0
268 hincol[icolz]++, kcez++; // expand the col
269 }
270 }
271 }
272
273# if PRESOLVE_DEBUG
274 printf(")\n") ;
275# endif
276
277 // delete the whole column
278 hincol[icoly] = 0;
279
280 return (false);
281}
282
283
284
285
286/*
287 *
288 * The col rep and row rep must be consistent.
289 */
290const CoinPresolveAction *tripleton_action::presolve(CoinPresolveMatrix *prob,
291 const CoinPresolveAction *next)
292{
293 double startTime = 0.0;
294 int startEmptyRows=0;
295 int startEmptyColumns = 0;
296 if (prob->tuning_) {
297 startTime = CoinCpuTime();
298 startEmptyRows = prob->countEmptyRows();
299 startEmptyColumns = prob->countEmptyCols();
300 }
301 double *colels = prob->colels_;
302 int *hrow = prob->hrow_;
303 CoinBigIndex *mcstrt = prob->mcstrt_;
304 int *hincol = prob->hincol_;
305 int ncols = prob->ncols_;
306
307 double *clo = prob->clo_;
308 double *cup = prob->cup_;
309
310 double *rowels = prob->rowels_;
311 int *hcol = prob->hcol_;
312 CoinBigIndex *mrstrt = prob->mrstrt_;
313 int *hinrow = prob->hinrow_;
314 int nrows = prob->nrows_;
315
316 double *rlo = prob->rlo_;
317 double *rup = prob->rup_;
318
319 presolvehlink *clink = prob->clink_;
320 presolvehlink *rlink = prob->rlink_;
321
322 const unsigned char *integerType = prob->integerType_;
323
324 double *cost = prob->cost_;
325
326 int numberLook = prob->numberRowsToDo_;
327 int iLook;
328 int * look = prob->rowsToDo_;
329 const double ztolzb = prob->ztolzb_;
330
331 action * actions = new action [nrows];
332# ifdef ZEROFAULT
333 // initialise alignment padding bytes
334 memset(actions,0,nrows*sizeof(action)) ;
335# endif
336 int nactions = 0;
337
338 int *zeros = prob->usefulColumnInt_; //new int[ncols];
339 char * mark = reinterpret_cast<char *>(zeros+ncols);
340 memset(mark,0,ncols);
341 int nzeros = 0;
342
343 // If rowstat exists then all do
344 unsigned char *rowstat = prob->rowstat_;
345 double *acts = prob->acts_;
346 // unsigned char * colstat = prob->colstat_;
347
348
349# if PRESOLVE_CONSISTENCY
350 presolve_links_ok(prob) ;
351# endif
352
353 // wasfor (int irow=0; irow<nrows; irow++)
354 for (iLook=0;iLook<numberLook;iLook++) {
355 int irow = look[iLook];
356 if (hinrow[irow] == 3 &&
357 fabs(rup[irow] - rlo[irow]) <= ZTOLDP) {
358 double rhs = rlo[irow];
359 CoinBigIndex krs = mrstrt[irow];
360 CoinBigIndex kre = krs + hinrow[irow];
361 int icolx, icoly, icolz;
362 double coeffx, coeffy, coeffz;
363 CoinBigIndex k;
364
365 /* locate first column */
366 for (k=krs; k<kre; k++) {
367 if (hincol[hcol[k]] > 0) {
368 break;
369 }
370 }
371 PRESOLVEASSERT(k<kre);
372 coeffx = rowels[k];
373 if (fabs(coeffx) < ZTOLDP2)
374 continue;
375 icolx = hcol[k];
376
377
378 /* locate second column */
379 for (k++; k<kre; k++) {
380 if (hincol[hcol[k]] > 0) {
381 break;
382 }
383 }
384 PRESOLVEASSERT(k<kre);
385 coeffy = rowels[k];
386 if (fabs(coeffy) < ZTOLDP2)
387 continue;
388 icoly = hcol[k];
389
390 /* locate third column */
391 for (k++; k<kre; k++) {
392 if (hincol[hcol[k]] > 0) {
393 break;
394 }
395 }
396 PRESOLVEASSERT(k<kre);
397 coeffz = rowels[k];
398 if (fabs(coeffz) < ZTOLDP2)
399 continue;
400 icolz = hcol[k];
401
402 // For now let's do obvious one
403 if (coeffx*coeffz>0.0) {
404 if(coeffx*coeffy>0.0)
405 continue;
406 } else if (coeffx*coeffy>0.0) {
407 int iTemp = icoly;
408 icoly=icolz;
409 icolz=iTemp;
410 double dTemp = coeffy;
411 coeffy=coeffz;
412 coeffz=dTemp;
413 } else {
414 int iTemp = icoly;
415 icoly=icolx;
416 icolx=iTemp;
417 double dTemp = coeffy;
418 coeffy=coeffx;
419 coeffx=dTemp;
420 }
421 // Not all same sign and y is odd one out
422 // don't bother with fixed variables
423 if (!(fabs(cup[icolx] - clo[icolx]) < ZTOLDP) &&
424 !(fabs(cup[icoly] - clo[icolx]) < ZTOLDP) &&
425 !(fabs(cup[icolz] - clo[icoly]) < ZTOLDP)) {
426 assert (coeffx*coeffz>0.0&&coeffx*coeffy<0.0);
427 // Only do if does not give implicit bounds on x and z
428 double cx = - coeffx/coeffy;
429 double cz = - coeffz/coeffy;
430 /* don't do if y integer for now */
431 if (integerType[icoly]) {
432#define PRESOLVE_DANGEROUS
433#ifndef PRESOLVE_DANGEROUS
434 continue;
435#else
436 if (!integerType[icolx]||!integerType[icolz])
437 continue;
438 if (cx!=floor(cx+0.5)||cz!=floor(cz+0.5))
439 continue;
440#endif
441 }
442 double rhsRatio = rhs/coeffy;
443 if (clo[icoly]>-1.0e30) {
444 if (clo[icolx]<-1.0e30||clo[icolz]<-1.0e30)
445 continue;
446 if (cx*clo[icolx]+cz*clo[icolz]+rhsRatio<clo[icoly]-ztolzb)
447 continue;
448 }
449 if (cup[icoly]<1.0e30) {
450 if (cup[icolx]>1.0e30||cup[icolz]>1.0e30)
451 continue;
452 if (cx*cup[icolx]+cz*cup[icolz]+rhsRatio>cup[icoly]+ztolzb)
453 continue;
454 }
455 CoinBigIndex krowx,krowy,krowz;
456 /* find this row in each of the columns and do counts */
457 bool singleton=false;
458 for (k=mcstrt[icoly]; k<mcstrt[icoly]+hincol[icoly]; k++) {
459 int jrow=hrow[k];
460 if (hinrow[jrow]==1)
461 singleton=true;
462 if (jrow == irow)
463 krowy=k;
464 else
465 prob->setRowUsed(jrow);
466 }
467 int nDuplicate=0;
468 for (k=mcstrt[icolx]; k<mcstrt[icolx]+hincol[icolx]; k++) {
469 int jrow=hrow[k];
470 if (jrow == irow)
471 krowx=k;
472 else if (prob->rowUsed(jrow))
473 nDuplicate++;
474 }
475 for (k=mcstrt[icolz]; k<mcstrt[icolz]+hincol[icolz]; k++) {
476 int jrow=hrow[k];
477 if (jrow == irow)
478 krowz=k;
479 else if (prob->rowUsed(jrow))
480 nDuplicate++;
481 }
482 int nAdded=hincol[icoly]-3-nDuplicate;
483 for (k=mcstrt[icoly]; k<mcstrt[icoly]+hincol[icoly]; k++) {
484 int jrow=hrow[k];
485 prob->unsetRowUsed(jrow);
486 }
487 // let singleton rows be taken care of first
488 if (singleton)
489 continue;
490 //if (nAdded<=1)
491 //printf("%d elements added, hincol %d , dups %d\n",nAdded,hincol[icoly],nDuplicate);
492 if (nAdded>2)
493 continue;
494
495 // it is possible that both x/z and y are singleton columns
496 // that can cause problems
497 if ((hincol[icolx] == 1 ||hincol[icolz] == 1) && hincol[icoly] == 1)
498 continue;
499
500 // common equations are of the form ax + by = 0, or x + y >= lo
501 {
502 action *s = &actions[nactions];
503 nactions++;
504 PRESOLVE_DETAIL_PRINT(printf("pre_tripleton %dR %dC %dC %dC E\n",
505 irow,icoly,icolx,icolz));
506
507 s->row = irow;
508 s->icolx = icolx;
509 s->icolz = icolz;
510
511 s->icoly = icoly;
512 s->cloy = clo[icoly];
513 s->cupy = cup[icoly];
514 s->costy = cost[icoly];
515
516 s->rlo = rlo[irow];
517 s->rup = rup[irow];
518
519 s->coeffx = coeffx;
520 s->coeffy = coeffy;
521 s->coeffz = coeffz;
522
523 s->ncoly = hincol[icoly];
524 s->colel = presolve_dupmajor(colels, hrow, hincol[icoly],
525 mcstrt[icoly]);
526 }
527
528 // costs
529 // the effect of maxmin cancels out
530 cost[icolx] += cost[icoly] * cx;
531 cost[icolz] += cost[icoly] * cz;
532
533 prob->change_bias(cost[icoly] * rhs / coeffy);
534 //if (cost[icoly]*rhs)
535 //printf("change %g col %d cost %g rhs %g coeff %g\n",cost[icoly]*rhs/coeffy,
536 // icoly,cost[icoly],rhs,coeffy);
537
538 if (rowstat) {
539 // update solution and basis
540 int numberBasic=0;
541 if (prob->columnIsBasic(icoly))
542 numberBasic++;
543 if (prob->rowIsBasic(irow))
544 numberBasic++;
545 if (numberBasic>1) {
546 if (!prob->columnIsBasic(icolx))
547 prob->setColumnStatus(icolx,CoinPrePostsolveMatrix::basic);
548 else
549 prob->setColumnStatus(icolz,CoinPrePostsolveMatrix::basic);
550 }
551 }
552
553 // Update next set of actions
554 {
555 prob->addCol(icolx);
556 int i,kcs,kce;
557 kcs = mcstrt[icoly];
558 kce = kcs + hincol[icoly];
559 for (i=kcs;i<kce;i++) {
560 int row = hrow[i];
561 prob->addRow(row);
562 }
563 kcs = mcstrt[icolx];
564 kce = kcs + hincol[icolx];
565 for (i=kcs;i<kce;i++) {
566 int row = hrow[i];
567 prob->addRow(row);
568 }
569 prob->addCol(icolz);
570 kcs = mcstrt[icolz];
571 kce = kcs + hincol[icolz];
572 for (i=kcs;i<kce;i++) {
573 int row = hrow[i];
574 prob->addRow(row);
575 }
576 }
577
578 /* transfer the colx factors to coly */
579 bool no_mem = elim_tripleton("ELIMT",
580 mcstrt, rlo, acts, rup, colels,
581 hrow, hcol, hinrow, hincol,
582 clink, ncols, rlink, nrows,
583 mrstrt, rowels,
584 cx,
585 cz,
586 rhs / coeffy,
587 irow, icolx, icoly,icolz);
588 if (no_mem)
589 throwCoinError("out of memory",
590 "tripleton_action::presolve");
591
592 // now remove irow from icolx and icolz in the col rep
593 // better if this were first.
594 presolve_delete_from_col(irow,icolx,mcstrt,hincol,hrow,colels) ;
595 presolve_delete_from_col(irow,icolz,mcstrt,hincol,hrow,colels) ;
596
597 // eliminate irow entirely from the row rep
598 hinrow[irow] = 0;
599
600 // eliminate irow entirely from the row rep
601 PRESOLVE_REMOVE_LINK(rlink, irow);
602
603 // eliminate coly entirely from the col rep
604 PRESOLVE_REMOVE_LINK(clink, icoly);
605 cost[icoly] = 0.0;
606
607 rlo[irow] = 0.0;
608 rup[irow] = 0.0;
609
610 if (!mark[icolx]) {
611 mark[icolx]=1;
612 zeros[nzeros++]=icolx;
613 }
614 if (!mark[icolz]) {
615 mark[icolz]=1;
616 zeros[nzeros++]=icolz;
617 }
618 }
619
620# if PRESOLVE_CONSISTENCY
621 presolve_links_ok(prob) ;
622 presolve_consistent(prob);
623# endif
624 }
625 }
626 if (nactions) {
627# if PRESOLVE_SUMMARY
628 printf("NTRIPLETONS: %d\n", nactions);
629# endif
630 action *actions1 = new action[nactions];
631 CoinMemcpyN(actions, nactions, actions1);
632
633 next = new tripleton_action(nactions, actions1, next);
634
635 if (nzeros) {
636 next = drop_zero_coefficients_action::presolve(prob, zeros, nzeros, next);
637 }
638 }
639
640 //delete[]zeros;
641 deleteAction(actions,action*);
642
643 if (prob->tuning_) {
644 double thisTime=CoinCpuTime();
645 int droppedRows = prob->countEmptyRows() - startEmptyRows ;
646 int droppedColumns = prob->countEmptyCols() - startEmptyColumns;
647 printf("CoinPresolveTripleton(8) - %d rows, %d columns dropped in time %g, total %g\n",
648 droppedRows,droppedColumns,thisTime-startTime,thisTime-prob->startTime_);
649 }
650 return (next);
651}
652
653
654void tripleton_action::postsolve(CoinPostsolveMatrix *prob) const
655{
656 const action *const actions = actions_;
657 const int nactions = nactions_;
658
659 double *colels = prob->colels_;
660 int *hrow = prob->hrow_;
661 CoinBigIndex *mcstrt = prob->mcstrt_;
662 int *hincol = prob->hincol_;
663 int *link = prob->link_;
664
665 double *clo = prob->clo_;
666 double *cup = prob->cup_;
667
668 double *rlo = prob->rlo_;
669 double *rup = prob->rup_;
670
671 double *dcost = prob->cost_;
672
673 double *sol = prob->sol_;
674 double *rcosts = prob->rcosts_;
675
676 double *acts = prob->acts_;
677 double *rowduals = prob->rowduals_;
678
679 unsigned char *colstat = prob->colstat_;
680 unsigned char *rowstat = prob->rowstat_;
681
682 const double maxmin = prob->maxmin_;
683
684# if PRESOLVE_DEBUG || PRESOLVE_CONSISTENCY
685 char *cdone = prob->cdone_;
686 char *rdone = prob->rdone_;
687# endif
688
689 CoinBigIndex &free_list = prob->free_list_;
690
691 const double ztolzb = prob->ztolzb_;
692 const double ztoldj = prob->ztoldj_;
693
694 // Space for accumulating two columns
695 int nrows = prob->nrows_;
696 int * index1 = new int[nrows];
697 double * element1 = new double[nrows];
698 memset(element1,0,nrows*sizeof(double));
699 int * index2 = new int[nrows];
700 double * element2 = new double[nrows];
701 memset(element2,0,nrows*sizeof(double));
702
703 for (const action *f = &actions[nactions-1]; actions<=f; f--) {
704 int irow = f->row;
705
706 // probably don't need this
707 double ylo0 = f->cloy;
708 double yup0 = f->cupy;
709
710 double coeffx = f->coeffx;
711 double coeffy = f->coeffy;
712 double coeffz = f->coeffz;
713 int jcolx = f->icolx;
714 int jcoly = f->icoly;
715 int jcolz = f->icolz;
716
717 // needed?
718 double rhs = f->rlo;
719
720 /* the column was in the reduced problem */
721 PRESOLVEASSERT(cdone[jcolx] && rdone[irow]==DROP_ROW&&cdone[jcolz]);
722 PRESOLVEASSERT(cdone[jcoly]==DROP_COL);
723
724 // probably don't need this
725 rlo[irow] = f->rlo;
726 rup[irow] = f->rup;
727
728 // probably don't need this
729 clo[jcoly] = ylo0;
730 cup[jcoly] = yup0;
731
732 dcost[jcoly] = f->costy;
733 dcost[jcolx] += f->costy*coeffx/coeffy;
734 dcost[jcolz] += f->costy*coeffz/coeffy;
735
736 // this is why we want coeffx < coeffy (55)
737 sol[jcoly] = (rhs - coeffx * sol[jcolx] - coeffz * sol[jcolz]) / coeffy;
738
739 // since this row is fixed
740 acts[irow] = rhs;
741
742 // acts[irow] always ok, since slack is fixed
743 if (rowstat)
744 prob->setRowStatus(irow,CoinPrePostsolveMatrix::atLowerBound);
745
746
747 // CLAIM:
748 // if the new pi value is chosen to keep the reduced cost
749 // of col x at its prior value, then the reduced cost of
750 // col y will be 0.
751
752 // also have to update row activities and bounds for rows affected by jcoly
753 //
754 // sol[jcolx] was found for coeffx that
755 // was += colels[kcoly] * coeff_factor;
756 // where coeff_factor == -coeffx / coeffy
757 //
758 // its contribution to activity was
759 // (colels[kcolx] + colels[kcoly] * coeff_factor) * sol[jcolx] (1)
760 //
761 // After adjustment, the two columns contribute:
762 // colels[kcoly] * sol[jcoly] + colels[kcolx] * sol[jcolx]
763 // == colels[kcoly] * ((rhs - coeffx * sol[jcolx]) / coeffy) + colels[kcolx] * sol[jcolx]
764 // == colels[kcoly] * rhs/coeffy + colels[kcoly] * coeff_factor * sol[jcolx] + colels[kcolx] * sol[jcolx]
765 // colels[kcoly] * rhs/coeffy + the expression (1)
766 //
767 // therefore, we must increase the row bounds by colels[kcoly] * rhs/coeffy,
768 // which is similar to the bias
769 double djy = maxmin * dcost[jcoly];
770 double djx = maxmin * dcost[jcolx];
771 double djz = maxmin * dcost[jcolz];
772 double bounds_factor = rhs/coeffy;
773 // need to reconstruct x and z
774 double multiplier1 = coeffx/coeffy;
775 double multiplier2 = coeffz/coeffy;
776 int * indy = reinterpret_cast<int *>(f->colel+f->ncoly);
777 int ystart = NO_LINK;
778 int nX=0,nZ=0;
779 int i,iRow;
780 for (i=0; i<f->ncoly; ++i) {
781 int iRow = indy[i];
782 double yValue = f->colel[i];
783 CoinBigIndex k = free_list;
784 assert(k >= 0 && k < prob->bulk0_) ;
785 free_list = link[free_list];
786 if (iRow != irow) {
787 // are these tests always true???
788
789 // undo elim_tripleton(1)
790 if (-PRESOLVE_INF < rlo[iRow])
791 rlo[iRow] += yValue * bounds_factor;
792
793 // undo elim_tripleton(2)
794 if (rup[iRow] < PRESOLVE_INF)
795 rup[iRow] += yValue * bounds_factor;
796
797 acts[iRow] += yValue * bounds_factor;
798
799 djy -= rowduals[iRow] * yValue;
800 }
801
802 hrow[k] = iRow;
803 PRESOLVEASSERT(rdone[hrow[k]] || hrow[k] == irow);
804 colels[k] = yValue;
805 link[k] = ystart;
806 ystart = k;
807 element1[iRow]=yValue*multiplier1;
808 index1[nX++]=iRow;
809 element2[iRow]=yValue*multiplier2;
810 index2[nZ++]=iRow;
811 }
812# if PRESOLVE_CONSISTENCY
813 presolve_check_free_list(prob) ;
814# endif
815 mcstrt[jcoly] = ystart;
816 hincol[jcoly] = f->ncoly;
817 // find the tail
818 CoinBigIndex k=mcstrt[jcolx];
819 CoinBigIndex last = NO_LINK;
820 int numberInColumn = hincol[jcolx];
821 int numberToDo=numberInColumn;
822 for (i=0; i<numberToDo; ++i) {
823 iRow = hrow[k];
824 assert (iRow>=0&&iRow<nrows);
825 double value = colels[k]+element1[iRow];
826 element1[iRow]=0.0;
827 if (fabs(value)>=1.0e-15) {
828 colels[k]=value;
829 last=k;
830 k = link[k];
831 if (iRow != irow)
832 djx -= rowduals[iRow] * value;
833 } else {
834 numberInColumn--;
835 // add to free list
836 int nextk = link[k];
837 link[k]=free_list;
838 free_list=k;
839 assert (k>=0);
840 k=nextk;
841 if (last!=NO_LINK)
842 link[last]=k;
843 else
844 mcstrt[jcolx]=k;
845 }
846 }
847 for (i=0;i<nX;i++) {
848 int iRow = index1[i];
849 double xValue = element1[iRow];
850 element1[iRow]=0.0;
851 if (fabs(xValue)>=1.0e-15) {
852 if (iRow != irow)
853 djx -= rowduals[iRow] * xValue;
854 numberInColumn++;
855 CoinBigIndex k = free_list;
856 assert(k >= 0 && k < prob->bulk0_) ;
857 free_list = link[free_list];
858 hrow[k] = iRow;
859 PRESOLVEASSERT(rdone[hrow[k]] || hrow[k] == irow);
860 colels[k] = xValue;
861 if (last!=NO_LINK)
862 link[last]=k;
863 else
864 mcstrt[jcolx]=k;
865 last = k;
866 }
867 }
868# if PRESOLVE_CONSISTENCY
869 presolve_check_free_list(prob) ;
870# endif
871 link[last]=NO_LINK;
872 assert(numberInColumn);
873 hincol[jcolx] = numberInColumn;
874 // find the tail
875 k=mcstrt[jcolz];
876 last = NO_LINK;
877 numberInColumn = hincol[jcolz];
878 numberToDo=numberInColumn;
879 for (i=0; i<numberToDo; ++i) {
880 iRow = hrow[k];
881 assert (iRow>=0&&iRow<nrows);
882 double value = colels[k]+element2[iRow];
883 element2[iRow]=0.0;
884 if (fabs(value)>=1.0e-15) {
885 colels[k]=value;
886 last=k;
887 k = link[k];
888 if (iRow != irow)
889 djz -= rowduals[iRow] * value;
890 } else {
891 numberInColumn--;
892 // add to free list
893 int nextk = link[k];
894 assert(free_list>=0);
895 link[k]=free_list;
896 free_list=k;
897 assert (k>=0);
898 k=nextk;
899 if (last!=NO_LINK)
900 link[last]=k;
901 else
902 mcstrt[jcolz]=k;
903 }
904 }
905 for (i=0;i<nZ;i++) {
906 int iRow = index2[i];
907 double zValue = element2[iRow];
908 element2[iRow]=0.0;
909 if (fabs(zValue)>=1.0e-15) {
910 if (iRow != irow)
911 djz -= rowduals[iRow] * zValue;
912 numberInColumn++;
913 CoinBigIndex k = free_list;
914 assert(k >= 0 && k < prob->bulk0_) ;
915 free_list = link[free_list];
916 hrow[k] = iRow;
917 PRESOLVEASSERT(rdone[hrow[k]] || hrow[k] == irow);
918 colels[k] = zValue;
919 if (last!=NO_LINK)
920 link[last]=k;
921 else
922 mcstrt[jcolz]=k;
923 last = k;
924 }
925 }
926# if PRESOLVE_CONSISTENCY
927 presolve_check_free_list(prob) ;
928# endif
929 link[last]=NO_LINK;
930 assert(numberInColumn);
931 hincol[jcolz] = numberInColumn;
932
933
934
935 // The only problem with keeping the reduced costs the way they were
936 // was that the variable's bound may have moved, requiring it
937 // to become basic.
938 //printf("djs x - %g (%g), y - %g (%g)\n",djx,coeffx,djy,coeffy);
939 if (colstat) {
940 if (prob->columnIsBasic(jcolx) ||
941 (fabs(clo[jcolx] - sol[jcolx]) < ztolzb && rcosts[jcolx] >= -ztoldj) ||
942 (fabs(cup[jcolx] - sol[jcolx]) < ztolzb && rcosts[jcolx] <= ztoldj) ||
943 (prob->getColumnStatus(jcolx) ==CoinPrePostsolveMatrix::isFree&&
944 fabs(rcosts[jcolx]) <= ztoldj)) {
945 // colx or y is fine as it is - make coly basic
946
947 prob->setColumnStatus(jcoly,CoinPrePostsolveMatrix::basic);
948 // this is the coefficient we need to force col y's reduced cost to 0.0;
949 // for example, this is obviously true if y is a singleton column
950 rowduals[irow] = djy / coeffy;
951 rcosts[jcolx] = djx - rowduals[irow] * coeffx;
952# if PRESOLVE_DEBUG
953 if (prob->columnIsBasic(jcolx)&&fabs(rcosts[jcolx])>1.0e-5)
954 printf("bad dj %d %g\n",jcolx,rcosts[jcolx]);
955# endif
956 rcosts[jcolz] = djz - rowduals[irow] * coeffz;
957 //if (prob->columnIsBasic(jcolz))
958 //assert (fabs(rcosts[jcolz])<1.0e-5);
959 rcosts[jcoly] = 0.0;
960 } else {
961 prob->setColumnStatus(jcolx,CoinPrePostsolveMatrix::basic);
962 prob->setColumnStatusUsingValue(jcoly);
963
964 // change rowduals[jcolx] enough to cancel out rcosts[jcolx]
965 rowduals[irow] = djx / coeffx;
966 rcosts[jcolx] = 0.0;
967 // change rowduals[jcolx] enough to cancel out rcosts[jcolx]
968 //rowduals[irow] = djz / coeffz;
969 //rcosts[jcolz] = 0.0;
970 rcosts[jcolz] = djz - rowduals[irow] * coeffz;
971 rcosts[jcoly] = djy - rowduals[irow] * coeffy;
972 }
973 } else {
974 // No status array
975 // this is the coefficient we need to force col y's reduced cost to 0.0;
976 // for example, this is obviously true if y is a singleton column
977 rowduals[irow] = djy / coeffy;
978 rcosts[jcoly] = 0.0;
979 }
980
981 // DEBUG CHECK
982# if PRESOLVE_DEBUG
983 {
984 CoinBigIndex k = mcstrt[jcolx];
985 int nx = hincol[jcolx];
986 double dj = maxmin * dcost[jcolx];
987
988 for (int i=0; i<nx; ++i) {
989 int row = hrow[k];
990 double coeff = colels[k];
991 k = link[k];
992
993 dj -= rowduals[row] * coeff;
994 }
995 if (! (fabs(rcosts[jcolx] - dj) < 100*ZTOLDP))
996 printf("BAD DOUBLE X DJ: %d %d %g %g\n",
997 irow, jcolx, rcosts[jcolx], dj);
998 rcosts[jcolx]=dj;
999 }
1000 {
1001 CoinBigIndex k = mcstrt[jcoly];
1002 int ny = hincol[jcoly];
1003 double dj = maxmin * dcost[jcoly];
1004
1005 for (int i=0; i<ny; ++i) {
1006 int row = hrow[k];
1007 double coeff = colels[k];
1008 k = link[k];
1009
1010 dj -= rowduals[row] * coeff;
1011 //printf("b %d coeff %g dual %g dj %g\n",
1012 // row,coeff,rowduals[row],dj);
1013 }
1014 if (! (fabs(rcosts[jcoly] - dj) < 100*ZTOLDP))
1015 printf("BAD DOUBLE Y DJ: %d %d %g %g\n",
1016 irow, jcoly, rcosts[jcoly], dj);
1017 rcosts[jcoly]=dj;
1018 //exit(0);
1019 }
1020# endif
1021
1022# if PRESOLVE_DEBUG || PRESOLVE_CONSISTENCY
1023 cdone[jcoly] = TRIPLETON;
1024 rdone[irow] = TRIPLETON;
1025# endif
1026 }
1027 delete [] index1;
1028 delete [] element1;
1029 delete [] index2;
1030 delete [] element2;
1031
1032# if PRESOLVE_CONSISTENCY
1033 presolve_check_threads(prob) ;
1034# endif
1035
1036}
1037
1038
1039tripleton_action::~tripleton_action()
1040{
1041 for (int i=nactions_-1; i>=0; i--) {
1042 delete[]actions_[i].colel;
1043 }
1044 deleteAction(actions_,action*);
1045}
1046
1047
1048
1049static double *tripleton_mult;
1050static int *tripleton_id;
1051void check_tripletons(const CoinPresolveAction * paction)
1052{
1053 const CoinPresolveAction * paction0 = paction;
1054
1055 if (paction) {
1056 check_tripletons(paction->next);
1057
1058 if (strcmp(paction0->name(), "tripleton_action") == 0) {
1059 const tripleton_action *daction = reinterpret_cast<const tripleton_action *>(paction0);
1060 for (int i=daction->nactions_-1; i>=0; --i) {
1061 int icolx = daction->actions_[i].icolx;
1062 int icoly = daction->actions_[i].icoly;
1063 double coeffx = daction->actions_[i].coeffx;
1064 double coeffy = daction->actions_[i].coeffy;
1065
1066 tripleton_mult[icoly] = -coeffx/coeffy;
1067 tripleton_id[icoly] = icolx;
1068 }
1069 }
1070 }
1071}
1072