All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
solverational.cpp
Go to the documentation of this file.
56 _solver.getBasis(_basisStatusRows.get_ptr(), _basisStatusCols.get_ptr(), _basisStatusRows.size(), _basisStatusCols.size());
59 // store objective, bounds, and sides of real LP in case they will be modified during iterative refinement
63 if( realParam(SoPlex::OBJLIMIT_LOWER) > -realParam(SoPlex::INFTY) || realParam(SoPlex::OBJLIMIT_UPPER) < realParam(SoPlex::INFTY) )
157 ///@todo this should be stored already earlier, possible switch use solRational above and solFeas here
292 if( _status == SPxSolver::OPTIMAL || _status == SPxSolver::INFEASIBLE || _status == SPxSolver::UNBOUNDED )
311 // restore objective, bounds, and sides of real LP in case they have been modified during iterative refinement
314 // since the real LP is loaded in the solver, we need to also pass the basis information to the solver if
330 void SoPlex::_performOptIRStable(SolRational& sol, bool acceptUnbounded, bool acceptInfeasible, int minRounds, bool& primalFeasible, bool& dualFeasible, bool& infeasible, bool& unbounded, bool& stopped, bool& error)
383 result = _solveRealStable(acceptUnbounded, acceptInfeasible, primalReal, dualReal, _basisStatusRows, _basisStatusCols, _hasBasis);
393 // the floating-point solve returns a Farkas ray if and only if the simplifier was not used, which is exactly
419 // store floating-point solution of original LP as current rational solution and ensure that solution vectors have
438 // it may happen that lower and upper are only equal in the real LP but different in the rational LP; we do
439 // not check this to avoid rational comparisons, but simply switch the basis status to the lower bound; this
456 // it may happen that left-hand and right-hand side are different in the rational, but equal in the real LP,
457 // leading to a fixed basis status; this is critical because rows with fixed basis status are ignored in the
458 // computation of the dual violation; to avoid rational comparisons we do not check this but simply switch to
469 // we assume that the objective function vector has less nonzeros than the reduced cost vector, and so multiplying
470 // with -1 first and subtracting the dual activity should be faster than adding the dual activity and negating
493 const int maxDimRational = numColsRational() > numRowsRational() ? numColsRational() : numRowsRational();
573 // if the activity is feasible, but too far from the bound, this violates complementary slackness; we
596 // if the activity is feasible, but too far from the bound, this violates complementary slackness; we
615 if( ((maximizing && basisStatusCol != SPxSolver::ON_LOWER) || (!maximizing && basisStatusCol != SPxSolver::ON_UPPER))
626 if( ((maximizing && basisStatusCol != SPxSolver::ON_UPPER) || (!maximizing && basisStatusCol != SPxSolver::ON_LOWER))
648 if( ((maximizing && basisStatusRow != SPxSolver::ON_LOWER) || (!maximizing && basisStatusRow != SPxSolver::ON_UPPER))
659 if( ((maximizing && basisStatusRow != SPxSolver::ON_UPPER) || (!maximizing && basisStatusRow != SPxSolver::ON_LOWER))
673 // output violations; the reduced cost violations for artificially introduced slack columns are actually violations of the dual multipliers
687 << std::setw(10) << rationalToString(boundsViolation > sideViolation ? boundsViolation : sideViolation, 2) << " & "
688 << std::setw(10) << rationalToString(redCostViolation > dualViolation ? redCostViolation : dualViolation, 2) << "\n");
702 MSG_INFO1( spxout, spxout << "Tolerances reached but minRounds forcing additional refinement rounds.\n" );
732 && lastStallRefinements >= intParam(SoPlex::RATFAC_MINSTALLS) && _hasBasis && factorSolNewBasis;
744 if( _reconstructSolutionRational(sol, _basisStatusRows, _basisStatusCols, stopped, error, maxViolation) )
752 MSG_DEBUG( spxout << "Next rational reconstruction after refinement " << nextRatrecRefinement << ".\n" );
790 // compute primal scaling factor; limit increase in scaling by tolerance used in floating point solve
904 // compute dual scaling factor; limit increase in scaling by tolerance used in floating point solve
967 // ensure that artificial slack columns are basic and inequality constraints are nonbasic; otherwise we may end
968 // up with dual violation on inequality constraints after removing the slack columns; do not change this in the
969 // floating-point solver, though, because the solver may require its original basis to detect optimality
993 MSG_DEBUG( spxout << "basis (status = " << _solver.basis().status() << ") desc before set:\n"; _solver.basis().desc().dump() );
995 MSG_DEBUG( spxout << "basis (status = " << _solver.basis().status() << ") desc after set:\n"; _solver.basis().desc().dump() );
998 MSG_DEBUG( spxout << "setting basis in solver " << (_hasBasis ? "successful" : "failed") << " (3)\n" );
1004 result = _solveRealStable(acceptUnbounded, acceptInfeasible, primalReal, dualReal, _basisStatusRows, _basisStatusCols, _hasBasis, primalScale > 1e20 || dualScale > 1e20);
1088 // it may happen that lower and upper are only equal in the real LP but different in the rational LP; we
1089 // do not check this to avoid rational comparisons, but simply switch the basis status to the lower
1184 if( ((maximizing && basisStatusCol != SPxSolver::ON_LOWER) || (!maximizing && basisStatusCol != SPxSolver::ON_UPPER))
1196 if( ((maximizing && basisStatusCol != SPxSolver::ON_UPPER) || (!maximizing && basisStatusCol != SPxSolver::ON_LOWER))
1222 if( ((maximizing && basisStatusRow != SPxSolver::ON_LOWER) || (!maximizing && basisStatusRow != SPxSolver::ON_UPPER))
1234 if( ((maximizing && basisStatusRow != SPxSolver::ON_UPPER) || (!maximizing && basisStatusRow != SPxSolver::ON_LOWER))
1258 if( debugRedCostViolation > _solver.opttol() || debugDualViolation > _solver.opttol() || debugBasicDualViolation > 1e-9 )
1278 // it may happen that left-hand and right-hand side are different in the rational, but equal in the real LP,
1279 // leading to a fixed basis status; this is critical because rows with fixed basis status are ignored in the
1280 // computation of the dual violation; to avoid rational comparisons we do not check this but simply switch
1306 // update or recompute reduced cost values depending on which looks faster; adding one to the length of the
1321 // we assume that the objective function vector has less nonzeros than the reduced cost vector, and so multiplying
1322 // with -1 first and subtracting the dual activity should be faster than adding the dual activity and negating
1331 MSG_INFO2( spxout, spxout << "Corrected " << numCorrectedPrimals << " primal variables and " << numCorrectedDuals << " dual values.\n" );
1369 void SoPlex::_performUnboundedIRStable(SolRational& sol, bool& hasUnboundedRay, bool& stopped, bool& error)
1386 _performOptIRStable(sol, false, false, 0, primalFeasible, dualFeasible, infeasible, unbounded, stopped, error);
1415 // because the right-hand side and all bounds (but tau's upper bound) are zero, tau should be approximately
1430 void SoPlex::_performFeasIRStable(SolRational& sol, bool& withDualFarkas, bool& stopped, bool& error)
1440 // if the problem has been found to be infeasible and an approximate Farkas proof is available, we compute a
1441 // scaled unit box around the origin that provably contains no feasible solution; this currently only works for
1460 _performOptIRStable(sol, false, false, 0, primalFeasible, dualFeasible, infeasible, unbounded, stopped, error);
1480 // else we should have either a refined Farkas proof or an approximate feasible solution to the original
1524 /// reduces matrix coefficient in absolute value by the lifting procedure of Thiele et al. 2013
1563 MSG_DEBUG( std::cout << " --> nonzero " << k << " has value " << rationalToString(value) << " in row " << colVector.index(k) << "\n" );
1583 _realLP->changeBounds(liftingColumnIndex, -realParam(SoPlex::INFTY), realParam(SoPlex::INFTY));
1630 MSG_DEBUG( std::cout << " --> nonzero " << k << " has value " << rationalToString(value) << " in row " << colVector.index(k) << "\n" );
1650 _realLP->changeBounds(liftingColumnIndex, -realParam(SoPlex::INFTY), realParam(SoPlex::INFTY));
1694 MSG_INFO1( spxout, spxout << "Added " << numColsRational() - _beforeLiftCols << " columns and "
1733 ///@todo if we know the mapping between original and lifting columns, we simply need to add the reduced cost of
1734 /// the lifting column to the reduced cost of the original column; this is not implemented now, because for
1763 MSG_INFO1( spxout, spxout << "Warning: lost basis during project phase because of nonbasic lifting column.\n" );
1772 MSG_INFO1( spxout, spxout << "Warning: lost basis during project phase because of basic lifting row.\n" );
1831 // representable, the user might have constructed the constraint in the real LP by rounding down the
1832 // left-hand side and rounding up the right-hand side; if the basis status is fixed, we need to adjust it
1867 /// introduces slack variables to transform inequality constraints into equations for both rational and real LP,
1939 MSG_INFO1( spxout, spxout << "Added " << _slackCols.num() << " slack columns to transform rows to equality form.\n" );
1992 assert(_basisStatusRows[row] != SPxSolver::BASIC || _basisStatusCols[col] != SPxSolver::BASIC);
2058 /// transforms LP to unboundedness problem by moving the objective function to the constraints, changing right-hand
2059 /// side and bounds to zero, and adding an auxiliary variable for the decrease in the objective function
2204 _hasBasis = (_basisStatusCols[numOrigCols] != SPxSolver::BASIC && _basisStatusRows[numOrigRows] == SPxSolver::BASIC);
2328 /// transforms LP to feasibility problem by removing the objective function, shifting variables, and homogenizing the
2367 // set objective coefficients to zero; shift primal space such as to guarantee that the zero solution is within
2570 LPColRational((intParam(SoPlex::OBJSENSE) == SoPlex::OBJSENSE_MAXIMIZE ? Rational::POSONE : Rational::NEGONE),
2680 bool shifted = (_lowerFinite(_colTypes[c]) && _feasLower[c] > 0) || (_upperFinite(_colTypes[c]) && _feasUpper[c] < 0);
2756 Given constraints of the form \f$ lhs <= Ax <= rhs \f$, a farkas proof y should satisfy \f$ y^T A = 0 \f$ and
2757 \f$ y_+^T lhs - y_-^T rhs > 0 \f$, where \f$ y_+, y_- \f$ denote the positive and negative parts of \f$ y \f$.
2758 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
2765 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
2766 bounds on \f$ x \f$ imply that all feasible \f$ x \f$ satisfy (*), and if not then compute bounds on \f$ x \f$ to
2775 \f$ B \f$ can be increased by iteratively including variable bounds smaller than \f$ B \f$. The speed of this
2776 method can be further improved by using interval arithmetic for all computations. For related information see
2798 // prepare ytransA and ytransb; since we want exact arithmetic, we set the zero threshold of the semi-sparse
2804 ///@todo this currently works only if all constraints are equations aggregate rows and sides using the multipliers of the Farkas ray
2817 // if we choose minus ytransb as vector of multipliers for the bound constraints on the variables, we obtain an
2818 // exactly feasible dual solution for the LP with zero objective function; we aggregate the bounds of the
2844 MSG_DEBUG( std::cout << "max ytransA*[x_l,x_u] = " << (isTempFinite ? rationalToString(temp) : "infinite") << "\n" );
2846 // ytransb - temp is the increase in the dual objective along the Farkas ray; if this is positive, the dual is
2858 // if ytransb is negative, try to make it zero by including a positive lower bound or a negative upper bound
2878 // if ytransb is still zero then the zero solution is inside the bounds and cannot be cut off by the Farkas
2882 MSG_INFO1( spxout, spxout << "Approximate Farkas proof to weak. Could not compute Farkas box. (1)\n" );
2892 // if the one norm is zero then ytransA is zero the Farkas proof should have been verified above
2895 // initialize variables in loop: size of Farkas box B, flag whether B has been increased, and number of current
2907 // if all nonzeros have been inspected once without increasing B, we abort; otherwise, we start another round
2921 // if the multiplier is positive we inspect the lower bound: if it is finite and within the Farkas box, we can
2938 MSG_INFO1( spxout, spxout << "Approximate Farkas proof to weak. Could not compute Farkas box. (2)\n" );
2956 // if the multiplier is negative we inspect the upper bound: if it is finite and within the Farkas box, we can
2971 MSG_INFO1( spxout, spxout << "Approximate Farkas proof to weak. Could not compute Farkas box. (2)\n" );
2992 // currently this bound cannot be used to increase B; we will check it again in the next round, because B might
3000 MSG_INFO1( spxout, spxout << "Computed Farkas box: provably no feasible solutions with components less than "
3008 SPxSolver::Status SoPlex::_solveRealForRational(bool fromscratch, VectorReal& primal, VectorReal& dual, DataArray< SPxSolver::VarStatus >& basisStatusRows, DataArray< SPxSolver::VarStatus >& basisStatusCols, bool& returnedBasis)
3051 simplificationStatus = _simplifier->simplify(_solver, realParam(SoPlex::EPSILON_ZERO), realParam(SoPlex::FPFEASTOL), realParam(SoPlex::FPOPTTOL), keepbounds);
3076 else if( simplificationStatus == SPxSimplifier::VANISHED || simplificationStatus == SPxSimplifier::OKAY )
3078 result = simplificationStatus == SPxSimplifier::VANISHED ? SPxSolver::OPTIMAL : _solver.status();
3085 // unsimplify if simplifier is active and LP is solved to optimality; this must be done here and not at solution
3090 assert(simplificationStatus == SPxSimplifier::VANISHED || simplificationStatus == SPxSimplifier::OKAY);
3121 _solver.getBasis(basisStatusRows.get_ptr(), basisStatusCols.get_ptr(), basisStatusRows.size(), basisStatusCols.size());
3125 _simplifier->unsimplify(tmpPrimal, tmpDual, tmpSlacks, tmpRedCost, basisStatusRows.get_ptr(), basisStatusCols.get_ptr());
3130 _simplifier->getBasis(basisStatusRows.get_ptr(), basisStatusCols.get_ptr(), basisStatusRows.size(), basisStatusCols.size());
3136 // if the original problem is not in the solver because of scaling, we also need to store the basis
3152 _solver.getBasis(basisStatusRows.get_ptr(), basisStatusCols.get_ptr(), basisStatusRows.size(), basisStatusCols.size());
3173 _solver.getBasis(basisStatusRows.get_ptr(), basisStatusCols.get_ptr(), basisStatusRows.size(), basisStatusCols.size());
3197 // if the original problem is not in the solver because of scaling, we also need to store the basis
3200 _solver.getBasis(basisStatusRows.get_ptr(), basisStatusCols.get_ptr(), basisStatusRows.size(), basisStatusCols.size());
3234 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)
3265 result = _solveRealForRational(fromScratch, primal, dual, basisStatusRows, basisStatusCols, returnedBasis);
3284 && (intParam(SoPlex::SIMPLIFIER) != SoPlex::SIMPLIFIER_OFF || intParam(SoPlex::SCALER) != SoPlex::SCALER_OFF) )
3442 void SoPlex::_factorizeColumnRational(SolRational& sol, DataArray< SPxSolver::VarStatus >& basisStatusRows, DataArray< SPxSolver::VarStatus >& basisStatusCols, bool& stopped, bool& error, bool& optimal)
3592 MSG_INFO2( spxout, spxout << "Rational factorization hit time limit while solving for primal.\n" );
3671 MSG_INFO2( spxout, spxout << "Rational factorization hit time limit while solving for dual.\n" );
3693 if( ((maximizing && basisStatusRows[i] != SPxSolver::ON_LOWER) || (!maximizing && basisStatusRows[i] != SPxSolver::ON_UPPER))
3703 << " and [lhs,rhs] = [" << rationalToString(lhsRational(i)) << "," << rationalToString(rhsRational(i)) << "]"
3709 if( ((maximizing && basisStatusRows[i] != SPxSolver::ON_UPPER) || (!maximizing && basisStatusRows[i] != SPxSolver::ON_LOWER))
3718 << " and [lhs,rhs] = [" << rationalToString(lhsRational(i)) << "," << rationalToString(rhsRational(i)) << "]"
3734 assert(basisStatusCols[i] != SPxSolver::BASIC || basicDual * colVectorRational(i) == objRational(i));
3743 if( ((maximizing && basisStatusCols[i] != SPxSolver::ON_LOWER) || (!maximizing && basisStatusCols[i] != SPxSolver::ON_UPPER))
3755 if( ((maximizing && basisStatusCols[i] != SPxSolver::ON_UPPER) || (!maximizing && basisStatusCols[i] != SPxSolver::ON_LOWER))
3840 bool SoPlex::_reconstructSolutionRational(SolRational& sol, DataArray< SPxSolver::VarStatus >& basisStatusRows, DataArray< SPxSolver::VarStatus >& basisStatusCols, bool& stopped, bool& error, const Rational& denomBoundSquared)
3880 // we want to notify the user whether the reconstructed solution is basic; otherwise, this would be redundant
3893 MSG_DEBUG( std::cout << "Lower bound of variable " << c << " violated by " << rationalToString(lowerRational(c) - _workSol._primal[c]) << "\n" );
3901 MSG_DEBUG( std::cout << "Upper bound of variable " << c << " violated by " << rationalToString(_workSol._primal[c] - upperRational(c)) << "\n" );
3916 // we want to notify the user whether the reconstructed solution is basic; otherwise, this would be redundant
3929 MSG_DEBUG( std::cout << "Lhs of row " << r << " violated by " << rationalToString(lhsRational(r) - _workSol._slacks[r]) << "\n" );
3937 MSG_DEBUG( std::cout << "Rhs of row " << r << " violated by " << rationalToString(_workSol._slacks[r] - rhsRational(r)) << "\n" );
3958 // check dual multipliers before reduced costs because this check is faster since it does not require the
4011 // compute reduced cost vector; we assume that the objective function vector has less nonzeros than the reduced
|