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);
(10) Event write_zero_model: "DIdxSet" sets "candidates.idx" to "nullptr". [details]
Also see events: [no_write_call][var_deref_model]
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 	            {
(23) Event var_deref_model: "index" dereferences null "candidates.idx". [details]
Also see events: [write_zero_model][no_write_call]
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