1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the class library */
4 /* SoPlex --- the Sequential object-oriented simPlex. */
5 /* */
6 /* Copyright (C) 1996-2022 Konrad-Zuse-Zentrum */
7 /* fuer Informationstechnik Berlin */
8 /* */
9 /* SoPlex is distributed under the terms of the ZIB Academic Licence. */
10 /* */
11 /* You should have received a copy of the ZIB Academic License */
12 /* along with SoPlex; see the file COPYING. If not email to soplex@zib.de. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15
16 #include <assert.h>
17 #include <iostream>
18
19 #include "soplex/spxdefines.h"
20 #include "soplex/rational.h"
21 #include "soplex/spxsolver.h"
22 #include "soplex/spxpricer.h"
23 #include "soplex/spxratiotester.h"
24 #include "soplex/spxdefaultrt.h"
25 #include "soplex/spxstarter.h"
26 #include "soplex/spxout.h"
27
28 #define MAXCYCLES 400
29 #define MAXSTALLS 10000
30 #define MAXSTALLRECOVERS 10
31 #define MAXREFACPIVOTS 10
32
33 namespace soplex
34 {
35
36 /**@todo check separately for ENTER and LEAVE algorithm */
37 template <class R>
38 bool SPxSolverBase<R>::precisionReached(R& newpricertol) const
39 {
40 R maxViolRedCost;
41 R sumViolRedCost;
42 R maxViolBounds;
43 R sumViolBounds;
44 R maxViolConst;
45 R sumViolConst;
46
47 qualRedCostViolation(maxViolRedCost, sumViolRedCost);
48 qualBoundViolation(maxViolBounds, sumViolBounds);
49 qualConstraintViolation(maxViolConst, sumViolConst);
50
51 // is the solution good enough ?
52 bool reached = maxViolRedCost < opttol() && maxViolBounds < feastol() && maxViolConst < feastol();
53
54 if(!reached)
55 {
56 newpricertol = thepricer->epsilon() / 10.0;
57
58 MSG_INFO3((*this->spxout), (*this->spxout) << "Precision not reached: Pricer tolerance = "
59 << thepricer->epsilon()
60 << " new tolerance = " << newpricertol
61 << std::endl
62 << " maxViolRedCost= " << maxViolRedCost
63 << " maxViolBounds= " << maxViolBounds
64 << " maxViolConst= " << maxViolConst
65 << std::endl
66 << " sumViolRedCost= " << sumViolRedCost
67 << " sumViolBounds= " << sumViolBounds
68 << " sumViolConst= " << sumViolConst
69 << std::endl;);
70 }
71
72 return reached;
73 }
74
75 template <class R>
76 void SPxSolverBase<R>::calculateProblemRanges()
77 {
78 // only collect absolute values
79 R minobj = R(infinity);
80 R maxobj = 0.0;
81 R minbound = R(infinity);
82 R maxbound = 0.0;
83 R minside = R(infinity);
84 R maxside = 0.0;
85
86 // get min and max absolute values of bounds and objective
87 for(int j = 0; j < this->nCols(); ++j)
88 {
89 R abslow = spxAbs(this->lower(j));
90 R absupp = spxAbs(this->lower(j));
91 R absobj = spxAbs(this->obj(j));
92
93 if(abslow < R(infinity))
94 {
95 minbound = MINIMUM(minbound, abslow);
96 maxbound = MAXIMUM(maxbound, abslow);
97 }
98
99 if(absupp < R(infinity))
100 {
101 minbound = MINIMUM(minbound, absupp);
102 maxbound = MAXIMUM(maxbound, absupp);
103 }
104
105 minobj = MINIMUM(minobj, absobj);
106 maxobj = MAXIMUM(maxobj, absobj);
107 }
108
109 // get min and max absoute values of sides
110 for(int i = 0; i < this->nRows(); ++i)
111 {
112 R abslhs = spxAbs(this->lhs(i));
113 R absrhs = spxAbs(this->rhs(i));
114
115 if(abslhs > R(infinity))
116 {
117 minside = MINIMUM(minside, abslhs);
118 maxside = MAXIMUM(maxside, abslhs);
119 }
120
121 if(absrhs < R(infinity))
122 {
123 minside = MINIMUM(minside, absrhs);
124 maxside = MAXIMUM(maxside, absrhs);
125 }
126 }
127
128 boundrange = maxbound - minbound;
129 siderange = maxside - minside;
130 objrange = maxobj - minobj;
131 }
132
133 template <class R>
134 typename SPxSolverBase<R>::Status SPxSolverBase<R>::solve(volatile bool* interrupt)
135 {
136
137 SPxId enterId;
138 int leaveNum;
139 int loopCount = 0;
140 R minShift = R(infinity);
141 int cycleCount = 0;
142 bool priced = false;
143 R lastDelta = 1;
144
145 /* allow clean up step only once */
146 recomputedVectors = false;
147
148 /* store the last (primal or dual) feasible objective value to recover/abort in case of stalling */
149 R stallRefValue;
150 R stallRefShift;
151 int stallRefIter;
152 int stallNumRecovers;
153
154 if(dim() <= 0 && coDim() <= 0) // no problem loaded
155 {
156 m_status = NO_PROBLEM;
157 throw SPxStatusException("XSOLVE01 No Problem loaded");
158 }
159
160 if(slinSolver() == 0) // linear system solver is required.
161 {
162 m_status = NO_SOLVER;
163 throw SPxStatusException("XSOLVE02 No Solver loaded");
164 }
165
166 if(thepricer == 0) // pricer is required.
167 {
168 m_status = NO_PRICER;
169 throw SPxStatusException("XSOLVE03 No Pricer loaded");
170 }
171
172 if(theratiotester == 0) // ratiotester is required.
173 {
174 m_status = NO_RATIOTESTER;
175 throw SPxStatusException("XSOLVE04 No RatioTester loaded");
176 }
177
178 theTime->reset();
179 theTime->start();
180
181 m_numCycle = 0;
182 this->iterCount = 0;
183 this->lastIterCount = 0;
184 this->iterDegenCheck = 0;
185
186 if(!isInitialized())
187 {
188 /*
189 if(SPxBasisBase<R>::status() <= NO_PROBLEM)
190 SPxBasisBase<R>::load(this);
191 */
192 /**@todo != REGULAR is not enough. Also OPTIMAL/DUAL/PRIMAL should
193 * be tested and acted accordingly.
194 */
195 if(thestarter != 0 && status() != REGULAR
196 && this->theLP->status() == NO_PROBLEM) // no basis and no starter.
197 thestarter->generate(*this); // generate start basis.
198
199 init();
200
201 // Inna/Tobi: init might fail, if the basis is singular
202 if(!isInitialized())
203 {
204 assert(SPxBasisBase<R>::status() == SPxBasisBase<R>::SINGULAR);
205 m_status = UNKNOWN;
206 return status();
207 }
208 }
209
210 //setType(type());
211
212 if(!this->matrixIsSetup)
213 SPxBasisBase<R>::load(this);
214
215 //factorized = false;
216
217 assert(thepricer->solver() == this);
218 assert(theratiotester->solver() == this);
219
220 // maybe this should be done in init() ?
221 thepricer->setType(type());
222 theratiotester->setType(type());
223
224 MSG_INFO3((*this->spxout),
225 (*this->spxout) << "starting value = " << value() << std::endl
226 << "starting shift = " << shift() << std::endl;
227 )
228
229 if(SPxBasisBase<R>::status() == SPxBasisBase<R>::OPTIMAL)
230 setBasisStatus(SPxBasisBase<R>::REGULAR);
231
232 m_status = RUNNING;
233 bool stop = terminate();
234 leaveCount = 0;
235 enterCount = 0;
236 primalCount = 0;
237 polishCount = 0;
238 boundflips = 0;
239 totalboundflips = 0;
240 enterCycles = 0;
241 leaveCycles = 0;
242 primalDegenSum = 0;
243 dualDegenSum = 0;
244
245 multSparseCalls = 0;
246 multFullCalls = 0;
247 multColwiseCalls = 0;
248 multUnsetupCalls = 0;
249
250 stallNumRecovers = 0;
251
252 /* if we run into a singular basis, we will retry from regulardesc with tighter tolerance in the ratio test */
253 typename SPxSolverBase<R>::Type tightenedtype = type();
254 bool tightened = false;
255
256 while(!stop)
257 {
258 const typename SPxBasisBase<R>::Desc regulardesc = this->desc();
259
260 // we need to reset these pointers to avoid unnecessary/wrong solves in leave() or enter()
261 solveVector2 = 0;
262 solveVector3 = 0;
263 coSolveVector2 = 0;
264 coSolveVector3 = 0;
265
266 updateViols.clear();
267 updateViolsCo.clear();
268
269 try
270 {
271
272 if(type() == ENTER)
273 {
274 forceRecompNonbasicValue();
275
276 int enterCycleCount = 0;
277 int enterFacPivotCount = 0;
278
279 instableEnterVal = 0;
280 instableEnterId = SPxId();
281 instableEnter = false;
282
283 stallRefIter = this->iteration() - 1;
284 stallRefShift = shift();
285 stallRefValue = value();
286
287 /* in the entering algorithm, entertol() should be maintained by the ratio test and leavetol() should be
288 * reached by the pricer
289 */
290 R maxpricertol = leavetol();
291 R minpricertol = 0.01 * maxpricertol;
292
293 thepricer->setEpsilon(maxpricertol);
294 priced = false;
295
296 // to avoid shifts we restrict tolerances in the ratio test
297 if(loopCount > 0)
298 {
299 lastDelta = (lastDelta < entertol()) ? lastDelta : entertol();
300 lastDelta *= 0.01;
301 theratiotester->setDelta(lastDelta);
302 assert(theratiotester->getDelta() > 0);
303 MSG_DEBUG(std::cout << "decreased delta for ratiotest to: " << theratiotester->getDelta() <<
304 std::endl;)
305 }
306 else
307 {
308 lastDelta = 1;
309 theratiotester->setDelta(entertol());
310 }
311
312 printDisplayLine(true);
313
314 do
315 {
316 printDisplayLine();
317
318 enterId = thepricer->selectEnter();
319
320 if(!enterId.isValid() && instableEnterId.isValid() && this->lastUpdate() == 0)
321 {
322 /* no entering variable was found, but because of valid instableEnterId we know
323 that this is due to the scaling of the test values. Thus, we use
324 instableEnterId and SPxFastRT<R>::selectEnter shall accept even an instable
325 leaving variable. */
326 MSG_INFO3((*this->spxout), (*this->spxout) << " --- trying instable enter iteration" << std::endl;)
327
328 enterId = instableEnterId;
329 instableEnter = true;
330 // we also need to reset the test() or coTest() value for getEnterVals()
331 assert(instableEnterVal < 0);
332
333 if(enterId.isSPxColId())
334 {
335 int idx = this->number(SPxColId(enterId));
336
337 if(rep() == COLUMN)
338 {
339 theTest[idx] = instableEnterVal;
340
341 if(sparsePricingEnterCo && isInfeasibleCo[idx] == SPxPricer<R>::NOT_VIOLATED)
342 {
343 infeasibilitiesCo.addIdx(idx);
344 isInfeasibleCo[idx] = SPxPricer<R>::VIOLATED;
345 }
346
347 if(hyperPricingEnter)
348 updateViolsCo.addIdx(idx);
349 }
350 else
351 {
352 theCoTest[idx] = instableEnterVal;
353
354 if(sparsePricingEnter && isInfeasible[idx] == SPxPricer<R>::NOT_VIOLATED)
355 {
356 infeasibilities.addIdx(idx);
357 isInfeasible[idx] = SPxPricer<R>::VIOLATED;
358 }
359
360 if(hyperPricingEnter)
361 updateViols.addIdx(idx);
362 }
363 }
364 else
365 {
366 int idx = this->number(SPxRowId(enterId));
367
368 if(rep() == COLUMN)
369 {
370 theCoTest[idx] = instableEnterVal;
371
372 if(sparsePricingEnter && isInfeasible[idx] == SPxPricer<R>::NOT_VIOLATED)
373 {
374 infeasibilities.addIdx(idx);
375 isInfeasible[idx] = SPxPricer<R>::VIOLATED;
376 }
377
378 if(hyperPricingEnter)
379 updateViols.addIdx(idx);
380 }
381 else
382 {
383 theTest[idx] = instableEnterVal;
384
385 if(sparsePricingEnterCo && isInfeasibleCo[idx] == SPxPricer<R>::NOT_VIOLATED)
386 {
387 infeasibilitiesCo.addIdx(idx);
388 isInfeasibleCo[idx] = SPxPricer<R>::VIOLATED;
389 }
390
391 if(hyperPricingEnter)
392 updateViolsCo.addIdx(idx);
393 }
394 }
395 }
396 else
397 {
398 instableEnter = false;
399 }
400
401 if(!enterId.isValid())
402 {
403 // we are not infeasible and have no shift
404 if(shift() <= epsilon()
405 && (SPxBasisBase<R>::status() == SPxBasisBase<R>::REGULAR
406 || SPxBasisBase<R>::status() == SPxBasisBase<R>::DUAL
407 || SPxBasisBase<R>::status() == SPxBasisBase<R>::PRIMAL))
408 {
409 R newpricertol = minpricertol;
410
411 // refactorize to eliminate accumulated errors from LU updates
412 if(this->lastUpdate() > 0)
413 factorize();
414
415 // recompute Fvec, Pvec and CoPvec to get a more precise solution and obj value
416 computeFrhs();
417 SPxBasisBase<R>::solve(*theFvec, *theFrhs);
418
419 computeEnterCoPrhs();
420 SPxBasisBase<R>::coSolve(*theCoPvec, *theCoPrhs);
421 computePvec();
422
423 forceRecompNonbasicValue();
424
425 MSG_INFO2((*this->spxout), (*this->spxout) << " --- checking feasibility and optimality\n")
426 computeCoTest();
427 computeTest();
428
429 // is the solution good enough ?
430 // max three times reduced
431 if((thepricer->epsilon() > minpricertol) && !precisionReached(newpricertol))
432 {
433 // no!
434 // we reduce the pricer tolerance. Note that if the pricer does not find a candiate
435 // with the reduced tolerance, we quit, regardless of the violations.
436 if(newpricertol < minpricertol)
437 newpricertol = minpricertol;
438
439 thepricer->setEpsilon(newpricertol);
440
441 MSG_INFO2((*this->spxout), (*this->spxout) << " --- setting pricer tolerance to "
442 << thepricer->epsilon()
443 << std::endl;)
444 }
445 }
446
447 // if the factorization is not fresh, we better refactorize and call the pricer again; however, this can
448 // create cycling, so it is performed only a limited number of times per ENTER round
449 if(this->lastUpdate() > 0 && enterFacPivotCount < MAXREFACPIVOTS)
450 {
451 MSG_INFO3((*this->spxout), (*this->spxout) << " --- solve(enter) triggers refactorization" <<
452 std::endl;)
453
454 factorize();
455
456 // if the factorization was found out to be singular, we have to quit
457 if(SPxBasisBase<R>::status() < SPxBasisBase<R>::REGULAR)
458 {
459 MSG_INFO1((*this->spxout), (*this->spxout) << "Something wrong with factorization, Basis status: "
460 << static_cast<int>(SPxBasisBase<R>::status()) << std::endl;)
461 stop = true;
462 break;
463 }
464
465 // call pricer again
466 enterId = thepricer->selectEnter();
467
468 // count how often the pricer has found something only after refactorizing
469 if(enterId.isValid())
470 enterFacPivotCount++;
471 }
472
473 if(!enterId.isValid())
474 {
475 priced = true;
476 break;
477 }
478 }
479
480 /* check if we have iterations left */
481 if(maxIters >= 0 && iterations() >= maxIters)
482 {
483 MSG_INFO2((*this->spxout), (*this->spxout) << " --- maximum number of iterations (" << maxIters
484 << ") reached" << std::endl;)
485 m_status = ABORT_ITER;
486 stop = true;
487 break;
488 }
489
490 if(interrupt != NULL && *interrupt)
491 {
492 MSG_INFO2((*this->spxout), (*this->spxout) << " --- aborted due to interrupt signal" << std::endl;)
493 m_status = ABORT_TIME;
494 stop = true;
495 break;
496 }
497
498 enter(enterId);
499 assert((testBounds(), 1));
500 thepricer->entered4(this->lastEntered(), this->lastIndex());
501 stop = terminate();
502 clearUpdateVecs();
503
504 /* if a successful pivot was performed or a nonbasic variable was flipped to its other bound, we reset the
505 * cycle counter
506 */
507 if(this->lastEntered().isValid())
508 enterCycleCount = 0;
509 else if(basis().status() != SPxBasisBase<R>::INFEASIBLE
510 && basis().status() != SPxBasisBase<R>::UNBOUNDED)
511 {
512 enterCycleCount++;
513
514 if(enterCycleCount > MAXCYCLES)
515 {
516 MSG_INFO2((*this->spxout), (*this->spxout) << " --- abort solving due to cycling in "
517 << "entering algorithm" << std::endl;);
518 m_status = ABORT_CYCLING;
519 stop = true;
520 }
521 }
522
523 /* only if the basis has really changed, we increase the iterations counter; this is not the case when only
524 * a nonbasic variable was flipped to its other bound
525 */
526 if(this->lastIndex() >= 0)
527 {
528 enterCount++;
529 assert(this->lastEntered().isValid());
530 }
531
532 /* check every MAXSTALLS iterations whether shift and objective value have not changed */
533 if((this->iteration() - stallRefIter) % MAXSTALLS == 0
534 && basis().status() != SPxBasisBase<R>::INFEASIBLE)
535 {
536 if(spxAbs(value() - stallRefValue) <= epsilon() && spxAbs(shift() - stallRefShift) <= epsilon())
537 {
538 if(stallNumRecovers < MAXSTALLRECOVERS)
539 {
540 /* try to recover by unshifting/switching algorithm up to MAXSTALLRECOVERS times (just a number picked) */
541 MSG_INFO3((*this->spxout), (*this->spxout) <<
542 " --- stalling detected - trying to recover by switching to LEAVING algorithm." << std::endl;)
543
544 ++stallNumRecovers;
545 break;
546 }
547 else
548 {
549 /* giving up */
550 MSG_INFO2((*this->spxout), (*this->spxout) <<
551 " --- abort solving due to stalling in entering algorithm." << std::endl;);
552
553 m_status = ABORT_CYCLING;
554 stop = true;
555 }
556 }
557 else
558 {
559 /* merely update reference values */
560 stallRefIter = this->iteration() - 1;
561 stallRefShift = shift();
562 stallRefValue = value();
563 }
564 }
565
566 //@ assert(isConsistent());
567 }
568 while(!stop);
569
570 MSG_INFO3((*this->spxout),
571 (*this->spxout) << " --- enter finished. iteration: " << this->iteration()
572 << ", value: " << value()
573 << ", shift: " << shift()
574 << ", epsilon: " << epsilon()
575 << ", feastol: " << feastol()
576 << ", opttol: " << opttol()
577 << std::endl
578 << "ISOLVE56 stop: " << stop
579 << ", basis status: " << static_cast<int>(SPxBasisBase<R>::status()) << " (" << static_cast<int>
580 (SPxBasisBase<R>::status()) << ")"
581 << ", solver status: " << static_cast<int>(m_status) << " (" << static_cast<int>
582 (m_status) << ")" << std::endl;
583 )
584
585 if(!stop)
586 {
587 /**@todo technically it would be ok to finish already when (priced && maxinfeas + shift() <= entertol()) is
588 * satisfied; maybe at least in the case when SoPlex keeps jumping back between ENTER and LEAVE always
589 * shifting (looping), we may relax this condition here;
590 * note also that unShift may increase shift() slightly due to roundoff errors
591 */
592 if(shift() <= epsilon())
593 {
594 // factorize();
595 unShift();
596
597 R maxinfeas = maxInfeas();
598
599 MSG_INFO3((*this->spxout),
600 (*this->spxout) << " --- maxInfeas: " << maxinfeas
601 << ", shift: " << shift()
602 << ", entertol: " << entertol() << std::endl;
603 )
604
605 if(priced && maxinfeas + shift() <= entertol())
606 {
607 setBasisStatus(SPxBasisBase<R>::OPTIMAL);
608 m_status = OPTIMAL;
609 break;
610 }
611 else if(loopCount > 2)
612 {
613 // calculate problem ranges if not done already
614 if(boundrange == 0.0 || siderange == 0.0 || objrange == 0.0)
615 calculateProblemRanges();
616
617 if(MAXIMUM(MAXIMUM(boundrange, siderange), objrange) >= 1e9)
618 {
619 SPxOut::setScientific(spxout->getCurrentStream(), 0);
620 MSG_INFO1((*this->spxout), (*this->spxout) <<
621 " --- termination despite violations (numerical difficulties,"
622 << " bound range = " << boundrange
623 << ", side range = " << siderange
624 << ", obj range = " << objrange
625 << ")" << std::endl;)
626 setBasisStatus(SPxBasisBase<R>::OPTIMAL);
627 m_status = OPTIMAL;
628 break;
629 }
630 else
631 {
632 m_status = ABORT_CYCLING;
633 throw SPxStatusException("XSOLVE14 Abort solving due to looping");
634 }
635 }
636
637 loopCount++;
638 }
639
640 setType(LEAVE);
641 init();
642 thepricer->setType(type());
643 theratiotester->setType(type());
644 }
645 }
646 else
647 {
648 assert(type() == LEAVE);
649
650 forceRecompNonbasicValue();
651
652 int leaveCycleCount = 0;
653 int leaveFacPivotCount = 0;
654
655 instableLeaveNum = -1;
656 instableLeave = false;
657 instableLeaveVal = 0;
658
659 stallRefIter = this->iteration() - 1;
660 stallRefShift = shift();
661 stallRefValue = value();
662
663 /* in the leaving algorithm, leavetol() should be maintained by the ratio test and entertol() should be reached
664 * by the pricer
665 */
666 R maxpricertol = entertol();
667 R minpricertol = 0.01 * maxpricertol;
668
669 thepricer->setEpsilon(maxpricertol);
670 priced = false;
671
672 // to avoid shifts we restrict tolerances in the ratio test
673 if(loopCount > 0)
674 {
675 lastDelta = (lastDelta < leavetol()) ? lastDelta : leavetol();
676 lastDelta *= 0.01;
677 theratiotester->setDelta(lastDelta);
678 assert(theratiotester->getDelta() > 0);
679 MSG_DEBUG(std::cout << "decreased delta for ratiotest to: " << theratiotester->getDelta() <<
680 std::endl;)
681 }
682 else
683 {
684 lastDelta = 1;
685 theratiotester->setDelta(leavetol());
686 }
687
688 printDisplayLine(true);
689
690 do
691 {
692 printDisplayLine();
693
694 leaveNum = thepricer->selectLeave();
695
696 if(leaveNum < 0 && instableLeaveNum >= 0 && this->lastUpdate() == 0)
697 {
698 /* no leaving variable was found, but because of instableLeaveNum >= 0 we know
699 that this is due to the scaling of theCoTest[...]. Thus, we use
700 instableLeaveNum and SPxFastRT<R>::selectEnter shall accept even an instable
701 entering variable. */
702 MSG_INFO3((*this->spxout),
703 (*this->spxout) << " --- trying instable leave iteration" << std::endl;
704 )
705
706 leaveNum = instableLeaveNum;
707 instableLeave = true;
708 // we also need to reset the fTest() value for getLeaveVals()
709 assert(instableLeaveVal < 0);
710 theCoTest[instableLeaveNum] = instableLeaveVal;
711
712 if(sparsePricingLeave)
713 {
714 if(isInfeasible[instableLeaveNum] == SPxPricer<R>::NOT_VIOLATED)
715 {
716 infeasibilities.addIdx(instableLeaveNum);
717 isInfeasible[instableLeaveNum] = SPxPricer<R>::VIOLATED;
718 }
719
720 if(hyperPricingLeave)
721 updateViols.addIdx(instableLeaveNum);
722 }
723 }
724 else
725 {
726 instableLeave = false;
727 }
728
729 if(leaveNum < 0)
730 {
731 // we are not infeasible and have no shift
732 if(shift() <= epsilon()
733 && (SPxBasisBase<R>::status() == SPxBasisBase<R>::REGULAR
734 || SPxBasisBase<R>::status() == SPxBasisBase<R>::DUAL
735 || SPxBasisBase<R>::status() == SPxBasisBase<R>::PRIMAL))
736 {
737 R newpricertol = minpricertol;
738
739 // refactorize to eliminate accumulated errors from LU updates
740 if(this->lastUpdate() > 0)
741 factorize();
742
743 // recompute Fvec, Pvec and CoPvec to get a more precise solution and obj value
744 computeFrhs();
745 SPxBasisBase<R>::solve(*theFvec, *theFrhs);
746
747 computeLeaveCoPrhs();
748 SPxBasisBase<R>::coSolve(*theCoPvec, *theCoPrhs);
749 computePvec();
750
751 forceRecompNonbasicValue();
752
753 MSG_INFO2((*this->spxout), (*this->spxout) << " --- checking feasibility and optimality\n")
754 computeFtest();
755
756 // is the solution good enough ?
757 // max three times reduced
758 if((thepricer->epsilon() > minpricertol) && !precisionReached(newpricertol))
759 {
760 // no
761 // we reduce the pricer tolerance. Note that if the pricer does not find a candiate
762 // with the reduced pricer tolerance, we quit, regardless of the violations.
763 if(newpricertol < minpricertol)
764 newpricertol = minpricertol;
765
766 thepricer->setEpsilon(newpricertol);
767
768 MSG_INFO2((*this->spxout), (*this->spxout) << " --- setting pricer tolerance to "
769 << thepricer->epsilon()
770 << std::endl;);
771 }
772 }
773
774 // if the factorization is not fresh, we better refactorize and call the pricer again; however, this can
775 // create cycling, so it is performed only a limited number of times per LEAVE round
776 if(this->lastUpdate() > 0 && leaveFacPivotCount < MAXREFACPIVOTS)
777 {
778 MSG_INFO3((*this->spxout), (*this->spxout) << " --- solve(leave) triggers refactorization" <<
779 std::endl;)
780
781 factorize();
782
783 // Inna/Tobi: if the factorization was found out to be singular, we have to quit
784 if(SPxBasisBase<R>::status() < SPxBasisBase<R>::REGULAR)
785 {
786 MSG_INFO1((*this->spxout), (*this->spxout) << "Something wrong with factorization, Basis status: "
787 << static_cast<int>(SPxBasisBase<R>::status()) << std::endl;)
788 stop = true;
789 break;
790 }
791
792 // call pricer again
793 leaveNum = thepricer->selectLeave();
794
795 // count how often the pricer has found something only after refactorizing
796 if(leaveNum >= 0)
797 leaveFacPivotCount++;
798 }
799
800 if(leaveNum < 0)
801 {
802 priced = true;
803 break;
804 }
805 }
806
807 /* check if we have iterations left */
808 if(maxIters >= 0 && iterations() >= maxIters)
809 {
810 MSG_INFO2((*this->spxout), (*this->spxout) << " --- maximum number of iterations (" << maxIters
811 << ") reached" << std::endl;)
812 m_status = ABORT_ITER;
813 stop = true;
814 break;
815 }
816
817 if(interrupt != NULL && *interrupt)
818 {
819 MSG_INFO2((*this->spxout), (*this->spxout) << " --- aborted due to interrupt signal" << std::endl;)
820 m_status = ABORT_TIME;
821 stop = true;
822 break;
823 }
824
825 leave(leaveNum);
826 assert((testBounds(), 1));
827 thepricer->left4(this->lastIndex(), this->lastLeft());
828 stop = terminate();
829 clearUpdateVecs();
830
831 /* if a successful pivot was performed or a nonbasic variable was flipped to its other bound, we reset the
832 * cycle counter
833 */
834 if(this->lastIndex() >= 0)
835 leaveCycleCount = 0;
836 else if(basis().status() != SPxBasisBase<R>::INFEASIBLE
837 && basis().status() != SPxBasisBase<R>::UNBOUNDED)
838 {
839 leaveCycleCount++;
840
841 if(leaveCycleCount > MAXCYCLES)
842 {
843 MSG_INFO2((*this->spxout), (*this->spxout) <<
844 " --- abort solving due to cycling in leaving algorithm" << std::endl;);
845 m_status = ABORT_CYCLING;
846 stop = true;
847 }
848 }
849
850 /* only if the basis has really changed, we increase the iterations counter; this is not the case when only
851 * a nonbasic variable was flipped to its other bound
852 */
853 if(this->lastEntered().isValid())
854 {
855 leaveCount++;
856 assert(this->lastIndex() >= 0);
857 }
858
859 /* check every MAXSTALLS iterations whether shift and objective value have not changed */
860 if((this->iteration() - stallRefIter) % MAXSTALLS == 0
861 && basis().status() != SPxBasisBase<R>::INFEASIBLE)
862 {
863 if(spxAbs(value() - stallRefValue) <= epsilon() && spxAbs(shift() - stallRefShift) <= epsilon())
864 {
865 if(stallNumRecovers < MAXSTALLRECOVERS)
866 {
867 /* try to recover by switching algorithm up to MAXSTALLRECOVERS times */
868 MSG_INFO3((*this->spxout), (*this->spxout) <<
869 " --- stalling detected - trying to recover by switching to ENTERING algorithm." << std::endl;)
870
871 ++stallNumRecovers;
872 break;
873 }
874 else
875 {
876 /* giving up */
877 MSG_INFO2((*this->spxout), (*this->spxout) <<
878 " --- abort solving due to stalling in leaving algorithm" << std::endl;);
879
880 m_status = ABORT_CYCLING;
881 stop = true;
882 }
883 }
884 else
885 {
886 /* merely update reference values */
887 stallRefIter = this->iteration() - 1;
888 stallRefShift = shift();
889 stallRefValue = value();
890 }
891 }
892
893 //@ assert(isConsistent());
894 }
895 while(!stop);
896
897 MSG_INFO3((*this->spxout),
898 (*this->spxout) << " --- leave finished. iteration: " << this->iteration()
899 << ", value: " << value()
900 << ", shift: " << shift()
901 << ", epsilon: " << epsilon()
902 << ", feastol: " << feastol()
903 << ", opttol: " << opttol()
904 << std::endl
905 << "ISOLVE57 stop: " << stop
906 << ", basis status: " << static_cast<int>(SPxBasisBase<R>::status()) << " (" << static_cast<int>
907 (SPxBasisBase<R>::status()) << ")"
908 << ", solver status: " << static_cast<int>(m_status) << " (" << static_cast<int>
909 (m_status) << ")" << std::endl;
910 )
911
912 if(!stop)
913 {
914 if(shift() < minShift)
915 {
916 minShift = shift();
917 cycleCount = 0;
918 }
919 else
920 {
921 cycleCount++;
922
923 if(cycleCount > MAXCYCLES)
924 {
925 m_status = ABORT_CYCLING;
926 throw SPxStatusException("XSOLVE13 Abort solving due to cycling");
927 }
928
929 MSG_INFO3((*this->spxout),
930 (*this->spxout) << " --- maxInfeas: " << maxInfeas()
931 << ", shift: " << shift()
932 << ", leavetol: " << leavetol()
933 << ", cycle count: " << cycleCount << std::endl;
934 )
935 }
936
937 /**@todo technically it would be ok to finish already when (priced && maxinfeas + shift() <= entertol()) is
938 * satisfied; maybe at least in the case when SoPlex keeps jumping back between ENTER and LEAVE always
939 * shifting (looping), we may relax this condition here;
940 * note also that unShift may increase shift() slightly due to roundoff errors
941 */
942 if(shift() <= epsilon())
943 {
944 cycleCount = 0;
945 // factorize();
946 unShift();
947
948 R maxinfeas = maxInfeas();
949
950 MSG_INFO3((*this->spxout),
951 (*this->spxout) << " --- maxInfeas: " << maxinfeas
952 << ", shift: " << shift()
953 << ", leavetol: " << leavetol() << std::endl;
954 )
955
956 // We stop if we are indeed optimal, or if we have already been
957 // two times at this place. In this case it seems futile to
958 // continue.
959 if(priced && maxinfeas + shift() <= leavetol())
960 {
961 setBasisStatus(SPxBasisBase<R>::OPTIMAL);
962 m_status = OPTIMAL;
963 break;
964 }
965 else if(loopCount > 2)
966 {
967 // calculate problem ranges if not done already
968 if(boundrange == 0.0 || siderange == 0.0 || objrange == 0.0)
969 calculateProblemRanges();
970
971 if(MAXIMUM(MAXIMUM(boundrange, siderange), objrange) >= 1e9)
972 {
973 SPxOut::setScientific(spxout->getCurrentStream(), 0);
974 MSG_INFO1((*this->spxout), (*this->spxout) <<
975 " --- termination despite violations (numerical difficulties,"
976 << " bound range = " << boundrange
977 << ", side range = " << siderange
978 << ", obj range = " << objrange
979 << ")" << std::endl;)
980 setBasisStatus(SPxBasisBase<R>::OPTIMAL);
981 m_status = OPTIMAL;
982 break;
983 }
984 else
985 {
986 m_status = ABORT_CYCLING;
987 throw SPxStatusException("XSOLVE14 Abort solving due to looping");
988 }
989 }
990
991 loopCount++;
992 }
993
994 setType(ENTER);
995 init();
996 thepricer->setType(type());
997 theratiotester->setType(type());
998 }
999 }
1000
1001 assert(m_status != SINGULAR);
1002
1003 }
1004 catch(const SPxException& E)
1005 {
1006 // if we stopped due to a singular basis, we reload the original basis and try again with tighter
1007 // tolerance (only once)
1008 if(m_status == SINGULAR && !tightened)
1009 {
1010 tightenedtype = type();
1011
1012 if(tightenedtype == ENTER)
1013 {
1014 m_entertol = 0.01 * m_entertol;
1015
1016 MSG_INFO2((*this->spxout), (*this->spxout) <<
1017 " --- basis singular: reloading basis and solving with tighter ratio test tolerance " << m_entertol
1018 << std::endl;)
1019 }
1020 else
1021 {
1022 m_leavetol = 0.01 * m_leavetol;
1023
1024 MSG_INFO2((*this->spxout), (*this->spxout) <<
1025 " --- basis singular: reloading basis and solving with tighter ratio test tolerance " << m_leavetol
1026 << std::endl;)
1027 }
1028
1029 // load original basis
1030 int niters = iterations();
1031 loadBasis(regulardesc);
1032
1033 // remember iteration count
1034 this->iterCount = niters;
1035
1036 // try initializing basis (might fail if starting basis was already singular)
1037 try
1038 {
1039 init();
1040 theratiotester->setType(type());
1041 }
1042 catch(const SPxException& Ex)
1043 {
1044 MSG_INFO2((*this->spxout), (*this->spxout) <<
1045 " --- reloaded basis singular, resetting original tolerances" << std::endl;)
1046
1047 if(tightenedtype == ENTER)
1048 m_entertol = 100.0 * m_entertol;
1049 else
1050 m_leavetol = 100.0 * m_leavetol;
1051
1052 theratiotester->setType(type());
1053
1054 throw Ex;
1055 }
1056
1057 // reset status and counters
1058 m_status = RUNNING;
1059 m_numCycle = 0;
1060 leaveCount = 0;
1061 enterCount = 0;
1062 stallNumRecovers = 0;
1063
1064 // continue
1065 stop = false;
1066 tightened = true;
1067 }
1068 // reset tolerance to its original value and pass on the exception
1069 else if(tightened)
1070 {
1071 if(tightenedtype == ENTER)
1072 m_entertol = 100.0 * m_entertol;
1073 else
1074 m_leavetol = 100.0 * m_leavetol;
1075
1076 theratiotester->setType(type());
1077
1078 throw E;
1079 }
1080 // pass on the exception
1081 else
1082 throw E;
1083 }
1084 }
1085
1086 // reset tolerance to its original value
1087 if(tightened)
1088 {
1089 if(tightenedtype == ENTER)
1090 m_entertol = 100.0 * m_entertol;
1091 else
1092 m_leavetol = 100.0 * m_leavetol;
1093
1094 theratiotester->setType(type());
1095 }
1096
1097 theTime->stop();
1098 theCumulativeTime += time();
1099
1100 if(m_status == RUNNING)
1101 {
1102 m_status = ERROR;
1103 throw SPxStatusException("XSOLVE05 Status is still RUNNING when it shouldn't be");
1104 }
1105
1106 MSG_INFO3((*this->spxout),
1107 (*this->spxout) << "Finished solving (status=" << static_cast<int>(status())
1108 << ", iters=" << this->iterCount
1109 << ", leave=" << leaveCount
1110 << ", enter=" << enterCount
1111 << ", flips=" << totalboundflips;
1112
1113 if(status() == OPTIMAL)
1114 (*this->spxout) << ", objValue=" << value();
1115 (*this->spxout) << ")" << std::endl;
1116 )
1117
1118 #ifdef ENABLE_ADDITIONAL_CHECKS
1119
1120 /* check if solution is really feasible */
1121 if(status() == OPTIMAL)
1122 {
1123 int c;
1124 R val;
1125 VectorBase<R> sol(this->nCols());
1126
1127 getPrimalSol(sol);
1128
1129 for(int row = 0; row < this->nRows(); ++row)
1130 {
1131 const SVectorBase<R>& rowvec = this->rowVector(row);
1132 val = 0.0;
1133
1134 for(c = 0; c < rowvec.size(); ++c)
1135 val += rowvec.value(c) * sol[rowvec.index(c)];
1136
1137 if(LT(val, this->lhs(row), feastol()) ||
1138 GT(val, this->rhs(row), feastol()))
1139 {
1140 // Minor rhs violations happen frequently, so print these
1141 // warnings only with verbose level INFO2 and higher.
1142 MSG_INFO2((*this->spxout), (*this->spxout) << "WSOLVE88 Warning! Constraint " << row
1143 << " is violated by solution" << std::endl
1144 << " lhs:" << this->lhs(row)
1145 << " <= val:" << val
1146 << " <= rhs:" << this->rhs(row) << std::endl;)
1147
1148 if(type() == LEAVE && isRowBasic(row))
1149 {
1150 // find basis variable
1151 for(c = 0; c < this->nRows(); ++c)
1152 if(basis().baseId(c).isSPxRowId()
1153 && (this->number(basis().baseId(c)) == row))
1154 break;
1155
1156 assert(c < this->nRows());
1157
1158 MSG_WARNING((*this->spxout), (*this->spxout) << "WSOLVE90 basis idx:" << c
1159 << " fVec:" << fVec()[c]
1160 << " fRhs:" << fRhs()[c]
1161 << " fTest:" << fTest()[c] << std::endl;)
1162 }
1163 }
1164 }
1165
1166 for(int col = 0; col < this->nCols(); ++col)
1167 {
1168 if(LT(sol[col], this->lower(col), feastol()) ||
1169 GT(sol[col], this->upper(col), feastol()))
1170 {
1171 // Minor bound violations happen frequently, so print these
1172 // warnings only with verbose level INFO2 and higher.
1173 MSG_INFO2((*this->spxout), (*this->spxout) << "WSOLVE91 Warning! Bound for column " << col
1174 << " is violated by solution" << std::endl
1175 << " lower:" << this->lower(col)
1176 << " <= val:" << sol[col]
1177 << " <= upper:" << this->upper(col) << std::endl;)
1178
1179 if(type() == LEAVE && isColBasic(col))
1180 {
1181 for(c = 0; c < this->nRows() ; ++c)
1182 if(basis().baseId(c).isSPxColId()
1183 && (this->number(basis().baseId(c)) == col))
1184 break;
1185
1186 assert(c < this->nRows());
1187 MSG_WARNING((*this->spxout), (*this->spxout) << "WSOLVE92 basis idx:" << c
1188 << " fVec:" << fVec()[c]
1189 << " fRhs:" << fRhs()[c]
1190 << " fTest:" << fTest()[c] << std::endl;)
1191 }
1192 }
1193 }
1194 }
1195
1196 #endif // ENABLE_ADDITIONAL_CHECKS
1197
1198
1199 primalCount = (rep() == SPxSolverBase<R>::COLUMN)
1200 ? enterCount
1201 : leaveCount;
1202
1203 printDisplayLine(true);
1204 performSolutionPolishing();
1205
1206 return status();
1207 }
1208
1209 template <class R>
1210 void SPxSolverBase<R>::performSolutionPolishing()
1211 {
1212 // catch rare case that the iteration limit is exactly reached at optimality
1213 bool stop = (maxIters >= 0 && iterations() >= maxIters && !isTimeLimitReached());
1214
1215 // only polish an already optimal basis
|
(1) Event cond_false: |
Condition "stop", taking false branch. |
|
(2) Event cond_false: |
Condition "this->polishObj == soplex::SPxSolverBase<double>::POLISH_OFF", taking false branch. |
|
(3) Event cond_false: |
Condition "this->status() != soplex::SPxSolverBase<double>::OPTIMAL", taking false branch. |
1216 if(stop || polishObj == POLISH_OFF || status() != OPTIMAL)
|
(4) Event if_end: |
End of if statement. |
1217 return;
1218
1219 int nSuccessfulPivots;
1220 const typename SPxBasisBase<R>::Desc& ds = this->desc();
1221 const typename SPxBasisBase<R>::Desc::Status* rowstatus = ds.rowStatus();
1222 const typename SPxBasisBase<R>::Desc::Status* colstatus = ds.colStatus();
1223 typename SPxBasisBase<R>::Desc::Status stat;
1224 SPxId polishId;
1225 bool success = false;
1226
|
(5) Event cond_true: |
Condition "this->spxout != NULL", taking true branch. |
|
(6) Event cond_true: |
Condition "soplex::SPxOut::INFO2 <= this->spxout->getVerbosity()", taking true branch. |
1227 MSG_INFO2((*this->spxout), (*this->spxout) << " --- perform solution polishing" << std::endl;)
1228
|
(7) Event cond_false: |
Condition "this->rep() == soplex::SPxSolverBase<double>::COLUMN", taking false branch. |
1229 if(rep() == COLUMN)
1230 {
1231 setType(ENTER); // use primal simplex to preserve feasibility
1232 init();
1233 #ifndef NDEBUG
1234 // allow a small relative deviation from the original values
1235 R alloweddeviation = entertol();
1236 R origval = value();
1237 R origshift = shift();
1238 #endif
1239 instableEnter = false;
1240 theratiotester->setType(type());
1241
1242 if(polishObj == POLISH_INTEGRALITY)
1243 {
1244 DIdxSet slackcandidates(this->nRows());
1245 DIdxSet continuousvars(this->nCols());
1246
1247 // collect nonbasic slack variables that could be made basic
1248 for(int i = 0; i < this->nRows(); ++i)
1249 {
1250 // only check nonbasic rows, skip equations
1251 if(rowstatus[i] == SPxBasisBase<R>::Desc::P_ON_LOWER
1252 || rowstatus[i] == SPxBasisBase<R>::Desc::P_ON_UPPER)
1253 {
1254 // only consider rows with zero dual multiplier to preserve optimality
1255 if(EQrel((*theCoPvec)[i], (R) 0))
1256 slackcandidates.addIdx(i);
1257 }
1258 }
1259
1260 // collect continuous variables that could be made basic
1261 if(integerVariables.size() == this->nCols())
1262 {
1263 for(int i = 0; i < this->nCols(); ++i)
1264 {
1265 // skip fixed variables
1266 if(colstatus[i] == SPxBasisBase<R>::Desc::P_ON_LOWER
1267 || colstatus[i] == SPxBasisBase<R>::Desc::P_ON_UPPER)
1268 {
1269 // only consider continuous variables with zero dual multiplier to preserve optimality
1270 if(EQrel(this->maxObj(i) - (*thePvec)[i], (R) 0) && integerVariables[i] == 0)
1271 continuousvars.addIdx(i);
1272 }
1273 }
1274 }
1275
1276 while(!stop)
1277 {
1278 nSuccessfulPivots = 0;
1279
1280 // identify nonbasic slack variables, i.e. rows, that may be moved into the basis
1281 for(int i = slackcandidates.size() - 1; i >= 0 && !stop; --i)
1282 {
1283 polishId = coId(slackcandidates.index(i));
1284 MSG_DEBUG(std::cout << "try pivoting: " << polishId << " stat: " << rowstatus[slackcandidates.index(
1285 i)];)
1286 success = enter(polishId, true);
1287 clearUpdateVecs();
1288 #ifndef NDEBUG
1289 assert(EQrel(value(), origval, alloweddeviation));
1290 assert(LErel(shift(), origshift, alloweddeviation));
1291 #endif
1292
1293 if(success)
1294 {
1295 MSG_DEBUG(std::cout << " -> success!";)
1296 ++nSuccessfulPivots;
1297 slackcandidates.remove(i);
1298
1299 if(maxIters >= 0 && iterations() >= maxIters)
1300 stop = true;
1301 }
1302
1303 MSG_DEBUG(std::cout << std::endl;)
1304
1305 if(isTimeLimitReached())
1306 stop = true;
1307 }
1308
1309 // identify nonbasic variables that may be moved into the basis
1310 for(int i = continuousvars.size() - 1; i >= 0 && !stop; --i)
1311 {
1312 polishId = id(continuousvars.index(i));
1313 MSG_DEBUG(std::cout << "try pivoting: " << polishId << " stat: " << colstatus[continuousvars.index(
1314 i)];)
1315 success = enter(polishId, true);
1316 clearUpdateVecs();
1317
1318 if(success)
1319 {
1320 MSG_DEBUG(std::cout << " -> success!";)
1321 ++nSuccessfulPivots;
1322 continuousvars.remove(i);
1323
1324 if(maxIters >= 0 && iterations() >= maxIters)
1325 stop = true;
1326 }
1327
1328 MSG_DEBUG(std::cout << std::endl;)
1329
1330 if(isTimeLimitReached())
1331 stop = true;
1332 }
1333
1334 // terminate if in the last round no more polishing steps were performed
1335 if(nSuccessfulPivots == 0)
1336 stop = true;
1337
1338 polishCount += nSuccessfulPivots;
1339 }
1340 }
1341 else
1342 {
1343 assert(polishObj == POLISH_FRACTIONALITY);
1344 DIdxSet candidates(dim());
1345
1346 // identify nonbasic variables, i.e. columns, that may be moved into the basis
1347 for(int i = 0; i < this->nCols() && !stop; ++i)
1348 {
1349 if(colstatus[i] == SPxBasisBase<R>::Desc::P_ON_LOWER
1350 || colstatus[i] == SPxBasisBase<R>::Desc::P_ON_UPPER)
1351 {
1352 // only consider variables with zero reduced costs to preserve optimality
1353 if(EQrel(this->maxObj(i) - (*thePvec)[i], (R) 0))
1354 candidates.addIdx(i);
1355 }
1356 }
1357
1358 while(!stop)
1359 {
1360 nSuccessfulPivots = 0;
1361
1362 for(int i = candidates.size() - 1; i >= 0 && !stop; --i)
1363 {
1364 polishId = id(candidates.index(i));
1365 MSG_DEBUG(std::cout << "try pivoting: " << polishId << " stat: " << colstatus[candidates.index(i)];)
1366 success = enter(polishId, true);
1367 clearUpdateVecs();
1368 #ifndef NDEBUG
1369 assert(EQrel(value(), origval, alloweddeviation));
1370 assert(LErel(shift(), origshift, alloweddeviation));
1371 #endif
1372
1373 if(success)
1374 {
1375 MSG_DEBUG(std::cout << " -> success!";)
1376 ++nSuccessfulPivots;
1377 candidates.remove(i);
1378
1379 if(maxIters >= 0 && iterations() >= maxIters)
1380 stop = true;
1381 }
1382
1383 MSG_DEBUG(std::cout << std::endl;)
1384
1385 if(isTimeLimitReached())
1386 stop = true;
1387 }
1388
1389 // terminate if in the last round no more polishing steps were performed
1390 if(nSuccessfulPivots == 0)
1391 stop = true;
1392
1393 polishCount += nSuccessfulPivots;
1394 }
1395 }
1396 }
1397 else
|
(8) Event else_branch: |
Reached else branch. |
1398 {
1399 setType(LEAVE); // use primal simplex to preserve feasibility
1400 init();
1401 #ifndef NDEBUG
1402 // allow a small relative deviation from the original values
1403 R alloweddeviation = leavetol();
1404 R origval = value();
1405 R origshift = shift();
1406 #endif
1407 instableLeave = false;
1408 theratiotester->setType(type());
1409 bool useIntegrality = false;
1410
|
(9) Event cond_true: |
Condition "this->integerVariables.size() == this->nCols()", taking true branch. |
1411 if(integerVariables.size() == this->nCols())
1412 useIntegrality = true;
1413
1414 // in ROW rep: pivot slack out of the basis
|
(10) Event cond_true: |
Condition "this->polishObj == soplex::SPxSolverBase<double>::POLISH_INTEGRALITY", taking true branch. |
1415 if(polishObj == POLISH_INTEGRALITY)
1416 {
1417 DIdxSet basiccandidates(dim());
1418
1419 // collect basic candidates that may be moved out of the basis
|
(12) Event cond_true: |
Condition "i < this->dim()", taking true branch. |
|
(20) Event loop_begin: |
Jumped back to beginning of loop. |
|
(21) Event cond_false: |
Condition "i < this->dim()", taking false branch. |
1420 for(int i = 0; i < dim(); ++i)
1421 {
1422 polishId = this->baseId(i);
1423
|
(13) Event cond_true: |
Condition "polishId.isSPxRowId()", taking true branch. |
1424 if(polishId.isSPxRowId())
|
(14) Event if_fallthrough: |
Falling through to end of if statement. |
1425 stat = ds.rowStatus(this->number(polishId));
1426 else
1427 {
1428 // skip (integer) variables
1429 if(!useIntegrality || integerVariables[this->number(SPxColId(polishId))] == 1)
1430 continue;
1431
1432 stat = ds.colStatus(this->number(polishId));
|
(15) Event if_end: |
End of if statement. |
1433 }
1434
|
(16) Event cond_true: |
Condition "stat == soplex::SPxBasisBase<double>::Desc::P_ON_LOWER", taking true branch. |
1435 if(stat == SPxBasisBase<R>::Desc::P_ON_LOWER || stat == SPxBasisBase<R>::Desc::P_ON_UPPER)
1436 {
|
(17) Event cond_true: |
Condition "soplex::EQrel((*this->theFvec)[i], (soplex::Real)0, soplex::Param::epsilon())", taking true branch. |
1437 if(EQrel((*theFvec)[i], (R) 0))
|
(18) Event no_write_call: |
Although "addIdx" does overwrite "basiccandidates->idx" on some paths, it also contains at least one feasible path which does not overwrite it. |
| Also see events: |
[write_zero_model][var_deref_model] |
1438 basiccandidates.addIdx(i);
1439 }
|
(19) Event loop: |
Jumping back to the beginning of the loop. |
|
(22) Event loop_end: |
Reached end of loop. |
1440 }
1441
|
(23) Event cond_true: |
Condition "!stop", taking true branch. |
1442 while(!stop)
1443 {
1444 nSuccessfulPivots = 0;
1445
|
(24) Event cond_true: |
Condition "i >= 0", taking true branch. |
|
(25) Event cond_true: |
Condition "!stop", taking true branch. |
1446 for(int i = basiccandidates.size() - 1; i >= 0 && !stop; --i)
1447 {
1448
1449 MSG_DEBUG(std::cout << "try pivoting: " << this->baseId(basiccandidates.index(i));)
1450 success = leave(basiccandidates.index(i), true);
1451 clearUpdateVecs();
1452 #ifndef NDEBUG
1453 assert(EQrel(value(), origval, alloweddeviation));
1454 assert(LErel(shift(), origshift, alloweddeviation));
1455 #endif
1456
1457 if(success)
1458 {
1459 MSG_DEBUG(std::cout << " -> success!";)
1460 ++nSuccessfulPivots;
1461 basiccandidates.remove(i);
1462
1463 if(maxIters >= 0 && iterations() >= maxIters)
1464 stop = true;
1465 }
1466
1467 MSG_DEBUG(std::cout << std::endl;)
1468
1469 if(isTimeLimitReached())
1470 stop = true;
1471 }
1472
1473 // terminate if in the last round no more polishing steps were performed
1474 if(nSuccessfulPivots == 0)
1475 stop = true;
1476
1477 polishCount += nSuccessfulPivots;
1478 }
1479 }
1480 else
1481 {
1482 assert(polishObj == POLISH_FRACTIONALITY);
1483 DIdxSet basiccandidates(dim());
1484
1485 // collect basic (integer) variables, that may be moved out of the basis
1486 for(int i = 0; i < dim(); ++i)
1487 {
1488 polishId = this->baseId(i);
1489
1490 if(polishId.isSPxRowId())
1491 continue;
1492 else
1493 {
1494 if(useIntegrality && integerVariables[this->number(SPxColId(polishId))] == 0)
1495 continue;
1496
1497 stat = ds.colStatus(i);
1498 }
1499
1500 if(stat == SPxBasisBase<R>::Desc::P_ON_LOWER || stat == SPxBasisBase<R>::Desc::P_ON_UPPER)
1501 {
1502 if(EQrel((*theFvec)[i], (R) 0))
1503 basiccandidates.addIdx(i);
1504 }
1505 }
1506
1507 while(!stop)
1508 {
1509 nSuccessfulPivots = 0;
1510
1511 for(int i = basiccandidates.size() - 1; i >= 0 && !stop; --i)
1512 {
1513 MSG_DEBUG(std::cout << "try pivoting: " << this->baseId(basiccandidates.index(i));)
1514 success = leave(basiccandidates.index(i), true);
1515 clearUpdateVecs();
1516 #ifndef NDEBUG
1517 assert(EQrel(value(), origval, alloweddeviation));
1518 assert(LErel(shift(), origshift, alloweddeviation));
1519 #endif
1520
1521 if(success)
1522 {
1523 MSG_DEBUG(std::cout << " -> success!";)
1524 ++nSuccessfulPivots;
1525 basiccandidates.remove(i);
1526
1527 if(maxIters >= 0 && iterations() >= maxIters)
1528 stop = true;
1529 }
1530
1531 MSG_DEBUG(std::cout << std::endl;)
1532
1533 if(isTimeLimitReached())
1534 stop = true;
1535 }
1536
1537 // terminate if in the last round no more polishing steps were performed
1538 if(nSuccessfulPivots == 0)
1539 stop = true;
1540
1541 polishCount += nSuccessfulPivots;
1542 }
1543 }
1544 }
1545
1546 MSG_INFO1((*this->spxout),
1547 (*this->spxout) << " --- finished solution polishing (" << polishCount << " pivots)" << std::endl;)
1548
1549 this->setStatus(SPxBasisBase<R>::OPTIMAL);
1550 }
1551
1552
1553 template <class R>
1554 void SPxSolverBase<R>::testVecs()
1555 {
1556
1557 assert(SPxBasisBase<R>::status() > SPxBasisBase<R>::SINGULAR);
1558
1559 VectorBase<R> tmp(dim());
1560
1561 tmp = *theCoPvec;
1562 this->multWithBase(tmp);
1563 tmp -= *theCoPrhs;
1564
1565 if(tmp.length() > leavetol())
1566 {
1567 MSG_INFO3((*this->spxout), (*this->spxout) << "ISOLVE93 " << this->iteration() <<
1568 ":\tcoP error = \t"
1569 << tmp.length() << std::endl;)
1570
1571 tmp.clear();
1572 SPxBasisBase<R>::coSolve(tmp, *theCoPrhs);
1573 this->multWithBase(tmp);
1574 tmp -= *theCoPrhs;
1575 MSG_INFO3((*this->spxout), (*this->spxout) << "ISOLVE94\t\t" << tmp.length() << std::endl;)
1576
1577 tmp.clear();
1578 SPxBasisBase<R>::coSolve(tmp, *theCoPrhs);
1579 tmp -= *theCoPvec;
1580 MSG_INFO3((*this->spxout), (*this->spxout) << "ISOLVE95\t\t" << tmp.length() << std::endl;)
1581 }
1582
1583 tmp = *theFvec;
1584 this->multBaseWith(tmp);
1585 tmp -= *theFrhs;
1586
1587 if(tmp.length() > entertol())
1588 {
1589 MSG_INFO3((*this->spxout), (*this->spxout) << "ISOLVE96 " << this->iteration() <<
1590 ":\t F error = \t"
1591 << tmp.length() << std::endl;)
1592
1593 tmp.clear();
1594 SPxBasisBase<R>::solve(tmp, *theFrhs);
1595 tmp -= *theFvec;
1596 MSG_INFO3((*this->spxout), (*this->spxout) << "ISOLVE97\t\t" << tmp.length() << std::endl;)
1597 }
1598
1599 if(type() == ENTER)
1600 {
1601 for(int i = 0; i < dim(); ++i)
1602 {
1603 if(theCoTest[i] < -leavetol() && isCoBasic(i))
1604 {
1605 /// @todo Error message "this shalt not be": shalt this be an assert (also below)?
1606 MSG_INFO1((*this->spxout), (*this->spxout) << "ESOLVE98 testVecs: theCoTest: this shalt not be!"
1607 << std::endl
1608 << " i=" << i
1609 << ", theCoTest[i]=" << theCoTest[i]
1610 << ", leavetol()=" << leavetol() << std::endl;)
1611 }
1612 }
1613
1614 for(int i = 0; i < coDim(); ++i)
1615 {
1616 if(theTest[i] < -leavetol() && isBasic(i))
1617 {
1618 MSG_INFO1((*this->spxout), (*this->spxout) << "ESOLVE99 testVecs: theTest: this shalt not be!"
1619 << std::endl
1620 << " i=" << i
1621 << ", theTest[i]=" << theTest[i]
1622 << ", leavetol()=" << leavetol() << std::endl;)
1623 }
1624 }
1625 }
1626 }
1627
1628
1629 /// print display line of flying table
1630 template <class R>
1631 void SPxSolverBase<R>::printDisplayLine(const bool force, const bool forceHead)
1632 {
1633 MSG_INFO1((*this->spxout),
1634
1635 if(forceHead || displayLine % (displayFreq * 30) == 0)
1636 {
1637 (*this->spxout)
1638 << "type | time | iters | facts | shift | viol sum | viol num | obj value ";
1639
1640 if(printBasisMetric >= 0)
1641 (*this->spxout) << " | basis metric";
1642
1643 (*this->spxout) << std::endl;
1644 }
1645 if((force || (displayLine % displayFreq == 0)) && !forceHead)
1646 {
1647 (type() == LEAVE)
1648 ? (*this->spxout) << " L |" : (*this->spxout) << " E |";
1649 (*this->spxout) << std::fixed << std::setw(7) << std::setprecision(1) << time() << " |";
1650 (*this->spxout) << std::scientific << std::setprecision(2);
1651 (*this->spxout) << std::setw(8) << this->iteration() << " | "
1652 << std::setw(5) << slinSolver()->getFactorCount() << " | "
1653 << shift() << " | "
1654 << MAXIMUM(0.0, m_pricingViol + m_pricingViolCo) << " | "
1655 << std::setw(8) << MAXIMUM(0, m_numViol) << " | "
1656 << std::setprecision(8) << value();
1657
1658 if(getStartingDecompBasis && rep() == SPxSolverBase<R>::ROW)
1659 (*this->spxout) << " (" << std::fixed << std::setprecision(2) << getDegeneracyLevel(fVec()) << ")";
1660
1661 if(printBasisMetric == 0)
1662 (*this->spxout) << " | " << std::scientific << std::setprecision(2) << getBasisMetric(0);
1663
1664 if(printBasisMetric == 1)
1665 (*this->spxout) << " | " << std::scientific << std::setprecision(2) << getBasisMetric(1);
1666
1667 if(printBasisMetric == 2)
1668 (*this->spxout) << " | " << std::scientific << std::setprecision(2) << getBasisMetric(2);
1669
1670 if(printBasisMetric == 3)
1671 (*this->spxout) << " | " << std::scientific << std::setprecision(2) <<
1672 basis().getEstimatedCondition();
1673
1674 (*this->spxout) << std::endl;
1675 }
1676 displayLine++;
1677 );
1678 }
1679
1680
1681 template <class R>
1682 bool SPxSolverBase<R>::terminate()
1683 {
1684 #ifdef ENABLE_ADDITIONAL_CHECKS
1685
1686 if(SPxBasisBase<R>::status() > SPxBasisBase<R>::SINGULAR)
1687 testVecs();
1688
1689 #endif
1690
1691 int redo = dim();
1692
1693 if(redo < 1000)
1694 redo = 1000;
1695
1696 if(this->iteration() > 10 && this->iteration() % redo == 0)
1697 {
1698 #ifdef ENABLE_ADDITIONAL_CHECKS
1699 VectorBase<R> cr(*theCoPrhs);
1700 VectorBase<R> fr(*theFrhs);
1701 #endif
1702
1703 if(type() == ENTER)
1704 computeEnterCoPrhs();
1705 else
1706 computeLeaveCoPrhs();
1707
1708 computeFrhs();
1709
1710 #ifdef ENABLE_ADDITIONAL_CHECKS
1711 cr -= *theCoPrhs;
1712 fr -= *theFrhs;
1713
1714 if(cr.length() > leavetol())
1715 MSG_WARNING((*this->spxout), (*this->spxout) << "WSOLVE50 unexpected change of coPrhs "
1716 << cr.length() << std::endl;)
1717 if(fr.length() > entertol())
1718 MSG_WARNING((*this->spxout), (*this->spxout) << "WSOLVE51 unexpected change of Frhs "
1719 << fr.length() << std::endl;)
1720 #endif
1721
1722 if(this->updateCount > 1)
1723 {
1724 MSG_INFO3((*this->spxout), (*this->spxout) << " --- terminate triggers refactorization"
1725 << std::endl;)
1726 factorize();
1727 }
1728
1729 SPxBasisBase<R>::coSolve(*theCoPvec, *theCoPrhs);
1730 SPxBasisBase<R>::solve(*theFvec, *theFrhs);
1731
1732 if(pricing() == FULL)
1733 {
1734 computePvec();
1735
1736 if(type() == ENTER)
1737 {
1738 computeCoTest();
1739 computeTest();
1740 }
1741 }
1742
1743 if(shift() > 0.0)
1744 unShift();
1745 }
1746
1747 // check time limit and objective limit only for non-terminal bases
1748 if(SPxBasisBase<R>::status() >= SPxBasisBase<R>::OPTIMAL ||
1749 SPxBasisBase<R>::status() <= SPxBasisBase<R>::SINGULAR)
1750 {
1751 m_status = UNKNOWN;
1752 return true;
1753 }
1754
1755 if(isTimeLimitReached())
1756 {
1757 MSG_INFO2((*this->spxout), (*this->spxout) << " --- timelimit (" << maxTime
1758 << ") reached" << std::endl;)
1759 m_status = ABORT_TIME;
1760 return true;
1761 }
1762
1763 // objLimit is set and we are running DUAL:
1764 // - objLimit is set if objLimit < R(infinity)
1765 // - DUAL is running if rep() * type() > 0 == DUAL (-1 == PRIMAL)
1766 //
1767 // In this case we have given a objective value limit, e.g, through a
1768 // MIP solver, and we want stop solving the LP if we figure out that the
1769 // optimal value of the current LP can not be better then this objective
1770 // limit. More precisely:
1771 // - MINIMIZATION Problem
1772 // We want stop the solving process if
1773 // objLimit <= current objective value of the DUAL LP
1774 // - MAXIMIZATION Problem
1775 // We want stop the solving process if
1776 // objLimit >= current objective value of the DUAL LP
1777 if(objLimit < R(infinity) && type() * rep() > 0)
1778 {
1779 // We have no bound shifts; therefore, we can trust the current
1780 // objective value.
1781 // It might be even possible to use this termination value in case of
1782 // bound violations (shifting) but in this case it is quite difficult
1783 // to determine if we already reached the limit.
1784 if(shift() < epsilon() && noViols(opttol() - shift()))
1785 {
1786 // SPxSense::MINIMIZE == -1, so we have sign = 1 on minimizing
1787 if(int(this->spxSense()) * value() <= int(this->spxSense()) * objLimit)
1788 {
1789 MSG_INFO2((*this->spxout), (*this->spxout) << " --- objective value limit (" << objLimit
1790 << ") reached" << std::endl;)
1791 MSG_DEBUG(
1792 (*this->spxout) << " --- objective value limit reached" << std::endl
1793 << " (value: " << value()
1794 << ", limit: " << objLimit << ")" << std::endl
1795 << " (spxSense: " << int(this->spxSense())
1796 << ", rep: " << int(rep())
1797 << ", type: " << int(type()) << ")" << std::endl;
1798 )
1799
1800 m_status = ABORT_VALUE;
1801 return true;
1802 }
1803 }
1804 }
1805
1806
1807
1808 if(getComputeDegeneracy() && this->iteration() > this->prevIteration())
1809 {
1810 VectorBase<R> degenvec(this->nCols());
1811
1812 if(rep() == ROW)
1813 {
1814 if(type() == ENTER) // dual simplex
1815 dualDegenSum += getDegeneracyLevel(fVec());
1816 else // primal simplex
1817 {
1818 getPrimalSol(degenvec);
1819 primalDegenSum += getDegeneracyLevel(degenvec);
1820 }
1821 }
1822 else
1823 {
1824 assert(rep() == COLUMN);
1825
1826 if(type() == LEAVE) // dual simplex
1827 dualDegenSum += getDegeneracyLevel(pVec());
1828 else
1829 {
1830 getPrimalSol(degenvec);
1831 primalDegenSum += getDegeneracyLevel(degenvec);
1832 }
1833 }
1834 }
1835
1836
1837 // the improved dual simplex requires a starting basis
1838 // if the flag getStartingDecompBasis is set to true the simplex will terminate when a dual basis is found
1839 if(getStartingDecompBasis)
1840 {
1841 R iterationFrac = 0.6;
1842
1843 if(type() == ENTER && SPxBasisBase<R>::status() == SPxBasisBase<R>::DUAL &&
1844 this->iteration() - this->lastDegenCheck() > getDegenCompOffset()/*iteration() % 10 == 0*/)
1845 {
1846 this->iterDegenCheck = this->iterCount;
1847
1848 if(SPxBasisBase<R>::status() >= SPxBasisBase<R>::OPTIMAL)
1849 {
1850 m_status = RUNNING;
1851 return true;
1852 }
1853
1854 R degeneracyLevel = 0;
1855 R degeneracyLB = 0.1;
1856 R degeneracyUB = 0.9;
1857 degeneracyLevel = getDegeneracyLevel(fVec());
1858
1859 if((degeneracyLevel < degeneracyUB && degeneracyLevel > degeneracyLB)
1860 && this->iteration() > this->nRows() * 0.2)
1861 {
1862 m_status = ABORT_DECOMP;
1863 return true;
1864 }
1865
1866 if(degeneracyLevel < degeneracyLB
1867 && this->iteration() > MINIMUM(getDecompIterationLimit(), int(this->nCols()*iterationFrac)))
1868 {
1869 setDecompIterationLimit(0);
1870 setDegenCompOffset(0);
1871 m_status = ABORT_EXDECOMP;
1872 return true;
1873 }
1874 }
1875 else if(type() == LEAVE
1876 && this->iteration() > MINIMUM(getDecompIterationLimit(), int(this->nCols()*iterationFrac)))
1877 {
1878 setDecompIterationLimit(0);
1879 setDegenCompOffset(0);
1880 m_status = ABORT_EXDECOMP;
1881 return true;
1882 }
1883 }
1884
1885 this->lastIterCount = this->iterCount;
1886
1887 return false;
1888 }
1889
1890 template <class R>
1891 typename SPxSolverBase<R>::Status SPxSolverBase<R>::getPrimalSol(VectorBase<R>& p_vector) const
1892 {
1893
1894 if(!isInitialized())
1895 {
1896 /* exit if presolving/simplifier cleared the problem */
1897 if(status() == NO_PROBLEM)
1898 return status();
1899
1900 throw SPxStatusException("XSOLVE06 Not Initialized");
1901 }
1902
1903 if(rep() == ROW)
1904 p_vector = coPvec();
1905 else
1906 {
1907 const typename SPxBasisBase<R>::Desc& ds = this->desc();
1908
1909 for(int i = 0; i < this->nCols(); ++i)
1910 {
1911 switch(ds.colStatus(i))
1912 {
1913 case SPxBasisBase<R>::Desc::P_ON_LOWER :
1914 p_vector[i] = SPxLPBase<R>::lower(i);
1915 break;
1916
1917 case SPxBasisBase<R>::Desc::P_ON_UPPER :
1918 case SPxBasisBase<R>::Desc::P_FIXED :
1919 p_vector[i] = SPxLPBase<R>::upper(i);
1920 break;
1921
1922 case SPxBasisBase<R>::Desc::P_FREE :
1923 p_vector[i] = 0;
1924 break;
1925
1926 case SPxBasisBase<R>::Desc::D_FREE :
1927 case SPxBasisBase<R>::Desc::D_ON_UPPER :
1928 case SPxBasisBase<R>::Desc::D_ON_LOWER :
1929 case SPxBasisBase<R>::Desc::D_ON_BOTH :
1930 case SPxBasisBase<R>::Desc::D_UNDEFINED :
1931 break;
1932
1933 default:
1934 throw SPxInternalCodeException("XSOLVE07 This should never happen.");
1935 }
1936 }
1937
1938 for(int j = 0; j < dim(); ++j)
1939 {
1940 if(this->baseId(j).isSPxColId())
1941 p_vector[ this->number(SPxColId(this->baseId(j))) ] = fVec()[j];
1942 }
1943 }
1944
1945 return status();
1946 }
1947
1948 template <class R>
1949 typename SPxSolverBase<R>::Status SPxSolverBase<R>::getDualSol(VectorBase<R>& p_vector) const
1950 {
1951
1952 assert(isInitialized());
1953
1954 if(!isInitialized())
1955 {
1956 /* exit if presolving/simplifier cleared the problem */
1957 if(status() == NO_PROBLEM)
1958 return status();
1959
1960 throw SPxStatusException("XSOLVE08 No Problem loaded");
1961 }
1962
1963 if(rep() == ROW)
1964 {
1965 int i;
1966 p_vector = this->maxRowObj();
1967
1968 for(i = this->nCols() - 1; i >= 0; --i)
1969 {
1970 if(this->baseId(i).isSPxRowId())
1971 p_vector[ this->number(SPxRowId(this->baseId(i))) ] = fVec()[i];
1972 }
1973 }
1974 else
1975 {
1976 const typename SPxBasisBase<R>::Desc& ds = this->desc();
1977
1978 for(int i = 0; i < this->nRows(); ++i)
1979 {
1980 switch(ds.rowStatus(i))
1981 {
1982 case SPxBasisBase<R>::Desc::D_FREE:
1983 case SPxBasisBase<R>::Desc::D_ON_UPPER:
1984 case SPxBasisBase<R>::Desc::D_ON_LOWER:
1985 case SPxBasisBase<R>::Desc::D_ON_BOTH:
1986 case SPxBasisBase<R>::Desc::D_UNDEFINED:
1987 p_vector[i] = 0;
1988 break;
1989
1990 default:
1991 p_vector[i] = (*theCoPvec)[i];
1992 }
1993 }
1994 }
1995
1996 p_vector *= Real(this->spxSense());
1997
1998 return status();
1999 }
2000
2001 template <class R>
2002 typename SPxSolverBase<R>::Status SPxSolverBase<R>::getRedCostSol(VectorBase<R>& p_vector) const
2003 {
2004
2005 assert(isInitialized());
2006
2007 if(!isInitialized())
2008 {
2009 throw SPxStatusException("XSOLVE09 No Problem loaded");
2010 // return NOT_INIT;
2011 }
2012
2013 if(rep() == ROW)
2014 {
2015 int i;
2016 p_vector.clear();
2017
2018 if(this->spxSense() == SPxLPBase<R>::MINIMIZE)
2019 {
2020 for(i = dim() - 1; i >= 0; --i)
2021 {
2022 if(this->baseId(i).isSPxColId())
2023 p_vector[ this->number(SPxColId(this->baseId(i))) ] = -fVec()[i];
2024 }
2025 }
2026 else
2027 {
2028 for(i = dim() - 1; i >= 0; --i)
2029 {
2030 if(this->baseId(i).isSPxColId())
2031 p_vector[ this->number(SPxColId(this->baseId(i))) ] = fVec()[i];
2032 }
2033 }
2034 }
2035 else
2036 {
2037 const typename SPxBasisBase<R>::Desc& ds = this->desc();
2038
2039 for(int i = 0; i < this->nCols(); ++i)
2040 {
2041 switch(ds.colStatus(i))
2042 {
2043 case SPxBasisBase<R>::Desc::D_FREE:
2044 case SPxBasisBase<R>::Desc::D_ON_UPPER:
2045 case SPxBasisBase<R>::Desc::D_ON_LOWER:
2046 case SPxBasisBase<R>::Desc::D_ON_BOTH:
2047 case SPxBasisBase<R>::Desc::D_UNDEFINED:
2048 p_vector[i] = 0;
2049 break;
2050
2051 default:
2052 p_vector[i] = this->maxObj()[i] - (*thePvec)[i];
2053 }
2054 }
2055
2056 if(this->spxSense() == SPxLPBase<R>::MINIMIZE)
2057 p_vector *= -1.0;
2058 }
2059
2060 return status();
2061 }
2062
2063 template <class R>
2064 typename SPxSolverBase<R>::Status SPxSolverBase<R>::getPrimalray(VectorBase<R>& p_vector) const
2065 {
2066
2067 assert(isInitialized());
2068
2069 if(!isInitialized())
2070 {
2071 throw SPxStatusException("XSOLVE10 No Problem loaded");
2072 // return NOT_INIT;
2073 }
2074
2075 assert(SPxBasisBase<R>::status() == SPxBasisBase<R>::UNBOUNDED);
2076 p_vector.clear();
2077 p_vector = primalRay;
2078
2079 return status();
2080 }
2081
2082 template <class R>
2083 typename SPxSolverBase<R>::Status SPxSolverBase<R>::getDualfarkas(VectorBase<R>& p_vector) const
2084 {
2085
2086 assert(isInitialized());
2087
2088 if(!isInitialized())
2089 {
2090 throw SPxStatusException("XSOLVE10 No Problem loaded");
2091 // return NOT_INIT;
2092 }
2093
2094 assert(SPxBasisBase<R>::status() == SPxBasisBase<R>::INFEASIBLE);
2095 p_vector.clear();
2096 p_vector = dualFarkas;
2097
2098 return status();
2099 }
2100
2101 template <class R>
2102 typename SPxSolverBase<R>::Status SPxSolverBase<R>::getSlacks(VectorBase<R>& p_vector) const
2103 {
2104
2105 assert(isInitialized());
2106
2107 if(!isInitialized())
2108 {
2109 throw SPxStatusException("XSOLVE11 No Problem loaded");
2110 // return NOT_INIT;
2111 }
2112
2113 if(rep() == COLUMN)
2114 {
2115 int i;
2116 const typename SPxBasisBase<R>::Desc& ds = this->desc();
2117
2118 for(i = this->nRows() - 1; i >= 0; --i)
2119 {
2120 switch(ds.rowStatus(i))
2121 {
2122 case SPxBasisBase<R>::Desc::P_ON_LOWER :
2123 p_vector[i] = this->lhs(i);
2124 break;
2125
2126 case SPxBasisBase<R>::Desc::P_ON_UPPER :
2127 case SPxBasisBase<R>::Desc::P_FIXED :
2128 p_vector[i] = this->rhs(i);
2129 break;
2130
2131 case SPxBasisBase<R>::Desc::P_FREE :
2132 p_vector[i] = 0;
2133 break;
2134
2135 case SPxBasisBase<R>::Desc::D_FREE :
2136 case SPxBasisBase<R>::Desc::D_ON_UPPER :
2137 case SPxBasisBase<R>::Desc::D_ON_LOWER :
2138 case SPxBasisBase<R>::Desc::D_ON_BOTH :
2139 case SPxBasisBase<R>::Desc::D_UNDEFINED :
2140 break;
2141
2142 default:
2143 throw SPxInternalCodeException("XSOLVE12 This should never happen.");
2144 }
2145 }
2146
2147 for(i = dim() - 1; i >= 0; --i)
2148 {
2149 if(this->baseId(i).isSPxRowId())
2150 p_vector[ this->number(SPxRowId(this->baseId(i))) ] = -(*theFvec)[i];
2151 }
2152 }
2153 else
2154 p_vector = pVec();
2155
2156 return status();
2157 }
2158
2159 template <class R>
2160 void SPxSolverBase<R>::setPrimal(VectorBase<R>& p_vector)
2161 {
2162
2163 if(!isInitialized())
2164 {
2165 throw SPxStatusException("XSOLVE20 Not Initialized");
2166 }
2167
2168 if(rep() == ROW)
2169 coPvec() = p_vector;
2170 else
2171 {
2172 for(int j = 0; j < dim(); ++j)
2173 {
2174 if(this->baseId(j).isSPxColId())
2175 fVec()[j] = p_vector[ this->number(SPxColId(this->baseId(j))) ];
2176 }
2177 }
2178 }
2179
2180 template <class R>
2181 void SPxSolverBase<R>::setDual(VectorBase<R>& p_vector)
2182 {
2183
2184 assert(isInitialized());
2185
2186 if(!isInitialized())
2187 {
2188 throw SPxStatusException("XSOLVE21 Not Initialized");
2189 }
2190
2191 if(rep() == ROW)
2192 {
2193 for(int i = this->nCols() - 1; i >= 0; --i)
2194 {
2195 if(this->baseId(i).isSPxRowId())
2196 {
2197 if(this->spxSense() == SPxLPBase<R>::MAXIMIZE)
2198 fVec()[i] = p_vector[ this->number(SPxRowId(this->baseId(i))) ];
2199 else
2200 fVec()[i] = -p_vector[ this->number(SPxRowId(this->baseId(i))) ];
2201 }
2202 }
2203 }
2204 else
2205 {
2206 coPvec() = p_vector;
2207
2208 if(this->spxSense() == SPxLPBase<R>::MINIMIZE)
2209 coPvec() *= -1.0;
2210 }
2211 }
2212
2213 template <class R>
2214 void SPxSolverBase<R>::setRedCost(VectorBase<R>& p_vector)
2215 {
2216
2217 assert(isInitialized());
2218
2219 if(!isInitialized())
2220 {
2221 throw SPxStatusException("XSOLVE22 Not Initialized");
2222 }
2223
2224 if(rep() == ROW)
2225 {
2226 for(int i = dim() - 1; i >= 0; --i)
2227 {
2228 if(this->baseId(i).isSPxColId())
2229 {
2230 if(this->spxSense() == SPxLPBase<R>::MINIMIZE)
2231 fVec()[i] = -p_vector[ this->number(SPxColId(this->baseId(i))) ];
2232 else
2233 fVec()[i] = p_vector[ this->number(SPxColId(this->baseId(i))) ];
2234 }
2235 }
2236 }
2237 else
2238 {
2239 pVec() = this->maxObj();
2240
2241 if(this->spxSense() == SPxLPBase<R>::MINIMIZE)
2242 pVec() += p_vector;
2243 else
2244 pVec() -= p_vector;
2245 }
2246 }
2247
2248 template <class R>
2249 void SPxSolverBase<R>::setSlacks(VectorBase<R>& p_vector)
2250 {
2251
2252 assert(isInitialized());
2253
2254 if(!isInitialized())
2255 {
2256 throw SPxStatusException("XSOLVE23 Not Initialized");
2257 }
2258
2259 if(rep() == COLUMN)
2260 {
2261 for(int i = dim() - 1; i >= 0; --i)
2262 {
2263 if(this->baseId(i).isSPxRowId())
2264 (*theFvec)[i] = -p_vector[ this->number(SPxRowId(this->baseId(i))) ];
2265 }
2266 }
2267 else
2268 pVec() = p_vector;
2269 }
2270
2271 template <class R>
2272 typename SPxSolverBase<R>::Status SPxSolverBase<R>::status() const
2273 {
2274 switch(m_status)
2275 {
2276 case UNKNOWN :
2277 switch(SPxBasisBase<R>::status())
2278 {
2279 case SPxBasisBase<R>::NO_PROBLEM :
2280 return NO_PROBLEM;
2281
2282 case SPxBasisBase<R>::SINGULAR :
2283 return SINGULAR;
2284
2285 case SPxBasisBase<R>::REGULAR :
2286 case SPxBasisBase<R>::DUAL :
2287 case SPxBasisBase<R>::PRIMAL :
2288 return UNKNOWN;
2289
2290 case SPxBasisBase<R>::OPTIMAL :
2291 return OPTIMAL;
2292
2293 case SPxBasisBase<R>::UNBOUNDED :
2294 return UNBOUNDED;
2295
2296 case SPxBasisBase<R>::INFEASIBLE :
2297 return INFEASIBLE;
2298
2299 default:
2300 return ERROR;
2301 }
2302
2303 case SINGULAR :
2304 return m_status;
2305
2306 case OPTIMAL :
2307 assert(SPxBasisBase<R>::status() == SPxBasisBase<R>::OPTIMAL);
2308
2309 /*lint -fallthrough*/
2310 case ABORT_EXDECOMP :
2311 case ABORT_DECOMP :
2312 case ABORT_CYCLING :
2313 case ABORT_TIME :
2314 case ABORT_ITER :
2315 case ABORT_VALUE :
2316 case RUNNING :
2317 case REGULAR :
2318 case NOT_INIT :
2319 case NO_SOLVER :
2320 case NO_PRICER :
2321 case NO_RATIOTESTER :
2322 case ERROR:
2323 return m_status;
2324
2325 default:
2326 return ERROR;
2327 }
2328 }
2329
2330 template <class R>
2331 typename SPxSolverBase<R>::Status SPxSolverBase<R>::getResult(
2332 R* p_value,
2333 VectorBase<R>* p_primal,
2334 VectorBase<R>* p_slacks,
2335 VectorBase<R>* p_dual,
2336 VectorBase<R>* reduCosts)
2337 {
2338 if(p_value)
2339 *p_value = this->value();
2340
2341 if(p_primal)
2342 getPrimalSol(*p_primal);
2343
2344 if(p_slacks)
2345 getSlacks(*p_slacks);
2346
2347 if(p_dual)
2348 getDualSol(*p_dual);
2349
2350 if(reduCosts)
2351 getRedCostSol(*reduCosts);
2352
2353 return status();
2354 }
2355 } // namespace soplex
2356