1 /* Copyright (c) 1997-2021
2 Ewgenij Gawrilow, Michael Joswig, and the polymake team
3 Technische Universität Berlin, Germany
4 https://polymake.org
5
6 This program is free software; you can redistribute it and/or modify it
7 under the terms of the GNU General Public License as published by the
8 Free Software Foundation; either version 2, or (at your option) any
9 later version: http://www.gnu.org/licenses/gpl.txt.
10
11 This program is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15 --------------------------------------------------------------------------------
16 */
17
18 #include "polymake/polytope/soplex_interface.h"
19 #include "polymake/Rational.h"
20
21 #if defined(__clang__)
22 #pragma clang diagnostic push
23 #pragma clang diagnostic ignored "-Wconversion"
24 #pragma clang diagnostic ignored "-Wzero-as-null-pointer-constant"
25 #elif defined(__GNUC__)
26 #pragma GCC diagnostic push
27 #pragma GCC diagnostic ignored "-Wconversion"
28 #pragma GCC diagnostic ignored "-Wzero-as-null-pointer-constant"
29 #endif
30
31 #define SOPLEX_WITH_GMP
32 #include "soplex.h"
33
34 #if defined(__clang__)
35 #pragma clang diagnostic pop
36 #elif defined(__GNUC__)
37 #pragma GCC diagnostic pop
38 #endif
39
40 namespace polymake { namespace polytope { namespace soplex_interface {
41
42 /** This method uses the exact version of Soplex to solve the given LP. */
43 LP_Solution<Rational>
solve(const Matrix<Rational> & Inequalities,const Matrix<Rational> & Equations,const Vector<Rational> & Objective,bool maximize,const Set<Int> & initial_basis) const44 Solver::solve(const Matrix<Rational>& Inequalities, const Matrix<Rational>& Equations,
45 const Vector<Rational>& Objective, bool maximize, const Set<Int>& initial_basis) const
46 {
47 const Int dim = Inequalities.cols(), num_constr = Inequalities.rows() + Equations.rows();
48 if (dim > std::numeric_limits<int>::max() || num_constr > std::numeric_limits<int>::max())
49 throw std::runtime_error("Problem is too big for soplex");
50
51 const int n = int(dim-1); // beware of leading constant term
52 const int m = int(num_constr);
53
54 // create soplex instance
55 soplex::SoPlex soplex;
56
57 // set parameters to activate rational solver
58 soplex.setRealParam(soplex::SoPlex::RealParam::FEASTOL, 0.0);
59 soplex.setRealParam(soplex::SoPlex::RealParam::OPTTOL, 0.0);
60 soplex.setIntParam(soplex::SoPlex::IntParam::SOLVEMODE, soplex::SoPlex::SOLVEMODE_RATIONAL);
61 soplex.setIntParam(soplex::SoPlex::IntParam::SYNCMODE, soplex::SoPlex::SYNCMODE_AUTO);
62 #if 0
63 soplex.setIntParam(soplex::SoPlex::IntParam::READMODE, soplex::SoPlex::READMODE_RATIONAL);
64 soplex.setIntParam(soplex::SoPlex::IntParam::CHECKMODE, soplex::SoPlex::CHECKMODE_RATIONAL);
65 #endif
66
67 // turn off output
68 soplex.setIntParam(soplex::SoPlex::IntParam::VERBOSITY, soplex::SoPlex::VERBOSITY_ERROR);
69 #if POLYMAKE_DEBUG
70 if (debug_print) {
71 soplex.setIntParam(soplex::SoPlex::IntParam::VERBOSITY, soplex::SoPlex::VERBOSITY_NORMAL);
72 }
73 #endif
74
75 // set objective sense
76 soplex.setIntParam(soplex::SoPlex::IntParam::OBJSENSE,
77 maximize ? soplex::SoPlex::OBJSENSE_MAXIMIZE : soplex::SoPlex::OBJSENSE_MINIMIZE);
78
79 // get value for infinity
80 const double realinfty = soplex.realParam(soplex::SoPlex::RealParam::INFTY);
81 const Rational plusinfinity(realinfty);
82 const Rational minusinfinity(-realinfty);
83
84 if (Equations.rows() > 0 && n != Equations.cols()-1) {
85 throw std::runtime_error("Dimensions do not fit.\n");
86 }
87
88 // create variables/columns
89 mpq_t* obj = new mpq_t [n];
90 mpq_t* lower = new mpq_t [n];
91 mpq_t* upper = new mpq_t [n];
92
93 for (int j = 0; j < n; ++j) {
94 mpq_init(obj[j]);
95 mpq_set(obj[j], Objective[j+1].get_rep());
96
97 mpq_init(lower[j]);
98 mpq_init(upper[j]);
99 mpq_set(lower[j], minusinfinity.get_rep());
100 mpq_set(upper[j], plusinfinity.get_rep());
101 }
102
103 // add columns
104 soplex.addColsRational(obj, lower, nullptr, nullptr, nullptr, nullptr, n, 0, upper);
105
106 // free space
107 for (int j = 0; j < n; ++j) {
108 mpq_clear(lower[j]);
109 mpq_clear(upper[j]);
110 mpq_clear(obj[j]);
111 }
112
113 delete [] upper;
114 delete [] lower;
115 delete [] obj;
116
117 // reserve space for rows
118 mpq_t* rowValues = new mpq_t [n];
119 int* rowIndices = new int [n];
120
121 for (int j = 0; j < n; ++j)
122 mpq_init(rowValues[j]);
123
124 mpq_t lhs;
125 mpq_t rhs;
126 mpq_init(rhs);
127 mpq_init(lhs);
128
129 // create rows - inequalities
130 for (int i = 0; i < Inequalities.rows(); ++i) {
131 int cnt = 0;
132 for (int j = 0; j < n; ++j) {
133 Rational val = Inequalities(i, j+1);
134 if (!is_zero(val)) {
135 mpq_set(rowValues[cnt], val.get_rep());
136 rowIndices[cnt++] = j;
137 }
138 }
139
140 Rational val = -Inequalities(i, 0);
141 mpq_set(rhs, plusinfinity.get_rep());
142 mpq_set(lhs, val.get_rep());
143
144 soplex.addRowRational(&lhs, rowValues, rowIndices, cnt, &rhs);
145 }
146
147 // create rows - equations
148 for (int i = 0; i < Equations.rows(); ++i) {
149 int cnt = 0;
150 for (int j = 0; j < n; ++j) {
151 Rational val = Equations(i, j+1);
152 if (!is_zero(val)) {
153 mpq_set(rowValues[cnt], val.get_rep());
154 rowIndices[cnt++] = j;
155 }
156 }
157
158 Rational val = -Equations(i, 0);
159 mpq_set(lhs, val.get_rep());
160 mpq_set(rhs, val.get_rep());
161
162 soplex.addRowRational(&lhs, rowValues, rowIndices, cnt, &rhs);
163 }
164 mpq_clear(lhs);
165 mpq_clear(rhs);
166
167 for (int j = 0; j < n; ++j)
168 mpq_clear(rowValues[j]);
169
170 delete [] rowValues;
171 delete [] rowIndices;
172
173 #if POLYMAKE_DEBUG
174 if (debug_print) {
175 std::cout << "Number of columns: " << soplex.numColsRational() << std::endl;
176 std::cout << "Number of rows: " << soplex.numRowsRational() << std::endl;
177 }
178 #endif
179
180 // add start basis
181 if (!initial_basis.empty()) {
182 soplex::SPxSolver::VarStatus* rows = new soplex::SPxSolver::VarStatus [m];
183 soplex::SPxSolver::VarStatus* cols = new soplex::SPxSolver::VarStatus [n];
184
185 // init basis
186 for (int j = 0; j < n; ++j)
187 cols[j] = soplex::SPxSolver::VarStatus::ZERO;
188
189 for (int i = 0; i < m; ++i)
190 rows[i] = soplex::SPxSolver::VarStatus::BASIC;
191
192 int count = 0;
193 for (const Int i : initial_basis) {
194 rows[i] = soplex::SPxSolver::VarStatus::ON_LOWER;
195 if (++count >= n)
196 break;
197 }
198 soplex.setBasis(rows, cols);
199
200 delete [] cols;
201 delete [] rows;
202 }
203
204 #if POLYMAKE_DEBUG
205 if (debug_print) {
206 soplex.writeFileRational("debug.lp", nullptr, nullptr, nullptr);
207 std::cout << "Wrote debug output of LP to 'debug.lp'." << std::endl;
208 }
209 #endif
210
211 LP_Solution<Rational> result;
212
213 // solves the LP
214 const soplex::SPxSolver::Status status = soplex.solve();
215
216 switch (status) {
217 case soplex::SPxSolver::ERROR:
218 case soplex::SPxSolver::NO_RATIOTESTER:
219 case soplex::SPxSolver::NO_PRICER:
220 case soplex::SPxSolver::NO_SOLVER:
221 case soplex::SPxSolver::NOT_INIT:
222 case soplex::SPxSolver::NO_PROBLEM:
223 case soplex::SPxSolver::UNKNOWN:
224 throw std::runtime_error("Error in setting up soplex.");
225
226 case soplex::SPxSolver::ABORT_TIME:
227 case soplex::SPxSolver::ABORT_ITER:
228 case soplex::SPxSolver::ABORT_VALUE:
229 case soplex::SPxSolver::REGULAR:
230 case soplex::SPxSolver::RUNNING:
231 throw std::runtime_error("Error in solving LP with soplex. This should not happen.");
232
233 case soplex::SPxSolver::SINGULAR:
234 case soplex::SPxSolver::ABORT_CYCLING:
235 throw std::runtime_error("Numerical problems while solving LP with soplex.");
236
237 case soplex::SPxSolver::OPTIMAL:
238 {
239 result.status = LP_status::valid;
240 // get solution
241 auto soplexsol = std::make_unique<mpq_t[]>(n+1);
242 for (int j = 0; j <= n; ++j)
243 mpq_init(soplexsol[j]);
244 mpq_set_si(soplexsol[0], 1, 1);
245 soplex.getPrimalRational(&soplexsol[1], n);
246
247 // convert to polymake rationals;
248 // add constant term from objective function
249 result.objective_value = Rational(std::move(soplex.objValueRational().getMpqRef())) + Objective[0];
250 result.solution = Vector<Rational>(n+1, enforce_movable_values(&soplexsol[0]));
251 break;
252 }
253 case soplex::SPxSolver::UNBOUNDED:
254 result.status = LP_status::unbounded;
255 break;
256
257 case soplex::SPxSolver::INFEASIBLE:
258 result.status = LP_status::infeasible;
259 break;
260
261 case soplex::SPxSolver::INForUNBD:
262 // not sure what to do ...
263 result.status = LP_status::infeasible;
264 break;
265
266 default:
267 throw std::runtime_error("Unknown error.");
268 }
269
270 return result;
271 }
272
273 } } }
274
275 // Local Variables:
276 // mode:C++
277 // c-basic-offset:3
278 // indent-tabs-mode:nil
279 // End:
280