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