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