Scippy

SoPlex

Sequential object-oriented simPlex

solverational.cpp
Go to the documentation of this file.
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the class library */
4 /* SoPlex --- the Sequential object-oriented simPlex. */
5 /* */
6 /* Copyright (C) 1996-2016 Konrad-Zuse-Zentrum */
7 /* fuer Informationstechnik Berlin */
8 /* */
9 /* SoPlex is distributed under the terms of the ZIB Academic Licence. */
10 /* */
11 /* You should have received a copy of the ZIB Academic License */
12 /* along with SoPlex; see the file COPYING. If not email to soplex@zib.de. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15 
16 #ifndef SOPLEX_LEGACY
17 #include <iostream>
18 #include <assert.h>
19 
20 #include "soplex.h"
21 #include "statistics.h"
22 #include "slufactor_rational.h"
23 #include "ratrecon.h"
24 
25 namespace soplex
26 {
27  /// solves rational LP
29  {
30  bool hasUnboundedRay = false;
31  bool infeasibilityNotCertified = false;
32  bool unboundednessNotCertified = false;
33 
34  // start timing
37 
38  // remember that last solve was rational
40 
41  // ensure that the solver has the original problem
42  if( !_isRealLPLoaded )
43  {
44  assert(_realLP != &_solver);
45 
48  _realLP = &_solver;
49  _isRealLPLoaded = true;
50  }
51  // during the rational solve, we always store basis information in the basis arrays
52  else if( _hasBasis )
53  {
54  _basisStatusRows.reSize(numRowsReal());
55  _basisStatusCols.reSize(numColsReal());
57  }
58 
59  // store objective, bounds, and sides of real LP in case they will be modified during iterative refinement
60  _storeLPReal();
61 
62  // deactivate objective limit in floating-point solver
64  {
65  MSG_INFO2( spxout, spxout << "Deactivating objective limit.\n" );
66  }
67 
69 
71 
72  // apply lifting to reduce range of nonzero matrix coefficients
74  _lift();
75 
76  // force column representation
77  ///@todo implement row objectives with row representation
78  int oldRepresentation = intParam(SoPlex::REPRESENTATION);
80 
81  // force ratio test (avoid bound flipping)
82  int oldRatiotester = intParam(SoPlex::RATIOTESTER);
84 
85  ///@todo implement handling of row objectives in Cplex interface
86 #ifdef SOPLEX_WITH_CPX
87  int oldEqtrans = boolParam(SoPlex::EQTRANS);
89 #endif
90 
91  // introduce slack variables to transform inequality constraints into equations
94 
95  _storedBasis = false;
96 
97  bool stoppedTime;
98  bool stoppedIter;
99  do
100  {
101  bool primalFeasible = false;
102  bool dualFeasible = false;
103  bool infeasible = false;
104  bool unbounded = false;
105  bool error = false;
106  stoppedTime = false;
107  stoppedIter = false;
108 
109  // solve problem with iterative refinement and recovery mechanism
110  _performOptIRStable(_solRational, !unboundednessNotCertified, !infeasibilityNotCertified, 0,
111  primalFeasible, dualFeasible, infeasible, unbounded, stoppedTime, stoppedIter, error);
112 
113  // case: an unrecoverable error occured
114  if( error )
115  {
117  break;
118  }
119  // case: stopped due to some limit
120  else if( stoppedTime )
121  {
123  break;
124  }
125  else if( stoppedIter )
126  {
128  break;
129  }
130  // case: unboundedness detected for the first time
131  else if( unbounded && !unboundednessNotCertified )
132  {
133  SolRational solUnbounded;
134 
135  _performUnboundedIRStable(solUnbounded, hasUnboundedRay, stoppedTime, stoppedIter, error);
136 
137  assert(!hasUnboundedRay || solUnbounded.hasPrimalRay());
138  assert(!solUnbounded.hasPrimalRay() || hasUnboundedRay);
139 
140  if( error )
141  {
142  MSG_INFO1( spxout, spxout << "Error while testing for unboundedness.\n" );
144  break;
145  }
146 
147  if( hasUnboundedRay )
148  {
149  MSG_INFO1( spxout, spxout << "Dual infeasible. Primal unbounded ray available.\n" );
150  }
151  else
152  {
153  MSG_INFO1( spxout, spxout << "Dual feasible. Rejecting primal unboundedness.\n" );
154  }
155 
156  unboundednessNotCertified = !hasUnboundedRay;
157 
158  if( stoppedTime )
159  {
161  break;
162  }
163  else if( stoppedIter )
164  {
166  break;
167  }
168 
169  _performFeasIRStable(_solRational, infeasible, stoppedTime, stoppedIter, error);
170 
171  ///@todo this should be stored already earlier, possible switch use solRational above and solFeas here
172  if( hasUnboundedRay )
173  {
174  _solRational._primalRay = solUnbounded._primalRay;
176  }
177 
178  if( error )
179  {
180  MSG_INFO1( spxout, spxout << "Error while testing for feasibility.\n" );
182  break;
183  }
184  else if( stoppedTime )
185  {
187  break;
188  }
189  else if( stoppedIter )
190  {
192  break;
193  }
194  else if( infeasible )
195  {
196  MSG_INFO1( spxout, spxout << "Primal infeasible. Dual Farkas ray available.\n" );
198  break;
199  }
200  else if( hasUnboundedRay )
201  {
202  MSG_INFO1( spxout, spxout << "Primal feasible and unbounded.\n" );
204  break;
205  }
206  else
207  {
208  MSG_INFO1( spxout, spxout << "Primal feasible and bounded.\n" );
209  continue;
210  }
211  }
212  // case: infeasibility detected
213  else if( infeasible && !infeasibilityNotCertified )
214  {
215  _storeBasis();
216 
217  _performFeasIRStable(_solRational, infeasible, stoppedTime, stoppedIter, error);
218 
219  if( error )
220  {
221  MSG_INFO1( spxout, spxout << "Error while testing for infeasibility.\n" );
223  _restoreBasis();
224  break;
225  }
226 
227  infeasibilityNotCertified = !infeasible;
228 
229  if( stoppedTime )
230  {
232  _restoreBasis();
233  break;
234  }
235  else if( stoppedIter )
236  {
238  _restoreBasis();
239  break;
240  }
241 
242  if( infeasible && boolParam(SoPlex::TESTDUALINF) )
243  {
244  SolRational solUnbounded;
245 
246  _performUnboundedIRStable(solUnbounded, hasUnboundedRay, stoppedTime, stoppedIter, error);
247 
248  assert(!hasUnboundedRay || solUnbounded.hasPrimalRay());
249  assert(!solUnbounded.hasPrimalRay() || hasUnboundedRay);
250 
251  if( error )
252  {
253  MSG_INFO1( spxout, spxout << "Error while testing for dual infeasibility.\n" );
255  _restoreBasis();
256  break;
257  }
258 
259  if( hasUnboundedRay )
260  {
261  MSG_INFO1( spxout, spxout << "Dual infeasible. Primal unbounded ray available.\n" );
262  _solRational._primalRay = solUnbounded._primalRay;
264  }
265  else if( solUnbounded._hasDual )
266  {
267  MSG_INFO1( spxout, spxout << "Dual feasible. Storing dual multipliers.\n" );
268  _solRational._dual = solUnbounded._dual;
269  _solRational._redCost = solUnbounded._redCost;
270  _solRational._hasDual = true;
271  }
272  else
273  {
274  assert(false);
275  MSG_INFO1( spxout, spxout << "Not dual infeasible.\n" );
276  }
277  }
278 
279  _restoreBasis();
280 
281  if( infeasible )
282  {
283  MSG_INFO1( spxout, spxout << "Primal infeasible. Dual Farkas ray available.\n" );
285  break;
286  }
287  else if( hasUnboundedRay )
288  {
289  MSG_INFO1( spxout, spxout << "Primal feasible and unbounded.\n" );
291  break;
292  }
293  else
294  {
295  MSG_INFO1( spxout, spxout << "Primal feasible. Optimizing again.\n" );
296  continue;
297  }
298  }
299  else if( primalFeasible && dualFeasible )
300  {
301  MSG_INFO1( spxout, spxout << "Solved to optimality.\n" );
303  break;
304  }
305  else
306  {
307  MSG_INFO1( spxout, spxout << "Terminating without success.\n" );
308  break;
309  }
310  }
311  while( !_isSolveStopped(stoppedTime, stoppedIter) );
312 
313  ///@todo set status to ABORT_VALUE if optimal solution exceeds objective limit
314 
316  _hasSolRational = true;
317 
318  // restore original problem
321 
322 #ifdef SOPLEX_WITH_CPX
323  setBoolParam(SoPlex::EQTRANS, oldEqtrans);
324 #endif
325 
326  // reset representation and ratio test
327  setIntParam(SoPlex::REPRESENTATION, oldRepresentation);
328  setIntParam(SoPlex::RATIOTESTER, oldRatiotester);
329 
330  // undo lifting
333 
334  // restore objective, bounds, and sides of real LP in case they have been modified during iterative refinement
335  _restoreLPReal();
336 
337  // since the real LP is loaded in the solver, we need to also pass the basis information to the solver if
338  // available
339  if( _hasBasis )
340  {
341  assert(_isRealLPLoaded);
342  _solver.setBasis(_basisStatusRows.get_const_ptr(), _basisStatusCols.get_const_ptr());
344  }
345 
346  // stop timing
348  }
349 
350 
351 
352  /// solves current problem with iterative refinement and recovery mechanism
354  SolRational& sol,
355  bool acceptUnbounded,
356  bool acceptInfeasible,
357  int minRounds,
358  bool& primalFeasible,
359  bool& dualFeasible,
360  bool& infeasible,
361  bool& unbounded,
362  bool& stoppedTime,
363  bool& stoppedIter,
364  bool& error)
365  {
366  // start rational solving timing
368 
369  primalFeasible = false;
370  dualFeasible = false;
371  infeasible = false;
372  unbounded = false;
373  stoppedTime = false;
374  stoppedIter = false;
375  error = false;
376 
377  // set working tolerances in floating-point solver
380 
381  // declare vectors and variables
383 
384  _modLower.reDim(numColsRational(), false);
385  _modUpper.reDim(numColsRational(), false);
386  _modLhs.reDim(numRowsRational(), false);
387  _modRhs.reDim(numRowsRational(), false);
388  _modObj.reDim(numColsRational(), false);
389 
390  DVectorReal primalReal(numColsRational());
391  DVectorReal dualReal(numRowsRational());
392 
393  Rational boundsViolation;
394  Rational sideViolation;
395  Rational redCostViolation;
396  Rational dualViolation;
397  Rational primalScale;
398  Rational dualScale;
399  Rational maxScale;
400 
401  // solve original LP
402  MSG_INFO1( spxout, spxout << "Initial floating-point solve . . .\n" );
403 
404  if( _hasBasis )
405  {
406  assert(_basisStatusRows.size() == numRowsRational());
407  assert(_basisStatusCols.size() == numColsRational());
408  _solver.setBasis(_basisStatusRows.get_const_ptr(), _basisStatusCols.get_const_ptr());
410  }
411 
412  for( int r = numRowsRational() - 1; r >= 0; r-- )
413  {
414  assert(_solver.maxRowObj(r) == 0.0);
415  }
416 
418  result = _solveRealStable(acceptUnbounded, acceptInfeasible, primalReal, dualReal, _basisStatusRows, _basisStatusCols, _hasBasis);
419 
420  // evaluate result
421  switch( result )
422  {
423  case SPxSolver::OPTIMAL:
424  MSG_INFO1( spxout, spxout << "Floating-point optimal.\n" );
425  break;
427  MSG_INFO1( spxout, spxout << "Floating-point infeasible.\n" );
428  // the floating-point solve returns a Farkas ray if and only if the simplifier was not used, which is exactly
429  // the case when a basis could be returned
430  if( _hasBasis )
431  {
432  sol._dualFarkas = dualReal;
433  sol._hasDualFarkas = true;
434  }
435  else
436  sol._hasDualFarkas = false;
437  infeasible = true;
438  return;
440  MSG_INFO1( spxout, spxout << "Floating-point unbounded.\n" );
441  unbounded = true;
442  return;
444  stoppedTime = true;
445  return;
447  stoppedIter = true;
448  return;
449  default:
450  error = true;
451  return;
452  }
453 
455 
456  // store floating-point solution of original LP as current rational solution and ensure that solution vectors have
457  // right dimension; ensure that solution is aligned with basis
458  sol._primal.reDim(numColsRational(), false);
459  sol._slacks.reDim(numRowsRational(), false);
460  sol._dual.reDim(numRowsRational(), false);
461  sol._redCost.reDim(numColsRational(), false);
462  sol._hasPrimal = true;
463  sol._hasDual = true;
464 
465  for( int c = numColsRational() - 1; c >= 0; c-- )
466  {
467  SPxSolver::VarStatus& basisStatusCol = _basisStatusCols[c];
468 
469  if( basisStatusCol == SPxSolver::ON_LOWER )
470  sol._primal[c] = lowerRational(c);
471  else if( basisStatusCol == SPxSolver::ON_UPPER )
472  sol._primal[c] = upperRational(c);
473  else if( basisStatusCol == SPxSolver::FIXED )
474  {
475  // it may happen that lower and upper are only equal in the real LP but different in the rational LP; we do
476  // not check this to avoid rational comparisons, but simply switch the basis status to the lower bound; this
477  // is necessary, because for fixed variables any reduced cost is feasible
478  sol._primal[c] = lowerRational(c);
479  basisStatusCol = SPxSolver::ON_LOWER;
480  }
481  else if( basisStatusCol == SPxSolver::ZERO )
482  sol._primal[c] = 0;
483  else
484  sol._primal[c] = primalReal[c];
485  }
487 
488  int dualSize = 0;
489  for( int r = numRowsRational() - 1; r >= 0; r-- )
490  {
491  SPxSolver::VarStatus& basisStatusRow = _basisStatusRows[r];
492 
493  // it may happen that left-hand and right-hand side are different in the rational, but equal in the real LP,
494  // leading to a fixed basis status; this is critical because rows with fixed basis status are ignored in the
495  // computation of the dual violation; to avoid rational comparisons we do not check this but simply switch to
496  // the left-hand side status
497  if( basisStatusRow == SPxSolver::FIXED )
498  basisStatusRow = SPxSolver::ON_LOWER;
499 
500  {
501  sol._dual[r] = dualReal[r];
502  if( dualReal[r] != 0.0 )
503  dualSize++;
504  }
505  }
506  // we assume that the objective function vector has less nonzeros than the reduced cost vector, and so multiplying
507  // with -1 first and subtracting the dual activity should be faster than adding the dual activity and negating
508  // afterwards
511 
512  // initial scaling factors are one
513  primalScale = Rational::POSONE;
514  dualScale = Rational::POSONE;
515 
516  // control progress
517  Rational maxViolation;
518  Rational bestViolation = _rationalPosInfty;
519  const Rational violationImprovementFactor = 0.9;
520  const Rational errorCorrectionFactor = 1.1;
521  Rational errorCorrection = 2;
522  int numFailedRefinements = 0;
523 
524  // store basis status in case solving modified problem failed
525  DataArray< SPxSolver::VarStatus > basisStatusRowsFirst;
526  DataArray< SPxSolver::VarStatus > basisStatusColsFirst;
527 
528  // refinement loop
529  const bool maximizing = (intParam(SoPlex::OBJSENSE) == SoPlex::OBJSENSE_MAXIMIZE);
530  const int maxDimRational = numColsRational() > numRowsRational() ? numColsRational() : numRowsRational();
531  SolRational factorSol;
532  bool factorSolNewBasis = true;
533  int lastStallRefinements = 0;
534  int nextRatrecRefinement = 0;
535  do
536  {
537  // decrement minRounds counter
538  minRounds--;
539 
540  MSG_DEBUG( std::cout << "Computing primal violations.\n" );
541 
542  // compute violation of bounds
543  boundsViolation = 0;
544  for( int c = numColsRational() - 1; c >= 0; c-- )
545  {
546  // lower bound
548  if( _lowerFinite(_colTypes[c]) )
549  {
550  if( lowerRational(c) == 0 )
551  {
552  _modLower[c] = sol._primal[c];
553  _modLower[c] *= -1;
554  if( _modLower[c] > boundsViolation )
555  boundsViolation = _modLower[c];
556  }
557  else
558  {
559  _modLower[c] = lowerRational(c);
560  _modLower[c] -= sol._primal[c];
561  if( _modLower[c] > boundsViolation )
562  boundsViolation = _modLower[c];
563  }
564  }
565 
566  // upper bound
568  if( _upperFinite(_colTypes[c]) )
569  {
570  if( upperRational(c) == 0 )
571  {
572  _modUpper[c] = sol._primal[c];
573  _modUpper[c] *= -1;
574  if( _modUpper[c] < -boundsViolation )
575  boundsViolation = -_modUpper[c];
576  }
577  else
578  {
579  _modUpper[c] = upperRational(c);
580  _modUpper[c] -= sol._primal[c];
581  if( _modUpper[c] < -boundsViolation )
582  boundsViolation = -_modUpper[c];
583  }
584  }
585  }
586 
587  // compute violation of sides
588  sideViolation = 0;
589  for( int r = numRowsRational() - 1; r >= 0; r-- )
590  {
591  const SPxSolver::VarStatus& basisStatusRow = _basisStatusRows[r];
592 
593  // left-hand side
594  assert((lhsRational(r) > _rationalNegInfty) == _lowerFinite(_rowTypes[r]));
595  if( _lowerFinite(_rowTypes[r]) )
596  {
597  if( lhsRational(r) == 0 )
598  {
599  _modLhs[r] = sol._slacks[r];
600  _modLhs[r] *= -1;
601  }
602  else
603  {
604  _modLhs[r] = lhsRational(r);
605  _modLhs[r] -= sol._slacks[r];
606  }
607 
608  if( _modLhs[r] > sideViolation )
609  sideViolation = _modLhs[r];
610  // if the activity is feasible, but too far from the bound, this violates complementary slackness; we
611  // count it as side violation here
612  else if( basisStatusRow == SPxSolver::ON_LOWER && _modLhs[r] < -sideViolation )
613  sideViolation = -_modLhs[r];
614  }
615 
616  // right-hand side
617  assert((rhsRational(r) < _rationalPosInfty) == _upperFinite(_rowTypes[r]));
618  if( _upperFinite(_rowTypes[r]) )
619  {
620  if( rhsRational(r) == 0 )
621  {
622  _modRhs[r] = sol._slacks[r];
623  _modRhs[r] *= -1;
624  }
625  else
626  {
627  _modRhs[r] = rhsRational(r);
628  _modRhs[r] -= sol._slacks[r];
629  }
630 
631  if( _modRhs[r] < -sideViolation )
632  sideViolation = -_modRhs[r];
633  // if the activity is feasible, but too far from the bound, this violates complementary slackness; we
634  // count it as side violation here
635  else if( basisStatusRow == SPxSolver::ON_UPPER && _modRhs[r] > sideViolation )
636  sideViolation = _modRhs[r];
637  }
638  }
639 
640  MSG_DEBUG( std::cout << "Computing dual violations.\n" );
641 
642  // compute reduced cost violation
643  redCostViolation = 0;
644  for( int c = numColsRational() - 1; c >= 0; c-- )
645  {
646  if( _colTypes[c] == RANGETYPE_FIXED )
647  continue;
648 
649  const SPxSolver::VarStatus& basisStatusCol = _basisStatusCols[c];
650  assert(basisStatusCol != SPxSolver::FIXED);
651 
652  if( ((maximizing && basisStatusCol != SPxSolver::ON_LOWER) || (!maximizing && basisStatusCol != SPxSolver::ON_UPPER))
653  && sol._redCost[c] < -redCostViolation )
654  {
655  MSG_DEBUG( std::cout << "basisStatusCol = " << basisStatusCol
656  << ", lower tight = " << bool(sol._primal[c] <= lowerRational(c))
657  << ", upper tight = " << bool(sol._primal[c] >= upperRational(c))
658  << ", sol._redCost[c] = " << rationalToString(sol._redCost[c])
659  << "\n" );
660  redCostViolation = -sol._redCost[c];
661  }
662 
663  if( ((maximizing && basisStatusCol != SPxSolver::ON_UPPER) || (!maximizing && basisStatusCol != SPxSolver::ON_LOWER))
664  && sol._redCost[c] > redCostViolation )
665  {
666  MSG_DEBUG( std::cout << "basisStatusCol = " << basisStatusCol
667  << ", lower tight = " << bool(sol._primal[c] <= lowerRational(c))
668  << ", upper tight = " << bool(sol._primal[c] >= upperRational(c))
669  << ", sol._redCost[c] = " << rationalToString(sol._redCost[c])
670  << "\n" );
671  redCostViolation = sol._redCost[c];
672  }
673  }
674 
675  // compute dual violation
676  dualViolation = 0;
677  for( int r = numRowsRational() - 1; r >= 0; r-- )
678  {
679  if( _rowTypes[r] == RANGETYPE_FIXED )
680  continue;
681 
682  const SPxSolver::VarStatus& basisStatusRow = _basisStatusRows[r];
683  assert(basisStatusRow != SPxSolver::FIXED);
684 
685  if( ((maximizing && basisStatusRow != SPxSolver::ON_LOWER) || (!maximizing && basisStatusRow != SPxSolver::ON_UPPER))
686  && sol._dual[r] < -dualViolation )
687  {
688  MSG_DEBUG( std::cout << "basisStatusRow = " << basisStatusRow
689  << ", lower tight = " << bool(sol._slacks[r] <= lhsRational(r))
690  << ", upper tight = " << bool(sol._slacks[r] >= rhsRational(r))
691  << ", sol._dual[r] = " << rationalToString(sol._dual[r])
692  << "\n" );
693  dualViolation = -sol._dual[r];
694  }
695 
696  if( ((maximizing && basisStatusRow != SPxSolver::ON_UPPER) || (!maximizing && basisStatusRow != SPxSolver::ON_LOWER))
697  && sol._dual[r] > dualViolation )
698  {
699  MSG_DEBUG( std::cout << "basisStatusRow = " << basisStatusRow
700  << ", lower tight = " << bool(sol._slacks[r] <= lhsRational(r))
701  << ", upper tight = " << bool(sol._slacks[r] >= rhsRational(r))
702  << ", sol._dual[r] = " << rationalToString(sol._dual[r])
703  << "\n" );
704  dualViolation = sol._dual[r];
705  }
706  }
707 
708  _modObj = sol._redCost;
709 
710  // output violations; the reduced cost violations for artificially introduced slack columns are actually violations of the dual multipliers
712  << "Max. bound violation = " << rationalToString(boundsViolation) << "\n"
713  << "Max. row violation = " << rationalToString(sideViolation) << "\n"
714  << "Max. reduced cost violation = " << rationalToString(redCostViolation) << "\n"
715  << "Max. dual violation = " << rationalToString(dualViolation) << "\n" );
716 
718  << std::fixed << std::setprecision(2) << std::setw(10)
719  << "Progress table: "
720  << std::setw(10) << _statistics->refinements << " & "
721  << std::setw(10) << _statistics->iterations << " & "
722  << std::setw(10) << _statistics->solvingTime->time() << " & "
723  << std::setw(10) << _statistics->rationalTime->time() << " & "
724  << std::setw(10) << rationalToString(boundsViolation > sideViolation ? boundsViolation : sideViolation, 2) << " & "
725  << std::setw(10) << rationalToString(redCostViolation > dualViolation ? redCostViolation : dualViolation, 2) << "\n");
726 
727  // terminate if tolerances are satisfied
728  primalFeasible = (boundsViolation <= _rationalFeastol && sideViolation <= _rationalFeastol);
729  dualFeasible = (redCostViolation <= _rationalOpttol && dualViolation <= _rationalOpttol);
730  if( primalFeasible && dualFeasible )
731  {
732  if( minRounds < 0 )
733  {
734  MSG_INFO1( spxout, spxout << "Tolerances reached.\n" );
735  break;
736  }
737  else
738  {
739  MSG_INFO1( spxout, spxout << "Tolerances reached but minRounds forcing additional refinement rounds.\n" );
740  }
741  }
742 
743  // terminate if some limit is reached
744  if( _isSolveStopped(stoppedTime, stoppedIter) )
745  break;
746 
747  // check progress
748  maxViolation = boundsViolation;
749  if( sideViolation > maxViolation )
750  maxViolation = sideViolation;
751  if( redCostViolation > maxViolation )
752  maxViolation = redCostViolation;
753  if( dualViolation > maxViolation )
754  maxViolation = dualViolation;
755  bestViolation *= violationImprovementFactor;
756  if( maxViolation > bestViolation )
757  {
758  MSG_INFO2( spxout, spxout << "Failed to reduce violation significantly.\n" );
759  numFailedRefinements++;
760  }
761  else
762  bestViolation = maxViolation;
763 
764  // decide whether to perform rational reconstruction and/or factorization
765  bool performRatfac = boolParam(SoPlex::RATFAC)
766  && lastStallRefinements >= intParam(SoPlex::RATFAC_MINSTALLS) && _hasBasis && factorSolNewBasis;
767  bool performRatrec = boolParam(SoPlex::RATREC)
768  && (_statistics->refinements >= nextRatrecRefinement || performRatfac);
769 
770  // attempt rational reconstruction
771  errorCorrection *= errorCorrectionFactor;
772  if( performRatrec && maxViolation > 0 )
773  {
774  MSG_INFO1( spxout, spxout << "Performing rational reconstruction . . .\n" );
775 
776  maxViolation *= errorCorrection; // only used for sign check later
777  maxViolation.invert();
779  {
780  MSG_INFO1( spxout, spxout << "Tolerances reached.\n" );
781  primalFeasible = true;
782  dualFeasible = true;
783  break;
784  }
785  nextRatrecRefinement = int(_statistics->refinements * realParam(SoPlex::RATREC_FREQ)) + 1;
786  MSG_DEBUG( spxout << "Next rational reconstruction after refinement " << nextRatrecRefinement << ".\n" );
787  }
788 
789  // solve basis systems exactly
790  if( performRatfac && maxViolation > 0 )
791  {
792  MSG_INFO1( spxout, spxout << "Performing rational factorization . . .\n" );
793 
794  bool optimal;
795  _factorizeColumnRational(sol, _basisStatusRows, _basisStatusCols, stoppedTime, stoppedIter, error, optimal);
796  factorSolNewBasis = false;
797 
798  if( stoppedTime )
799  {
800  MSG_INFO1( spxout, spxout << "Stopped rational factorization.\n" );
801  }
802  else if( error )
803  {
804  // message was already printed; reset error flag and continue without factorization
805  error = false;
806  }
807  else if( optimal )
808  {
809  MSG_INFO1( spxout, spxout << "Tolerances reached.\n" );
810  primalFeasible = true;
811  dualFeasible = true;
812  break;
813  }
814  else if( boolParam(SoPlex::RATFACJUMP) )
815  {
816  MSG_INFO1( spxout, spxout << "Jumping to exact basic solution.\n" );
817  minRounds++;
818  continue;
819  }
820  }
821 
822  // start refinement
823 
824  // compute primal scaling factor; limit increase in scaling by tolerance used in floating point solve
825  maxScale = primalScale;
826  maxScale *= _rationalMaxscaleincr;
827 
828  primalScale = boundsViolation > sideViolation ? boundsViolation : sideViolation;
829  if( primalScale < redCostViolation )
830  primalScale = redCostViolation;
831  assert(primalScale >= 0);
832 
833  if( primalScale > 0 )
834  {
835  primalScale.invert();
836  if( primalScale > maxScale )
837  primalScale = maxScale;
838  }
839  else
840  primalScale = maxScale;
841 
843  primalScale.powRound();
844 
845  // apply scaled bounds
846  if( primalScale <= 1 )
847  {
848  if( primalScale < 1 )
849  primalScale = 1;
850  for( int c = numColsRational() - 1; c >= 0; c-- )
851  {
852  if( _lowerFinite(_colTypes[c]) )
853  {
854  if( _modLower[c] <= _rationalNegInfty )
856  else
858  }
859  if( _upperFinite(_colTypes[c]) )
860  {
861  if( _modUpper[c] >= _rationalPosInfty )
863  else
865  }
866  }
867  }
868  else
869  {
870  MSG_INFO2( spxout, spxout << "Scaling primal by " << rationalToString(primalScale) << ".\n" );
871 
872  for( int c = numColsRational() - 1; c >= 0; c-- )
873  {
874  if( _lowerFinite(_colTypes[c]) )
875  {
876  _modLower[c] *= primalScale;
877  if( _modLower[c] <= _rationalNegInfty )
879  else
881  }
882  if( _upperFinite(_colTypes[c]) )
883  {
884  _modUpper[c] *= primalScale;
885  if( _modUpper[c] >= _rationalPosInfty )
887  else
889  }
890  }
891  }
892 
893  // apply scaled sides
894  assert(primalScale >= 1);
895  if( primalScale == 1 )
896  {
897  for( int r = numRowsRational() - 1; r >= 0; r-- )
898  {
899  if( _lowerFinite(_rowTypes[r]) )
900  {
901  if( _modLhs[r] <= _rationalNegInfty )
903  else
904  _solver.changeLhs(r, Real(_modLhs[r]));
905  }
906  if( _upperFinite(_rowTypes[r]) )
907  {
908  if( _modRhs[r] >= _rationalPosInfty )
910  else
911  _solver.changeRhs(r, Real(_modRhs[r]));
912  }
913  }
914  }
915  else
916  {
917  for( int r = numRowsRational() - 1; r >= 0; r-- )
918  {
919  if( _lowerFinite(_rowTypes[r]) )
920  {
921  _modLhs[r] *= primalScale;
922  if( _modLhs[r] <= _rationalNegInfty )
924  else
925  _solver.changeLhs(r, Real(_modLhs[r]));
926  }
927  if( _upperFinite(_rowTypes[r]) )
928  {
929  _modRhs[r] *= primalScale;
930  if( _modRhs[r] >= _rationalPosInfty )
932  else
933  _solver.changeRhs(r, Real(_modRhs[r]));
934  }
935  }
936  }
937 
938  // compute dual scaling factor; limit increase in scaling by tolerance used in floating point solve
939  maxScale = dualScale;
940  maxScale *= _rationalMaxscaleincr;
941 
942  dualScale = redCostViolation > dualViolation ? redCostViolation : dualViolation;
943  assert(dualScale >= 0);
944 
945  if( dualScale > 0 )
946  {
947  dualScale.invert();
948  if( dualScale > maxScale )
949  dualScale = maxScale;
950  }
951  else
952  dualScale = maxScale;
953 
955  dualScale.powRound();
956 
957  if( dualScale > primalScale )
958  dualScale = primalScale;
959 
960  if( dualScale < 1 )
961  dualScale = 1;
962  else
963  {
964  MSG_INFO2( spxout, spxout << "Scaling dual by " << rationalToString(dualScale) << ".\n" );
965 
966  // perform dual scaling
967  ///@todo remove _modObj and use dualScale * sol._redCost directly
968  _modObj *= dualScale;
969  }
970 
971  // apply scaled objective function
972  for( int c = numColsRational() - 1; c >= 0; c-- )
973  {
974  if( _modObj[c] >= _rationalPosInfty )
976  else if( _modObj[c] <= _rationalNegInfty )
978  else
979  _solver.changeObj(c, Real(_modObj[c]));
980  }
981  for( int r = numRowsRational() - 1; r >= 0; r-- )
982  {
983  Rational newRowObj;
984  if( _rowTypes[r] == RANGETYPE_FIXED )
985  _solver.changeRowObj(r, 0.0);
986  else
987  {
988  newRowObj = sol._dual[r];
989  newRowObj *= dualScale;
990  if( newRowObj >= _rationalPosInfty )
992  else if( newRowObj <= _rationalNegInfty )
994  else
995  _solver.changeRowObj(r, -Real(newRowObj));
996  }
997  }
998 
999  MSG_INFO1( spxout, spxout << "Refined floating-point solve . . .\n" );
1000 
1001  // ensure that artificial slack columns are basic and inequality constraints are nonbasic; otherwise we may end
1002  // up with dual violation on inequality constraints after removing the slack columns; do not change this in the
1003  // floating-point solver, though, because the solver may require its original basis to detect optimality
1004  if( _slackCols.num() > 0 && _hasBasis )
1005  {
1006  int numOrigCols = numColsRational() - _slackCols.num();
1007  assert(_slackCols.num() <= 0 || boolParam(SoPlex::EQTRANS));
1008  for( int i = 0; i < _slackCols.num(); i++ )
1009  {
1010  int row = _slackCols.colVector(i).index(0);
1011  int col = numOrigCols + i;
1012 
1013  assert(row >= 0);
1014  assert(row < numRowsRational());
1015 
1017  {
1018  _basisStatusRows[row] = _basisStatusCols[col];
1020  }
1021  }
1022  }
1023 
1024  // load basis
1026  {
1027  MSG_DEBUG( spxout << "basis (status = " << _solver.basis().status() << ") desc before set:\n"; _solver.basis().desc().dump() );
1028  _solver.setBasis(_basisStatusRows.get_const_ptr(), _basisStatusCols.get_const_ptr());
1029  MSG_DEBUG( spxout << "basis (status = " << _solver.basis().status() << ") desc after set:\n"; _solver.basis().desc().dump() );
1030 
1032  MSG_DEBUG( spxout << "setting basis in solver " << (_hasBasis ? "successful" : "failed") << " (3)\n" );
1033  }
1034 
1035  // solve modified problem
1036  int prevIterations = _statistics->iterations;
1038  result = _solveRealStable(acceptUnbounded, acceptInfeasible, primalReal, dualReal, _basisStatusRows, _basisStatusCols, _hasBasis, primalScale > 1e20 || dualScale > 1e20);
1039 
1040  // count refinements and remember whether we moved to a new basis
1042  if( _statistics->iterations <= prevIterations )
1043  {
1044  lastStallRefinements++;
1046  }
1047  else
1048  {
1049  factorSolNewBasis = true;
1050  lastStallRefinements = 0;
1052  }
1053 
1054  // evaluate result; if modified problem was not solved to optimality, stop refinement
1055  switch( result )
1056  {
1057  case SPxSolver::OPTIMAL:
1058  MSG_INFO1( spxout, spxout << "Floating-point optimal.\n" );
1059  break;
1060  case SPxSolver::INFEASIBLE:
1061  MSG_INFO1( spxout, spxout << "Floating-point infeasible.\n" );
1062  sol._dualFarkas = dualReal;
1063  sol._hasDualFarkas = true;
1064  infeasible = true;
1066  return;
1067  case SPxSolver::UNBOUNDED:
1068  MSG_INFO1( spxout, spxout << "Floating-point unbounded.\n" );
1069  unbounded = true;
1071  return;
1072  case SPxSolver::ABORT_TIME:
1073  stoppedTime = true;
1074  return;
1075  case SPxSolver::ABORT_ITER:
1076  stoppedIter = true;
1078  return;
1079  default:
1080  error = true;
1082  return;
1083  }
1084 
1086 
1087  // correct primal solution and align with basis
1088  MSG_DEBUG( std::cout << "Correcting primal solution.\n" );
1089 
1090  int primalSize = 0;
1091  Rational primalScaleInverse = primalScale;
1092  primalScaleInverse.invert();
1094  for( int c = numColsRational() - 1; c >= 0; c-- )
1095  {
1096  // force values of nonbasic variables to bounds
1097  SPxSolver::VarStatus& basisStatusCol = _basisStatusCols[c];
1098 
1099  if( basisStatusCol == SPxSolver::ON_LOWER )
1100  {
1101  if( sol._primal[c] != lowerRational(c) )
1102  {
1103  int i = _primalDualDiff.size();
1105  _primalDualDiff.add(c);
1107  _primalDualDiff.value(i) -= sol._primal[c];
1108  sol._primal[c] = lowerRational(c);
1109  }
1110  }
1111  else if( basisStatusCol == SPxSolver::ON_UPPER )
1112  {
1113  if( sol._primal[c] != upperRational(c) )
1114  {
1115  int i = _primalDualDiff.size();
1117  _primalDualDiff.add(c);
1119  _primalDualDiff.value(i) -= sol._primal[c];
1120  sol._primal[c] = upperRational(c);
1121  }
1122  }
1123  else if( basisStatusCol == SPxSolver::FIXED )
1124  {
1125  // it may happen that lower and upper are only equal in the real LP but different in the rational LP; we
1126  // do not check this to avoid rational comparisons, but simply switch the basis status to the lower
1127  // bound; this is necessary, because for fixed variables any reduced cost is feasible
1128  basisStatusCol = SPxSolver::ON_LOWER;
1129  if( sol._primal[c] != lowerRational(c) )
1130  {
1131  int i = _primalDualDiff.size();
1133  _primalDualDiff.add(c);
1135  _primalDualDiff.value(i) -= sol._primal[c];
1136  sol._primal[c] = lowerRational(c);
1137  }
1138  }
1139  else if( basisStatusCol == SPxSolver::ZERO )
1140  {
1141  if( sol._primal[c] != 0 )
1142  {
1143  int i = _primalDualDiff.size();
1145  _primalDualDiff.add(c);
1146  _primalDualDiff.value(i) = sol._primal[c];
1147  _primalDualDiff.value(i) *= -1;
1148  sol._primal[c] = 0;
1149  }
1150  }
1151  else
1152  {
1153  if( primalReal[c] == 1.0 )
1154  {
1155  int i = _primalDualDiff.size();
1157  _primalDualDiff.add(c);
1158  _primalDualDiff.value(i) = primalScaleInverse;
1159  sol._primal[c] += _primalDualDiff.value(i);
1160  }
1161  else if( primalReal[c] == -1.0 )
1162  {
1163  int i = _primalDualDiff.size();
1165  _primalDualDiff.add(c);
1166  _primalDualDiff.value(i) = primalScaleInverse;
1167  _primalDualDiff.value(i) *= -1;
1168  sol._primal[c] += _primalDualDiff.value(i);
1169  }
1170  else if( primalReal[c] != 0.0 )
1171  {
1172  int i = _primalDualDiff.size();
1174  _primalDualDiff.add(c);
1175  _primalDualDiff.value(i) = primalReal[c];
1176  _primalDualDiff.value(i) *= primalScaleInverse;
1177  sol._primal[c] += _primalDualDiff.value(i);
1178  }
1179  }
1180 
1181  if( sol._primal[c] != 0 )
1182  primalSize++;
1183  }
1184 
1185  // update or recompute slacks depending on which looks faster
1186  if( _primalDualDiff.size() < primalSize )
1187  {
1189 #ifndef NDEBUG
1190  {
1191  DVectorRational activity(numRowsRational());
1192  _rationalLP->computePrimalActivity(sol._primal, activity);
1193  assert(sol._slacks == activity);
1194  }
1195 #endif
1196  }
1197  else
1199  const int numCorrectedPrimals = _primalDualDiff.size();
1200 
1201  // correct dual solution and align with basis
1202  MSG_DEBUG( std::cout << "Correcting dual solution.\n" );
1203 
1204 #ifndef NDEBUG
1205  {
1206  // compute reduced cost violation
1207  DVectorRational debugRedCost(numColsRational());
1208  debugRedCost = DVectorRational(_realLP->maxObj());
1209  debugRedCost *= -1;
1210  _rationalLP->subDualActivity(DVectorRational(dualReal), debugRedCost);
1211 
1212  Rational debugRedCostViolation = 0;
1213  for( int c = numColsRational() - 1; c >= 0; c-- )
1214  {
1215  if( _colTypes[c] == RANGETYPE_FIXED )
1216  continue;
1217 
1218  const SPxSolver::VarStatus& basisStatusCol = _basisStatusCols[c];
1219  assert(basisStatusCol != SPxSolver::FIXED);
1220 
1221  if( ((maximizing && basisStatusCol != SPxSolver::ON_LOWER) || (!maximizing && basisStatusCol != SPxSolver::ON_UPPER))
1222  && debugRedCost[c] < -debugRedCostViolation )
1223  {
1224  MSG_DEBUG( std::cout << "basisStatusCol = " << basisStatusCol
1225  << ", lower tight = " << bool(sol._primal[c] <= lowerRational(c))
1226  << ", upper tight = " << bool(sol._primal[c] >= upperRational(c))
1227  << ", obj[c] = " << _realLP->obj(c)
1228  << ", debugRedCost[c] = " << rationalToString(debugRedCost[c])
1229  << "\n" );
1230  debugRedCostViolation = -debugRedCost[c];
1231  }
1232 
1233  if( ((maximizing && basisStatusCol != SPxSolver::ON_UPPER) || (!maximizing && basisStatusCol != SPxSolver::ON_LOWER))
1234  && debugRedCost[c] > debugRedCostViolation )
1235  {
1236  MSG_DEBUG( std::cout << "basisStatusCol = " << basisStatusCol
1237  << ", lower tight = " << bool(sol._primal[c] <= lowerRational(c))
1238  << ", upper tight = " << bool(sol._primal[c] >= upperRational(c))
1239  << ", obj[c] = " << _realLP->obj(c)
1240  << ", debugRedCost[c] = " << rationalToString(debugRedCost[c])
1241  << "\n" );
1242  debugRedCostViolation = debugRedCost[c];
1243  }
1244  }
1245 
1246  // compute dual violation
1247  Rational debugDualViolation = 0;
1248  Rational debugBasicDualViolation = 0;
1249  for( int r = numRowsRational() - 1; r >= 0; r-- )
1250  {
1251  if( _rowTypes[r] == RANGETYPE_FIXED )
1252  continue;
1253 
1254  const SPxSolver::VarStatus& basisStatusRow = _basisStatusRows[r];
1255  assert(basisStatusRow != SPxSolver::FIXED);
1256 
1257  Rational val = (-dualScale * sol._dual[r]) - Rational(dualReal[r]);
1258 
1259  if( ((maximizing && basisStatusRow != SPxSolver::ON_LOWER) || (!maximizing && basisStatusRow != SPxSolver::ON_UPPER))
1260  && val > debugDualViolation )
1261  {
1262  MSG_DEBUG( std::cout << "basisStatusRow = " << basisStatusRow
1263  << ", lower tight = " << bool(sol._slacks[r] <= lhsRational(r))
1264  << ", upper tight = " << bool(sol._slacks[r] >= rhsRational(r))
1265  << ", dualReal[r] = " << rationalToString(val)
1266  << ", dualReal[r] = " << dualReal[r]
1267  << "\n" );
1268  debugDualViolation = val;
1269  }
1270 
1271  if( ((maximizing && basisStatusRow != SPxSolver::ON_UPPER) || (!maximizing && basisStatusRow != SPxSolver::ON_LOWER))
1272  && val < -debugDualViolation )
1273  {
1274  MSG_DEBUG( std::cout << "basisStatusRow = " << basisStatusRow
1275  << ", lower tight = " << bool(sol._slacks[r] <= lhsRational(r))
1276  << ", upper tight = " << bool(sol._slacks[r] >= rhsRational(r))
1277  << ", dualReal[r] = " << rationalToString(val)
1278  << ", dualReal[r] = " << dualReal[r]
1279  << "\n" );
1280  debugDualViolation = -val;
1281  }
1282 
1283  if( basisStatusRow == SPxSolver::BASIC && spxAbs(val) > debugBasicDualViolation )
1284  {
1285  MSG_DEBUG( std::cout << "basisStatusRow = " << basisStatusRow
1286  << ", lower tight = " << bool(sol._slacks[r] <= lhsRational(r))
1287  << ", upper tight = " << bool(sol._slacks[r] >= rhsRational(r))
1288  << ", dualReal[r] = " << rationalToString(val)
1289  << ", dualReal[r] = " << dualReal[r]
1290  << "\n" );
1291  debugBasicDualViolation = spxAbs(val);
1292  }
1293  }
1294 
1295  if( debugRedCostViolation > _solver.opttol() || debugDualViolation > _solver.opttol() || debugBasicDualViolation > 1e-9 )
1296  {
1297  MSG_WARNING( spxout, spxout << "Warning: floating-point dual solution with violation "
1298  << rationalToString(debugRedCostViolation) << " / "
1299  << rationalToString(debugDualViolation) << " / "
1300  << rationalToString(debugBasicDualViolation)
1301  << " (red. cost, dual, basic).\n" );
1302  }
1303  }
1304 #endif
1305 
1306  Rational dualScaleInverseNeg = dualScale;
1307  dualScaleInverseNeg.invert();
1308  dualScaleInverseNeg *= -1;
1310  dualSize = 0;
1311  for( int r = numRowsRational() - 1; r >= 0; r-- )
1312  {
1313  SPxSolver::VarStatus& basisStatusRow = _basisStatusRows[r];
1314 
1315  // it may happen that left-hand and right-hand side are different in the rational, but equal in the real LP,
1316  // leading to a fixed basis status; this is critical because rows with fixed basis status are ignored in the
1317  // computation of the dual violation; to avoid rational comparisons we do not check this but simply switch
1318  // to the left-hand side status
1319  if( basisStatusRow == SPxSolver::FIXED )
1320  basisStatusRow = SPxSolver::ON_LOWER;
1321 
1322  {
1323  if( dualReal[r] != 0 )
1324  {
1325  int i = _primalDualDiff.size();
1327  _primalDualDiff.add(r);
1328  _primalDualDiff.value(i) = dualReal[r];
1329  _primalDualDiff.value(i) *= dualScaleInverseNeg;
1330  sol._dual[r] -= _primalDualDiff.value(i);
1331 
1332  dualSize++;
1333  }
1334  else
1335  {
1336  // we do not check whether the dual value is nonzero, because it probably is; this gives us an
1337  // overestimation of the number of nonzeros in the dual solution
1338  dualSize++;
1339  }
1340  }
1341  }
1342 
1343  // update or recompute reduced cost values depending on which looks faster; adding one to the length of the
1344  // dual vector accounts for the objective function vector
1345  if( _primalDualDiff.size() < dualSize + 1 )
1346  {
1348 #ifndef NDEBUG
1349  {
1350  DVectorRational activity(_rationalLP->maxObj());
1351  activity *= -1;
1352  _rationalLP->subDualActivity(sol._dual, activity);
1353  }
1354 #endif
1355  }
1356  else
1357  {
1358  // we assume that the objective function vector has less nonzeros than the reduced cost vector, and so multiplying
1359  // with -1 first and subtracting the dual activity should be faster than adding the dual activity and negating
1360  // afterwards
1361  _rationalLP->getObj(sol._redCost);
1363  }
1364  const int numCorrectedDuals = _primalDualDiff.size();
1365 
1366  if( numCorrectedPrimals + numCorrectedDuals > 0 )
1367  {
1368  MSG_INFO2( spxout, spxout << "Corrected " << numCorrectedPrimals << " primal variables and " << numCorrectedDuals << " dual values.\n" );
1369  }
1370  }
1371  while( true );
1372 
1373  // correct basis status for restricted inequalities
1374  if( _hasBasis )
1375  {
1376  for( int r = numRowsRational() - 1; r >= 0; r-- )
1377  {
1378  assert((lhsRational(r) == rhsRational(r)) == (_rowTypes[r] == RANGETYPE_FIXED));
1380  _basisStatusRows[r] = (maximizing == (sol._dual[r] < 0))
1383  }
1384  }
1385 
1386  // compute objective function values
1387  assert(sol._hasPrimal == sol._hasDual);
1388  if( sol._hasPrimal )
1389  {
1390  sol._primalObjVal = sol._primal * _rationalLP->maxObj();
1392  sol._primalObjVal *= -1;
1393  sol._dualObjVal = sol._primalObjVal;
1394  }
1395 
1396  // set objective coefficients for all rows to zero
1398 
1399  // stop rational solving time
1401  }
1402 
1403 
1404 
1405  /// performs iterative refinement on the auxiliary problem for testing unboundedness
1407  SolRational& sol,
1408  bool& hasUnboundedRay,
1409  bool& stopped,
1410  bool& stoppedIter,
1411  bool& error)
1412  {
1413  bool primalFeasible;
1414  bool dualFeasible;
1415  bool infeasible;
1416  bool unbounded;
1417 
1418  // move objective function to constraints and adjust sides and bounds
1420 
1421  // invalidate solution
1422  sol.invalidate();
1423 
1424  // remember current number of refinements
1425  int oldRefinements = _statistics->refinements;
1426 
1427  // perform iterative refinement
1428  _performOptIRStable(sol, false, false, 0, primalFeasible, dualFeasible, infeasible, unbounded, stopped, stoppedIter, error);
1429 
1430  // update unbounded refinement counter
1431  _statistics->unbdRefinements += _statistics->refinements - oldRefinements;
1432 
1433  // stopped due to some limit
1434  if( stopped )
1435  {
1436  sol.invalidate();
1437  hasUnboundedRay = false;
1438  error = false;
1439  }
1440  // the unbounded problem should always be solved to optimality
1441  else if( error || unbounded || infeasible || !primalFeasible || !dualFeasible )
1442  {
1443  sol.invalidate();
1444  hasUnboundedRay = false;
1445  stopped = false;
1446  error = true;
1447  }
1448  else
1449  {
1450  const Rational& tau = sol._primal[numColsRational() - 1];
1451 
1452  MSG_DEBUG( std::cout << "tau = " << tau << " (roughly " << rationalToString(tau) << ")\n" );
1453 
1454  assert(tau <= 1.0 + 2.0 * realParam(SoPlex::FEASTOL));
1455  assert(tau >= -realParam(SoPlex::FEASTOL));
1456 
1457  // because the right-hand side and all bounds (but tau's upper bound) are zero, tau should be approximately
1458  // zero if basic; otherwise at its upper bound 1
1459  error = !(tau >= Rational::POSONE || tau <= _rationalFeastol);
1460  assert(!error);
1461 
1462  hasUnboundedRay = (tau >= 1);
1463  }
1464 
1465  // restore problem
1466  _untransformUnbounded(sol, hasUnboundedRay);
1467  }
1468 
1469 
1470 
1471  /// performs iterative refinement on the auxiliary problem for testing feasibility
1473  SolRational& sol,
1474  bool& withDualFarkas,
1475  bool& stopped,
1476  bool& stoppedIter,
1477  bool& error)
1478  {
1479  bool primalFeasible;
1480  bool dualFeasible;
1481  bool infeasible;
1482  bool unbounded;
1483  bool success = false;
1484  error = false;
1485 
1486 #if 0
1487  // if the problem has been found to be infeasible and an approximate Farkas proof is available, we compute a
1488  // scaled unit box around the origin that provably contains no feasible solution; this currently only works for
1489  // equality form
1490  ///@todo check whether approximate Farkas proof can be used
1492  ///@todo if approx Farkas proof is good enough then exit without doing any transformation
1493 #endif
1494 
1495  // remove objective function, shift, homogenize
1497 
1498  // invalidate solution
1499  sol.invalidate();
1500 
1501  do
1502  {
1503  // remember current number of refinements
1504  int oldRefinements = _statistics->refinements;
1505 
1506  // perform iterative refinement
1507  _performOptIRStable(sol, false, false, 0, primalFeasible, dualFeasible, infeasible, unbounded, stopped, stoppedIter, error);
1508 
1509  // update feasible refinement counter
1510  _statistics->feasRefinements += _statistics->refinements - oldRefinements;
1511 
1512  // stopped due to some limit
1513  if( stopped )
1514  {
1515  sol.invalidate();
1516  withDualFarkas = false;
1517  error = false;
1518  }
1519  // the feasibility problem should always be solved to optimality
1520  else if( error || unbounded || infeasible || !primalFeasible || !dualFeasible )
1521  {
1522  sol.invalidate();
1523  withDualFarkas = false;
1524  stopped = false;
1525  error = true;
1526  }
1527  // else we should have either a refined Farkas proof or an approximate feasible solution to the original
1528  else
1529  {
1530  const Rational& tau = sol._primal[numColsRational() - 1];
1531 
1532  MSG_DEBUG( std::cout << "tau = " << tau << " (roughly " << rationalToString(tau) << ")\n" );
1533 
1534  assert(tau >= -realParam(SoPlex::FEASTOL));
1535  assert(tau <= 1.0 + realParam(SoPlex::FEASTOL));
1536 
1537  error = (tau < -_rationalFeastol || tau > Rational::POSONE + _rationalFeastol);
1538  withDualFarkas = (tau < Rational::POSONE);
1539 
1540  if( withDualFarkas )
1541  {
1544 
1545 #if 0
1546  // check if we can compute sufficiently large Farkas box
1548 #endif
1549  if( true )//@todo check if computeInfeasBox found a sufficient box
1550  {
1551 
1552  success = true;
1553  sol._hasPrimal = false;
1554  }
1555  }
1556  else
1557  {
1558  sol._hasDual = false;
1559  success = true; //successfully found approximate feasible solution
1560  }
1561  }
1562  }
1563  while(!error && !success && !stopped);
1564 
1565  // restore problem
1566  _untransformFeasibility(sol, withDualFarkas);
1567  }
1568 
1569 
1570 
1571  /// reduces matrix coefficient in absolute value by the lifting procedure of Thiele et al. 2013
1573  {
1574  MSG_DEBUG( std::cout << "Reducing matrix coefficients by lifting.\n" );
1575 
1576  // start timing
1578 
1579  MSG_DEBUG( _realLP->writeFile("beforeLift.lp", 0, 0, 0) );
1580 
1581  // remember unlifted state
1584 
1585  // allocate vector memory
1586  DSVectorRational colVector;
1587  SVectorRational::Element liftingRowMem[2];
1588  SVectorRational liftingRowVector(2, liftingRowMem);
1589 
1590  // search each column for large nonzeros entries
1591  const Rational maxValue = realParam(SoPlex::LIFTMAXVAL);
1592 
1593  for( int i = 0; i < numColsRational(); i++ )
1594  {
1595  MSG_DEBUG( std::cout << "in lifting: examining column " << i << "\n" );
1596 
1597  // get column vector
1598  colVector = colVectorRational(i);
1599 
1600  bool addedLiftingRow = false;
1601  int liftingColumnIndex = -1;
1602 
1603  // go through nonzero entries of the column
1604  for( int k = colVector.size() - 1; k >= 0; k-- )
1605  {
1606  const Rational& value = colVector.value(k);
1607 
1608  if( spxAbs(value) > maxValue )
1609  {
1610  MSG_DEBUG( std::cout << " --> nonzero " << k << " has value " << rationalToString(value) << " in row " << colVector.index(k) << "\n" );
1611 
1612  // add new column equal to maxValue times original column
1613  if( !addedLiftingRow )
1614  {
1615  MSG_DEBUG( std::cout << " --> adding lifting row\n" );
1616 
1617  assert(liftingRowVector.size() == 0);
1618 
1619  liftingColumnIndex = numColsRational();
1620  liftingRowVector.add(i, maxValue);
1621  liftingRowVector.add(liftingColumnIndex, -1);
1622 
1623  _rationalLP->addRow(LPRowRational(0, liftingRowVector, 0));
1624  _realLP->addRow(LPRowReal(0.0, DSVectorReal(liftingRowVector), 0.0));
1625 
1626  assert(liftingColumnIndex == numColsRational() - 1);
1627  assert(liftingColumnIndex == numColsReal() - 1);
1628 
1631 
1632  liftingRowVector.clear();
1633  addedLiftingRow = true;
1634  }
1635 
1636  // get row index
1637  int rowIndex = colVector.index(k);
1638  assert(rowIndex >= 0);
1639  assert(rowIndex < _beforeLiftRows);
1640  assert(liftingColumnIndex == numColsRational() - 1);
1641 
1642  MSG_DEBUG( std::cout << " --> changing matrix\n" );
1643 
1644  // remove nonzero from original column
1645  _rationalLP->changeElement(rowIndex, i, 0);
1646  _realLP->changeElement(rowIndex, i, 0.0);
1647 
1648  // add nonzero divided by maxValue to new column
1649  Rational newValue(value);
1650  newValue /= maxValue;
1651  _rationalLP->changeElement(rowIndex, liftingColumnIndex, newValue);
1652  _realLP->changeElement(rowIndex, liftingColumnIndex, Real(newValue));
1653  }
1654  }
1655  }
1656 
1657  // search each column for small nonzeros entries
1658  const Rational minValue = realParam(SoPlex::LIFTMINVAL);
1659 
1660  for( int i = 0; i < numColsRational(); i++ )
1661  {
1662  MSG_DEBUG( std::cout << "in lifting: examining column " << i << "\n" );
1663 
1664  // get column vector
1665  colVector = colVectorRational(i);
1666 
1667  bool addedLiftingRow = false;
1668  int liftingColumnIndex = -1;
1669 
1670  // go through nonzero entries of the column
1671  for( int k = colVector.size() - 1; k >= 0; k-- )
1672  {
1673  const Rational& value = colVector.value(k);
1674 
1675  if( spxAbs(value) < minValue )
1676  {
1677  MSG_DEBUG( std::cout << " --> nonzero " << k << " has value " << rationalToString(value) << " in row " << colVector.index(k) << "\n" );
1678 
1679  // add new column equal to maxValue times original column
1680  if( !addedLiftingRow )
1681  {
1682  MSG_DEBUG( std::cout << " --> adding lifting row\n" );
1683 
1684  assert(liftingRowVector.size() == 0);
1685 
1686  liftingColumnIndex = numColsRational();
1687  liftingRowVector.add(i, minValue);
1688  liftingRowVector.add(liftingColumnIndex, -1);
1689 
1690  _rationalLP->addRow(LPRowRational(0, liftingRowVector, 0));
1691  _realLP->addRow(LPRowReal(0.0, DSVectorReal(liftingRowVector), 0.0));
1692 
1693  assert(liftingColumnIndex == numColsRational() - 1);
1694  assert(liftingColumnIndex == numColsReal() - 1);
1695 
1698 
1699  liftingRowVector.clear();
1700  addedLiftingRow = true;
1701  }
1702 
1703  // get row index
1704  int rowIndex = colVector.index(k);
1705  assert(rowIndex >= 0);
1706  assert(rowIndex < _beforeLiftRows);
1707  assert(liftingColumnIndex == numColsRational() - 1);
1708 
1709  MSG_DEBUG( std::cout << " --> changing matrix\n" );
1710 
1711  // remove nonzero from original column
1712  _rationalLP->changeElement(rowIndex, i, 0);
1713  _realLP->changeElement(rowIndex, i, 0.0);
1714 
1715  // add nonzero divided by maxValue to new column
1716  Rational newValue(value);
1717  newValue /= minValue;
1718  _rationalLP->changeElement(rowIndex, liftingColumnIndex, newValue);
1719  _realLP->changeElement(rowIndex, liftingColumnIndex, Real(newValue));
1720  }
1721  }
1722  }
1723 
1724  // adjust basis
1725  if( _hasBasis )
1726  {
1727  assert(numColsRational() >= _beforeLiftCols);
1728  assert(numRowsRational() >= _beforeLiftRows);
1729 
1730  _basisStatusCols.append(numColsRational() - _beforeLiftCols, SPxSolver::BASIC);
1732  }
1733 
1734  MSG_DEBUG( _realLP->writeFile("afterLift.lp", 0, 0, 0) );
1735 
1736  // stop timing
1738 
1739  if( numColsRational() > _beforeLiftCols || numRowsRational() > _beforeLiftRows )
1740  {
1741  MSG_INFO1( spxout, spxout << "Added " << numColsRational() - _beforeLiftCols << " columns and "
1742  << numRowsRational() - _beforeLiftRows << " rows to reduce large matrix coefficients\n." );
1743  }
1744  }
1745 
1746 
1747 
1748  /// undoes lifting
1750  {
1751  // start timing
1753 
1754  // print LP if in debug mode
1755  MSG_DEBUG( _realLP->writeFile("beforeProject.lp", 0, 0, 0) );
1756 
1757  assert(numColsRational() >= _beforeLiftCols);
1758  assert(numRowsRational() >= _beforeLiftRows);
1759 
1760  // shrink rational LP to original size
1763 
1764  // shrink real LP to original size
1767 
1768  // adjust solution
1769  if( sol.hasPrimal() )
1770  {
1773  }
1774 
1775  if( sol.hasPrimalRay() )
1776  {
1778  }
1779 
1780  ///@todo if we know the mapping between original and lifting columns, we simply need to add the reduced cost of
1781  /// the lifting column to the reduced cost of the original column; this is not implemented now, because for
1782  /// optimal solutions the reduced costs of the lifting columns are zero
1783  const Rational maxValue = realParam(SoPlex::LIFTMAXVAL);
1784 
1785  for( int i = _beforeLiftCols; i < numColsRational() && sol._hasDual; i++ )
1786  {
1787  if( spxAbs(maxValue * sol._redCost[i]) > _rationalOpttol )
1788  {
1789  MSG_INFO1( spxout, spxout << "Warning: lost dual solution during project phase.\n" );
1790  sol._hasDual = false;
1791  }
1792  }
1793 
1794  if( sol.hasDual() )
1795  {
1798  }
1799 
1800  if( sol.hasDualFarkas() )
1801  {
1803  }
1804 
1805  // adjust basis
1806  for( int i = _beforeLiftCols; i < numColsRational() && _hasBasis; i++ )
1807  {
1809  {
1810  MSG_INFO1( spxout, spxout << "Warning: lost basis during project phase because of nonbasic lifting column.\n" );
1811  _hasBasis = false;
1812  }
1813  }
1814 
1815  for( int i = _beforeLiftRows; i < numRowsRational() && _hasBasis; i++ )
1816  {
1818  {
1819  MSG_INFO1( spxout, spxout << "Warning: lost basis during project phase because of basic lifting row.\n" );
1820  _hasBasis = false;
1821  }
1822  }
1823 
1824  if( _hasBasis )
1825  {
1828  }
1829 
1830  // print LP if in debug mode
1831  MSG_DEBUG( _realLP->writeFile("afterProject.lp", 0, 0, 0) );
1832 
1833  // stop timing
1835  }
1836 
1837 
1838 
1839  /// stores objective, bounds, and sides of real LP
1841  {
1842 #ifndef SOPLEX_MANUAL_ALT
1844  {
1846  return;
1847  }
1848 #endif
1849 
1850  _manualLower = _realLP->lower();
1851  _manualUpper = _realLP->upper();
1852  _manualLhs = _realLP->lhs();
1853  _manualRhs = _realLP->rhs();
1856  }
1857 
1858 
1859 
1860  /// restores objective, bounds, and sides of real LP
1862  {
1864  {
1865 #ifndef SOPLEX_MANUAL_ALT
1867 #else
1873 #endif
1874 
1875  if( _hasBasis )
1876  {
1877  // in manual sync mode, if the right-hand side of an equality constraint is not floating-point
1878  // representable, the user might have constructed the constraint in the real LP by rounding down the
1879  // left-hand side and rounding up the right-hand side; if the basis status is fixed, we need to adjust it
1880  for( int i = 0; i < _solver.nRows(); i++ )
1881  {
1882  if( _basisStatusRows[i] == SPxSolver::FIXED && _solver.lhs(i) != _solver.rhs(i) )
1883  {
1884  assert(_solver.rhs(i) == spxNextafter(_solver.lhs(i), infinity));
1888  {
1890  }
1891  else
1892  {
1894  }
1895  }
1896  }
1897 
1898  _solver.setBasis(_basisStatusRows.get_const_ptr(), _basisStatusCols.get_const_ptr());
1900  }
1901  }
1902  else
1903  {
1909  }
1910  }
1911 
1912 
1913 
1914  /// introduces slack variables to transform inequality constraints into equations for both rational and real LP,
1915  /// which should be in sync
1917  {
1918  MSG_DEBUG( std::cout << "Transforming rows to equation form.\n" );
1919 
1920  // start timing
1922 
1923  MSG_DEBUG( _realLP->writeFile("beforeTransEqu.lp", 0, 0, 0) );
1924 
1925  // clear array of slack columns
1926  _slackCols.clear();
1927 
1928  // add artificial slack variables to convert inequality to equality constraints
1929  for( int i = 0; i < numRowsRational(); i++ )
1930  {
1931  assert((lhsRational(i) == rhsRational(i)) == (_rowTypes[i] == RANGETYPE_FIXED));
1932  if( _rowTypes[i] != RANGETYPE_FIXED )
1933  {
1935  if( _rationalLP->lhs(i) != 0 )
1937  if( _rationalLP->rhs(i) != 0 )
1939  assert(_rationalLP->lhs(i) == 0);
1940  assert(_rationalLP->rhs(i) == 0);
1941  _realLP->changeRange(i, 0.0, 0.0);
1944  }
1945  }
1946 
1949 
1950  // adjust basis
1951  if( _hasBasis )
1952  {
1953  for( int i = 0; i < _slackCols.num(); i++ )
1954  {
1955  int row = _slackCols.colVector(i).index(0);
1956 
1957  assert(row >= 0);
1958  assert(row < numRowsRational());
1959 
1960  switch( _basisStatusRows[row] )
1961  {
1962  case SPxSolver::ON_LOWER:
1964  break;
1965  case SPxSolver::ON_UPPER:
1967  break;
1968  case SPxSolver::BASIC:
1969  case SPxSolver::FIXED:
1970  default:
1971  _basisStatusCols.append(_basisStatusRows[row]);
1972  break;
1973  }
1974 
1976  }
1977  }
1978 
1979  MSG_DEBUG( _realLP->writeFile("afterTransEqu.lp", 0, 0, 0) );
1980 
1981  // stop timing
1983 
1984  if( _slackCols.num() > 0 )
1985  {
1986  MSG_INFO1( spxout, spxout << "Added " << _slackCols.num() << " slack columns to transform rows to equality form.\n" );
1987  }
1988  }
1989 
1990 
1991 
1992  /// restores original problem
1994  {
1995  // start timing
1997 
1998  // print LP if in debug mode
1999  MSG_DEBUG( _realLP->writeFile("beforeUntransEqu.lp", 0, 0, 0) );
2000 
2001  int numCols = numColsRational();
2002  int numOrigCols = numColsRational() - _slackCols.num();
2003 
2004  // adjust solution
2005  if( sol.hasPrimal() )
2006  {
2007  for( int i = 0; i < _slackCols.num(); i++ )
2008  {
2009  int col = numOrigCols + i;
2010  int row = _slackCols.colVector(i).index(0);
2011 
2012  assert(row >= 0);
2013  assert(row < numRowsRational());
2014 
2015  sol._slacks[row] -= sol._primal[col];
2016  }
2017 
2018  sol._primal.reDim(numOrigCols);
2019  }
2020 
2021  if( sol.hasPrimalRay() )
2022  {
2023  sol._primalRay.reDim(numOrigCols);
2024  }
2025 
2026  // adjust basis
2027  if( _hasBasis )
2028  {
2029  for( int i = 0; i < _slackCols.num(); i++ )
2030  {
2031  int col = numOrigCols + i;
2032  int row = _slackCols.colVector(i).index(0);
2033 
2034  assert(row >= 0);
2035  assert(row < numRowsRational());
2036  assert(_basisStatusRows[row] != SPxSolver::UNDEFINED);
2037  assert(_basisStatusRows[row] != SPxSolver::ZERO || lhsRational(row) == 0);
2038  assert(_basisStatusRows[row] != SPxSolver::ZERO || rhsRational(row) == 0);
2040 
2041  MSG_DEBUG( std::cout << "slack column " << col << " for row " << row
2042  << ": col status=" << _basisStatusCols[col]
2043  << ", row status=" << _basisStatusRows[row]
2044  << ", redcost=" << rationalToString(sol._redCost[col])
2045  << ", dual=" << rationalToString(sol._dual[row]) << "\n" );
2046 
2047  if( _basisStatusRows[row] != SPxSolver::BASIC )
2048  {
2049  switch( _basisStatusCols[col] )
2050  {
2051  case SPxSolver::ON_LOWER:
2053  break;
2054  case SPxSolver::ON_UPPER:
2056  break;
2057  case SPxSolver::BASIC:
2058  case SPxSolver::FIXED:
2059  default:
2060  _basisStatusRows[row] = _basisStatusCols[col];
2061  break;
2062  }
2063  }
2064  }
2065 
2066  _basisStatusCols.reSize(numOrigCols);
2067  }
2068 
2069  // not earlier because of debug message
2070  if( sol.hasDual() )
2071  {
2072  sol._redCost.reDim(numOrigCols);
2073  }
2074 
2075  // restore sides and remove slack columns
2076  for( int i = 0; i < _slackCols.num(); i++ )
2077  {
2078  int col = numOrigCols + i;
2079  int row = _slackCols.colVector(i).index(0);
2080 
2081  if( upperRational(col) != 0 )
2082  _rationalLP->changeLhs(row, -upperRational(col));
2083  if( lowerRational(col) != 0 )
2084  _rationalLP->changeRhs(row, -lowerRational(col));
2085  assert(_rationalLP->lhs(row) == -upperRational(col));
2086  assert(_rationalLP->rhs(row) == -lowerRational(col));
2087  _rowTypes[row] = _switchRangeType(_colTypes[col]);
2088  }
2089 
2090  _rationalLP->removeColRange(numOrigCols, numCols - 1);
2091  _realLP->removeColRange(numOrigCols, numCols - 1);
2092  _colTypes.reSize(numOrigCols);
2093 
2094  // objective, bounds, and sides of real LP are restored only after _solveRational()
2095 
2096  // print LP if in debug mode
2097  MSG_DEBUG( _realLP->writeFile("afterUntransEqu.lp", 0, 0, 0) );
2098 
2099  // stop timing
2101  }
2102 
2103 
2104 
2105  /// transforms LP to unboundedness problem by moving the objective function to the constraints, changing right-hand
2106  /// side and bounds to zero, and adding an auxiliary variable for the decrease in the objective function
2108  {
2109  MSG_INFO1( spxout, spxout << "Setting up LP to compute primal unbounded ray.\n" );
2110 
2111  // start timing
2113 
2114  // print LP if in debug mode
2115  MSG_DEBUG( _realLP->writeFile("beforeTransUnbounded.lp", 0, 0, 0) );
2116 
2117  // store bounds
2120  for( int c = numColsRational() - 1; c >= 0; c-- )
2121  {
2122  if( _lowerFinite(_colTypes[c]) )
2124  if( _upperFinite(_colTypes[c]) )
2126  }
2127 
2128  // store sides
2131  for( int r = numRowsRational() - 1; r >= 0; r-- )
2132  {
2133  if( _lowerFinite(_rowTypes[r]) )
2134  _unboundedLhs[r] = lhsRational(r);
2135  if( _upperFinite(_rowTypes[r]) )
2136  _unboundedRhs[r] = rhsRational(r);
2137  }
2138 
2139  // make right-hand side zero
2140  for( int r = numRowsRational() - 1; r >= 0; r-- )
2141  {
2142  assert((lhsRational(r) > _rationalNegInfty) == _lowerFinite(_rowTypes[r]));
2143  if( _lowerFinite(_rowTypes[r]) )
2144  {
2145  _rationalLP->changeLhs(r, 0);
2146  _realLP->changeLhs(r, 0.0);
2147  }
2148  else if( _realLP->lhs(r) > -realParam(SoPlex::INFTY) )
2150 
2151  assert((rhsRational(r) < _rationalPosInfty) == _upperFinite(_rowTypes[r]));
2152  if( _upperFinite(_rowTypes[r]) )
2153  {
2154  _rationalLP->changeRhs(r, 0);
2155  _realLP->changeRhs(r, 0.0);
2156  }
2157  else if( _realLP->rhs(r) < realParam(SoPlex::INFTY) )
2159  }
2160 
2161  // transform objective function to constraint and add auxiliary variable
2162  int numOrigCols = numColsRational();
2163  DSVectorRational obj(numOrigCols + 1);
2164  ///@todo implement this without copying the objective function
2165  obj = _rationalLP->maxObj();
2166  obj.add(numOrigCols, -1);
2167  _rationalLP->addRow(LPRowRational(0, obj, 0));
2168  _realLP->addRow(LPRowReal(0, DSVectorReal(obj), 0));
2170 
2171  assert(numColsRational() == numOrigCols + 1);
2172 
2173  // set objective coefficient and bounds for auxiliary variable
2174  _rationalLP->changeMaxObj(numOrigCols, 1);
2175  _realLP->changeMaxObj(numOrigCols, 1.0);
2176 
2177  _rationalLP->changeBounds(numOrigCols, _rationalNegInfty, 1);
2178  _realLP->changeBounds(numOrigCols, -realParam(SoPlex::INFTY), 1.0);
2180 
2181  // set objective coefficients to zero and adjust bounds for problem variables
2182  for( int c = numColsRational() - 2; c >= 0; c-- )
2183  {
2184  _rationalLP->changeObj(c, 0);
2185  _realLP->changeObj(c, 0.0);
2186 
2187  assert((lowerRational(c) > _rationalNegInfty) == _lowerFinite(_colTypes[c]));
2188  if( _lowerFinite(_colTypes[c]) )
2189  {
2190  _rationalLP->changeLower(c, 0);
2191  _realLP->changeLower(c, 0.0);
2192  }
2193  else if( _realLP->lower(c) > -realParam(SoPlex::INFTY) )
2195 
2196  assert((upperRational(c) < _rationalPosInfty) == _upperFinite(_colTypes[c]));
2197  if( _upperFinite(_colTypes[c]) )
2198  {
2199  _rationalLP->changeUpper(c, 0);
2200  _realLP->changeUpper(c, 0.0);
2201  }
2202  else if( _realLP->upper(c) < realParam(SoPlex::INFTY) )
2204  }
2205 
2206  // adjust basis
2207  if( _hasBasis )
2208  {
2211  }
2212 
2213  // print LP if in debug mode
2214  MSG_DEBUG( _realLP->writeFile("afterTransUnbounded.lp", 0, 0, 0) );
2215 
2216  // stop timing
2218  }
2219 
2220 
2221 
2222  /// undoes transformation to unboundedness problem
2223  void SoPlex::_untransformUnbounded(SolRational& sol, bool unbounded)
2224  {
2225  // start timing
2227 
2228  // print LP if in debug mode
2229  MSG_DEBUG( _realLP->writeFile("beforeUntransUnbounded.lp", 0, 0, 0) );
2230 
2231  int numOrigCols = numColsRational() - 1;
2232  int numOrigRows = numRowsRational() - 1;
2233  const Rational& tau = sol._primal[numOrigCols];
2234 
2235  // adjust solution and basis
2236  if( unbounded )
2237  {
2238  assert(tau >= Rational::POSONE);
2239 
2240  sol._hasPrimal = false;
2241  sol._hasPrimalRay = true;
2242  sol._hasDual = false;
2243  sol._hasDualFarkas = false;
2244 
2245  if( tau != 1 )
2246  sol._primal /= tau;
2247 
2248  sol._primalRay = sol._primal;
2249  sol._primalRay.reDim(numOrigCols);
2250 
2251  _hasBasis = (_basisStatusCols[numOrigCols] != SPxSolver::BASIC && _basisStatusRows[numOrigRows] == SPxSolver::BASIC);
2252  _basisStatusCols.reSize(numOrigCols);
2253  _basisStatusRows.reSize(numOrigRows);
2254  }
2255  else if( boolParam(SoPlex::TESTDUALINF) && tau < _rationalFeastol )
2256  {
2257  const Rational& alpha = sol._dual[numOrigRows];
2258 
2259  assert(sol._hasDual);
2260  assert(alpha <= _rationalFeastol - Rational::POSONE);
2261 
2262  sol._hasPrimal = false;
2263  sol._hasPrimalRay = false;
2264  sol._hasDualFarkas = false;
2265 
2266  if( alpha != -1 )
2267  {
2268  sol._dual /= -alpha;
2269  sol._redCost /= -alpha;
2270  }
2271  sol._dual.reDim(numOrigRows);
2272  sol._redCost.reDim(numOrigCols);
2273  }
2274  else
2275  {
2276  sol.invalidate();
2277  _hasBasis = false;
2278  _basisStatusCols.reSize(numOrigCols);
2279  _basisStatusCols.reSize(numOrigRows);
2280  }
2281 
2282  // recover objective function
2283  const SVectorRational& objRowVector = _rationalLP->rowVector(numOrigRows);
2284  for( int i = objRowVector.size() - 1; i >= 0; i-- )
2285  {
2286  _rationalLP->changeMaxObj(objRowVector.index(i), objRowVector.value(i));
2287  _realLP->changeMaxObj(objRowVector.index(i), Real(objRowVector.value(i)));
2288  }
2289 
2290  // remove objective function constraint and auxiliary variable
2291  _rationalLP->removeRow(numOrigRows);
2292  _realLP->removeRow(numOrigRows);
2293  _rowTypes.reSize(numOrigRows);
2294 
2295  _rationalLP->removeCol(numOrigCols);
2296  _realLP->removeCol(numOrigCols);
2297  _colTypes.reSize(numOrigCols);
2298 
2299  // restore objective, sides and bounds
2300  for( int r = numRowsRational() - 1; r >= 0; r-- )
2301  {
2302  if( _lowerFinite(_rowTypes[r]) )
2303  {
2306  }
2307  if( _upperFinite(_rowTypes[r]) )
2308  {
2311  }
2312  assert((lhsRational(r) > _rationalNegInfty) == _lowerFinite(_rowTypes[r]));
2313  assert((rhsRational(r) < _rationalPosInfty) == _upperFinite(_rowTypes[r]));
2314  assert((lhsReal(r) > -realParam(SoPlex::INFTY)) == _lowerFinite(_rowTypes[r]));
2315  assert((rhsReal(r) < realParam(SoPlex::INFTY)) == _upperFinite(_rowTypes[r]));
2316  }
2317  for( int c = numColsRational() - 1; c >= 0; c-- )
2318  {
2319  if( _lowerFinite(_colTypes[c]) )
2320  {
2323  }
2324  if( _upperFinite(_colTypes[c]) )
2325  {
2328  }
2329  assert((lowerRational(c) > _rationalNegInfty) == _lowerFinite(_colTypes[c]));
2330  assert((upperRational(c) < _rationalPosInfty) == _upperFinite(_colTypes[c]));
2331  assert((lowerReal(c) > -realParam(SoPlex::INFTY)) == _lowerFinite(_colTypes[c]));
2332  assert((upperReal(c) < realParam(SoPlex::INFTY)) == _upperFinite(_colTypes[c]));
2333  }
2334 
2335  // print LP if in debug mode
2336  MSG_DEBUG( _realLP->writeFile("afterUntransUnbounded.lp", 0, 0, 0) );
2337 
2338  // stop timing
2340  }
2341 
2342 
2343 
2344  /// store basis
2346  {
2347  assert(!_storedBasis);
2348 
2349  if( _hasBasis )
2350  {
2351  _storedBasis = true;
2354  }
2355  else
2356  _storedBasis = false;
2357  }
2358 
2359 
2360 
2361  /// restore basis
2363  {
2364  if( _storedBasis )
2365  {
2366  _hasBasis = true;
2369  _storedBasis = false;
2370  }
2371  }
2372 
2373 
2374 
2375  /// transforms LP to feasibility problem by removing the objective function, shifting variables, and homogenizing the
2376  /// right-hand side
2378  {
2379  MSG_INFO1( spxout, spxout << "Setting up LP to test for feasibility.\n" );
2380 
2381  // start timing
2383 
2384  // print LP if in debug mode
2385  MSG_DEBUG( _realLP->writeFile("beforeTransFeas.lp", 0, 0, 0) );
2386 
2387  // store objective function
2389  for( int c = numColsRational() - 1; c >= 0; c-- )
2390  _feasObj[c] = _rationalLP->maxObj(c);
2391 
2392  // store bounds
2395  for( int c = numColsRational() - 1; c >= 0; c-- )
2396  {
2397  if( _lowerFinite(_colTypes[c]) )
2398  _feasLower[c] = lowerRational(c);
2399  if( _upperFinite(_colTypes[c]) )
2400  _feasUpper[c] = upperRational(c);
2401  }
2402 
2403  // store sides
2406  for( int r = numRowsRational() - 1; r >= 0; r-- )
2407  {
2408  if( _lowerFinite(_rowTypes[r]) )
2409  _feasLhs[r] = lhsRational(r);
2410  if( _upperFinite(_rowTypes[r]) )
2411  _feasRhs[r] = rhsRational(r);
2412  }
2413 
2414  // set objective coefficients to zero; shift primal space such as to guarantee that the zero solution is within
2415  // the bounds
2416  Rational shiftValue;
2417  Rational shiftValue2;
2418  for( int c = numColsRational() - 1; c >= 0; c-- )
2419  {
2420  _rationalLP->changeMaxObj(c, 0);
2421  _realLP->changeMaxObj(c, 0.0);
2422 
2423  if( lowerRational(c) > 0 )
2424  {
2425  const SVectorRational& colVector = colVectorRational(c);
2426 
2427  for( int i = 0; i < colVector.size(); i++ )
2428  {
2429  shiftValue = colVector.value(i);
2430  shiftValue *= lowerRational(c);
2431  int r = colVector.index(i);
2432 
2433  assert((lhsRational(r) > _rationalNegInfty) == _lowerFinite(_rowTypes[r]));
2434  assert((rhsRational(r) < _rationalPosInfty) == _upperFinite(_rowTypes[r]));
2435 
2436  if( _lowerFinite(_rowTypes[r]) && _upperFinite(_rowTypes[r]) )
2437  {
2438  shiftValue2 = lhsRational(r);
2439  shiftValue2 -= shiftValue;
2440  _rationalLP->changeLhs(r, shiftValue2);
2441  _realLP->changeLhs(r, Real(shiftValue2));
2442 
2443  shiftValue -= rhsRational(r);
2444  shiftValue *= -1;
2445  _rationalLP->changeRhs(r, shiftValue);
2446  _realLP->changeRhs(r, Real(shiftValue));
2447  }
2448  else if( _lowerFinite(_rowTypes[r]) )
2449  {
2450  shiftValue -= lhsRational(r);
2451  shiftValue *= -1;
2452  _rationalLP->changeLhs(r, shiftValue);
2453  _realLP->changeLhs(r, Real(shiftValue));
2454  }
2455  else if( _upperFinite(_rowTypes[r]) )
2456  {
2457  shiftValue -= rhsRational(r);
2458  shiftValue *= -1;
2459  _rationalLP->changeRhs(r, shiftValue);
2460  _realLP->changeRhs(r, Real(shiftValue));
2461  }
2462  }
2463 
2464  assert((upperRational(c) < _rationalPosInfty) == _upperFinite(_colTypes[c]));
2465  if( _upperFinite(_colTypes[c]) )
2466  {
2468  _realLP->changeBounds(c, 0.0, Real(upperRational(c)));
2469  }
2470  else if( _realLP->upper(c) < realParam(SoPlex::INFTY) )
2471  {
2472  _rationalLP->changeLower(c, 0);
2474  }
2475  else
2476  {
2477  _rationalLP->changeLower(c, 0);
2478  _realLP->changeLower(c, 0.0);
2479  }
2480  }
2481  else if( upperRational(c) < 0 )
2482  {
2483  const SVectorRational& colVector = colVectorRational(c);
2484 
2485  for( int i = 0; i < colVector.size(); i++ )
2486  {
2487  shiftValue = colVector.value(i);
2488  shiftValue *= upperRational(c);
2489  int r = colVector.index(i);
2490 
2491  assert((lhsRational(r) > _rationalNegInfty) == _lowerFinite(_rowTypes[r]));
2492  assert((rhsRational(r) < _rationalPosInfty) == _upperFinite(_rowTypes[r]));
2493 
2494  if( _lowerFinite(_rowTypes[r]) && _upperFinite(_rowTypes[r]) )
2495  {
2496  shiftValue2 = lhsRational(r);
2497  shiftValue2 -= shiftValue;
2498  _rationalLP->changeLhs(r, shiftValue2);
2499  _realLP->changeLhs(r, Real(shiftValue2));
2500 
2501  shiftValue -= rhsRational(r);
2502  shiftValue *= -1;
2503  _rationalLP->changeRhs(r, shiftValue);
2504  _realLP->changeRhs(r, Real(shiftValue));
2505  }
2506  else if( _lowerFinite(_rowTypes[r]) )
2507  {
2508  shiftValue -= lhsRational(r);
2509  shiftValue *= -1;
2510  _rationalLP->changeLhs(r, shiftValue);
2511  _realLP->changeLhs(r, Real(shiftValue));
2512  }
2513  else if( _upperFinite(_rowTypes[r]) )
2514  {
2515  shiftValue -= rhsRational(r);
2516  shiftValue *= -1;
2517  _rationalLP->changeRhs(r, shiftValue);
2518  _realLP->changeRhs(r, Real(shiftValue));
2519  }
2520  }
2521 
2522  assert((lowerRational(c) > _rationalNegInfty) == _lowerFinite(_colTypes[c]));
2523  if( _lowerFinite(_colTypes[c]) )
2524  {
2526  _realLP->changeBounds(c, Real(lowerRational(c)), 0.0);
2527  }
2528  else if( _realLP->lower(c) > -realParam(SoPlex::INFTY) )
2529  {
2530  _rationalLP->changeUpper(c, 0);
2532  }
2533  else
2534  {
2535  _rationalLP->changeUpper(c, 0);
2536  _realLP->changeUpper(c, 0.0);
2537  }
2538  }
2539  else
2540  {
2541  if( _lowerFinite(_colTypes[c]) )
2543  else if( _realLP->lower(c) > -realParam(SoPlex::INFTY) )
2545  if( _upperFinite(_colTypes[c]) )
2547  else if( _realLP->upper(c) < realParam(SoPlex::INFTY) )
2549  }
2550 
2551  assert(lowerReal(c) <= upperReal(c));
2552  }
2553 
2554  // homogenize sides
2555  _tauColVector.clear();
2556  for( int r = numRowsRational() - 1; r >= 0; r-- )
2557  {
2558  if( lhsRational(r) > 0 )
2559  {
2560  _tauColVector.add(r, lhsRational(r));
2561  assert((rhsRational(r) < _rationalPosInfty) == _upperFinite(_rowTypes[r]));
2562  if( _upperFinite(_rowTypes[r]) )
2563  {
2565  _realLP->changeRange(r, 0.0, Real(rhsRational(r)));
2566  }
2567  else
2568  {
2569  _rationalLP->changeLhs(r, 0);
2570  _realLP->changeLhs(r, 0.0);
2571  if( _realLP->rhs(r) < realParam(SoPlex::INFTY) )
2573  }
2574  }
2575  else if( rhsRational(r) < 0 )
2576  {
2577  _tauColVector.add(r, rhsRational(r));
2578  assert((lhsRational(r) > _rationalNegInfty) == _lowerFinite(_rowTypes[r]));
2579  if( _lowerFinite(_rowTypes[r]) )
2580  {
2582  _realLP->changeRange(r, Real(lhsRational(r)), 0.0);
2583  }
2584  else
2585  {
2586  _rationalLP->changeRhs(r, 0);
2587  _realLP->changeRhs(r, 0.0);
2588  if( _realLP->lhs(r) > -realParam(SoPlex::INFTY) )
2590  }
2591  }
2592  else
2593  {
2594  if( _lowerFinite(_rowTypes[r]) )
2595  _realLP->changeLhs(r, Real(lhsRational(r)));
2596  else if( _realLP->lhs(r) > -realParam(SoPlex::INFTY) )
2598  if( _upperFinite(_rowTypes[r]) )
2599  _realLP->changeRhs(r, Real(rhsRational(r)));
2600  else if( _realLP->rhs(r) < realParam(SoPlex::INFTY) )
2602  }
2603 
2604  assert(rhsReal(r) <= rhsReal(r));
2605  }
2606 
2607  ///@todo exploit this case by returning without LP solving
2608  if( _tauColVector.size() == 0 )
2609  {
2610  MSG_INFO3( spxout, spxout << "LP is trivially feasible.\n" );
2611  }
2612 
2613  // add artificial column
2614  SPxColId id;
2615  _tauColVector *= -1;
2616  _rationalLP->addCol(id,
2618  _tauColVector, 1, 0));
2619  _realLP->addCol(id,
2621  DSVectorReal(_tauColVector), 1.0, 0.0));
2623 
2624  // adjust basis
2625  if( _hasBasis )
2626  {
2628  }
2629 
2630  // print LP if in debug mode
2631  MSG_DEBUG( _realLP->writeFile("afterTransFeas.lp", 0, 0, 0) );
2632 
2633  // stop timing
2635  }
2636 
2637 
2638 
2639  /// undoes transformation to feasibility problem
2640  void SoPlex::_untransformFeasibility(SolRational& sol, bool infeasible)
2641  {
2642  // start timing
2644 
2645  // print LP if in debug mode
2646  MSG_DEBUG( _realLP->writeFile("beforeUntransFeas.lp", 0, 0, 0) );
2647 
2648  int numOrigCols = numColsRational() - 1;
2649 
2650  // adjust solution and basis
2651  if( infeasible )
2652  {
2653  assert(sol._hasDual);
2654  assert(sol._primal[numOrigCols] < 1);
2655 
2656  sol._hasPrimal = false;
2657  sol._hasPrimalRay = false;
2658  sol._hasDual = false;
2659  sol._hasDualFarkas = true;
2660 
2661  sol._dualFarkas = sol._dual;
2662 
2663  _hasBasis = false;
2664  _basisStatusCols.reSize(numOrigCols);
2665  }
2666  else if( sol._hasPrimal )
2667  {
2668  assert(sol._primal[numOrigCols] >= 1);
2669 
2670  sol._hasPrimalRay = false;
2671  sol._hasDual = false;
2672  sol._hasDualFarkas = false;
2673 
2674  if( sol._primal[numOrigCols] != 1 )
2675  {
2676  sol._slacks /= sol._primal[numOrigCols];
2677  for( int i = 0; i < numOrigCols; i++ )
2678  sol._primal[i] /= sol._primal[numOrigCols];
2679  sol._primal[numOrigCols] = 1;
2680  }
2681 
2682  sol._primal.reDim(numOrigCols);
2683  sol._slacks -= _rationalLP->colVector(numOrigCols);
2684 
2685  _hasBasis = (_basisStatusCols[numOrigCols] != SPxSolver::BASIC);
2686  _basisStatusCols.reSize(numOrigCols);
2687  }
2688  else
2689  {
2690  _hasBasis = false;
2691  _basisStatusCols.reSize(numOrigCols);
2692  }
2693 
2694  // restore right-hand side
2695  for( int r = numRowsRational() - 1; r >= 0; r-- )
2696  {
2698  || _feasLhs[r] - lhsRational(r) == _feasRhs[r] - rhsRational(r));
2699 
2700  if( _lowerFinite(_rowTypes[r]) )
2701  {
2702  _rationalLP->changeLhs(r, _feasLhs[r]);
2703  _realLP->changeLhs(r, Real(_feasLhs[r]));
2704  }
2705  else if( _realLP->lhs(r) > -realParam(SoPlex::INFTY) )
2707  assert(_lowerFinite(_rowTypes[r]) == (lhsRational(r) > _rationalNegInfty));
2708  assert(_lowerFinite(_rowTypes[r]) == (lhsReal(r) > -realParam(SoPlex::INFTY)));
2709 
2710  if( _upperFinite(_rowTypes[r]) )
2711  {
2712  _rationalLP->changeRhs(r, _feasRhs[r]);
2713  _realLP->changeRhs(r, Real(_feasRhs[r]));
2714  }
2715  else if( _realLP->rhs(r) < realParam(SoPlex::INFTY) )
2717  assert(_upperFinite(_rowTypes[r]) == (rhsRational(r) < _rationalPosInfty));
2718  assert(_upperFinite(_rowTypes[r]) == (rhsReal(r) < realParam(SoPlex::INFTY)));
2719 
2720  assert(lhsReal(r) <= rhsReal(r));
2721  }
2722 
2723  // unshift primal space and restore objective coefficients
2724  Rational shiftValue;
2725  for( int c = numOrigCols - 1; c >= 0; c-- )
2726  {
2727  bool shifted = (_lowerFinite(_colTypes[c]) && _feasLower[c] > 0) || (_upperFinite(_colTypes[c]) && _feasUpper[c] < 0);
2728  assert(shifted || !_lowerFinite(_colTypes[c]) || _feasLower[c] == lowerRational(c));
2729  assert(shifted || !_upperFinite(_colTypes[c]) || _feasUpper[c] == upperRational(c));
2731  || _feasLower[c] - lowerRational(c) == _feasUpper[c] - upperRational(c));
2732 
2733  if( shifted )
2734  {
2735  if( _lowerFinite(_colTypes[c]) )
2736  {
2737  shiftValue = _feasLower[c];
2738  shiftValue -= lowerRational(c);
2739  }
2740  else if( _upperFinite(_colTypes[c]) )
2741  {
2742  shiftValue = _feasUpper[c];
2743  shiftValue -= upperRational(c);
2744  }
2745  if( sol._hasPrimal )
2746  {
2747  sol._primal[c] += shiftValue;
2748  sol._slacks.multAdd(shiftValue, _rationalLP->colVector(c));
2749  }
2750  }
2751 
2752  if( _lowerFinite(_colTypes[c]) )
2753  {
2754  if( shifted )
2757  }
2758  else if( _realLP->lower(c) > -realParam(SoPlex::INFTY) )
2760  assert(_lowerFinite(_colTypes[c]) == (lowerRational(c) > -_rationalPosInfty));
2761  assert(_lowerFinite(_colTypes[c]) == (lowerReal(c) > -realParam(SoPlex::INFTY)));
2762 
2763  if( _upperFinite(_colTypes[c]) )
2764  {
2765  if( shifted )
2768  }
2769  else if( _realLP->upper(c) < realParam(SoPlex::INFTY) )
2771  assert(_upperFinite(_colTypes[c]) == (upperRational(c) < _rationalPosInfty));
2772  assert(_upperFinite(_colTypes[c]) == (upperReal(c) < realParam(SoPlex::INFTY)));
2773 
2775  _realLP->changeMaxObj(c, Real(_feasObj[c]));
2776 
2777  assert(lowerReal(c) <= upperReal(c));
2778  }
2779 
2780  // remove last column
2781  _rationalLP->removeCol(numOrigCols);
2782  _realLP->removeCol(numOrigCols);
2783  _colTypes.reSize(numOrigCols);
2784 
2785  // print LP if in debug mode
2786  MSG_DEBUG( _realLP->writeFile("afterUntransFeas.lp", 0, 0, 0) );
2787 
2788  // stop timing
2790 
2791 #ifndef NDEBUG
2792  if( sol._hasPrimal )
2793  {
2794  DVectorRational activity(numRowsRational());
2795  _rationalLP->computePrimalActivity(sol._primal, activity);
2796  assert(sol._slacks == activity);
2797  }
2798 #endif
2799  }
2800 
2801  /** computes radius of infeasibility box implied by an approximate Farkas' proof
2802 
2803  Given constraints of the form \f$ lhs <= Ax <= rhs \f$, a farkas proof y should satisfy \f$ y^T A = 0 \f$ and
2804  \f$ y_+^T lhs - y_-^T rhs > 0 \f$, where \f$ y_+, y_- \f$ denote the positive and negative parts of \f$ y \f$.
2805  If \f$ y \f$ is approximate, it may not satisfy \f$ y^T A = 0 \f$ exactly, but the proof is still valid as long
2806  as the following holds for all potentially feasible \f$ x \f$:
2807 
2808  \f[
2809  y^T Ax < (y_+^T lhs - y_-^T rhs) (*)
2810  \f]
2811 
2812  we may therefore calculate \f$ y^T A \f$ and \f$ y_+^T lhs - y_-^T rhs \f$ exactly and check if the upper and lower
2813  bounds on \f$ x \f$ imply that all feasible \f$ x \f$ satisfy (*), and if not then compute bounds on \f$ x \f$ to
2814  guarantee (*). The simplest way to do this is to compute
2815 
2816  \f[
2817  B = (y_+^T lhs - y_-^T rhs) / \sum_i(|(y^T A)_i|)
2818  \f]
2819 
2820  noting that if every component of \f$ x \f$ has \f$ |x_i| < B \f$, then (*) holds.
2821 
2822  \f$ B \f$ can be increased by iteratively including variable bounds smaller than \f$ B \f$. The speed of this
2823  method can be further improved by using interval arithmetic for all computations. For related information see
2824  Sec. 4 of Neumaier and Shcherbina, Mathematical Programming A, 2004.
2825 
2826  Set transformed to true if this method is called after _transformFeasibility().
2827  */
2828  void SoPlex::_computeInfeasBox(SolRational& sol, bool transformed)
2829  {
2830  assert(sol.hasDualFarkas());
2831 
2832  const VectorRational& lower = transformed ? _feasLower : lowerRational();
2833  const VectorRational& upper = transformed ? _feasUpper : upperRational();
2834  const VectorRational& lhs = transformed ? _feasLhs : lhsRational();
2835  const VectorRational& rhs = transformed ? _feasRhs : rhsRational();
2836  const VectorRational& y = sol._dualFarkas;
2837 
2838  const int numRows = numRowsRational();
2839  const int numCols = transformed ? numColsRational() - 1 : numColsRational();
2840 
2841  SSVectorRational ytransA(numColsRational());
2842  Rational ytransb;
2843  Rational temp;
2844 
2845  // prepare ytransA and ytransb; since we want exact arithmetic, we set the zero threshold of the semi-sparse
2846  // vector to zero
2847  ytransA.setEpsilon(0);
2848  ytransA.clear();
2849  ytransb = 0;
2850 
2851  ///@todo this currently works only if all constraints are equations aggregate rows and sides using the multipliers of the Farkas ray
2852  for( int r = 0; r < numRows; r++ )
2853  {
2854  ytransA += y[r] * _rationalLP->rowVector(r);
2855  ytransb += y[r] * (y[r] > 0 ? lhs[r] : rhs[r]);
2856  }
2857 
2858  // if we work on the feasibility problem, we ignore the last column
2859  if( transformed )
2860  ytransA.reDim(numCols);
2861 
2862  MSG_DEBUG( std::cout << "ytransb = " << rationalToString(ytransb) << "\n" );
2863 
2864  // if we choose minus ytransb as vector of multipliers for the bound constraints on the variables, we obtain an
2865  // exactly feasible dual solution for the LP with zero objective function; we aggregate the bounds of the
2866  // variables accordingly and store its negation in temp
2867  temp = 0;
2868  bool isTempFinite = true;
2869  for( int c = 0; c < numCols && isTempFinite; c++ )
2870  {
2871  const Rational& minusRedCost = ytransA[c];
2872 
2873  if( minusRedCost > 0 )
2874  {
2875  assert((upper[c] < _rationalPosInfty) == _upperFinite(_colTypes[c]));
2876  if( _upperFinite(_colTypes[c]) )
2877  temp.addProduct(minusRedCost, upper[c]);
2878  else
2879  isTempFinite = false;
2880  }
2881  else if( minusRedCost < 0 )
2882  {
2883  assert((lower[c] > _rationalNegInfty) == _lowerFinite(_colTypes[c]));
2884  if( _lowerFinite(_colTypes[c]) )
2885  temp.addProduct(minusRedCost, lower[c]);
2886  else
2887  isTempFinite = false;
2888  }
2889  }
2890 
2891  MSG_DEBUG( std::cout << "max ytransA*[x_l,x_u] = " << (isTempFinite ? rationalToString(temp) : "infinite") << "\n" );
2892 
2893  // ytransb - temp is the increase in the dual objective along the Farkas ray; if this is positive, the dual is
2894  // unbounded and certifies primal infeasibility
2895  if( isTempFinite && temp < ytransb )
2896  {
2897  MSG_INFO1( spxout, spxout << "Farkas infeasibility proof verified exactly. (1)\n" );
2898  return;
2899  }
2900 
2901  // ensure that array of nonzero elements in ytransA is available
2902  assert(ytransA.isSetup());
2903  ytransA.setup();
2904 
2905  // if ytransb is negative, try to make it zero by including a positive lower bound or a negative upper bound
2906  if( ytransb < 0 )
2907  {
2908  for( int c = 0; c < numCols; c++ )
2909  {
2910  if( lower[c] > 0 )
2911  {
2912  ytransA.setValue(c, ytransA[c] - ytransb / lower[c]);
2913  ytransb = 0;
2914  break;
2915  }
2916  else if( upper[c] < 0 )
2917  {
2918  ytransA.setValue(c, ytransA[c] - ytransb / upper[c]);
2919  ytransb = 0;
2920  break;
2921  }
2922  }
2923  }
2924 
2925  // if ytransb is still zero then the zero solution is inside the bounds and cannot be cut off by the Farkas
2926  // constraint; in this case, we cannot compute a Farkas box
2927  if( ytransb < 0 )
2928  {
2929  MSG_INFO1( spxout, spxout << "Approximate Farkas proof to weak. Could not compute Farkas box. (1)\n" );
2930  return;
2931  }
2932 
2933  // compute the one norm of ytransA
2934  temp = 0;
2935  const int size = ytransA.size();
2936  for( int n = 0; n < size; n++ )
2937  temp += spxAbs(ytransA.value(n));
2938 
2939  // if the one norm is zero then ytransA is zero the Farkas proof should have been verified above
2940  assert(temp != 0);
2941 
2942  // initialize variables in loop: size of Farkas box B, flag whether B has been increased, and number of current
2943  // nonzero in ytransA
2944  Rational B = ytransb / temp;
2945  bool success = false;
2946  int n = 0;
2947 
2948  // loop through nonzeros of ytransA
2949  MSG_DEBUG( std::cout << "B = " << rationalToString(B) << "\n" );
2950  assert(ytransb >= 0);
2951 
2952  while( true )
2953  {
2954  // if all nonzeros have been inspected once without increasing B, we abort; otherwise, we start another round
2955  if( n >= ytransA.size() )
2956  {
2957  if( !success )
2958  break;
2959 
2960  success = false;
2961  n = 0;
2962  }
2963 
2964  // get Farkas multiplier of the bound constraint as minus the value in ytransA
2965  const Rational& minusRedCost = ytransA.value(n);
2966  int colIdx = ytransA.index(n);
2967 
2968  // if the multiplier is positive we inspect the lower bound: if it is finite and within the Farkas box, we can
2969  // increase B by including it in the Farkas proof
2970  assert((upper[colIdx] < _rationalPosInfty) == _upperFinite(_colTypes[colIdx]));
2971  assert((lower[colIdx] > _rationalNegInfty) == _lowerFinite(_colTypes[colIdx]));
2972  if( minusRedCost < 0 && lower[colIdx] > -B && _lowerFinite(_colTypes[colIdx]) )
2973  {
2974  ytransA.clearNum(n);
2975  ytransb.subProduct(minusRedCost, lower[colIdx]);
2976  temp += minusRedCost;
2977 
2978  assert(ytransb >= 0);
2979  assert(temp >= 0);
2980  assert(temp == 0 || ytransb / temp > B);
2981 
2982  // if ytransA and ytransb are zero, we have 0^T x >= 0 and cannot compute a Farkas box
2983  if( temp == 0 && ytransb == 0 )
2984  {
2985  MSG_INFO1( spxout, spxout << "Approximate Farkas proof to weak. Could not compute Farkas box. (2)\n" );
2986  return;
2987  }
2988  // if ytransb is positive and ytransA is zero, we have 0^T x > 0, proving infeasibility
2989  else if( temp == 0 )
2990  {
2991  assert(ytransb > 0);
2992  MSG_INFO1( spxout, spxout << "Farkas infeasibility proof verified exactly. (2)\n" );
2993  return;
2994  }
2995  else
2996  {
2997  B = ytransb / temp;
2998  MSG_DEBUG( std::cout << "B = " << rationalToString(B) << "\n" );
2999  }
3000 
3001  success = true;
3002  }
3003  // if the multiplier is negative we inspect the upper bound: if it is finite and within the Farkas box, we can
3004  // increase B by including it in the Farkas proof
3005  else if( minusRedCost > 0 && upper[colIdx] < B && _upperFinite(_colTypes[colIdx]) )
3006  {
3007  ytransA.clearNum(n);
3008  ytransb.subProduct(minusRedCost, upper[colIdx]);
3009  temp -= minusRedCost;
3010 
3011  assert(ytransb >= 0);
3012  assert(temp >= 0);
3013  assert(temp == 0 || ytransb / temp > B);
3014 
3015  // if ytransA and ytransb are zero, we have 0^T x >= 0 and cannot compute a Farkas box
3016  if( temp == 0 && ytransb == 0 )
3017  {
3018  MSG_INFO1( spxout, spxout << "Approximate Farkas proof to weak. Could not compute Farkas box. (2)\n" );
3019  return;
3020  }
3021  // if ytransb is positive and ytransA is zero, we have 0^T x > 0, proving infeasibility
3022  else if( temp == 0 )
3023  {
3024  assert(ytransb > 0);
3025  MSG_INFO1( spxout, spxout << "Farkas infeasibility proof verified exactly. (2)\n" );
3026  return;
3027  }
3028  else
3029  {
3030  B = ytransb / temp;
3031  MSG_DEBUG( std::cout << "B = " << rationalToString(B) << "\n" );
3032  }
3033 
3034  success = true;
3035  }
3036  // the multiplier is zero, we can ignore the bound constraints on this variable
3037  else if( minusRedCost == 0 )
3038  ytransA.clearNum(n);
3039  // currently this bound cannot be used to increase B; we will check it again in the next round, because B might
3040  // have increased by then
3041  else
3042  n++;
3043  }
3044 
3045  if( B > 0 )
3046  {
3047  MSG_INFO1( spxout, spxout << "Computed Farkas box: provably no feasible solutions with components less than "
3048  << rationalToString(B) << " in absolute value.\n" );
3049  }
3050  }
3051 
3052 
3053 
3054  /// solves real LP during iterative refinement
3055  SPxSolver::Status SoPlex::_solveRealForRational(bool fromscratch, VectorReal& primal, VectorReal& dual, DataArray< SPxSolver::VarStatus >& basisStatusRows, DataArray< SPxSolver::VarStatus >& basisStatusCols, bool& returnedBasis)
3056  {
3057  assert(_isConsistent());
3058 
3059  assert(_solver.nRows() == numRowsRational());
3060  assert(_solver.nCols() == numColsRational());
3061  assert(primal.dim() == numColsRational());
3062  assert(dual.dim() == numRowsRational());
3063 
3065 
3066 #ifndef SOPLEX_MANUAL_ALT
3067  if( fromscratch || !_hasBasis )
3069  else
3071 #else
3073 #endif
3074 
3075  // start timing
3077 
3078  // if preprocessing is applied, we need to restore the original LP at the end
3079  SPxLPRational* rationalLP = 0;
3080  if( _simplifier != 0 || _scaler != 0 )
3081  {
3082  spx_alloc(rationalLP);
3083  rationalLP = new (rationalLP) SPxLPRational(_solver);
3084  }
3085 
3086  // stop timing
3088 
3089  returnedBasis = false;
3090  try
3091  {
3092  // apply problem simplification
3093  SPxSimplifier::Result simplificationStatus = SPxSimplifier::OKAY;
3094  if( _simplifier != 0 )
3095  {
3096  // do not remove bounds of boxed variables or sides of ranged rows if bound flipping is used
3099  }
3100 
3101  // apply scaling after the simplification
3102  if( _scaler != 0 && simplificationStatus == SPxSimplifier::OKAY )
3103  _scaler->scale(_solver);
3104 
3105  // run the simplex method if problem has not been solved by the simplifier
3106  if( simplificationStatus == SPxSimplifier::OKAY )
3107  {
3108  MSG_INFO1( spxout, spxout << std::endl );
3109 
3111 
3112  MSG_INFO1( spxout, spxout << std::endl );
3113  }
3114 
3115  ///@todo move to private helper methods
3116  // evaluate status flag
3117  if( simplificationStatus == SPxSimplifier::INFEASIBLE )
3118  result = SPxSolver::INFEASIBLE;
3119  else if( simplificationStatus == SPxSimplifier::DUAL_INFEASIBLE )
3120  result = SPxSolver::INForUNBD;
3121  else if( simplificationStatus == SPxSimplifier::UNBOUNDED )
3122  result = SPxSolver::UNBOUNDED;
3123  else if( simplificationStatus == SPxSimplifier::VANISHED || simplificationStatus == SPxSimplifier::OKAY )
3124  {
3125  result = simplificationStatus == SPxSimplifier::VANISHED ? SPxSolver::OPTIMAL : _solver.status();
3126 
3127  ///@todo move to private helper methods
3128  // process result
3129  switch( result )
3130  {
3131  case SPxSolver::OPTIMAL:
3132  // unsimplify if simplifier is active and LP is solved to optimality; this must be done here and not at solution
3133  // query, because we want to have the basis for the original problem
3134  if( _simplifier != 0 )
3135  {
3136  assert(!_simplifier->isUnsimplified());
3137  assert(simplificationStatus == SPxSimplifier::VANISHED || simplificationStatus == SPxSimplifier::OKAY);
3138 
3139  bool vanished = simplificationStatus == SPxSimplifier::VANISHED;
3140 
3141  // get solution vectors for transformed problem
3142  DVectorReal tmpPrimal(vanished ? 0 : _solver.nCols());
3143  DVectorReal tmpSlacks(vanished ? 0 : _solver.nRows());
3144  DVectorReal tmpDual(vanished ? 0 : _solver.nRows());
3145  DVectorReal tmpRedCost(vanished ? 0 : _solver.nCols());
3146 
3147  if( !vanished )
3148  {
3149  assert(_solver.status() == SPxSolver::OPTIMAL);
3150 
3151  _solver.getPrimal(tmpPrimal);
3152  _solver.getSlacks(tmpSlacks);
3153  _solver.getDual(tmpDual);
3154  _solver.getRedCost(tmpRedCost);
3155 
3156  // unscale vectors
3157  if( _scaler != 0 )
3158  {
3159  _scaler->unscalePrimal(tmpPrimal);
3160  _scaler->unscaleSlacks(tmpSlacks);
3161  _scaler->unscaleDual(tmpDual);
3162  _scaler->unscaleRedCost(tmpRedCost);
3163  }
3164 
3165  // get basis of transformed problem
3166  basisStatusRows.reSize(_solver.nRows());
3167  basisStatusCols.reSize(_solver.nCols());
3168  _solver.getBasis(basisStatusRows.get_ptr(), basisStatusCols.get_ptr(), basisStatusRows.size(), basisStatusCols.size());
3169  }
3170 
3171  ///@todo catch exception
3172  _simplifier->unsimplify(tmpPrimal, tmpDual, tmpSlacks, tmpRedCost, basisStatusRows.get_ptr(), basisStatusCols.get_ptr());
3173 
3174  // store basis for original problem
3175  basisStatusRows.reSize(numRowsRational());
3176  basisStatusCols.reSize(numColsRational());
3177  _simplifier->getBasis(basisStatusRows.get_ptr(), basisStatusCols.get_ptr(), basisStatusRows.size(), basisStatusCols.size());
3178  returnedBasis = true;
3179 
3180  primal = _simplifier->unsimplifiedPrimal();
3181  dual = _simplifier->unsimplifiedDual();
3182  }
3183  // if the original problem is not in the solver because of scaling, we also need to store the basis
3184  else
3185  {
3186  _solver.getPrimal(primal);
3187  _solver.getDual(dual);
3188 
3189  // unscale vectors
3190  if( _scaler != 0 )
3191  {
3192  _scaler->unscalePrimal(primal);
3193  _scaler->unscaleDual(dual);
3194  }
3195 
3196  // get basis of transformed problem
3197  basisStatusRows.reSize(_solver.nRows());
3198  basisStatusCols.reSize(_solver.nCols());
3199  _solver.getBasis(basisStatusRows.get_ptr(), basisStatusCols.get_ptr(), basisStatusRows.size(), basisStatusCols.size());
3200  returnedBasis = true;
3201  }
3202  break;
3203 
3206  {
3207  _solver.getPrimal(primal);
3208  _solver.getDual(dual);
3209 
3210  // unscale vectors
3211  if( _scaler != 0 )
3212  {
3213  _scaler->unscalePrimal(primal);
3214  _scaler->unscaleDual(dual);
3215  }
3216 
3217  // get basis of transformed problem
3218  basisStatusRows.reSize(_solver.nRows());
3219  basisStatusCols.reSize(_solver.nCols());
3220  _solver.getBasis(basisStatusRows.get_ptr(), basisStatusCols.get_ptr(), basisStatusRows.size(), basisStatusCols.size());
3221  returnedBasis = true;
3222  result = SPxSolver::OPTIMAL;
3223  }
3224  break;
3225  case SPxSolver::ABORT_TIME:
3226  case SPxSolver::ABORT_ITER:
3228  case SPxSolver::REGULAR:
3229  case SPxSolver::RUNNING:
3230  case SPxSolver::UNBOUNDED:
3231  break;
3232  case SPxSolver::INFEASIBLE:
3233  // if simplifier is active we cannot return a Farkas ray currently
3234  if( _simplifier != 0 )
3235  break;
3236 
3237  // return Farkas ray as dual solution
3238  _solver.getDualfarkas(dual);
3239 
3240  // unscale vectors
3241  if( _scaler != 0 )
3242  _scaler->unscaleDual(dual);
3243 
3244  // if the original problem is not in the solver because of scaling, we also need to store the basis
3245  basisStatusRows.reSize(_solver.nRows());
3246  basisStatusCols.reSize(_solver.nCols());
3247  _solver.getBasis(basisStatusRows.get_ptr(), basisStatusCols.get_ptr(), basisStatusRows.size(), basisStatusCols.size());
3248  returnedBasis = true;
3249  break;
3250 
3251  case SPxSolver::INForUNBD:
3252  case SPxSolver::SINGULAR:
3253  default:
3254  break;
3255  }
3256  }
3257  }
3258  catch( ... )
3259  {
3260  MSG_INFO1( spxout, spxout << "Exception thrown during floating-point solve.\n" );
3261  result = SPxSolver::ERROR;
3262  }
3263 
3264  // restore original LP if necessary
3265  if( _simplifier != 0 || _scaler != 0 )
3266  {
3267  assert(rationalLP != 0);
3268  _solver.loadLP((SPxLPReal)(*rationalLP));
3269  rationalLP->~SPxLPRational();
3270  spx_free(rationalLP);
3271  if( _hasBasis )
3272  _solver.setBasis(basisStatusRows.get_ptr(), basisStatusCols.get_ptr());
3273  }
3274 
3275  return result;
3276  }
3277 
3278 
3279 
3280  /// solves real LP with recovery mechanism
3281  SPxSolver::Status SoPlex::_solveRealStable(bool acceptUnbounded, bool acceptInfeasible, VectorReal& primal, VectorReal& dual, DataArray< SPxSolver::VarStatus >& basisStatusRows, DataArray< SPxSolver::VarStatus >& basisStatusCols, bool& returnedBasis, const bool forceNoSimplifier)
3282  {
3284 
3285  bool fromScratch = false;
3286  bool solved = false;
3287  bool solvedFromScratch = false;
3288  bool initialSolve = true;
3289  bool increasedMarkowitz = false;
3290  bool relaxedTolerances = false;
3291  bool tightenedTolerances = false;
3292  bool switchedScaler = false;
3293  bool switchedSimplifier = false;
3294  bool switchedRatiotester = false;
3295  bool switchedPricer = false;
3296  bool turnedoffPre = false;
3297 
3298  Real markowitz = _slufactor.markowitz();
3299  int ratiotester = intParam(SoPlex::RATIOTESTER);
3300  int pricer = intParam(SoPlex::PRICER);
3301  int simplifier = intParam(SoPlex::SIMPLIFIER);
3302  int scaler = intParam(SoPlex::SCALER);
3303  int type = intParam(SoPlex::ALGORITHM);
3304 
3305  if( forceNoSimplifier )
3307 
3308  while( true )
3309  {
3310  assert(!increasedMarkowitz || GE(_slufactor.markowitz(), 0.9));
3311 
3312  result = _solveRealForRational(fromScratch, primal, dual, basisStatusRows, basisStatusCols, returnedBasis);
3313 
3314  solved = (result == SPxSolver::OPTIMAL)
3315  || (result == SPxSolver::INFEASIBLE && acceptInfeasible)
3316  || (result == SPxSolver::UNBOUNDED && acceptUnbounded);
3317 
3318  if( solved )
3319  break;
3320 
3321 // if( _isSolveStopped() )
3322 // break;
3323 
3324  if( initialSolve )
3325  {
3326  MSG_INFO1( spxout, spxout << "Numerical troubles during floating-point solve." << std::endl );
3327  initialSolve = false;
3328  }
3329 
3330  if( !turnedoffPre
3332  {
3333  MSG_INFO1( spxout, spxout << "Turning off preprocessing." << std::endl );
3334 
3335  turnedoffPre = true;
3336 
3339 
3340  fromScratch = true;
3341  _solver.reLoad();
3342  solvedFromScratch = true;
3343  continue;
3344  }
3345 
3346  setIntParam(SoPlex::SCALER, scaler);
3347  setIntParam(SoPlex::SIMPLIFIER, simplifier);
3348 
3349  if( !increasedMarkowitz )
3350  {
3351  MSG_INFO1( spxout, spxout << "Increasing Markowitz threshold." << std::endl );
3352 
3353  _slufactor.setMarkowitz(0.9);
3354  increasedMarkowitz = true;
3355  try
3356  {
3357  _solver.factorize();
3358  continue;
3359  }
3360  catch( ... )
3361  {
3362  MSG_DEBUG( std::cout << std::endl << "Factorization failed." << std::endl );
3363  }
3364  }
3365 
3366  if( !solvedFromScratch )
3367  {
3368  MSG_INFO1( spxout, spxout << "Solving from scratch." << std::endl );
3369 
3370  fromScratch = true;
3371  _solver.reLoad();
3372 
3373  solvedFromScratch = true;
3374  continue;
3375  }
3376 
3377  setIntParam(SoPlex::RATIOTESTER, ratiotester);
3378  setIntParam(SoPlex::PRICER, pricer);
3379 
3380  if( !switchedScaler )
3381  {
3382  MSG_INFO1( spxout, spxout << "Switching scaling." << std::endl );
3383 
3384  if( scaler == int(SoPlex::SCALER_OFF) )
3386  else
3388 
3389  fromScratch = true;
3390  _solver.reLoad();
3391 
3392  solvedFromScratch = true;
3393  switchedScaler = true;
3394  continue;
3395  }
3396 
3397  if( !switchedSimplifier && !forceNoSimplifier )
3398  {
3399  MSG_INFO1( spxout, spxout << "Switching simplification." << std::endl );
3400 
3401  if( simplifier == int(SoPlex::SIMPLIFIER_OFF) )
3403  else
3405 
3406  fromScratch = true;
3407  _solver.reLoad();
3408 
3409  solvedFromScratch = true;
3410  switchedSimplifier = true;
3411  continue;
3412  }
3413 
3415 
3416  if( !relaxedTolerances )
3417  {
3418  MSG_INFO1( spxout, spxout << "Relaxing tolerances." << std::endl );
3419 
3421  _solver.setDelta((_solver.feastol() * 1e3 > 1e-3) ? 1e-3 : (_solver.feastol() * 1e3));
3422  relaxedTolerances = _solver.feastol() >= 1e-3;
3423  solvedFromScratch = false;
3424  continue;
3425  }
3426 
3427  if( !tightenedTolerances && result != SPxSolver::INFEASIBLE )
3428  {
3429  MSG_INFO1( spxout, spxout << "Tightening tolerances." << std::endl );
3430 
3432  _solver.setDelta(_solver.feastol() * 1e-3 < 1e-9 ? 1e-9 : _solver.feastol() * 1e-3);
3433  tightenedTolerances = _solver.feastol() <= 1e-9;
3434  solvedFromScratch = false;
3435  continue;
3436  }
3437 
3439 
3440  if( !switchedRatiotester )
3441  {
3442  MSG_INFO1( spxout, spxout << "Switching ratio test." << std::endl );
3443 
3447  else
3449  switchedRatiotester = true;
3450  solvedFromScratch = false;
3451  continue;
3452  }
3453 
3454  if( !switchedPricer )
3455  {
3456  MSG_INFO1( spxout, spxout << "Switching pricer." << std::endl );
3457 
3459  if( _solver.pricer() != (SPxPricer*)&_pricerDevex )
3461  else
3463  switchedPricer = true;
3464  solvedFromScratch = false;
3465  continue;
3466  }
3467 
3468  MSG_INFO1( spxout, spxout << "Giving up." << std::endl );
3469 
3470  break;
3471  }
3472 
3475  _slufactor.setMarkowitz(markowitz);
3476 
3477  setIntParam(SoPlex::RATIOTESTER, ratiotester);
3478  setIntParam(SoPlex::PRICER, pricer);
3479  setIntParam(SoPlex::SIMPLIFIER, simplifier);
3480  setIntParam(SoPlex::SCALER, scaler);
3482 
3483  return result;
3484  }
3485 
3486 
3487 
3488  /// factorizes rational basis matrix in column representation
3489  void SoPlex::_factorizeColumnRational(SolRational& sol, DataArray< SPxSolver::VarStatus >& basisStatusRows, DataArray< SPxSolver::VarStatus >& basisStatusCols, bool& stoppedTime, bool& stoppedIter, bool& error, bool& optimal)
3490  {
3491  // start rational solving time
3493 
3494  stoppedTime = false;
3495  stoppedIter = false;
3496  error = false;
3497  optimal = false;
3498 
3499  SLUFactorRational linsolver;
3500  const bool maximizing = (intParam(SoPlex::OBJSENSE) == SoPlex::OBJSENSE_MAXIMIZE);
3501  const int matrixdim = numRowsRational();
3502  DataArray< const SVectorRational* > matrix(matrixdim);
3503  int numBasicRows;
3504 
3505  _workSol._primal.reDim(matrixdim);
3506  _workSol._slacks.reDim(matrixdim);
3507  _workSol._dual.reDim(matrixdim);
3508  _workSol._redCost.reDim(numColsRational() > matrixdim ? numColsRational() : matrixdim);
3509 
3510  DVectorRational& basicPrimalRhs = _workSol._slacks;
3511  DVectorRational& basicDualRhs = _workSol._redCost;
3512  DVectorRational& basicPrimal = _workSol._primal;
3513  DVectorRational& basicDual = _workSol._dual;
3514 
3515  Rational violation;
3516  Rational primalViolation;
3517  Rational dualViolation;
3518  bool primalFeasible = false;
3519  bool dualFeasible = false;
3520 
3521  assert(basisStatusCols.size() == numColsRational());
3522  assert(basisStatusRows.size() == numRowsRational());
3523 
3524  int j = 0;
3525  for( int i = 0; i < basisStatusRows.size(); i++ )
3526  {
3527  if( basisStatusRows[i] == SPxSolver::BASIC && j < matrixdim )
3528  {
3529  basicPrimalRhs[i] = 0;
3530  basicDualRhs[j] = 0;
3531  matrix[j] = _unitVectorRational(i);
3532  j++;
3533  }
3534  else if( basisStatusRows[i] == SPxSolver::ON_LOWER )
3535  basicPrimalRhs[i] = lhsRational(i);
3536  else if( basisStatusRows[i] == SPxSolver::ON_UPPER )
3537  basicPrimalRhs[i] = rhsRational(i);
3538  else if( basisStatusRows[i] == SPxSolver::ZERO )
3539  basicPrimalRhs[i] = 0;
3540  else if( basisStatusRows[i] == SPxSolver::FIXED )
3541  {
3542  assert(lhsRational(i) == rhsRational(i));
3543  basicPrimalRhs[i] = lhsRational(i);
3544  }
3545  else if( basisStatusRows[i] == SPxSolver::UNDEFINED )
3546  {
3547  MSG_ERROR( std::cerr << "Undefined basis status of row in rational factorization.\n" );
3548  error = true;
3549  goto TERMINATE;
3550  }
3551  else
3552  {
3553  assert(basisStatusRows[i] == SPxSolver::BASIC);
3554  MSG_ERROR( std::cerr << "Too many basic rows in rational factorization.\n" );
3555  error = true;
3556  goto TERMINATE;
3557  }
3558  }
3559  numBasicRows = j;
3560 
3561  for( int i = 0; i < basisStatusCols.size(); i++ )
3562  {
3563  if( basisStatusCols[i] == SPxSolver::BASIC && j < matrixdim )
3564  {
3565  basicDualRhs[j] = objRational(i);
3566  matrix[j] = &colVectorRational(i);
3567  j++;
3568  }
3569  else if( basisStatusCols[i] == SPxSolver::ON_LOWER )
3570  basicPrimalRhs.multAdd(-lowerRational(i), colVectorRational(i));
3571  else if( basisStatusCols[i] == SPxSolver::ON_UPPER )
3572  basicPrimalRhs.multAdd(-upperRational(i), colVectorRational(i));
3573  else if( basisStatusCols[i] == SPxSolver::ZERO )
3574  {}
3575  else if( basisStatusCols[i] == SPxSolver::FIXED )
3576  {
3577  assert(lowerRational(i) == upperRational(i));
3578  basicPrimalRhs.multAdd(-lowerRational(i), colVectorRational(i));
3579  }
3580  else if( basisStatusCols[i] == SPxSolver::UNDEFINED )
3581  {
3582  MSG_ERROR( std::cerr << "Undefined basis status of column in rational factorization.\n" );
3583  error = true;
3584  goto TERMINATE;
3585  }
3586  else
3587  {
3588  assert(basisStatusCols[i] == SPxSolver::BASIC);
3589  MSG_ERROR( std::cerr << "Too many basic columns in rational factorization.\n" );
3590  error = true;
3591  goto TERMINATE;
3592  }
3593  }
3594 
3595  if( j != matrixdim )
3596  {
3597  MSG_ERROR( std::cerr << "Too few basic entries in rational factorization.\n" );
3598  error = true;
3599  goto TERMINATE;
3600  }
3601 
3602  // load rational basis matrix
3605  else
3606  linsolver.setTimeLimit(-1.0);
3607  linsolver.load(matrix.get_ptr(), matrixdim);
3608 
3609  // record statistics
3612  linsolver.resetCounters();
3613 
3614  if( linsolver.status() == SLinSolverRational::TIME )
3615  {
3616  MSG_INFO2( spxout, spxout << "Rational factorization hit time limit.\n" );
3617  stoppedTime = true;
3618  return;
3619  }
3620  else if( linsolver.status() != SLinSolverRational::OK )
3621  {
3622  MSG_ERROR( std::cerr << "Error performing rational LU factorization.\n" );
3623  error = true;
3624  return;
3625  }
3626 
3627  // solve for primal solution
3630  else
3631  linsolver.setTimeLimit(-1.0);
3632  linsolver.solveRight(basicPrimal, basicPrimalRhs);
3633 
3634  // record statistics
3636  linsolver.resetCounters();
3637 
3638  if( _isSolveStopped(stoppedTime, stoppedIter) )
3639  {
3640  MSG_INFO2( spxout, spxout << "Rational factorization hit time limit while solving for primal.\n" );
3641  return;
3642  }
3643 
3644  // check bound violation on basic rows and columns
3645  j = 0;
3646  primalViolation = 0;
3647  primalFeasible = true;
3648  for( int i = 0; i < basisStatusRows.size(); i++ )
3649  {
3650  if( basisStatusRows[i] == SPxSolver::BASIC )
3651  {
3652  assert(j < matrixdim);
3653  violation = lhsRational(i);
3654  violation += basicPrimal[j];
3655  if( violation > primalViolation )
3656  {
3657  primalFeasible = false;
3658  primalViolation = violation;
3659  }
3660  violation = rhsRational(i);
3661  violation += basicPrimal[j];
3662  violation *= -1;
3663  if( violation > primalViolation )
3664  {
3665  primalFeasible = false;
3666  primalViolation = violation;
3667  }
3668  j++;
3669  }
3670  }
3671  for( int i = 0; i < basisStatusCols.size(); i++ )
3672  {
3673  if( basisStatusCols[i] == SPxSolver::BASIC )
3674  {
3675  assert(j < matrixdim);
3676  if( basicPrimal[j] < lowerRational(i) )
3677  {
3678  violation = lowerRational(i);
3679  violation -= basicPrimal[j];
3680  if( violation > primalViolation )
3681  {
3682  primalFeasible = false;
3683  primalViolation = violation;
3684  }
3685  }
3686  if( basicPrimal[j] > upperRational(i) )
3687  {
3688  violation = basicPrimal[j];
3689  violation -= upperRational(i);
3690  if( violation > primalViolation )
3691  {
3692  primalFeasible = false;
3693  primalViolation = violation;
3694  }
3695  }
3696  j++;
3697  }
3698  }
3699 
3700  if( !primalFeasible )
3701  {
3702  MSG_INFO1( spxout, spxout << "Rational solution primal infeasible.\n" );
3703  }
3704 
3705  // solve for dual solution
3708  else
3709  linsolver.setTimeLimit(-1.0);
3710  linsolver.solveLeft(basicDual, basicDualRhs);
3711 
3712  // record statistics
3714  linsolver.resetCounters();
3715 
3716  if( _isSolveStopped(stoppedTime, stoppedIter) )
3717  {
3718  MSG_INFO2( spxout, spxout << "Rational factorization hit time limit while solving for dual.\n" );
3719  return;
3720  }
3721 
3722  // check dual violation on nonbasic rows
3723  dualViolation = 0;
3724  dualFeasible = true;
3725  for( int i = 0; i < basisStatusRows.size(); i++ )
3726  {
3727  if( _rowTypes[i] == RANGETYPE_FIXED
3728  && (basisStatusRows[i] == SPxSolver::ON_LOWER || basisStatusRows[i] == SPxSolver::ON_UPPER) )
3729  {
3730  assert(lhsRational(i) == rhsRational(i));
3731  basisStatusRows[i] = SPxSolver::FIXED;
3732  }
3733 
3734  assert(basisStatusRows[i] != SPxSolver::BASIC || basicDual[i] == 0);
3735  if( basisStatusRows[i] == SPxSolver::BASIC || basisStatusRows[i] == SPxSolver::FIXED )
3736  continue;
3737  else if( basicDual[i] < 0 )
3738  {
3739  if( ((maximizing && basisStatusRows[i] != SPxSolver::ON_LOWER) || (!maximizing && basisStatusRows[i] != SPxSolver::ON_UPPER))
3740  && (basisStatusRows[i] != SPxSolver::ZERO || rhsRational(i) != 0) )
3741  {
3742  dualFeasible = false;
3743  violation = -basicDual[i];
3744  if( violation > dualViolation )
3745  dualViolation = violation;
3746  MSG_DEBUG( spxout << "negative dual multliplier for row " << i
3747  << " with dual " << rationalToString(basicDual[i])
3748  << " and status " << basisStatusRows[i]
3749  << " and [lhs,rhs] = [" << rationalToString(lhsRational(i)) << "," << rationalToString(rhsRational(i)) << "]"
3750  << "\n" );
3751  }
3752  }
3753  else if( basicDual[i] > 0 )
3754  {
3755  if( ((maximizing && basisStatusRows[i] != SPxSolver::ON_UPPER) || (!maximizing && basisStatusRows[i] != SPxSolver::ON_LOWER))
3756  && (basisStatusRows[i] != SPxSolver::ZERO || lhsRational(i) == 0) )
3757  {
3758  dualFeasible = false;
3759  if( basicDual[i] > dualViolation )
3760  dualViolation = basicDual[i];
3761  MSG_DEBUG( spxout << "positive dual multliplier for row " << i
3762  << " with dual " << rationalToString(basicDual[i])
3763  << " and status " << basisStatusRows[i]
3764  << " and [lhs,rhs] = [" << rationalToString(lhsRational(i)) << "," << rationalToString(rhsRational(i)) << "]"
3765  << "\n" );
3766  }
3767  }
3768  }
3769 
3770  // check reduced cost violation on nonbasic columns
3771  for( int i = 0; i < basisStatusCols.size(); i++ )
3772  {
3773  if( _colTypes[i] == RANGETYPE_FIXED
3774  && (basisStatusCols[i] == SPxSolver::ON_LOWER || basisStatusCols[i] == SPxSolver::ON_UPPER) )
3775  {
3776  assert(lowerRational(i) == upperRational(i));
3777  basisStatusCols[i] = SPxSolver::FIXED;
3778  }
3779 
3780  assert(basisStatusCols[i] != SPxSolver::BASIC || basicDual * colVectorRational(i) == objRational(i));
3781  if( basisStatusCols[i] == SPxSolver::BASIC || basisStatusCols[i] == SPxSolver::FIXED )
3782  continue;
3783  else
3784  {
3785  _workSol._redCost[i] = basicDual * colVectorRational(i);
3786  _workSol._redCost[i] -= objRational(i);
3787  if( _workSol._redCost[i] > 0 )
3788  {
3789  if( ((maximizing && basisStatusCols[i] != SPxSolver::ON_LOWER) || (!maximizing && basisStatusCols[i] != SPxSolver::ON_UPPER))
3790  && (basisStatusCols[i] != SPxSolver::ZERO || upperRational(i) != 0) )
3791  {
3792  dualFeasible = false;
3793  if( _workSol._redCost[i] > dualViolation )
3794  dualViolation = _workSol._redCost[i];
3795  }
3796  _workSol._redCost[i] *= -1;
3797  }
3798  else if( _workSol._redCost[i] < 0 )
3799  {
3800  _workSol._redCost[i] *= -1;
3801  if( ((maximizing && basisStatusCols[i] != SPxSolver::ON_UPPER) || (!maximizing && basisStatusCols[i] != SPxSolver::ON_LOWER))
3802  && (basisStatusCols[i] != SPxSolver::ZERO || lowerRational(i) != 0) )
3803  {
3804  dualFeasible = false;
3805  if( _workSol._redCost[i] > dualViolation )
3806  dualViolation = _workSol._redCost[i];
3807  }
3808  }
3809  else
3810  _workSol._redCost[i] *= -1;
3811  }
3812  }
3813 
3814  if( !dualFeasible )
3815  {
3816  MSG_INFO1( spxout, spxout << "Rational solution dual infeasible.\n" );
3817  }
3818 
3819  // store solution
3820  optimal = primalFeasible && dualFeasible;
3821  if( optimal || boolParam(SoPlex::RATFACJUMP) )
3822  {
3823  _hasBasis = true;
3824  if( &basisStatusRows != &_basisStatusRows )
3825  _basisStatusRows = basisStatusRows;
3826  if( &basisStatusCols != &_basisStatusCols )
3827  _basisStatusCols = basisStatusCols;
3828 
3829  sol._primal.reDim(numColsRational());
3830  j = numBasicRows;
3831  for( int i = 0; i < basisStatusCols.size(); i++ )
3832  {
3833  if( basisStatusCols[i] == SPxSolver::BASIC )
3834  {
3835  assert(j < matrixdim);
3836  sol._primal[i] = basicPrimal[j];
3837  j++;
3838  }
3839  else if( basisStatusCols[i] == SPxSolver::ON_LOWER )
3840  sol._primal[i] = lowerRational(i);
3841  else if( basisStatusCols[i] == SPxSolver::ON_UPPER )
3842  sol._primal[i] = upperRational(i);
3843  else if( basisStatusCols[i] == SPxSolver::ZERO )
3844  sol._primal[i] = 0;
3845  else if( basisStatusCols[i] == SPxSolver::FIXED )
3846  {
3847  assert(lowerRational(i) == upperRational(i));
3848  sol._primal[i] = lowerRational(i);
3849  }
3850  else
3851  {
3852  assert(basisStatusCols[i] == SPxSolver::UNDEFINED);
3853  MSG_ERROR( std::cerr << "Undefined basis status of column in rational factorization.\n" );
3854  error = true;
3855  goto TERMINATE;
3856  }
3857  }
3858  sol._slacks.reDim(numRowsRational());
3860  sol._hasPrimal = true;
3861 
3862  sol._dual = basicDual;
3863  for( int i = 0; i < numColsRational(); i++ )
3864  {
3865  if( basisStatusCols[i] == SPxSolver::BASIC )
3866  sol._redCost[i] = 0;
3867  else if( basisStatusCols[i] == SPxSolver::FIXED )
3868  {
3869  sol._redCost[i] = basicDual * colVectorRational(i);
3870  sol._redCost[i] -= objRational(i);
3871  sol._redCost[i] *= -1;
3872  }
3873  else
3874  sol._redCost[i] = _workSol._redCost[i];
3875  }
3876  sol._hasDual = true;
3877  }
3878 
3879  TERMINATE:
3880  // stop rational solving time
3882  return;
3883  }
3884 
3885  /// attempts rational reconstruction of primal-dual solution
3887  {
3888  bool success;
3889  bool isSolBasic;
3890  DIdxSet basicIndices(numColsRational());
3891 
3892  success = false;
3893  isSolBasic = true;
3894 
3895  if( !sol.hasPrimal() || !sol.hasDual() )
3896  return success;
3897 
3898  // start timing and increment statistics counter
3901 
3902  // reconstruct primal vector
3903  _workSol._primal = sol._primal;
3904  for( int j = 0; j < numColsRational(); ++j )
3905  {
3906  if( basisStatusCols[j] == SPxSolver::BASIC )
3907  basicIndices.addIdx(j);
3908  }
3909 
3910  success = reconstructVector(_workSol._primal, denomBoundSquared, &basicIndices);
3911 
3912  if( !success )
3913  {
3914  MSG_INFO1( spxout, spxout << "Rational reconstruction of primal solution failed.\n" );
3916  return success;
3917  }
3918 
3919  MSG_DEBUG( spxout << "Rational reconstruction of primal solution successful.\n" );
3920 
3921  // check violation of bounds
3922  for( int c = numColsRational() - 1; c >= 0; c-- )
3923  {
3924  // we want to notify the user whether the reconstructed solution is basic; otherwise, this would be redundant
3925  SPxSolver::VarStatus& basisStatusCol = _basisStatusCols[c];
3926  if( (basisStatusCol == SPxSolver::FIXED && _workSol._primal[c] != lowerRational(c))
3927  || (basisStatusCol == SPxSolver::ON_LOWER && _workSol._primal[c] != lowerRational(c))
3928  || (basisStatusCol == SPxSolver::ON_UPPER && _workSol._primal[c] != upperRational(c))
3929  || (basisStatusCol == SPxSolver::ZERO && _workSol._primal[c] != 0)
3930  || (basisStatusCol == SPxSolver::UNDEFINED) )
3931  {
3932  isSolBasic = false;
3933  }
3934 
3935  if( _lowerFinite(_colTypes[c]) && _workSol._primal[c] < lowerRational(c) )
3936  {
3937  MSG_DEBUG( std::cout << "Lower bound of variable " << c << " violated by " << rationalToString(lowerRational(c) - _workSol._primal[c]) << "\n" );
3938  MSG_INFO1( spxout, spxout << "Reconstructed solution primal infeasible (1).\n" );
3940  return false;
3941  }
3942 
3943  if( _upperFinite(_colTypes[c]) && _workSol._primal[c] > upperRational(c) )
3944  {
3945  MSG_DEBUG( std::cout << "Upper bound of variable " << c << " violated by " << rationalToString(_workSol._primal[c] - upperRational(c)) << "\n" );
3946  MSG_INFO1( spxout, spxout << "Reconstructed solution primal infeasible (2).\n" );
3948  return false;
3949  }
3950  }
3951 
3952  // compute slacks
3953  ///@todo we should compute them one by one so we can abort when encountering an infeasibility
3956 
3957  // check violation of sides
3958  for( int r = numRowsRational() - 1; r >= 0; r-- )
3959  {
3960  // we want to notify the user whether the reconstructed solution is basic; otherwise, this would be redundant
3961  SPxSolver::VarStatus& basisStatusRow = _basisStatusRows[r];
3962  if( (basisStatusRow == SPxSolver::FIXED && _workSol._slacks[r] != lhsRational(r))
3963  || (basisStatusRow == SPxSolver::ON_LOWER && _workSol._slacks[r] != lhsRational(r))
3964  || (basisStatusRow == SPxSolver::ON_UPPER && _workSol._slacks[r] != rhsRational(r))
3965  || (basisStatusRow == SPxSolver::ZERO && _workSol._slacks[r] != 0)
3966  || (basisStatusRow == SPxSolver::UNDEFINED) )
3967  {
3968  isSolBasic = false;
3969  }
3970 
3971  if( _lowerFinite(_rowTypes[r]) && _workSol._slacks[r] < lhsRational(r) )
3972  {
3973  MSG_DEBUG( std::cout << "Lhs of row " << r << " violated by " << rationalToString(lhsRational(r) - _workSol._slacks[r]) << "\n" );
3974  MSG_INFO1( spxout, spxout << "Reconstructed solution primal infeasible (3).\n" );
3976  return false;
3977  }
3978 
3979  if( _upperFinite(_rowTypes[r]) && _workSol._slacks[r] > rhsRational(r) )
3980  {
3981  MSG_DEBUG( std::cout << "Rhs of row " << r << " violated by " << rationalToString(_workSol._slacks[r] - rhsRational(r)) << "\n" );
3982  MSG_INFO1( spxout, spxout << "Reconstructed solution primal infeasible (4).\n" );
3984  return false;
3985  }
3986  }
3987 
3988  // reconstruct dual vector
3989  _workSol._dual = sol._dual;
3990 
3991  success = reconstructVector(_workSol._dual, denomBoundSquared);
3992 
3993  if( !success )
3994  {
3995  MSG_INFO1( spxout, spxout << "Rational reconstruction of dual solution failed.\n" );
3997  return success;
3998  }
3999 
4000  MSG_DEBUG( spxout << "Rational reconstruction of dual vector successful.\n" );
4001 
4002  // check dual multipliers before reduced costs because this check is faster since it does not require the
4003  // computation of reduced costs
4004  const bool maximizing = (intParam(SoPlex::OBJSENSE) == SoPlex::OBJSENSE_MAXIMIZE);
4005  for( int r = numRowsRational() - 1; r >= 0; r-- )
4006  {
4007  int sig = sign(_workSol._dual[r]);
4008 
4009  if( (!maximizing && sig > 0) || (maximizing && sig < 0) )
4010  {
4011  if( !_lowerFinite(_rowTypes[r]) || _workSol._slacks[r] > lhsRational(r) )
4012  {
4013  MSG_DEBUG( std::cout << "complementary slackness violated by row " << r
4014  << " with dual " << rationalToString(_workSol._dual[r])
4015  << " and slack " << rationalToString(_workSol._slacks[r])
4016  << " not at lhs " << rationalToString(lhsRational(r))
4017  << "\n" );
4018  MSG_INFO1( spxout, spxout << "Reconstructed solution dual infeasible (1).\n" );
4020  return false;
4021  }
4022 
4024  {
4026  isSolBasic = false;
4027  else
4029  }
4030  }
4031  else if( (!maximizing && sig < 0) || (maximizing && sig > 0) )
4032  {
4033  if( !_upperFinite(_rowTypes[r]) || _workSol._slacks[r] < rhsRational(r) )
4034  {
4035  MSG_DEBUG( std::cout << "complementary slackness violated by row " << r
4036  << " with dual " << rationalToString(_workSol._dual[r])
4037  << " and slack " << rationalToString(_workSol._slacks[r])
4038  << " not at rhs " << rationalToString(rhsRational(r))
4039  << "\n" );
4040  MSG_INFO1( spxout, spxout << "Reconstructed solution dual infeasible (2).\n" );
4042  return false;
4043  }
4044 
4046  {
4048  isSolBasic = false;
4049  else
4051  }
4052  }
4053  }
4054 
4055  // compute reduced cost vector; we assume that the objective function vector has less nonzeros than the reduced
4056  // cost vector, and so multiplying with -1 first and subtracting the dual activity should be faster than adding
4057  // the dual activity and negating afterwards
4058  ///@todo we should compute them one by one so we can abort when encountering an infeasibility
4062 
4063  // check reduced cost violation
4064  for( int c = numColsRational() - 1; c >= 0; c-- )
4065  {
4066  int sig = sign(_workSol._redCost[c]);
4067 
4068  if( (!maximizing && sig > 0) || (maximizing && sig < 0) )
4069  {
4070  if( !_lowerFinite(_colTypes[c]) || _workSol._primal[c] > lowerRational(c) )
4071  {
4072  MSG_DEBUG( std::cout << "complementary slackness violated by column " << c
4073  << " with reduced cost " << rationalToString(_workSol._redCost[c])
4074  << " and value " << rationalToString(_workSol._primal[c])
4075  << " not at lower bound " << rationalToString(lowerRational(c))
4076  << "\n" );
4077  MSG_INFO1( spxout, spxout << "Reconstructed solution dual infeasible (3).\n" );
4079  return false;
4080  }
4081 
4083  {
4085  isSolBasic = false;
4086  else
4088  }
4089  }
4090  else if( (!maximizing && sig < 0) || (maximizing && sig > 0) )
4091  {
4092  if( !_upperFinite(_colTypes[c]) || _workSol._primal[c] < upperRational(c) )
4093  {
4094  MSG_DEBUG( std::cout << "complementary slackness violated by column " << c
4095  << " with reduced cost " << rationalToString(_workSol._redCost[c])
4096  << " and value " << rationalToString(_workSol._primal[c])
4097  << " not at upper bound " << rationalToString(upperRational(c))
4098  << "\n" );
4099  MSG_INFO1( spxout, spxout << "Reconstructed solution dual infeasible (4).\n" );
4101  return false;
4102  }
4103 
4105  {
4107  isSolBasic = false;
4108  else
4110  }
4111  }
4112  }
4113 
4114  // update solution
4115  sol._primal = _workSol._primal;
4116  sol._slacks = _workSol._slacks;
4117  sol._dual = _workSol._dual;
4118  sol._redCost = _workSol._redCost;
4119 
4120  if( !isSolBasic )
4121  {
4122  MSG_WARNING( spxout, spxout << "Warning: Reconstructed solution not basic.\n" );
4123  _hasBasis = false;
4124  }
4125 
4126  // stop timing
4128 
4129  return success;
4130  }
4131 } // namespace soplex
4132 #endif
Rational spxAbs(const Rational &r)
Absolute.
Definition: rational.cpp:3925
int _lastSolveMode
Definition: soplex.h:1429
virtual void subDualActivity(const VectorBase< R > &dual, VectorBase< R > &activity) const
Updates "dual" activity of the columns for a given dual vector, i.e., y^T A; activity does not need t...
Definition: spxlpbase.h:1727
DVectorBase< R > _slacks
Definition: solbase.h:223
textbook ratio test without stabilization
Definition: soplex.h:1054
R obj(int i) const
Returns objective value of column i.
Definition: spxlpbase.h:362
int numColsReal() const
returns number of columns
Definition: soplex.cpp:779
const VectorReal & lowerReal() const
returns lower bound vector
Definition: soplex.cpp:932
virtual Status getPrimal(Vector &vector) const
get solution vector for primal variables.
Definition: spxsolve.cpp:1219
int unbdRefinements
number of refinement steps during undboundedness test
Definition: statistics.h:113
virtual void removeRow(int i)
Removes i &#39;th row.
Definition: spxlpbase.h:813
Rational _rationalMaxscaleincr
Definition: soplex.h:1325
free variable fixed to zero.
Definition: spxsolver.h:182
DVectorRational _unboundedRhs
Definition: soplex.h:1380
Safe arrays of data objects.Class DataArray provides safe arrays of Data Objects. For general C++ obj...
Definition: dataarray.h:63
virtual void scale(SPxLP &lp)=0
scale SPxLP.
bool GE(Real a, Real b, Real eps=Param::epsilon())
returns true iff a >= b + eps
Definition: spxdefines.h:401
SolRational _workSol
Definition: soplex.h:1436
void reDim(int newdim, const bool setZero=true)
Resets DVectorBase&#39;s dimension to newdim.
Definition: dvectorbase.h:249
Rational & addProduct(const Rational &r, const Rational &s)
add product of two rationals
Definition: rational.cpp:3372
void _transformEquality()
introduces slack variables to transform inequality constraints into equations for both rational and r...
DVectorBase< R > _dualFarkas
Definition: solbase.h:227
Timer * syncTime
time for synchronization between real and rational LP (included in solving time)
Definition: statistics.h:92
Basis is not known to be dual nor primal feasible.
Definition: spxbasis.h:94
int feasRefinements
number of refinement steps during infeasibility test
Definition: statistics.h:112
Implementation of Sparse Linear Solver with Rational precision.This class implements a SLinSolverRati...
virtual void getBasis(SPxSolver::VarStatus[], SPxSolver::VarStatus[], const int rowsSize=-1, const int colsSize=-1) const =0
get optimal basis.
virtual void unsimplify(const Vector &, const Vector &, const Vector &, const Vector &, const SPxSolver::VarStatus[], const SPxSolver::VarStatus[])
reconstructs an optimal solution for the unsimplified LP.
const VectorBase< R > & maxObj() const
Returns objective vector for maximization problem.
Definition: spxlpbase.h:384
The time limit has been hit.
upper limit on objective value
Definition: soplex.h:1172
virtual void changeRhs(const VectorBase< R > &newRhs)
Changes right hand side vector for constraints to newRhs.
Definition: spxlpbase.h:1371
void _transformUnbounded()
transforms LP to unboundedness problem by moving the objective function to the constraints, changing right-hand side and bounds to zero, and adding an auxiliary variable for the decrease in the objective function
general zero tolerance
Definition: soplex.h:1151
LPRowBase< Rational > LPRowRational
Definition: lprow.h:29
Result
Result of the simplification.
Definition: spxsimplifier.h:81
type of ratio test
Definition: soplex.h:891
bool setBoolParam(const BoolParam param, const bool value, const bool quiet=false, const bool init=false)
sets boolean parameter value; returns true on success
Definition: soplex.cpp:4896
virtual void unscaleSlacks(Vector &s) const
unscale dense slack vector given in s.
Definition: spxscaler.cpp:279
the problem was so much simplified that it vanished
Definition: spxsimplifier.h:87
continue iterative refinement with exact basic solution if not optimal?
Definition: soplex.h:833
void setBasis(const VarStatus rows[], const VarStatus cols[])
set the lp solver&#39;s basis.
Definition: spxsolver.cpp:1795
virtual Result simplify(SPxLP &lp, Real eps, Real delta)=0
simplify SPxLP lp with identical primal and dual feasibility tolerance.
DVectorRational _feasLhs
Definition: soplex.h:1383
virtual void changeObj(const VectorBase< R > &newObj)
Changes objective vector to newObj.
Definition: spxlpbase.h:1182
type of computational form, i.e., column or row representation
Definition: soplex.h:852
int stallRefinements
number of refinement steps without pivots
Definition: statistics.h:110
T * get_ptr()
get a C pointer to the data.
Definition: dataarray.h:110
virtual void removeRowRange(int start, int end, int perm[]=0)
Removes rows from start to end (including both).
Definition: spxlpbase.h:883
bool _isConsistent() const
checks consistency
Definition: soplex.cpp:6344
int num() const
Returns the number of LPColBases currently in LPColSetBase.
Definition: lpcolsetbase.h:80
SPxLPRational * _rationalLP
Definition: soplex.h:1374
SPxStatus status() const
returns current SPxStatus.
Definition: spxbasis.h:419
static const Rational ZERO
rational zero
Definition: rational.h:63
SPxSolver::Status _solveRealStable(bool acceptUnbounded, bool acceptInfeasible, VectorReal &primal, VectorReal &dual, DataArray< SPxSolver::VarStatus > &basisStatusRows, DataArray< SPxSolver::VarStatus > &basisStatusCols, bool &returnedBasis, const bool forceNoSimplifier=false)
solves real LP with recovery mechanism
bool hasPrimal() const
is a primal feasible solution available?
Definition: solbase.h:51
virtual Status getRedCost(Vector &vector) const
get vector of reduced costs.
Definition: spxsolve.cpp:1299
virtual bool isUnsimplified() const
specifies whether an optimal solution has already been unsimplified.
time limit in seconds (INFTY if unlimited)
Definition: soplex.h:1166
const VectorRational & rhsRational() const
returns right-hand side vector
Definition: soplex.cpp:1084
virtual Real time() const =0
Rational _rationalFeastol
Definition: soplex.h:1323
DVectorBase< R > _primalRay
Definition: solbase.h:224
modified Harris ratio test
Definition: soplex.h:1060
DVectorBase< R > _redCost
Definition: solbase.h:226
virtual void clearRowObjs()
Definition: spxsolver.h:837
primal feasibility tolerance
Definition: soplex.h:1145
void setType(Type tp)
set LEAVE or ENTER algorithm.
Definition: spxsolver.cpp:169
int numRowsReal() const
message handler
Definition: soplex.cpp:770
bound flipping ratio test for long steps in the dual simplex
Definition: soplex.h:1063
DVectorRational _feasRhs
Definition: soplex.h:1384
DVectorRational _feasLower
Definition: soplex.h:1385
solve() aborted due to iteration limit.
Definition: spxsolver.h:199
SPxScaler * _scaler
Definition: soplex.h:1356
const SVectorBase< R > & colVector(int i) const
Returns column vector of column i.
Definition: spxlpbase.h:341
virtual void addCols(const LPColSetBase< R > &pset)
Definition: spxlpbase.h:702
mode for synchronizing real and rational LP
Definition: soplex.h:894
void setValue(int i, R x)
Sets i &#39;th element to x.
Definition: ssvectorbase.h:218
virtual void changeLhs(const VectorBase< R > &newLhs)
Changes left hand side vector for constraints to newLhs.
Definition: spxlpbase.h:1342
int size() const
Number of used indices.
Definition: svectorbase.h:152
void _transformFeasibility()
transforms LP to feasibility problem by removing the objective function, shifting variables...
bool hasDualFarkas() const
is a dual farkas ray available?
Definition: solbase.h:114
Abstract ratio test base class.Class SPxRatioTester is the virtual base class for computing the ratio...
should cycling solutions be accepted during iterative refinement?
Definition: soplex.h:824
void _disableSimplifierAndScaler()
disables simplifier and scaler
Definition: soplex.cpp:7060
const VectorBase< R > & lower() const
Returns lower bound vector.
Definition: spxlpbase.h:420
DataArray< SPxSolver::VarStatus > _storedBasisStatusCols
Definition: soplex.h:1394
void add(const LPColBase< R > &pcol)
Definition: lpcolsetbase.h:253
unsigned int _hasPrimalRay
Definition: solbase.h:233
Implementation of Sparse Linear Solver with Rational precision.
void _performFeasIRStable(SolRational &sol, bool &withDualFarkas, bool &stopped, bool &stoppedIter, bool &error)
performs iterative refinement on the auxiliary problem for testing feasibility
virtual void changeRhs(const Vector &newRhs)
void _storeBasis()
store basis
Rational & subProduct(const Rational &r, const Rational &s)
subtract product of two rationals
Definition: rational.cpp:3381
Rational & invert()
inversion
Definition: rational.cpp:3408
bool _upperFinite(const RangeType &rangeType) const
checks whether RangeType corresponds to finite upper bound
Definition: soplex.cpp:6450
int refinements
number of refinement steps
Definition: statistics.h:109
variable fixed to identical bounds.
Definition: spxsolver.h:181
DVectorReal _manualObj
Definition: soplex.h:1365
objective sense
Definition: soplex.h:849
virtual void unscaleDual(Vector &pi) const
unscale dense dual solution vector given in pi.
Definition: spxscaler.cpp:296
int numColsRational() const
returns number of columns
Definition: soplex.cpp:1021
virtual void removeCol(int i)
Removes i &#39;th column.
Definition: spxlpbase.h:910
LP has been proven to be primal infeasible.
Definition: spxsolver.h:208
R & value(int n)
Reference to value of n &#39;th nonzero.
Definition: svectorbase.h:254
SPxSolver::Status _solveRealForRational(bool fromscratch, VectorReal &primal, VectorReal &dual, DataArray< SPxSolver::VarStatus > &basisStatusRows, DataArray< SPxSolver::VarStatus > &basisStatusCols, bool &returnedBasis)
solves real LP during iterative refinement
void _untransformUnbounded(SolRational &sol, bool unbounded)
undoes transformation to unboundedness problem
bool setIntParam(const IntParam param, const int value, const bool quiet=false, const bool init=false)
sets integer parameter value; returns true on success
Definition: soplex.cpp:4939
Ids for LP columns.Class SPxColId provides DataKeys for the column indices of an SPxLP.
Definition: spxid.h:36
void setOpttol(Real d)
set parameter opttol.
Definition: spxsolver.cpp:916
minimum number of stalling refinements since last pivot to trigger rational factorization ...
Definition: soplex.h:912
int sign(const Rational &r)
Sign function; returns 1 if r > 0, 0 if r = 0, and -1 if r < 0.
Definition: rational.cpp:3938
Wrapper for GMP type mpq_class.We wrap mpq_class so that we can replace it by a double type if GMP is...
Definition: rational.h:45
automatic choice
Definition: soplex.h:990
virtual void changeRowObj(const Vector &newObj)
const VectorBase< R > & upper() const
Returns upper bound vector.
Definition: spxlpbase.h:402
should a rational factorization be performed after iterative refinement?
Definition: soplex.h:821
const VectorReal & upperReal() const
returns upper bound vector
Definition: soplex.cpp:914
Real opttol() const
allowed optimality, i.e., dual feasibility tolerance.
Definition: spxsolver.h:697
LPColBase< Real > LPColReal
Definition: lpcol.h:30
Rational _rationalPosInfty
Definition: soplex.h:1321
static const Rational NEGONE
rational minus one
Definition: rational.h:65
Rational _rationalOpttol
Definition: soplex.h:1324
bool _isSolveStopped(bool &stoppedTime, bool &stoppedIter) const
should solving process be stopped?
Definition: soplex.cpp:6364
Real realParam(const RealParam param) const
returns real parameter value
Definition: soplex.cpp:4866
steepest edge pricer with exact initialization of norms
Definition: soplex.h:1047
virtual void changeLower(const Vector &newLower)
virtual void changeUpper(const VectorBase< R > &newUpper)
Changes vector of upper bounds to newUpper.
Definition: spxlpbase.h:1278
void add(const SVectorBase< S > &vec)
Append nonzeros of sv.
Definition: dsvectorbase.h:216
DVectorRational _modLhs
Definition: soplex.h:1389
DVectorReal _manualLower
Definition: soplex.h:1361
virtual void start()=0
start timer, resume accounting user, system and real time.
Entering Simplex.
Definition: spxsolver.h:134
virtual Real stop()=0
stop timer, return accounted user time.
void spx_alloc(T &p, int n=1)
Allocate memory.
Definition: spxalloc.h:48
simplification could be done
Definition: spxsimplifier.h:83
const UnitVectorRational * _unitVectorRational(const int i)
returns pointer to a constant unit vector available until destruction of the SoPlex class ...
Definition: soplex.cpp:7376
virtual void addCol(const LPColBase< R > &col)
Definition: spxlpbase.h:650
const SPxRatioTester * ratiotester() const
return loaded SPxRatioTester.
Definition: spxsolver.h:1636
bool _lowerFinite(const RangeType &rangeType) const
checks whether RangeType corresponds to finite lower bound
Definition: soplex.cpp:6442
LP is primal infeasible or unbounded.
Definition: spxsolver.h:209
Rational reconstruction of solution vector.
Leaving Simplex.
Definition: spxsolver.h:143
virtual void changeLower(const VectorBase< R > &newLower)
Changes vector of lower bounds to newLower.
Definition: spxlpbase.h:1247
virtual void changeRange(const VectorBase< R > &newLhs, const VectorBase< R > &newRhs)
Changes left and right hand side vectors.
Definition: spxlpbase.h:1392
DataArray< SPxSolver::VarStatus > _basisStatusCols
Definition: soplex.h:1432
DVectorRational _unboundedLower
Definition: soplex.h:1377
bool reconstructVector(VectorRational &input, const Rational &denomBoundSquared, const DIdxSet *indexSet)
Definition: ratrecon.cpp:228
nothing known about basis status (possibly due to a singular basis in transformed problem) ...
Definition: spxsolver.h:184
bool boolParam(const BoolParam param) const
returns boolean parameter value
Definition: soplex.cpp:4846
DVectorReal _manualRhs
Definition: soplex.h:1364
DataArray< RangeType > _colTypes
Definition: soplex.h:1419
type of scaler
Definition: soplex.h:882
double Real
SOPLEX_DEBUG.
Definition: spxdefines.h:200
#define MSG_INFO2(spxout, x)
Prints out message x if the verbosity level is at least SPxOut::INFO2.
Definition: spxdefines.h:115
bool _isRealLPLoaded
Definition: soplex.h:1359
virtual Status getDual(Vector &vector) const
get current solution vector for dual variables.
Definition: spxsolve.cpp:1268
void _solveRealLPAndRecordStatistics()
call floating-point solver and update statistics on iterations etc.
Definition: soplex.cpp:7109
DVectorBase< R > _primal
Definition: solbase.h:222
int index(int n) const
Returns index of the n &#39;th nonzero element.
Definition: ssvectorbase.h:174
#define MSG_DEBUG(x)
Definition: spxdefines.h:127
void setMarkowitz(Real m)
sets minimum Markowitz threshold.
Definition: slufactor.h:128
lower threshold in lifting (nonzero matrix coefficients with smaller absolute value will be reformula...
Definition: soplex.h:1184
int & index(int n)
Reference to index of n &#39;th nonzero.
Definition: svectorbase.h:236
const VectorReal & rhsReal() const
returns right-hand side vector
Definition: soplex.cpp:842
Timer * rationalTime
time for rational LP solving (included in solving time)
Definition: statistics.h:94
int _beforeLiftCols
Definition: soplex.h:1398
void _performOptIRStable(SolRational &sol, bool acceptUnbounded, bool acceptInfeasible, int minRounds, bool &primalFeasible, bool &dualFeasible, bool &infeasible, bool &unbounded, bool &stopped, bool &stoppedIter, bool &error)
solves current problem with iterative refinement and recovery mechanism
dual simplex algorithm, i.e., leaving for column and entering for row representation ...
Definition: soplex.h:948
Status getBasis(VarStatus rows[], VarStatus cols[], const int rowsSize=-1, const int colsSize=-1) const
get current basis, and return solver status.
Definition: spxsolver.cpp:1722
LP has been solved to optimality.
Definition: spxsolver.h:206
int numRowsRational() const
returns number of rows
Definition: soplex.cpp:1012
void setTimeLimit(const Real limit)
set time limit on factorization
const VectorRational & lowerRational() const
returns lower bound vector
Definition: soplex.cpp:1174
bool _hasBasis
Definition: soplex.h:1438
unsigned int _hasDualFarkas
Definition: solbase.h:235
int nCols() const
Returns number of columns in LP.
Definition: spxlpbase.h:133
void append(const T &t)
append element t.
Definition: dataarray.h:121
LPRowBase< Real > LPRowReal
Definition: lprow.h:28
void _restoreLPReal()
restores objective, bounds, and sides of real LP
DVectorBase< Rational > DVectorRational
Definition: dvector.h:30
void _project(SolRational &sol)
undoes lifting
const SPxBasis & basis() const
Return current basis.
Definition: spxsolver.h:1616
Class for collecting statistical information.
void clearNum(int n)
Sets n &#39;th nonzero element to 0 (index n must exist).
Definition: ssvectorbase.h:258
solve() aborted due to time limit.
Definition: spxsolver.h:198
void _performUnboundedIRStable(SolRational &sol, bool &hasUnboundedRay, bool &stopped, bool &stoppedIter, bool &error)
performs iterative refinement on the auxiliary problem for testing unboundedness
DVectorRational _feasObj
Definition: soplex.h:1382
an error occured.
Definition: spxsolver.h:192
void addIdx(int i)
adds index i to the index set
Definition: didxset.h:75
round scaling factors for iterative refinement to powers of two?
Definition: soplex.h:830
int size() const
return nr. of elements.
Definition: dataarray.h:211
void _computeInfeasBox(SolRational &sol, bool transformed)
Statistics * _statistics
statistics since last call to solveReal() or solveRational()
Definition: soplex.h:1308
Real getSolveTime() const
time spent in solves
void _solveRational()
solves rational LP
LPColSetRational _slackCols
Definition: soplex.h:1376
type of simplifier
Definition: soplex.h:879
virtual void changeLhs(const Vector &newLhs)
RangeType _switchRangeType(const RangeType &rangeType) const
switches RANGETYPE_LOWER to RANGETYPE_UPPER and vice versa
Definition: soplex.cpp:6429
void solveRight(VectorRational &x, const VectorRational &b)
Solves .
const SPxPricer * pricer() const
return loaded SPxPricer.
Definition: spxsolver.h:1626
DVectorRational _feasUpper
Definition: soplex.h:1386
#define MSG_ERROR(x)
Prints out message x if the verbosity level is at least SPxOut::ERROR.
Definition: spxdefines.h:109
void clear()
Clears vector.
Definition: ssvectorbase.h:572
int getFactorCount() const
number of factorizations performed
Abstract pricer base class.Class SPxPricer is a pure virtual class defining the interface for pricer ...
Definition: spxpricer.h:46
bool _reconstructSolutionRational(SolRational &sol, DataArray< SPxSolver::VarStatus > &basisStatusRows, DataArray< SPxSolver::VarStatus > &basisStatusCols, const Rational &denomBoundSquared)
attempts rational reconstruction of primal-dual solution
virtual void unscalePrimal(Vector &x) const
unscale dense primal solution vector given in x.
Definition: spxscaler.cpp:262
SLUFactor _slufactor
Definition: soplex.h:1334
variable set to its upper bound.
Definition: spxsolver.h:179
void setEpsilon(R eps)
Changes the non-zero epsilon, invalidating the setup. */.
Definition: ssvectorbase.h:110
Dynamic index set.Class DIdxSet provides dynamic IdxSet in the sense, that no restrictions are posed ...
Definition: didxset.h:42
void getObj(VectorBase< R > &pobj) const
Gets objective vector.
Definition: spxlpbase.h:353
#define MSG_INFO1(spxout, x)
Prints out message x if the verbosity level is at least SPxOut::INFO1.
Definition: spxdefines.h:113
int dim() const
Dimension of vector.
Definition: vectorbase.h:174
Real luSolveTimeRational
time for solving linear systems in rational precision
Definition: statistics.h:100
Type type() const
return current Type.
Definition: spxsolver.h:412
primal infeasibility was detected
Definition: spxsimplifier.h:84
DVectorRational _modUpper
Definition: soplex.h:1388
user sync of real and rational LP
Definition: soplex.h:1076
void _ensureDSVectorRationalMemory(DSVectorRational &vec, const int newmax) const
extends sparse vector to hold newmax entries if and only if it holds no more free entries ...
Definition: soplex.cpp:6301
SPxLPBase< Rational > SPxLPRational
Definition: spxlp.h:37
SPxOut spxout
Definition: soplex.h:1215
infinity threshold
Definition: soplex.h:1163
algorithm is running
Definition: spxsolver.h:204
SPxDefaultRT _ratiotesterTextbook
Definition: soplex.h:1349
Status load(const SVectorRational *vec[], int dim)
virtual void changeBounds(const VectorBase< R > &newLower, const VectorBase< R > &newUpper)
Changes variable bounds to newLower and newUpper.
Definition: spxlpbase.h:1309
virtual const Vector & unsimplifiedPrimal()=0
returns a reference to the unsimplified primal solution.
variable set to its lower bound.
Definition: spxsolver.h:180
LPColBase< Rational > LPColRational
Definition: lpcol.h:31
Preconfigured SoPlex LP solver.
const SVectorRational & colVectorRational(int i) const
returns vector of column i
Definition: soplex.cpp:1147
static const Rational POSONE
rational plus one
Definition: rational.h:64
void add(int i, const R &v)
Append one nonzero (i,v).
Definition: svectorbase.h:272
Real markowitz()
returns Markowitz threshold.
Definition: slufactor.h:141
void reDim(int newdim)
Resets dimension to newdim.
Definition: ssvectorbase.h:538
bool isSetup() const
Returns setup status.
Definition: ssvectorbase.h:120
int nRows() const
Returns number of rows in LP.
Definition: spxlpbase.h:127
int _beforeLiftRows
Definition: soplex.h:1397
virtual void unscaleRedCost(Vector &r) const
unscale dense reduced cost vector given in r.
Definition: spxscaler.cpp:313
lower and upper bound finite, but different
Definition: soplex.h:1413
int luFactorizationsRational
number of basis matrix factorizations in rational precision
Definition: statistics.h:107
DVectorReal _manualLhs
Definition: soplex.h:1363
Timer * preprocessingTime
preprocessing time
Definition: statistics.h:90
equilibrium scaling on rows and columns
Definition: soplex.h:1003
bool _storedBasis
Definition: soplex.h:1396
void _storeLPReal()
stores objective, bounds, and sides of real LP
Rational _rationalNegInfty
Definition: soplex.h:1322
const VectorReal & lhsReal() const
returns left-hand side vector
Definition: soplex.cpp:860
VectorBase< R > & multAdd(const S &x, const VectorBase< T > &vec)
Addition of scaled vector.
Definition: vectorbase.h:369
lower bound equals upper bound
Definition: soplex.h:1416
column representation Ax - s = 0, lower <= x <= upper, lhs <= s <= rhs
Definition: soplex.h:935
Timer * reconstructionTime
time for rational reconstructions
Definition: statistics.h:95
const VectorBase< R > & lhs() const
Returns left hand side vector.
Definition: spxlpbase.h:242
Everything should be within this namespace.
DVectorReal _manualUpper
Definition: soplex.h:1362
working tolerance for feasibility in floating-point solver during iterative refinement ...
Definition: soplex.h:1175
SPxLPReal * _realLP
Definition: soplex.h:1354
SPxDevexPR _pricerDevex
Definition: soplex.h:1346
void clear()
Remove all indices.
Definition: svectorbase.h:396
DataArray< RangeType > _rowTypes
Definition: soplex.h:1420
DVectorRational _modObj
Definition: soplex.h:1391
bool hasPrimalRay() const
is a primal unbounded ray available?
Definition: solbase.h:75
apply rational reconstruction after each iterative refinement?
Definition: soplex.h:827
solve() aborted due to detection of cycling.
Definition: spxsolver.h:197
type of pricer
Definition: soplex.h:888
Saving LPs in a form suitable for SoPlex.Class SPxLPBase provides the data structures required for sa...
Definition: spxlpbase.h:76
DSVectorBase< Real > DSVectorReal
Definition: dsvector.h:29
primal unboundedness was detected
Definition: spxsimplifier.h:86
working tolerance for optimality in floating-point solver during iterative refinement ...
Definition: soplex.h:1178
#define MSG_WARNING(spxout, x)
Prints out message x if the verbosity level is at least SPxOut::WARNING.
Definition: spxdefines.h:111
Real spxNextafter(Real x, Real y)
Definition: spxdefines.h:332
int iterations
number of iterations/pivots
Definition: statistics.h:101
Real feastol() const
allowed primal feasibility tolerance.
Definition: spxsolver.h:689
DVectorRational _unboundedLhs
Definition: soplex.h:1379
Real getFactorTime() const
time spent in factorizations
DSVectorRational _primalDualDiff
Definition: soplex.h:1392
void setup()
Initializes nonzero indices for elements with absolute values above epsilon and sets all other elemen...
Definition: ssvectorbase.h:132
DVectorRational _unboundedUpper
Definition: soplex.h:1378
SPxSolver _solver
Definition: soplex.h:1333
nothing known on loaded problem.
Definition: spxsolver.h:205
variable is basic.
Definition: spxsolver.h:183
type of algorithm, i.e., primal or dual
Definition: soplex.h:855
virtual Status getSlacks(Vector &vector) const
get vector of slack variables.
Definition: spxsolve.cpp:1378
const VectorRational & lhsRational() const
returns left-hand side vector
Definition: soplex.cpp:1102
Status status() const
Status of solution process.
Definition: spxsolve.cpp:1536
DataArray< SPxSolver::VarStatus > _basisStatusRows
Definition: soplex.h:1431
virtual void changeElement(int i, int j, const R &val)
Changes LP element (i, j) to val.
Definition: spxlpbase.h:1532
virtual void changeMaxObj(const VectorBase< R > &newObj)
Changes objective vector to newObj.
Definition: spxlpbase.h:1218
Rational objRational(int i) const
returns objective value of column i
Definition: soplex.cpp:1211
virtual const Vector & unsimplifiedDual()=0
returns a reference to the unsimplified dual solution.
should dual infeasibility be tested in order to try to return a dual solution even if primal infeasib...
Definition: soplex.h:818
void _factorizeColumnRational(SolRational &sol, DataArray< SPxSolver::VarStatus > &basisStatusRows, DataArray< SPxSolver::VarStatus > &basisStatusCols, bool &stoppedTime, bool &stoppedIter, bool &error, bool &optimal)
factorizes rational basis matrix in column representation
unsigned int _hasPrimal
Definition: solbase.h:232
Timer * solvingTime
solving time
Definition: statistics.h:89
void dump() const
Prints out status.
Definition: spxdesc.cpp:113
int intParam(const IntParam param) const
returns integer parameter value
Definition: soplex.cpp:4856
should lifting be used to reduce range of nonzero matrix coefficients?
Definition: soplex.h:812
const Real infinity
Definition: spxdefines.cpp:26
virtual void removeColRange(int start, int end, int perm[]=0)
Removes columns from start to end (including both).
Definition: spxlpbase.h:980
const VectorBase< R > & rhs() const
Returns right hand side vector.
Definition: spxlpbase.h:224
lower limit on objective value
Definition: soplex.h:1169
DSVectorRational _tauColVector
Definition: soplex.h:1381
DVectorRational _modLower
Definition: soplex.h:1387
const SVectorBase< R > & rowVector(int i) const
Gets row vector of row i.
Definition: spxlpbase.h:212
virtual Status getDualfarkas(Vector &vector) const
get dual farkas proof of infeasibility.
Definition: spxsolve.cpp:1360
void setFeastol(Real d)
set parameter feastol.
Definition: spxsolver.cpp:904
virtual void addPrimalActivity(const SVectorBase< R > &primal, VectorBase< R > &activity) const
Updates activity of the rows for a given primal vector; activity does not need to be zero...
Definition: spxlpbase.h:1655
virtual void loadLP(const SPxLP &LP)
copy LP.
Definition: spxsolver.cpp:68
unsigned int _hasDual
Definition: solbase.h:234
virtual void addRow(const LPRowBase< R > &row)
Definition: spxlpbase.h:488
const Desc & desc() const
Definition: spxbasis.h:457
Timer * transformTime
time for transforming LPs (included in solving time)
Definition: statistics.h:93
void invalidate()
invalidate solution
Definition: solbase.h:213
void _restoreBasis()
restore basis
SolRational _solRational
Definition: soplex.h:1435
std::string rationalToString(const Rational &r, const int precision)
convert rational number to string
Definition: rational.cpp:3468
virtual void addDualActivity(const SVectorBase< R > &dual, VectorBase< R > &activity) const
Updates "dual" activity of the columns for a given dual vector, i.e., y^T A; activity does not need t...
Definition: spxlpbase.h:1709
geometric frequency at which to apply rational reconstruction
Definition: soplex.h:1196
#define MSG_INFO3(spxout, x)
Prints out message x if the verbosity level is at least SPxOut::INFO3.
Definition: spxdefines.h:117
bool _hasSolRational
Definition: soplex.h:1440
void _enableSimplifierAndScaler()
enables simplifier and scaler according to current parameters
Definition: soplex.cpp:7017
primal simplex algorithm, i.e., entering for column and leaving for row representation ...
Definition: soplex.h:945
DVectorRational _modRhs
Definition: soplex.h:1390
should LP be transformed to equality form before a rational solve?
Definition: soplex.h:815
R value(int n) const
Returns value of the n &#39;th nonzero element.
Definition: ssvectorbase.h:182
DVectorBase< R > _dual
Definition: solbase.h:225
virtual void setTerminationValue(Real value=infinity)
set objective limit.
Definition: spxsolver.cpp:1542
void clear()
Removes all LPColBases from the set.
Definition: lpcolsetbase.h:485
void solveLeft(VectorRational &x, const VectorRational &b)
Solves .
const VectorRational & upperRational() const
returns upper bound vector
Definition: soplex.cpp:1156
virtual void writeFile(const char *filename, const NameSet *rowNames=0, const NameSet *colNames=0, const DIdxSet *p_intvars=0) const
Write loaded LP to filename.
Definition: spxlpbase.h:1080
virtual void reLoad()
reload LP.
Definition: spxsolver.cpp:55
upper threshold in lifting (nonzero matrix coefficients with larger absolute value will be reformulat...
Definition: soplex.h:1187
DataArray< SPxSolver::VarStatus > _storedBasisStatusRows
Definition: soplex.h:1393
SPxLPReal _manualRealLP
Definition: soplex.h:1366
void resetCounters()
reset timers and counters
bool hasDual() const
is a dual solution available?
Definition: solbase.h:90
void reSize(int newsize)
reset size to newsize.
Definition: dataarray.h:223
Real luFactorizationTimeRational
time for factorizing bases matrices in rational precision
Definition: statistics.h:99
void _lift()
reduces matrix coefficient in absolute value by the lifting procedure of Thiele et al...
force iterative refinement
Definition: soplex.h:1099
SPxSolver::Status _status
Definition: soplex.h:1428
virtual void computePrimalActivity(const VectorBase< R > &primal, VectorBase< R > &activity) const
Computes activity of the rows for a given primal vector; activity does not need to be zero...
Definition: spxlpbase.h:1619
int rationalReconstructions
number of rational reconstructions performed
Definition: statistics.h:108
solve() aborted due to objective limit.
Definition: spxsolver.h:200
void spx_free(T &p)
Release memory.
Definition: spxalloc.h:109
SPxSimplifier * _simplifier
Definition: soplex.h:1355
const VectorBase< R > & maxRowObj() const
Definition: spxlpbase.h:286
void _untransformEquality(SolRational &sol)
undoes transformation to equality form
int pivotRefinements
number of refinement steps until final basis is reached
Definition: statistics.h:111
int size() const
Returns the number of nonzeros.
Definition: ssvectorbase.h:199
Basis is singular, numerical troubles?
Definition: spxsolver.h:201
const SVectorBase< R > & colVector(int i) const
Returns colVector of i &#39;th LPColBase in LPColSetBase.
Definition: lpcolsetbase.h:206
virtual void factorize()
Factorize basis matrix.
Definition: spxsolver.cpp:525
upper bound is finite, lower bound is infinite
Definition: soplex.h:1410
LP has a usable Basis (maybe LP is changed).
Definition: spxsolver.h:203
virtual void changeUpper(const Vector &newUpper)
dual infeasibility was detected
Definition: spxsimplifier.h:85
void _untransformFeasibility(SolRational &sol, bool infeasible)
undoes transformation to feasibility problem
LP has been proven to be primal unbounded.
Definition: spxsolver.h:207
No Problem has been loaded to the basis.
Definition: spxbasis.h:92
void setDelta(Real d)
set parameter delta, i.e., set feastol and opttol to same value.
Definition: spxsolver.cpp:928
virtual void changeObj(const Vector &newObj)
Rational & powRound()
round up to next power of two
Definition: rational.cpp:3417