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 <iostream>
26   	#include <assert.h>
27   	
28   	#include "soplex.h"
29   	#include "soplex/statistics.h"
30   	
31   	#define SOPLEX_ALLOWED_UNSCALE_PERCENTAGE    0.1
32   	#define SOPLEX_MIN_OPT_CALLS_WITH_SCALING     10
33   	
34   	namespace soplex
35   	{
36   	/// solves real LP
37   	template <class R>
38   	void SoPlexBase<R>::_optimize(volatile bool* interrupt)
39   	{
40   	   assert(_realLP != 0);
41   	   assert(_realLP == &_solver);
42   	
43   	   _solReal.invalidate();
44   	   ++_optimizeCalls;
45   	
46   	   // start timing
47   	   _statistics->solvingTime->start();
48   	
49   	   if(boolParam(SoPlexBase<R>::PERSISTENTSCALING))
50   	   {
51   	      // scale original problem; overwriting _realLP
52   	      if(_scaler && !_realLP->isScaled() && _reapplyPersistentScaling())
53   	      {
54   	#ifdef SOPLEX_DEBUG
55   	         SPxLPBase<R>* origLP = 0;
56   	         spx_alloc(origLP);
57   	         origLP = new(origLP) SPxLPBase<R>(*_realLP);
58   	#endif
59   	         _scaler->scale(*_realLP, true);
60   	         _isRealLPScaled = _realLP->isScaled(); // a scaler might decide not to apply scaling
61   	         _solver.invalidateBasis();
62   	#ifdef SOPLEX_DEBUG
63   	         _checkScaling(origLP);
64   	#endif
65   	      }
66   	      // unscale previously scaled problem, overwriting _realLP
67   	      else if(!_scaler && _realLP->isScaled())
68   	      {
69   	         _solver.unscaleLPandReloadBasis();
70   	         _isRealLPScaled = false;
71   	         ++_unscaleCalls;
72   	      }
73   	   }
74   	
75   	   // remember that last solve was in floating-point
76   	   _lastSolveMode = SOLVEMODE_REAL;
77   	
78   	   // solve and store solution; if we have a starting basis, do not apply preprocessing; if we are solving from
79   	   // scratch, apply preprocessing according to parameter settings
80   	   if(!_hasBasis && realParam(SoPlexBase<R>::OBJLIMIT_LOWER) == -realParam(SoPlexBase<R>::INFTY)
81   	         && realParam(SoPlexBase<R>::OBJLIMIT_UPPER) == realParam(SoPlexBase<R>::INFTY))
82   	      _preprocessAndSolveReal(true, interrupt);
83   	   else
84   	      _preprocessAndSolveReal(false, interrupt);
85   	
86   	   _statistics->finalBasisCondition = _solver.getBasisMetric(0);
87   	
88   	   // stop timing
89   	   _statistics->solvingTime->stop();
90   	}
91   	
92   	
93   	
94   	/// check whether persistent scaling is supposed to be reapplied again after unscaling
95   	template <class R>
96   	bool SoPlexBase<R>::_reapplyPersistentScaling() const
97   	{
98   	   if((_unscaleCalls > _optimizeCalls * SOPLEX_ALLOWED_UNSCALE_PERCENTAGE)
99   	         && _optimizeCalls > SOPLEX_MIN_OPT_CALLS_WITH_SCALING)
100  	      return false;
101  	   else
102  	      return true;
103  	}
104  	
105  	
106  	
107  	/// checks result of the solving process and solves again without preprocessing if necessary
108  	template <class R>
109  	void SoPlexBase<R>::_evaluateSolutionReal(typename SPxSimplifier<R>::Result simplificationStatus)
110  	{
111  	   // if the simplifier detected infeasibility or unboundedness we optimize again
112  	   // just to get the proof (primal or dual ray)
113  	   // todo get infeasibility proof from simplifier
114  	   switch(simplificationStatus)
115  	   {
116  	   case SPxSimplifier<R>::INFEASIBLE:
117  	   case SPxSimplifier<R>::DUAL_INFEASIBLE:
118  	   case SPxSimplifier<R>::UNBOUNDED:
119  	      _hasBasis = false;
120  	
121  	      if(boolParam(SoPlexBase<R>::ENSURERAY))
122  	      {
123  	         SPX_MSG_INFO1(spxout, spxout <<
124  	                       "simplifier detected infeasibility or unboundedness - solve again without simplifying" << std::endl;
125  	                      )
126  	         _preprocessAndSolveReal(false);
127  	      }
128  	      else
129  	      {
130  	         if(simplificationStatus == SPxSimplifier<R>::INFEASIBLE)
131  	            _status = SPxSolverBase<R>::INFEASIBLE;
132  	         else if(simplificationStatus == SPxSimplifier<R>::UNBOUNDED)
133  	            _status = SPxSolverBase<R>::UNBOUNDED;
134  	         else
135  	            _status = SPxSolverBase<R>::INForUNBD;
136  	
137  	         // load original LP to restore clean solver state
138  	         _loadRealLP(false);
139  	      }
140  	
141  	      return;
142  	
143  	   case SPxSimplifier<R>::VANISHED:
144  	      _status = SPxSolverBase<R>::OPTIMAL;
145  	      _storeSolutionRealFromPresol();
146  	      return;
147  	
148  	   case SPxSimplifier<R>::OKAY:
149  	      _status = _solver.status();
150  	   }
151  	
152  	   // process result
153  	   switch(_status)
154  	   {
155  	   case SPxSolverBase<R>::OPTIMAL:
156  	      _storeSolutionReal(!_isRealLPLoaded || _isRealLPScaled);
157  	
158  	      // apply polishing on original problem
159  	      if(_applyPolishing)
160  	      {
161  	         int polishing = intParam(SoPlexBase<R>::SOLUTION_POLISHING);
162  	         setIntParam(SoPlexBase<R>::SOLUTION_POLISHING, polishing);
163  	         _preprocessAndSolveReal(false);
164  	      }
165  	
166  	      break;
167  	
168  	   case SPxSolverBase<R>::UNBOUNDED:
169  	   case SPxSolverBase<R>::INFEASIBLE:
170  	   case SPxSolverBase<R>::INForUNBD:
171  	
172  	      // in case of infeasibility or unboundedness, we currently can not unsimplify, but have to solve the original LP again
173  	      if(!_isRealLPLoaded && boolParam(SoPlexBase<R>::ENSURERAY))
174  	      {
175  	         SPX_MSG_INFO1(spxout, spxout << " --- loading original problem" << std::endl;)
176  	         _solver.changeObjOffset(realParam(SoPlexBase<R>::OBJ_OFFSET));
177  	         // we cannot do more to remove violations
178  	         _resolveWithoutPreprocessing(simplificationStatus);
179  	      }
180  	      else
181  	      {
182  	         _storeSolutionReal(false);
183  	      }
184  	
185  	      break;
186  	
187  	   case SPxSolverBase<R>::SINGULAR:
188  	
189  	      // if preprocessing was applied, try to run again without to avoid singularity
190  	      if(!_isRealLPLoaded)
191  	      {
192  	         SPX_MSG_INFO1(spxout, spxout <<
193  	                       "encountered singularity - trying to solve again without simplifying" <<
194  	                       std::endl;)
195  	         _preprocessAndSolveReal(false);
196  	         return;
197  	      }
198  	
199  	      _hasBasis = false;
200  	      break;
201  	
202  	   case SPxSolverBase<R>::ABORT_VALUE:
203  	
204  	      // If we aborted the solve for some reason and there is still a shift, ensure that the basis status is correct
205  	      if(_solver.shift() > _solver.epsilon())
206  	         _solver.setBasisStatus(SPxBasisBase<R>::REGULAR);
207  	
208  	      _storeSolutionReal(true);
209  	      break;
210  	
211  	   case SPxSolverBase<R>::ABORT_CYCLING:
212  	
213  	      // if preprocessing or scaling was applied, try to run again without to
214  	      // avoid cycling
215  	      if(!_isRealLPLoaded || _isRealLPScaled)
216  	      {
217  	         SPX_MSG_INFO1(spxout, spxout << "encountered cycling - trying to solve again without simplifying" <<
218  	                       std::endl;)
219  	         // store and unsimplify sub-optimal solution and basis, may trigger re-solve
220  	         _storeSolutionReal(true);
221  	         return;
222  	      }
223  	
224  	      if(_solReal.isPrimalFeasible() || _solReal.isDualFeasible())
225  	         _status = SPxSolverBase<R>::OPTIMAL_UNSCALED_VIOLATIONS;
226  	
227  	   // FALLTHROUGH
228  	   case SPxSolverBase<R>::ABORT_TIME:
229  	   case SPxSolverBase<R>::ABORT_ITER:
230  	   case SPxSolverBase<R>::REGULAR:
231  	   case SPxSolverBase<R>::RUNNING:
232  	
233  	      // If we aborted the solve for some reason and there is still a shift, ensure that the basis status is correct
234  	      if(_solver.shift() > _solver.epsilon())
235  	         _solver.setBasisStatus(SPxBasisBase<R>::REGULAR);
236  	
237  	      _storeSolutionReal(false);
238  	      break;
239  	
240  	   default:
241  	      _hasBasis = false;
242  	      break;
243  	   }
244  	}
245  	
246  	
247  	
248  	/// solves real LP with/without preprocessing
249  	template <class R>
250  	void SoPlexBase<R>::_preprocessAndSolveReal(bool applySimplifier, volatile bool* interrupt)
251  	{
252  	   _solver.changeObjOffset(realParam(SoPlexBase<R>::OBJ_OFFSET));
253  	   _statistics->preprocessingTime->start();
254  	
255  	   _applyPolishing = false;
256  	
257  	   if(applySimplifier)
258  	      _enableSimplifierAndScaler();
259  	   else
260  	      _disableSimplifierAndScaler();
261  	
262  	   // create a copy of the LP when simplifying or when using internal scaling, i.e. w/o persistent scaling
263  	   bool copyLP = (_simplifier != 0 || (_scaler && !_isRealLPScaled));
264  	
265  	   _solver.setTerminationValue(intParam(SoPlexBase<R>::OBJSENSE) == SoPlexBase<R>::OBJSENSE_MINIMIZE
266  	                               ? realParam(SoPlexBase<R>::OBJLIMIT_UPPER) : realParam(SoPlexBase<R>::OBJLIMIT_LOWER));
267  	
268  	   if(_isRealLPLoaded)
269  	   {
270  	      assert(_realLP == &_solver);
271  	
272  	      // preprocessing is always applied to the LP in the solver; hence we have to create a copy of the original LP
273  	      // if simplifier is turned on
274  	      if(copyLP)
275  	      {
276  	         _realLP = nullptr;
277  	         spx_alloc(_realLP);
278  	         _realLP = new(_realLP) SPxLPBase<R>(_solver);
279  	         _isRealLPLoaded = false;
280  	      }
281  	   }
282  	   else
283  	   {
284  	      assert(_realLP != &_solver);
285  	
286  	      // load real LP and basis if available
287  	      if(_hasBasis)
288  	      {
289  	         assert(_basisStatusRows.size() == numRows());
290  	         assert(_basisStatusCols.size() == this->numCols());
291  	
292  	         _solver.loadLP(*_realLP, false);
293  	         _solver.setBasis(_basisStatusRows.get_const_ptr(), _basisStatusCols.get_const_ptr());
294  	      }
295  	      // load real LP and set up slack basis
296  	      else
297  	         _solver.loadLP(*_realLP, true);
298  	
299  	      // if there is no simplifier, then the original and the transformed problem are identical and it is more
300  	      // memory-efficient to keep only the problem in the solver
301  	      if(!copyLP)
302  	      {
303  	         _realLP->~SPxLPBase<R>();
304  	         spx_free(_realLP);
305  	         _realLP = &_solver;
306  	         _isRealLPLoaded = true;
307  	      }
308  	   }
309  	
310  	   // assert that we have two problems if and only if we apply the simplifier
311  	   assert(_realLP == &_solver || copyLP);
312  	   assert(_realLP != &_solver || !copyLP);
313  	
314  	   // apply problem simplification
315  	   typename SPxSimplifier<R>::Result simplificationStatus = SPxSimplifier<R>::OKAY;
316  	
317  	   if(_simplifier)
318  	   {
319  	      assert(!_isRealLPLoaded);
320  	      // do not remove bounds of boxed variables or sides of ranged rows if bound flipping is used; also respect row-boundflip parameter
321  	      bool keepbounds = intParam(SoPlexBase<R>::RATIOTESTER) == SoPlexBase<R>::RATIOTESTER_BOUNDFLIPPING;
322  	
323  	      if(intParam(SoPlexBase<R>::REPRESENTATION) == SoPlexBase<R>::REPRESENTATION_ROW
324  	            || (intParam(SoPlexBase<R>::REPRESENTATION) == SoPlexBase<R>::REPRESENTATION_AUTO
325  	                && (_solver.nCols() + 1) * realParam(SoPlexBase<R>::REPRESENTATION_SWITCH) < (_solver.nRows() + 1)))
326  	         keepbounds &= boolParam(SoPlexBase<R>::ROWBOUNDFLIPS);
327  	
328  	      Real remainingTime = _solver.getMaxTime() - _solver.time();
329  	      simplificationStatus = _simplifier->simplify(_solver, remainingTime, keepbounds,
330  	                             _solver.random.getSeed());
331  	      _solver.changeObjOffset(_simplifier->getObjoffset() + realParam(SoPlexBase<R>::OBJ_OFFSET));
332  	      _solver.setScalingInfo(false);
333  	      _applyPolishing = true;
334  	
335  	      _solver.setSolutionPolishing(SPxSolverBase<R>::POLISH_OFF);
336  	   }
337  	
338  	   _statistics->preprocessingTime->stop();
339  	
340  	   // run the simplex method if problem has not been solved by the simplifier
341  	   if(simplificationStatus == SPxSimplifier<R>::OKAY)
342  	   {
343  	      if(_scaler && !_solver.isScaled())
344  	      {
345  	         _scaler->scale(_solver, false);
346  	         _solver.invalidateBasis();
347  	      }
348  	
349  	      _solveRealLPAndRecordStatistics(interrupt);
350  	   }
351  	
352  	   _evaluateSolutionReal(simplificationStatus);
353  	}
354  	
355  	
356  	
357  	/// loads original problem into solver and solves again after it has been solved to infeasibility or unboundedness with preprocessing
358  	template <class R>
359  	void SoPlexBase<R>::_resolveWithoutPreprocessing(typename SPxSimplifier<R>::Result
360  	      simplificationStatus)
361  	{
362  	   assert(!_isRealLPLoaded || _scaler != nullptr);
363  	   assert(_simplifier != 0 || _scaler != nullptr);
364  	   assert(_status == SPxSolverBase<R>::INFEASIBLE || _status == SPxSolverBase<R>::INForUNBD
365  	          || _status == SPxSolverBase<R>::UNBOUNDED);
366  	
367  	   // if simplifier was active, then we unsimplify to get the basis
368  	   if(_simplifier)
369  	   {
370  	      assert(!_simplifier->isUnsimplified());
371  	      assert(simplificationStatus == SPxSimplifier<R>::OKAY);
372  	
373  	      // get temporary solution vectors for transformed problem
374  	      VectorBase<R> primal(_solver.nCols());
375  	      VectorBase<R> slacks(_solver.nRows());
376  	      VectorBase<R> dual(_solver.nRows());
377  	      VectorBase<R> redCost(_solver.nCols());
378  	
379  	      _basisStatusRows.reSize(numRows());
380  	      _basisStatusCols.reSize(numCols());
381  	      assert(_basisStatusRows.size() >= _solver.nRows());
382  	      assert(_basisStatusCols.size() >= _solver.nCols());
383  	
384  	      // get solution data from transformed problem
385  	      _solver.getPrimalSol(primal);
386  	      _solver.getSlacks(slacks);
387  	      _solver.getDualSol(dual);
388  	      _solver.getRedCostSol(redCost);
389  	
390  	      // unscale vectors
391  	      if(_scaler && _solver.isScaled())
392  	      {
393  	         _scaler->unscalePrimal(_solver, primal);
394  	         _scaler->unscaleSlacks(_solver, slacks);
395  	         _scaler->unscaleDual(_solver, dual);
396  	         _scaler->unscaleRedCost(_solver, redCost);
397  	      }
398  	
399  	      // get basis of transformed problem
400  	      _solver.getBasis(_basisStatusRows.get_ptr(), _basisStatusCols.get_ptr(), _basisStatusRows.size(),
401  	                       _basisStatusCols.size());
402  	
403  	      try
404  	      {
405  	         _simplifier->unsimplify(primal, dual, slacks, redCost, _basisStatusRows.get_ptr(),
406  	                                 _basisStatusCols.get_ptr(), false);
407  	         _simplifier->getBasis(_basisStatusRows.get_ptr(), _basisStatusCols.get_ptr(),
408  	                               _basisStatusRows.size(), _basisStatusCols.size());
409  	         _hasBasis = true;
410  	      }
411  	      catch(const SPxException& E)
412  	      {
413  	         SPX_MSG_INFO1(spxout, spxout << "Caught exception <" << E.what() <<
414  	                       "> during unsimplification. Resolving without simplifier and scaler.\n");
415  	         _hasBasis = false;
416  	      }
417  	   }
418  	   // if the original problem is not in the solver because of scaling, we also need to store the basis
419  	   else if(_scaler != nullptr)
420  	   {
421  	      _basisStatusRows.reSize(numRows());
422  	      _basisStatusCols.reSize(numCols());
423  	      assert(_basisStatusRows.size() == _solver.nRows());
424  	      assert(_basisStatusCols.size() == _solver.nCols());
425  	
426  	      _solver.getBasis(_basisStatusRows.get_ptr(), _basisStatusCols.get_ptr(), _basisStatusRows.size(),
427  	                       _basisStatusCols.size());
428  	      _hasBasis = true;
429  	   }
430  	
431  	   // resolve the original problem
432  	   _preprocessAndSolveReal(false);
433  	   return;
434  	}
435  	
436  	
437  	
438  	/// verify computed solution based on status and resolve if claimed primal or dual feasibility is not fulfilled
439  	template <class R>
440  	void SoPlexBase<R>::_verifySolutionReal()
441  	{
442  	   assert(_hasSolReal);
443  	
444  	   SPX_MSG_INFO1(spxout, spxout << " --- verifying computed solution" << std::endl;)
445  	
446  	   R sumviol = 0;
447  	   R boundviol = 0;
448  	   R rowviol = 0;
449  	   R dualviol = 0;
450  	   R redcostviol = 0;
451  	
452  	   (void) getBoundViolation(boundviol, sumviol);
453  	   (void) getRowViolation(rowviol, sumviol);
454  	   (void) getDualViolation(dualviol, sumviol);
455  	   (void) getRedCostViolation(redcostviol, sumviol);
456  	
457  	   if(boundviol >= _solver.tolerances()->floatingPointFeastol()
458  	         || rowviol >= _solver.tolerances()->floatingPointFeastol()
459  	         || dualviol >= _solver.tolerances()->floatingPointOpttol()
460  	         || redcostviol >= _solver.tolerances()->floatingPointOpttol())
461  	   {
462  	      assert(&_solver == _realLP);
463  	      assert(_isRealLPLoaded);
464  	      SPX_MSG_INFO3(spxout, spxout << "bound violation: " << boundviol
465  	                    << ", row violation: " << rowviol
466  	                    << ", dual violation: " << dualviol
467  	                    << ", redcost violation: " << redcostviol << std::endl;)
468  	      SPX_MSG_INFO1(spxout, spxout <<
469  	                    " --- detected violations in original problem space -- solve again without presolving/scaling" <<
470  	                    std::endl;)
471  	
472  	      if(_isRealLPScaled)
473  	      {
474  	         _solver.unscaleLPandReloadBasis();
475  	         _isRealLPScaled = false;
476  	         ++_unscaleCalls;
477  	      }
478  	
479  	      _preprocessAndSolveReal(false);
480  	   }
481  	}
482  	
483  	/// verify computed solution based on status and resolve if claimed primal or dual feasibility is not fulfilled
484  	template <class R>
485  	void SoPlexBase<R>::_verifyObjLimitReal()
486  	{
487  	   SPX_MSG_INFO1(spxout, spxout << " --- verifying objective limit" << std::endl;)
488  	
489  	   R sumviol = 0;
490  	   R dualviol = 0;
491  	   R redcostviol = 0;
492  	
493  	   bool dualfeasible = true;
494  	
495  	   if(!getDualViolation(dualviol, sumviol))
496  	      dualfeasible = false;
497  	
498  	   if(!getRedCostViolation(redcostviol, sumviol))
499  	      dualfeasible = false;
500  	
501  	   if(!dualfeasible || dualviol >= _solver.tolerances()->floatingPointOpttol()
502  	         || redcostviol >= _solver.tolerances()->floatingPointOpttol())
503  	   {
504  	      assert(&_solver == _realLP);
505  	      assert(_isRealLPLoaded);
506  	      SPX_MSG_INFO3(spxout, spxout << ", dual violation: " << dualviol
507  	                    << ", redcost violation: " << redcostviol << std::endl;)
508  	      SPX_MSG_INFO1(spxout, spxout <<
509  	                    " --- detected violations in original problem space -- solve again without presolving/scaling" <<
510  	                    std::endl;)
511  	
512  	      if(_isRealLPScaled)
513  	      {
514  	         _solver.unscaleLPandReloadBasis();
515  	         _isRealLPScaled = false;
516  	         ++_unscaleCalls;
517  	      }
518  	
519  	      _preprocessAndSolveReal(false);
520  	   }
521  	}
522  	
523  	
524  	/// stores solution data from the solver, possibly after applying unscaling and unsimplifying
525  	template <class R>
526  	void SoPlexBase<R>::_storeSolutionReal(bool verify)
527  	{
528  	   // prepare storage for basis (enough to fit the original basis)
529  	   _basisStatusRows.reSize(numRows());
530  	   _basisStatusCols.reSize(numCols());
531  	
532  	   // prepare storage for the solution data (only in transformed space due to unscaling), w/o setting it to zero
533  	   _solReal._primal.reDim(_solver.nCols(), false);
534  	   _solReal._slacks.reDim(_solver.nRows(), false);
535  	   _solReal._dual.reDim(_solver.nRows(), false);
536  	   _solReal._redCost.reDim(_solver.nCols(), false);
537  	
538  	   // check primal status consistency and query solution status
539  	   assert(_solver.basis().status() != SPxBasisBase<R>::PRIMAL || status() != SPxSolverBase<R>::ERROR);
540  	   assert(_solver.basis().status() != SPxBasisBase<R>::PRIMAL
541  	          || status() != SPxSolverBase<R>::NO_RATIOTESTER);
542  	   assert(_solver.basis().status() != SPxBasisBase<R>::PRIMAL
543  	          || status() != SPxSolverBase<R>::NO_PRICER);
544  	   assert(_solver.basis().status() != SPxBasisBase<R>::PRIMAL
545  	          || status() != SPxSolverBase<R>::NO_SOLVER);
546  	   assert(_solver.basis().status() != SPxBasisBase<R>::PRIMAL
547  	          || status() != SPxSolverBase<R>::NOT_INIT);
548  	   assert(_solver.basis().status() != SPxBasisBase<R>::PRIMAL
549  	          || status() != SPxSolverBase<R>::SINGULAR);
550  	   assert(_solver.basis().status() != SPxBasisBase<R>::PRIMAL
551  	          || status() != SPxSolverBase<R>::NO_PROBLEM);
552  	   assert(_solver.basis().status() != SPxBasisBase<R>::PRIMAL
553  	          || status() != SPxSolverBase<R>::UNBOUNDED);
554  	   assert(_solver.basis().status() != SPxBasisBase<R>::PRIMAL
555  	          || status() != SPxSolverBase<R>::INFEASIBLE);
556  	   assert(_solver.basis().status() != SPxBasisBase<R>::UNBOUNDED
557  	          || status() == SPxSolverBase<R>::UNBOUNDED);
558  	   assert(_solver.basis().status() == SPxBasisBase<R>::UNBOUNDED
559  	          || _solver.basis().status() == SPxBasisBase<R>::NO_PROBLEM
560  	          || status() != SPxSolverBase<R>::UNBOUNDED);
561  	
562  	   _solReal._isPrimalFeasible = (status() == SPxSolverBase<R>::OPTIMAL
563  	                                 || ((_solver.basis().status() == SPxBasisBase<R>::PRIMAL
564  	                                      || _solver.basis().status() == SPxBasisBase<R>::UNBOUNDED)
565  	                                     && _solver.shift() < 10.0 * realParam(SoPlexBase<R>::EPSILON_ZERO)));
566  	
567  	   _solReal._hasPrimalRay = (status() == SPxSolverBase<R>::UNBOUNDED && _isRealLPLoaded);
568  	
569  	   // check dual status consistency and query solution status
570  	   assert(_solver.basis().status() != SPxBasisBase<R>::DUAL || status() != SPxSolverBase<R>::ERROR);
571  	   assert(_solver.basis().status() != SPxBasisBase<R>::DUAL
572  	          || status() != SPxSolverBase<R>::NO_RATIOTESTER);
573  	   assert(_solver.basis().status() != SPxBasisBase<R>::DUAL
574  	          || status() != SPxSolverBase<R>::NO_PRICER);
575  	   assert(_solver.basis().status() != SPxBasisBase<R>::DUAL
576  	          || status() != SPxSolverBase<R>::NO_SOLVER);
577  	   assert(_solver.basis().status() != SPxBasisBase<R>::DUAL || status() != SPxSolverBase<R>::NOT_INIT);
578  	   assert(_solver.basis().status() != SPxBasisBase<R>::DUAL || status() != SPxSolverBase<R>::SINGULAR);
579  	   assert(_solver.basis().status() != SPxBasisBase<R>::DUAL
580  	          || status() != SPxSolverBase<R>::NO_PROBLEM);
581  	   assert(_solver.basis().status() != SPxBasisBase<R>::DUAL
582  	          || status() != SPxSolverBase<R>::UNBOUNDED);
583  	   assert(_solver.basis().status() != SPxBasisBase<R>::DUAL
584  	          || status() != SPxSolverBase<R>::INFEASIBLE);
585  	   assert(_solver.basis().status() != SPxBasisBase<R>::INFEASIBLE
586  	          || status() == SPxSolverBase<R>::INFEASIBLE);
587  	   assert(_solver.basis().status() == SPxBasisBase<R>::INFEASIBLE
588  	          || _solver.basis().status() == SPxBasisBase<R>::NO_PROBLEM
589  	          || status() != SPxSolverBase<R>::INFEASIBLE);
590  	
591  	   _solReal._isDualFeasible = (status() == SPxSolverBase<R>::OPTIMAL
592  	                               || ((_solver.basis().status() == SPxBasisBase<R>::DUAL
593  	                                    || _solver.basis().status() == SPxBasisBase<R>::INFEASIBLE)
594  	                                   && _solver.shift() < 10.0 * realParam(SoPlexBase<R>::EPSILON_ZERO)));
595  	
596  	   _solReal._hasDualFarkas = (status() == SPxSolverBase<R>::INFEASIBLE && _isRealLPLoaded);
597  	
598  	   // get infeasibility or unboundedness proof if available
599  	   if(_solReal._hasPrimalRay)
600  	   {
601  	      _solReal._primalRay.reDim(_solver.nCols(), false);
602  	      _solver.getPrimalray(_solReal._primalRay);
603  	   }
604  	
605  	   if(_solReal._hasDualFarkas)
606  	   {
607  	      _solReal._dualFarkas.reDim(_solver.nRows(), false);
608  	      _solver.getDualfarkas(_solReal._dualFarkas);
609  	   }
610  	
611  	   // get solution data from the solver; independent of solution status
612  	   _solver.getBasis(_basisStatusRows.get_ptr(), _basisStatusCols.get_ptr(),
613  	                    _basisStatusRows.size(), _basisStatusCols.size());
614  	
615  	   _solver.getPrimalSol(_solReal._primal);
616  	   _solver.getSlacks(_solReal._slacks);
617  	   _solver.getDualSol(_solReal._dual);
618  	   _solver.getRedCostSol(_solReal._redCost);
619  	
620  	   _hasBasis = true;
621  	
622  	   // get primal and/or dual objective function value depending on status
623  	   _solver.forceRecompNonbasicValue();
624  	   _solReal._objVal = _solver.objValue();
625  	
626  	   // infeasible solutions shall also be stored and be accessible
627  	   _hasSolReal = true;
628  	
629  	   // unscale vectors
630  	   if(_solver.isScaled() && !_isRealLPLoaded)
631  	      _unscaleSolutionReal(_solver, false);
632  	
633  	   // get unsimplified solution data from simplifier
634  	   if(_simplifier)
635  	   {
636  	      assert(!_simplifier->isUnsimplified());
637  	      assert(_simplifier->result() == SPxSimplifier<R>::OKAY);
638  	      assert(_realLP != &_solver);
639  	
640  	      typename SPxBasisBase<R>::SPxStatus simplifiedBasisStatus = _solver.getBasisStatus();
641  	
642  	      try
643  	      {
644  	         // pass solution data of transformed problem to simplifier
645  	         _simplifier->unsimplify(_solReal._primal, _solReal._dual, _solReal._slacks, _solReal._redCost,
646  	                                 _basisStatusRows.get_ptr(), _basisStatusCols.get_ptr(), status() == SPxSolverBase<R>::OPTIMAL);
647  	      }
648  	      catch(const SPxException& E)
649  	      {
650  	         SPX_MSG_INFO1(spxout, spxout << "Caught exception <" << E.what() <<
651  	                       "> during unsimplification. Resolving without simplifier and scaler.\n");
652  	         _hasBasis = false;
653  	         _preprocessAndSolveReal(false);
654  	         return;
655  	      }
656  	
657  	      // copy unsimplified solution data from simplifier (size and dimension is adapted during copy)
658  	      _solReal._primal  = _simplifier->unsimplifiedPrimal();
659  	      _solReal._slacks  = _simplifier->unsimplifiedSlacks();
660  	      _solReal._dual    = _simplifier->unsimplifiedDual();
661  	      _solReal._redCost = _simplifier->unsimplifiedRedCost();
662  	
663  	      // overwrite the transformed basis with the unsimplified one
664  	      _simplifier->getBasis(_basisStatusRows.get_ptr(), _basisStatusCols.get_ptr(),
665  	                            _basisStatusRows.size(), _basisStatusCols.size());
666  	
667  	      // load original problem but don't setup a slack basis
668  	      _loadRealLP(false);
669  	
670  	      // since we presolved the problem, we should not allow access to the dual norms
671  	      _solver.weightsAreSetup = false;
672  	
673  	      assert(_realLP == &_solver);
674  	
675  	      // reset basis status
676  	      _solver.setBasisStatus(simplifiedBasisStatus);
677  	
678  	      _solver.setBasis(_basisStatusRows.get_const_ptr(), _basisStatusCols.get_const_ptr());
679  	      // load unsimplified basis into solver
680  	      assert(_basisStatusRows.size() == numRows());
681  	      assert(_basisStatusCols.size() == this->numCols());
682  	      _hasBasis = true;
683  	   }
684  	
685  	   // load realLP into the solver again (internal scaling was applied)
686  	   else if(_realLP != &_solver)
687  	   {
688  	      assert(_solver.isScaled());
689  	      _loadRealLP(false);
690  	   }
691  	
692  	   // unscale stored solution (removes persistent scaling)
693  	   if(_isRealLPScaled)
694  	      _unscaleSolutionReal(*_realLP, true);
695  	
696  	   // check solution for violations and solve again if necessary
697  	   if(verify)
698  	   {
699  	      if(_status == SPxSolverBase<R>::ABORT_VALUE)
700  	         _verifyObjLimitReal();
701  	      else
702  	         _verifySolutionReal();
703  	   }
704  	
705  	   assert(_solver.nCols() == this->numCols());
706  	   assert(_solver.nRows() == numRows());
707  	}
708  	
709  	
710  	
711  	template <class R>
712  	void SoPlexBase<R>::_storeSolutionRealFromPresol()
713  	{
714  	   assert(_simplifier);
715  	   assert(_simplifier->result() == SPxSimplifier<R>::VANISHED);
716  	
717  	   // prepare storage for basis (enough to fit the original basis)
718  	   _basisStatusRows.reSize(numRows());
719  	   _basisStatusCols.reSize(numCols());
720  	
721  	   // prepare storage for the solution data and initialize it to zero
722  	   _solReal._primal.reDim(numCols(), true);
723  	   _solReal._slacks.reDim(numRows(), true);
724  	   _solReal._dual.reDim(numRows(), true);
725  	   _solReal._redCost.reDim(numCols(), true);
726  	
727  	   // load original LP and setup slack basis for unsimplifying
728  	   _loadRealLP(true);
729  	
730  	   // store slack basis
731  	   _solver.getBasis(_basisStatusRows.get_ptr(), _basisStatusCols.get_ptr(),
732  	                    _basisStatusRows.size(), _basisStatusCols.size());
733  	
734  	   assert(!_simplifier->isUnsimplified());
735  	
736  	   try
737  	   {
738  	      // unsimplify basis and solution data
739  	      _simplifier->unsimplify(_solReal._primal, _solReal._dual, _solReal._slacks, _solReal._redCost,
740  	                              _basisStatusRows.get_ptr(), _basisStatusCols.get_ptr());
741  	
742  	   }
743  	   catch(const SPxException& E)
744  	   {
745  	      SPX_MSG_INFO1(spxout, spxout << "Caught exception <" << E.what() <<
746  	                    "> during unsimplification. Resolving without simplifier and scaler.\n");
747  	      _preprocessAndSolveReal(false);
748  	      return;
749  	   }
750  	
751  	   // copy unsimplified solution data from simplifier
752  	   _solReal._primal  = _simplifier->unsimplifiedPrimal();
753  	   _solReal._slacks  = _simplifier->unsimplifiedSlacks();
754  	   _solReal._dual    = _simplifier->unsimplifiedDual();
755  	   _solReal._redCost = _simplifier->unsimplifiedRedCost();
756  	
757  	   // unscale stored solution (removes persistent scaling)
758  	   if(_isRealLPScaled)
759  	      _unscaleSolutionReal(*_realLP, true);
760  	
761  	   // compute the original objective function value
762  	   StableSum<R> objVal(realParam(SoPlexBase<R>::OBJ_OFFSET));
763  	
764  	   for(int i = 0; i < numCols(); ++i)
765  	      objVal += _solReal._primal[i] * objReal(i);
766  	
767  	   _solReal._objVal = R(objVal);
768  	
769  	   // store the unsimplified basis
770  	   _simplifier->getBasis(_basisStatusRows.get_ptr(), _basisStatusCols.get_ptr(),
771  	                         _basisStatusRows.size(), _basisStatusCols.size());
772  	   _hasBasis = true;
773  	   _hasSolReal = true;
774  	   _solReal._isPrimalFeasible = true;
775  	   _solReal._isDualFeasible = true;
776  	   _solver.setBasisStatus(SPxBasisBase<R>::OPTIMAL);
777  	
778  	   // check solution for violations and solve again if necessary
779  	   _verifySolutionReal();
780  	}
781  	
782  	
783  	
784  	/// load original LP and possibly setup a slack basis
785  	template <class R>
786  	void SoPlexBase<R>::_loadRealLP(bool initBasis)
787  	{
788  	   _solver.loadLP(*_realLP, initBasis);
789  	   _isRealLPLoaded = true;
790  	   _realLP->~SPxLPBase<R>();
791  	   spx_free(_realLP);
792  	   _realLP = &_solver;
793  	
794  	   if(initBasis)
795  	      _solver.init();
796  	}
797  	
798  	
799  	
800  	/// unscales stored solution to remove internal or external scaling of LP
801  	template <class R>
802  	void SoPlexBase<R>::_unscaleSolutionReal(SPxLPBase<R>& LP, bool persistent)
803  	{
804  	   SPX_MSG_INFO1(spxout, spxout << " --- unscaling " << (persistent ? "external" : "internal") <<
805  	                 " solution" << std::endl;)
806  	   assert(_scaler);
807  	   assert(!persistent || (boolParam(SoPlexBase<R>::PERSISTENTSCALING) && _isRealLPScaled));
808  	   _scaler->unscalePrimal(LP, _solReal._primal);
809  	   _scaler->unscaleSlacks(LP, _solReal._slacks);
810  	   _scaler->unscaleDual(LP, _solReal._dual);
811  	   _scaler->unscaleRedCost(LP, _solReal._redCost);
812  	
813  	   if(_solReal.hasPrimalRay())
814  	      _scaler->unscalePrimalray(LP, _solReal._primalRay);
815  	
816  	   if(_solReal.hasDualFarkas())
817  	      _scaler->unscaleDualray(LP, _solReal._dualFarkas);
818  	}
819  	} // namespace soplex
820