1/* $Id: CoinPresolveSingleton.cpp 1463 2011-07-30 16:31:31Z 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 "CoinHelperFunctions.hpp"
10#include "CoinPresolveMatrix.hpp"
11
12#include "CoinPresolveEmpty.hpp" // for DROP_COL/DROP_ROW
13#include "CoinPresolveFixed.hpp"
14#include "CoinPresolveSingleton.hpp"
15#include "CoinPresolvePsdebug.hpp"
16#include "CoinMessage.hpp"
17#include "CoinFinite.hpp"
18
19// #define PRESOLVE_DEBUG 1
20// #define PRESOLVE_SUMMARY 1
21
22/*
23 * Transfers singleton row bound information to the corresponding column bounds.
24 * What I refer to as a row singleton would be called a doubleton
25 * in the paper, since my terminology doesn't refer to the slacks.
26 * In terms of the paper, we transfer the bounds of the slack onto
27 * the variable (vii) and then "substitute" the slack out of the problem
28 * (which is a noop).
29 */
30const CoinPresolveAction *
31slack_doubleton_action::presolve(CoinPresolveMatrix *prob,
32 const CoinPresolveAction *next,
33 bool & notFinished)
34{
35 double startTime = 0.0;
36 int startEmptyRows=0;
37 int startEmptyColumns = 0;
38 if (prob->tuning_) {
39 startTime = CoinCpuTime();
40 startEmptyRows = prob->countEmptyRows();
41 startEmptyColumns = prob->countEmptyCols();
42 }
43 double *colels = prob->colels_;
44 int *hrow = prob->hrow_;
45 CoinBigIndex *mcstrt = prob->mcstrt_;
46 int *hincol = prob->hincol_;
47 //int ncols = prob->ncols_;
48
49 double *clo = prob->clo_;
50 double *cup = prob->cup_;
51
52 double *rowels = prob->rowels_;
53 const int *hcol = prob->hcol_;
54 const CoinBigIndex *mrstrt = prob->mrstrt_;
55 int *hinrow = prob->hinrow_;
56 // int nrows = prob->nrows_;
57
58 double *rlo = prob->rlo_;
59 double *rup = prob->rup_;
60
61 // If rowstat exists then all do
62 unsigned char *rowstat = prob->rowstat_;
63 double *acts = prob->acts_;
64 double * sol = prob->sol_;
65 // unsigned char * colstat = prob->colstat_;
66
67# if PRESOLVE_DEBUG
68 std::cout << "Entering slack_doubleton_action::presolve." << std::endl ;
69 presolve_check_sol(prob) ;
70 presolve_check_nbasic(prob) ;
71# endif
72
73 const unsigned char *integerType = prob->integerType_;
74
75 const double ztolzb = prob->ztolzb_;
76
77 int numberLook = prob->numberRowsToDo_;
78 int iLook;
79 int * look = prob->rowsToDo_;
80 bool fixInfeasibility = (prob->presolveOptions_&16384)!=0;
81
82 action * actions = new action[numberLook];
83 int nactions = 0;
84 notFinished = false;
85
86 int * fixed_cols = prob->usefulColumnInt_; //new int [ncols];
87 int nfixed_cols = 0;
88
89 // this loop can apparently create new singleton rows;
90 // I haven't bothered to detect this.
91 // wasfor (int irow=0; irow<nrows; irow++)
92 for (iLook=0;iLook<numberLook;iLook++) {
93 int irow = look[iLook];
94 if (hinrow[irow] == 1) {
95 int jcol = hcol[mrstrt[irow]];
96 double coeff = rowels[mrstrt[irow]];
97 double lo = rlo[irow];
98 double up = rup[irow];
99 double acoeff = fabs(coeff);
100 // const bool singleton_col = (hincol[jcol] == 1);
101
102 if (acoeff < ZTOLDP2)
103 continue;
104
105 // don't bother with fixed cols
106 if (fabs(cup[jcol] - clo[jcol]) < ztolzb)
107 continue;
108
109 {
110 // put column on stack of things to do next time
111 prob->addCol(jcol);
112 PRESOLVE_DETAIL_PRINT(printf("pre_singleton %dC %dR E\n",jcol,irow));
113 action *s = &actions[nactions];
114 nactions++;
115
116 s->col = jcol;
117 s->clo = clo[jcol];
118 s->cup = cup[jcol];
119
120 s->row = irow;
121 s->rlo = rlo[irow];
122 s->rup = rup[irow];
123
124 s->coeff = coeff;
125#if 0
126 if (prob->tuning_) {
127 // really for gcc 4.6 compiler bug
128 printf("jcol %d %g %g irow %d %g %g coeff %g\n",
129 jcol,clo[jcol],cup[jcol],irow,rlo[irow],rup[irow],coeff);
130 }
131#endif
132 }
133
134 if (coeff < 0.0) {
135 CoinSwap(lo, up);
136 lo = -lo;
137 up = -up;
138 }
139
140 if (lo <= -PRESOLVE_INF)
141 lo = -PRESOLVE_INF;
142 else {
143 lo /= acoeff;
144 if (lo <= -PRESOLVE_INF)
145 lo = -PRESOLVE_INF;
146 }
147
148 if (up > PRESOLVE_INF)
149 up = PRESOLVE_INF;
150 else {
151 up /= acoeff;
152 if (up > PRESOLVE_INF)
153 up = PRESOLVE_INF;
154 }
155
156 if (clo[jcol] < lo) {
157 // If integer be careful
158 if (integerType[jcol]) {
159 //#define COIN_DEVELOP
160#ifdef COIN_DEVELOP
161 double lo2=lo;
162#endif
163 if (fabs(lo-floor(lo+0.5))<0.000001)
164 lo=floor(lo+0.5);
165#ifdef COIN_DEVELOP
166 if (lo!=lo2&&fabs(lo-floor(lo+0.5))<0.01)
167 printf("first lo %g second %g orig %g\n",lo2,lo,clo[jcol]);
168#endif
169 if (clo[jcol] < lo)
170 clo[jcol] = lo;
171 } else {
172 clo[jcol] = lo;
173 }
174 }
175
176 if (cup[jcol] > up) {
177 // If integer be careful
178 if (integerType[jcol]) {
179#ifdef COIN_DEVELOP
180 double up2=up;
181#endif
182 if (fabs(up-floor(up+0.5))<0.000001)
183 up=floor(up+0.5);
184#ifdef COIN_DEVELOP
185 if (up!=up2&&fabs(up-floor(up+0.5))<0.01)
186 printf("first up %g second %g orig %g\n",up2,up,cup[jcol]);
187#endif
188 if (cup[jcol] > up)
189 cup[jcol] = up;
190 } else {
191 cup[jcol] = up;
192 }
193 }
194 if (fabs(cup[jcol] - clo[jcol]) < ZTOLDP) {
195 fixed_cols[nfixed_cols++] = jcol;
196 }
197
198 if (lo > up) {
199 if (lo <= up + prob->feasibilityTolerance_||fixInfeasibility) {
200 // If close to integer then go there
201 double nearest = floor(lo+0.5);
202 if (fabs(nearest-lo)<2.0*prob->feasibilityTolerance_) {
203 lo = nearest;
204 up = nearest;
205 } else {
206 lo = up;
207 }
208 clo[jcol] = lo;
209 cup[jcol] = up;
210 } else {
211 prob->status_ |= 1;
212 prob->messageHandler()->message(COIN_PRESOLVE_COLINFEAS,
213 prob->messages())
214 <<jcol
215 <<lo
216 <<up
217 <<CoinMessageEol;
218 break;
219 }
220 }
221
222# if PRESOLVE_DEBUG
223 printf("SINGLETON R-%d C-%d\n", irow, jcol);
224# endif
225
226 // eliminate the row entirely from the row rep
227 hinrow[irow] = 0;
228 PRESOLVE_REMOVE_LINK(prob->rlink_,irow) ;
229
230 // just to make things squeeky
231 rlo[irow] = 0.0;
232 rup[irow] = 0.0;
233
234 if (rowstat&&sol) {
235 // update solution and basis
236 int basisChoice=0;
237 int numberBasic=0;
238 double movement = 0 ;
239 if (prob->columnIsBasic(jcol)) {
240 numberBasic++;
241 basisChoice=2; // move to row to keep consistent
242 }
243 if (prob->rowIsBasic(irow))
244 numberBasic++;
245 if (sol[jcol] <= clo[jcol]+ztolzb) {
246 movement = clo[jcol]-sol[jcol] ;
247 sol[jcol] = clo[jcol];
248 prob->setColumnStatus(jcol,CoinPrePostsolveMatrix::atLowerBound);
249 } else if (sol[jcol] >= cup[jcol]-ztolzb) {
250 movement = cup[jcol]-sol[jcol] ;
251 sol[jcol] = cup[jcol];
252 prob->setColumnStatus(jcol,CoinPrePostsolveMatrix::atUpperBound);
253 } else {
254 basisChoice=1;
255 }
256 if (numberBasic>1||basisChoice==1)
257 prob->setColumnStatus(jcol,CoinPrePostsolveMatrix::basic);
258 else if (basisChoice==2)
259 prob->setRowStatus(irow,CoinPrePostsolveMatrix::basic);
260 if (movement) {
261 CoinBigIndex k;
262 for (k = mcstrt[jcol] ; k < mcstrt[jcol]+hincol[jcol] ; k++) {
263 int row = hrow[k];
264 if (hinrow[row])
265 acts[row] += movement*colels[k];
266 }
267 }
268 }
269
270/*
271 Remove the row from this col in the col rep. It can happen that this will
272 empty the column, in which case we can delink it.
273*/
274 presolve_delete_from_col(irow,jcol,mcstrt,hincol,hrow,colels) ;
275 if (hincol[jcol] == 0)
276 { PRESOLVE_REMOVE_LINK(prob->clink_,jcol) ; }
277
278 if (nactions >= numberLook) {
279 notFinished=true;
280 break;
281 }
282 }
283 }
284
285 if (nactions) {
286# if PRESOLVE_SUMMARY
287 printf("SINGLETON ROWS: %d\n", nactions);
288# endif
289 action *save_actions = new action[nactions];
290 CoinMemcpyN(actions, nactions, save_actions);
291 next = new slack_doubleton_action(nactions, save_actions, next);
292
293 if (nfixed_cols)
294 next = make_fixed_action::presolve(prob, fixed_cols, nfixed_cols,
295 true, // arbitrary
296 next);
297 }
298 //delete [] fixed_cols;
299 if (prob->tuning_) {
300 double thisTime=CoinCpuTime();
301 int droppedRows = prob->countEmptyRows() - startEmptyRows ;
302 int droppedColumns = prob->countEmptyCols() - startEmptyColumns;
303 printf("CoinPresolveSingleton(2) - %d rows, %d columns dropped in time %g, total %g\n",
304 droppedRows,droppedColumns,thisTime-startTime,thisTime-prob->startTime_);
305 }
306
307# if PRESOLVE_DEBUG
308 presolve_check_sol(prob) ;
309 presolve_check_nbasic(prob) ;
310 std::cout << "Leaving slack_doubleton_action::presolve." << std::endl ;
311# endif
312 delete [] actions;
313
314 return (next);
315}
316
317void slack_doubleton_action::postsolve(CoinPostsolveMatrix *prob) const
318{
319 const action *const actions = actions_;
320 const int nactions = nactions_;
321
322 double *colels = prob->colels_;
323 int *hrow = prob->hrow_;
324 CoinBigIndex *mcstrt = prob->mcstrt_;
325 int *hincol = prob->hincol_;
326 int *link = prob->link_;
327 // int ncols = prob->ncols_;
328
329 double *clo = prob->clo_;
330 double *cup = prob->cup_;
331
332 double *rlo = prob->rlo_;
333 double *rup = prob->rup_;
334
335 double *sol = prob->sol_;
336 double *rcosts = prob->rcosts_;
337
338 double *acts = prob->acts_;
339 double *rowduals = prob->rowduals_;
340
341 unsigned char *colstat = prob->colstat_;
342 // unsigned char *rowstat = prob->rowstat_;
343
344# if PRESOLVE_DEBUG
345 char *rdone = prob->rdone_;
346
347 std::cout << "Entering slack_doubleton_action::postsolve." << std::endl ;
348 presolve_check_sol(prob) ;
349 presolve_check_nbasic(prob) ;
350# endif
351 CoinBigIndex &free_list = prob->free_list_;
352
353 const double ztolzb = prob->ztolzb_;
354
355 for (const action *f = &actions[nactions-1]; actions<=f; f--) {
356 int irow = f->row;
357 double lo0 = f->clo;
358 double up0 = f->cup;
359 double coeff = f->coeff;
360 int jcol = f->col;
361
362 rlo[irow] = f->rlo;
363 rup[irow] = f->rup;
364
365 clo[jcol] = lo0;
366 cup[jcol] = up0;
367
368 acts[irow] = coeff * sol[jcol];
369
370 // add new element
371 {
372 CoinBigIndex k = free_list;
373 assert(k >= 0 && k < prob->bulk0_) ;
374 free_list = link[free_list];
375 hrow[k] = irow;
376 colels[k] = coeff;
377 link[k] = mcstrt[jcol];
378 mcstrt[jcol] = k;
379 }
380 hincol[jcol]++; // right?
381
382 /*
383 * Have to compute rowstat and rowduals, since we are adding the row.
384 * that satisfy complentarity slackness.
385 *
386 * We may also have to modify colstat and rcosts if bounds
387 * have been relaxed.
388 */
389 if (!colstat) {
390 // ????
391 rowduals[irow] = 0.0;
392 } else {
393 if (prob->columnIsBasic(jcol)) {
394 /* variable is already basic, so slack in this row must be basic */
395 prob->setRowStatus(irow,CoinPrePostsolveMatrix::basic);
396
397 rowduals[irow] = 0.0;
398 /* hence no reduced costs change */
399 /* since the column was basic, it doesn't matter if it is now
400 away from its bounds. */
401 /* the slack is basic and its reduced cost is 0 */
402 } else if ((fabs(sol[jcol] - lo0) <= ztolzb &&
403 rcosts[jcol] >= 0) ||
404
405 (fabs(sol[jcol] - up0) <= ztolzb &&
406 rcosts[jcol] <= 0)) {
407 /* up against its bound but the rcost is still ok */
408 prob->setRowStatus(irow,CoinPrePostsolveMatrix::basic);
409
410 rowduals[irow] = 0.0;
411 /* hence no reduced costs change */
412 } else if (! (fabs(sol[jcol] - lo0) <= ztolzb) &&
413 ! (fabs(sol[jcol] - up0) <= ztolzb)) {
414 /* variable not marked basic,
415 * so was at a bound in the reduced problem,
416 * but its bounds were tighter in the reduced problem,
417 * so now it has to be made basic.
418 */
419 prob->setColumnStatus(jcol,CoinPrePostsolveMatrix::basic);
420 prob->setRowStatusUsingValue(irow);
421
422 /* choose a value for rowduals[irow] that forces rcosts[jcol] to 0.0 */
423 /* new reduced cost = 0 = old reduced cost - new dual * coeff */
424 rowduals[irow] = rcosts[jcol] / coeff;
425 rcosts[jcol] = 0.0;
426
427 /* this value is provably of the right sign for the slack */
428 /* SHOULD SHOW THIS */
429 } else {
430 /* the solution is at a bound, but the rcost is wrong */
431
432 prob->setColumnStatus(jcol,CoinPrePostsolveMatrix::basic);
433 prob->setRowStatusUsingValue(irow);
434
435 /* choose a value for rowduals[irow] that forcesd rcosts[jcol] to 0.0 */
436 rowduals[irow] = rcosts[jcol] / coeff;
437 rcosts[jcol] = 0.0;
438
439 /* this value is provably of the right sign for the slack */
440 /*rcosts[irow] = -rowduals[irow];*/
441 }
442 }
443
444# if PRESOLVE_DEBUG
445 rdone[irow] = SLACK_DOUBLETON;
446# endif
447 }
448
449# if PRESOLVE_CONSISTENCY
450 presolve_check_threads(prob) ;
451 presolve_check_sol(prob) ;
452 presolve_check_nbasic(prob) ;
453 std::cout << "Leaving slack_doubleton_action::postsolve." << std::endl ;
454# endif
455
456 return ;
457}
458/*
459 If we have a variable with one entry and no cost then we can
460 transform the row from E to G etc.
461 If there is a row objective region then we may be able to do
462 this even with a cost.
463*/
464const CoinPresolveAction *
465slack_singleton_action::presolve(CoinPresolveMatrix *prob,
466 const CoinPresolveAction *next,
467 double * rowObjective)
468{
469 double startTime = 0.0;
470 int startEmptyRows=0;
471 int startEmptyColumns = 0;
472 if (prob->tuning_) {
473 startTime = CoinCpuTime();
474 startEmptyRows = prob->countEmptyRows();
475 startEmptyColumns = prob->countEmptyCols();
476 }
477 double *colels = prob->colels_;
478 int *hrow = prob->hrow_;
479 CoinBigIndex *mcstrt = prob->mcstrt_;
480 int *hincol = prob->hincol_;
481 //int ncols = prob->ncols_;
482
483 double *clo = prob->clo_;
484 double *cup = prob->cup_;
485
486 double *rowels = prob->rowels_;
487 int *hcol = prob->hcol_;
488 CoinBigIndex *mrstrt = prob->mrstrt_;
489 int *hinrow = prob->hinrow_;
490 int nrows = prob->nrows_;
491
492 double *rlo = prob->rlo_;
493 double *rup = prob->rup_;
494
495 // If rowstat exists then all do
496 unsigned char *rowstat = prob->rowstat_;
497 double *acts = prob->acts_;
498 double * sol = prob->sol_;
499 //unsigned char * colstat = prob->colstat_;
500
501 const unsigned char *integerType = prob->integerType_;
502
503 const double ztolzb = prob->ztolzb_;
504 double *dcost = prob->cost_;
505 //const double maxmin = prob->maxmin_;
506
507# if PRESOLVE_DEBUG
508 std::cout << "Entering slack_singleton_action::presolve." << std::endl ;
509 presolve_check_sol(prob) ;
510 presolve_check_nbasic(prob) ;
511# endif
512
513 int numberLook = prob->numberColsToDo_;
514 int iLook;
515 int * look = prob->colsToDo_;
516 // Make sure we allocate at least one action
517 int maxActions = CoinMin(numberLook,nrows/10)+1;
518 action * actions = new action[maxActions];
519 int nactions = 0;
520 int * fixed_cols = new int [numberLook];
521 int nfixed_cols=0;
522 int nWithCosts=0;
523 double costOffset=0.0;
524 for (iLook=0;iLook<numberLook;iLook++) {
525 int iCol = look[iLook];
526 if (dcost[iCol])
527 continue;
528 if (hincol[iCol] == 1) {
529 int iRow=hrow[mcstrt[iCol]];
530 double coeff = colels[mcstrt[iCol]];
531 double acoeff = fabs(coeff);
532 if (acoeff < ZTOLDP2)
533 continue;
534 // don't bother with fixed cols
535 if (fabs(cup[iCol] - clo[iCol]) < ztolzb)
536 continue;
537 if (integerType&&integerType[iCol]) {
538 // only possible if everything else integer and unit coefficient
539 // check everything else a bit later
540 if (acoeff!=1.0)
541 continue;
542 double currentLower = rlo[iRow];
543 double currentUpper = rup[iRow];
544 if (coeff==1.0&&currentLower==1.0&&currentUpper==1.0) {
545 // leave if integer slack on sum x == 1
546 bool allInt=true;
547 for (CoinBigIndex j=mrstrt[iRow];
548 j<mrstrt[iRow]+hinrow[iRow];j++) {
549 int iColumn = hcol[j];
550 double value = fabs(rowels[j]);
551 if (!integerType[iColumn]||value!=1.0) {
552 allInt=false;
553 break;
554 }
555 }
556 if (allInt)
557 continue; // leave as may help search
558 }
559 }
560 if (!prob->colProhibited(iCol)) {
561 double currentLower = rlo[iRow];
562 double currentUpper = rup[iRow];
563 if (!rowObjective) {
564 if (dcost[iCol])
565 continue;
566 } else if ((dcost[iCol]&&currentLower!=currentUpper)||rowObjective[iRow]) {
567 continue;
568 }
569 double newLower=currentLower;
570 double newUpper=currentUpper;
571 if (coeff<0.0) {
572 if (currentUpper>1.0e20||cup[iCol]>1.0e20) {
573 newUpper=COIN_DBL_MAX;
574 } else {
575 newUpper -= coeff*cup[iCol];
576 if (newUpper>1.0e20)
577 newUpper=COIN_DBL_MAX;
578 }
579 if (currentLower<-1.0e20||clo[iCol]<-1.0e20) {
580 newLower=-COIN_DBL_MAX;
581 } else {
582 newLower -= coeff*clo[iCol];
583 if (newLower<-1.0e20)
584 newLower=-COIN_DBL_MAX;
585 }
586 } else {
587 if (currentUpper>1.0e20||clo[iCol]<-1.0e20) {
588 newUpper=COIN_DBL_MAX;
589 } else {
590 newUpper -= coeff*clo[iCol];
591 if (newUpper>1.0e20)
592 newUpper=COIN_DBL_MAX;
593 }
594 if (currentLower<-1.0e20||cup[iCol]>1.0e20) {
595 newLower=-COIN_DBL_MAX;
596 } else {
597 newLower -= coeff*cup[iCol];
598 if (newLower<-1.0e20)
599 newLower=-COIN_DBL_MAX;
600 }
601 }
602 if (integerType&&integerType[iCol]) {
603 // only possible if everything else integer
604 if (newLower>-1.0e30) {
605 if (newLower!=floor(newLower+0.5))
606 continue;
607 }
608 if (newUpper<1.0e30) {
609 if (newUpper!=floor(newUpper+0.5))
610 continue;
611 }
612 bool allInt=true;
613 for (CoinBigIndex j=mrstrt[iRow];
614 j<mrstrt[iRow]+hinrow[iRow];j++) {
615 int iColumn = hcol[j];
616 double value = fabs(rowels[j]);
617 if (!integerType[iColumn]||value!=floor(value+0.5)) {
618 allInt=false;
619 break;
620 }
621 }
622 if (!allInt)
623 continue; // no good
624 }
625 if (nactions>=maxActions) {
626 maxActions += CoinMin(numberLook-iLook,maxActions);
627 action * temp = new action[maxActions];
628 memcpy(temp,actions,nactions*sizeof(action));
629 // changed as 4.6 compiler bug! CoinMemcpyN(actions,nactions,temp);
630 delete [] actions;
631 actions=temp;
632 }
633
634 action *s = &actions[nactions];
635 nactions++;
636
637 s->col = iCol;
638 s->clo = clo[iCol];
639 s->cup = cup[iCol];
640
641 s->row = iRow;
642 s->rlo = rlo[iRow];
643 s->rup = rup[iRow];
644
645 s->coeff = coeff;
646
647 presolve_delete_from_row(iRow,iCol,mrstrt,hinrow,hcol,rowels) ;
648 if (!hinrow[iRow])
649 PRESOLVE_REMOVE_LINK(prob->rlink_,iRow) ;
650 // put row on stack of things to do next time
651 prob->addRow(iRow);
652#ifdef PRINTCOST
653 if (rowObjective&&dcost[iCol]) {
654 printf("Singleton %d had coeff of %g in row %d - bounds %g %g - cost %g\n",
655 iCol,coeff,iRow,clo[iCol],cup[iCol],dcost[iCol]);
656 printf("Row bounds were %g %g now %g %g\n",
657 rlo[iRow],rup[iRow],newLower,newUpper);
658 }
659#endif
660 // Row may be redundant but let someone else do that
661 rlo[iRow]=newLower;
662 rup[iRow]=newUpper;
663 if (rowstat&&sol) {
664 // update solution and basis
665 if ((sol[iCol] < cup[iCol]-ztolzb&&
666 sol[iCol] > clo[iCol]+ztolzb)||prob->columnIsBasic(iCol))
667 prob->setRowStatus(iRow,CoinPrePostsolveMatrix::basic);
668 prob->setColumnStatusUsingValue(iCol);
669 }
670 // Force column to zero
671 clo[iCol]=0.0;
672 cup[iCol]=0.0;
673 if (rowObjective&&dcost[iCol]) {
674 rowObjective[iRow]=-dcost[iCol]/coeff;
675 nWithCosts++;
676 // adjust offset
677 costOffset += currentLower*rowObjective[iRow];
678 prob->dobias_ -= currentLower*rowObjective[iRow];
679 }
680 if (sol) {
681 double movement;
682 if (fabs(sol[iCol]-clo[iCol])<fabs(sol[iCol]-cup[iCol])) {
683 movement = clo[iCol]-sol[iCol] ;
684 sol[iCol]=clo[iCol];
685 } else {
686 movement = cup[iCol]-sol[iCol] ;
687 sol[iCol]=cup[iCol];
688 }
689 if (movement)
690 acts[iRow] += movement*coeff;
691 }
692 /*
693 Remove the row from this col in the col rep.and delink it.
694 */
695 presolve_delete_from_col(iRow,iCol,mcstrt,hincol,hrow,colels) ;
696 assert (hincol[iCol] == 0);
697 PRESOLVE_REMOVE_LINK(prob->clink_,iCol) ;
698 //clo[iCol] = 0.0;
699 //cup[iCol] = 0.0;
700 fixed_cols[nfixed_cols++] = iCol;
701 //presolve_consistent(prob);
702 }
703 }
704 }
705
706 if (nactions) {
707# if PRESOLVE_SUMMARY
708 printf("SINGLETON COLS: %d\n", nactions);
709# endif
710#ifdef COIN_DEVELOP
711 printf("%d singletons, %d with costs - offset %g\n",nactions,
712 nWithCosts, costOffset);
713#endif
714 action *save_actions = new action[nactions];
715 CoinMemcpyN(actions, nactions, save_actions);
716 next = new slack_singleton_action(nactions, save_actions, next);
717
718 if (nfixed_cols)
719 next = make_fixed_action::presolve(prob, fixed_cols, nfixed_cols,
720 true, // arbitrary
721 next);
722 }
723 delete [] actions;
724 delete [] fixed_cols;
725 if (prob->tuning_) {
726 double thisTime=CoinCpuTime();
727 int droppedRows = prob->countEmptyRows() - startEmptyRows ;
728 int droppedColumns = prob->countEmptyCols() - startEmptyColumns;
729 printf("CoinPresolveSingleton(3) - %d rows, %d columns dropped in time %g, total %g\n",
730 droppedRows,droppedColumns,thisTime-startTime,thisTime-prob->startTime_);
731 }
732
733# if PRESOLVE_DEBUG
734 presolve_check_sol(prob) ;
735 presolve_check_nbasic(prob) ;
736 std::cout << "Leaving slack_singleton_action::presolve." << std::endl ;
737# endif
738
739 return (next);
740}
741
742void slack_singleton_action::postsolve(CoinPostsolveMatrix *prob) const
743{
744 const action *const actions = actions_;
745 const int nactions = nactions_;
746
747 double *colels = prob->colels_;
748 int *hrow = prob->hrow_;
749 CoinBigIndex *mcstrt = prob->mcstrt_;
750 int *hincol = prob->hincol_;
751 int *link = prob->link_;
752 // int ncols = prob->ncols_;
753
754 //double *rowels = prob->rowels_;
755 //int *hcol = prob->hcol_;
756 //CoinBigIndex *mrstrt = prob->mrstrt_;
757 //int *hinrow = prob->hinrow_;
758
759 double *clo = prob->clo_;
760 double *cup = prob->cup_;
761
762 double *rlo = prob->rlo_;
763 double *rup = prob->rup_;
764
765 double *sol = prob->sol_;
766 double *rcosts = prob->rcosts_;
767
768 double *acts = prob->acts_;
769 double *rowduals = prob->rowduals_;
770 double *dcost = prob->cost_;
771 //const double maxmin = prob->maxmin_;
772
773 unsigned char *colstat = prob->colstat_;
774 // unsigned char *rowstat = prob->rowstat_;
775
776# if PRESOLVE_DEBUG
777 char *rdone = prob->rdone_;
778
779 std::cout << "Entering slack_singleton_action::postsolve." << std::endl ;
780 presolve_check_sol(prob) ;
781 presolve_check_nbasic(prob) ;
782# endif
783
784 CoinBigIndex &free_list = prob->free_list_;
785
786 const double ztolzb = prob->ztolzb_;
787#ifdef CHECK_ONE_ROW
788 {
789 double act=0.0;
790 for (int i=0;i<prob->ncols_;i++) {
791 double solV = sol[i];
792 assert (solV>=clo[i]-ztolzb&&solV<=cup[i]+ztolzb);
793 int j=mcstrt[i];
794 for (int k=0;k<hincol[i];k++) {
795 if (hrow[j]==CHECK_ONE_ROW) {
796 act += colels[j]*solV;
797 }
798 j=link[j];
799 }
800 }
801 assert (act>=rlo[CHECK_ONE_ROW]-ztolzb&&act<=rup[CHECK_ONE_ROW]+ztolzb);
802 printf("start %g %g %g %g\n",rlo[CHECK_ONE_ROW],act,acts[CHECK_ONE_ROW],rup[CHECK_ONE_ROW]);
803 }
804#endif
805 for (const action *f = &actions[nactions-1]; actions<=f; f--) {
806 int iRow = f->row;
807 double lo0 = f->clo;
808 double up0 = f->cup;
809 double coeff = f->coeff;
810 int iCol = f->col;
811 assert (!hincol[iCol]);
812#ifdef CHECK_ONE_ROW
813 if (iRow==CHECK_ONE_ROW)
814 printf("Col %d coeff %g old bounds %g,%g new %g,%g - new rhs %g,%g - act %g\n",
815 iCol,coeff,clo[iCol],cup[iCol],lo0,up0,f->rlo,f->rup,acts[CHECK_ONE_ROW]);
816#endif
817 rlo[iRow] = f->rlo;
818 rup[iRow] = f->rup;
819
820 clo[iCol] = lo0;
821 cup[iCol] = up0;
822 double movement=0.0;
823 // acts was without coefficient - adjust
824 acts[iRow] += coeff*sol[iCol];
825 if (acts[iRow]<rlo[iRow]-ztolzb)
826 movement = rlo[iRow]-acts[iRow];
827 else if (acts[iRow]>rup[iRow]+ztolzb)
828 movement = rup[iRow]-acts[iRow];
829 double cMove = movement/coeff;
830 sol[iCol] += cMove;
831 acts[iRow] += movement;
832 if (!dcost[iCol]) {
833 // and to get column feasible
834 cMove=0.0;
835 if (sol[iCol]>cup[iCol]+ztolzb)
836 cMove = cup[iCol]-sol[iCol];
837 else if (sol[iCol]<clo[iCol]-ztolzb)
838 cMove = clo[iCol]-sol[iCol];
839 sol[iCol] += cMove;
840 acts[iRow] += cMove*coeff;
841 /*
842 * Have to compute status. At most one can be basic. It's possible that
843 both are nonbasic and nonbasic status must change.
844 */
845 if (colstat) {
846 int numberBasic =0;
847 if (prob->columnIsBasic(iCol))
848 numberBasic++;
849 if (prob->rowIsBasic(iRow))
850 numberBasic++;
851#ifdef COIN_DEVELOP
852 if (numberBasic>1)
853 printf("odd in singleton\n");
854#endif
855 if (sol[iCol]>clo[iCol]+ztolzb&&sol[iCol]<cup[iCol]-ztolzb) {
856 prob->setColumnStatus(iCol,CoinPrePostsolveMatrix::basic);
857 prob->setRowStatusUsingValue(iRow);
858 } else if (acts[iRow]>rlo[iRow]+ztolzb&&acts[iRow]<rup[iRow]-ztolzb) {
859 prob->setRowStatus(iRow,CoinPrePostsolveMatrix::basic);
860 prob->setColumnStatusUsingValue(iCol);
861 } else if (numberBasic) {
862 prob->setRowStatus(iRow,CoinPrePostsolveMatrix::basic);
863 prob->setColumnStatusUsingValue(iCol);
864 } else {
865 prob->setRowStatusUsingValue(iRow);
866 prob->setColumnStatusUsingValue(iCol);
867 }
868 }
869# if PRESOLVE_DEBUG
870 printf("SLKSING: %d = %g restored %d lb = %g ub = %g.\n",
871 iCol,sol[iCol],prob->getColumnStatus(iCol),clo[iCol],cup[iCol]) ;
872# endif
873 } else {
874 // must have been equality row
875 assert (rlo[iRow]==rup[iRow]);
876 double cost = rcosts[iCol];
877 // adjust for coefficient
878 cost -= rowduals[iRow]*coeff;
879 bool basic=true;
880 if (fabs(sol[iCol]-cup[iCol])<ztolzb&&cost<-1.0e-6)
881 basic=false;
882 else if (fabs(sol[iCol]-clo[iCol])<ztolzb&&cost>1.0e-6)
883 basic=false;
884 //printf("Singleton %d had coeff of %g in row %d (dual %g) - bounds %g %g - cost %g, (dj %g)\n",
885 // iCol,coeff,iRow,rowduals[iRow],clo[iCol],cup[iCol],dcost[iCol],rcosts[iCol]);
886 //if (prob->columnIsBasic(iCol))
887 //printf("column basic! ");
888 //if (prob->rowIsBasic(iRow))
889 //printf("row basic ");
890 //printf("- make column basic %s\n",basic ? "yes" : "no");
891 if (basic&&!prob->rowIsBasic(iRow)) {
892#ifdef PRINTCOST
893 printf("Singleton %d had coeff of %g in row %d (dual %g) - bounds %g %g - cost %g, (dj %g - new %g)\n",
894 iCol,coeff,iRow,rowduals[iRow],clo[iCol],cup[iCol],dcost[iCol],rcosts[iCol],cost);
895#endif
896#ifdef COIN_DEVELOP
897 if (prob->columnIsBasic(iCol))
898 printf("column basic!\n");
899#endif
900 basic=false;
901 }
902 if (fabs(rowduals[iRow])>1.0e-6&&prob->rowIsBasic(iRow))
903 basic=true;
904 if (basic) {
905 // Make basic have zero reduced cost
906 rowduals[iRow] = rcosts[iCol] / coeff;
907 rcosts[iCol] = 0.0;
908 } else {
909 rcosts[iCol]=cost;
910 //rowduals[iRow]=0.0;
911 }
912 if (colstat) {
913 if (basic) {
914 if (!prob->rowIsBasic(iRow)) {
915#if 0
916 // find column in row
917 int jCol=-1;
918 //for (CoinBigIndex j=mrstrt[iRow];j<mrstrt
919 for (int k=0;k<prob->ncols0_;k++) {
920 CoinBigIndex j=mcstrt[k];
921 for (int i=0;i<hincol[k];i++) {
922 if (hrow[k]==iRow) {
923 break;
924 }
925 k=link[k];
926 }
927 }
928#endif
929 } else {
930 prob->setColumnStatus(iCol,CoinPrePostsolveMatrix::basic);
931 }
932 prob->setRowStatusUsingValue(iRow);
933 } else {
934 //prob->setRowStatus(iRow,CoinPrePostsolveMatrix::basic);
935 prob->setColumnStatusUsingValue(iCol);
936 }
937 }
938 }
939#if 0
940 int nb=0;
941 int kk;
942 for (kk=0;kk<prob->nrows_;kk++)
943 if (prob->rowIsBasic(kk))
944 nb++;
945 for (kk=0;kk<prob->ncols_;kk++)
946 if (prob->columnIsBasic(kk))
947 nb++;
948 assert (nb==prob->nrows_);
949#endif
950 // add new element
951 {
952 CoinBigIndex k = free_list;
953 assert(k >= 0 && k < prob->bulk0_) ;
954 free_list = link[free_list];
955 hrow[k] = iRow;
956 colels[k] = coeff;
957 link[k] = mcstrt[iCol];
958 mcstrt[iCol] = k;
959 }
960 hincol[iCol]++; // right?
961#ifdef CHECK_ONE_ROW
962 {
963 double act=0.0;
964 for (int i=0;i<prob->ncols_;i++) {
965 double solV = sol[i];
966 assert (solV>=clo[i]-ztolzb&&solV<=cup[i]+ztolzb);
967 int j=mcstrt[i];
968 for (int k=0;k<hincol[i];k++) {
969 if (hrow[j]==CHECK_ONE_ROW) {
970 //printf("c %d el %g sol %g old act %g new %g\n",
971 // i,colels[j],solV,act, act+colels[j]*solV);
972 act += colels[j]*solV;
973 }
974 j=link[j];
975 }
976 }
977 assert (act>=rlo[CHECK_ONE_ROW]-ztolzb&&act<=rup[CHECK_ONE_ROW]+ztolzb);
978 printf("rhs now %g %g %g %g\n",rlo[CHECK_ONE_ROW],act,acts[CHECK_ONE_ROW],rup[CHECK_ONE_ROW]);
979 }
980#endif
981
982# if PRESOLVE_DEBUG
983 rdone[iRow] = SLACK_SINGLETON;
984# endif
985 }
986
987# if PRESOLVE_DEBUG
988 presolve_check_sol(prob) ;
989 presolve_check_nbasic(prob) ;
990 std::cout << "Leaving slack_singleton_action::postsolve." << std::endl ;
991# endif
992
993 return ;
994}
995