1 // Copyright 2019 RTE
2 // Licensed under the Apache License, Version 2.0 (the "License");
3 // you may not use this file except in compliance with the License.
4 // You may obtain a copy of the License at
5 //
6 // http://www.apache.org/licenses/LICENSE-2.0
7 //
8 // Unless required by applicable law or agreed to in writing, software
9 // distributed under the License is distributed on an "AS IS" BASIS,
10 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
11 // See the License for the specific language governing permissions and
12 // limitations under the License.
13
14 // Initial version of this code was provided by RTE
15
16 #if defined(USE_XPRESS)
17
18 #include <algorithm>
19 #include <limits>
20 #include <memory>
21 #include <string>
22
23 #include "absl/strings/str_format.h"
24 #include "ortools/base/integral_types.h"
25 #include "ortools/base/logging.h"
26 #include "ortools/base/timer.h"
27 #include "ortools/linear_solver/linear_solver.h"
28
29 extern "C" {
30 #include "xprs.h"
31 }
32
33 #define XPRS_INTEGER 'I'
34 #define XPRS_CONTINUOUS 'C'
35 #define STRINGIFY2(X) #X
36 #define STRINGIFY(X) STRINGIFY2(X)
37
printError(const XPRSprob & mLp,int line)38 void printError(const XPRSprob& mLp, int line) {
39 char errmsg[512];
40 XPRSgetlasterror(mLp, errmsg);
41 VLOG(0) << absl::StrFormat("Function line %d did not execute correctly: %s\n",
42 line, errmsg);
43 exit(0);
44 }
45
XPRSgetnumcols(const XPRSprob & mLp)46 int XPRSgetnumcols(const XPRSprob& mLp) {
47 int nCols = 0;
48 XPRSgetintattrib(mLp, XPRS_COLS, &nCols);
49 return nCols;
50 }
51
XPRSgetnumrows(const XPRSprob & mLp)52 int XPRSgetnumrows(const XPRSprob& mLp) {
53 int nRows = 0;
54 XPRSgetintattrib(mLp, XPRS_ROWS, &nRows);
55 return nRows;
56 }
57
XPRSgetitcnt(const XPRSprob & mLp)58 int XPRSgetitcnt(const XPRSprob& mLp) {
59 int nIters = 0;
60 XPRSgetintattrib(mLp, XPRS_SIMPLEXITER, &nIters);
61 return nIters;
62 }
63
XPRSgetnodecnt(const XPRSprob & mLp)64 int XPRSgetnodecnt(const XPRSprob& mLp) {
65 int nNodes = 0;
66 XPRSgetintattrib(mLp, XPRS_NODES, &nNodes);
67 return nNodes;
68 }
69
XPRSsetobjoffset(const XPRSprob & mLp,double value)70 int XPRSsetobjoffset(const XPRSprob& mLp, double value) {
71 XPRSsetdblcontrol(mLp, XPRS_OBJRHS, value);
72 return 0;
73 }
74
75 enum XPRS_BASIS_STATUS {
76 XPRS_AT_LOWER = 0,
77 XPRS_BASIC = 1,
78 XPRS_AT_UPPER = 2,
79 XPRS_FREE_SUPER = 3
80 };
81
82 // In case we need to return a double but don't have a value for that
83 // we just return a NaN.
84 #if !defined(CPX_NAN)
85 #define XPRS_NAN std::numeric_limits<double>::quiet_NaN()
86 #endif
87
88 // The argument to this macro is the invocation of a XPRS function that
89 // returns a status. If the function returns non-zero the macro aborts
90 // the program with an appropriate error message.
91 #define CHECK_STATUS(s) \
92 do { \
93 int const status_ = s; \
94 CHECK_EQ(0, status_); \
95 } while (0)
96
97 namespace operations_research {
98
99 using std::unique_ptr;
100
101 // For a model that is extracted to an instance of this class there is a
102 // 1:1 corresponence between MPVariable instances and XPRESS columns: the
103 // index of an extracted variable is the column index in the XPRESS model.
104 // Similar for instances of MPConstraint: the index of the constraint in
105 // the model is the row index in the XPRESS model.
106 class XpressInterface : public MPSolverInterface {
107 public:
108 // NOTE: 'mip' specifies the type of the problem (either continuous or
109 // mixed integer. This type is fixed for the lifetime of the
110 // instance. There are no dynamic changes to the model type.
111 explicit XpressInterface(MPSolver* const solver, bool mip);
112 ~XpressInterface();
113
114 // Sets the optimization direction (min/max).
115 virtual void SetOptimizationDirection(bool maximize);
116
117 // ----- Solve -----
118 // Solve the problem using the parameter values specified.
119 virtual MPSolver::ResultStatus Solve(MPSolverParameters const& param);
120
121 // ----- Model modifications and extraction -----
122 // Resets extracted model
123 virtual void Reset();
124
125 virtual void SetVariableBounds(int var_index, double lb, double ub);
126 virtual void SetVariableInteger(int var_index, bool integer);
127 virtual void SetConstraintBounds(int row_index, double lb, double ub);
128
129 virtual void AddRowConstraint(MPConstraint* const ct);
130 virtual void AddVariable(MPVariable* const var);
131 virtual void SetCoefficient(MPConstraint* const constraint,
132 MPVariable const* const variable,
133 double new_value, double old_value);
134
135 // Clear a constraint from all its terms.
136 virtual void ClearConstraint(MPConstraint* const constraint);
137 // Change a coefficient in the linear objective
138 virtual void SetObjectiveCoefficient(MPVariable const* const variable,
139 double coefficient);
140 // Change the constant term in the linear objective.
141 virtual void SetObjectiveOffset(double value);
142 // Clear the objective from all its terms.
143 virtual void ClearObjective();
144
145 // ------ Query statistics on the solution and the solve ------
146 // Number of simplex iterations
147 virtual int64_t iterations() const;
148 // Number of branch-and-bound nodes. Only available for discrete problems.
149 virtual int64_t nodes() const;
150
151 // Returns the basis status of a row.
152 virtual MPSolver::BasisStatus row_status(int constraint_index) const;
153 // Returns the basis status of a column.
154 virtual MPSolver::BasisStatus column_status(int variable_index) const;
155
156 // ----- Misc -----
157
158 // Query problem type.
159 // Remember that problem type is a static property that is set
160 // in the constructor and never changed.
IsContinuous() const161 virtual bool IsContinuous() const { return IsLP(); }
IsLP() const162 virtual bool IsLP() const { return !mMip; }
IsMIP() const163 virtual bool IsMIP() const { return mMip; }
164
165 virtual void ExtractNewVariables();
166 virtual void ExtractNewConstraints();
167 virtual void ExtractObjective();
168
169 virtual std::string SolverVersion() const;
170
underlying_solver()171 virtual void* underlying_solver() { return reinterpret_cast<void*>(mLp); }
172
ComputeExactConditionNumber() const173 virtual double ComputeExactConditionNumber() const {
174 if (!IsContinuous()) {
175 LOG(DFATAL) << "ComputeExactConditionNumber not implemented for"
176 << " XPRESS_MIXED_INTEGER_PROGRAMMING";
177 return 0.0;
178 }
179
180 // TODO(user,user): Not yet working.
181 LOG(DFATAL) << "ComputeExactConditionNumber not implemented for"
182 << " XPRESS_LINEAR_PROGRAMMING";
183 return 0.0;
184 }
185
186 protected:
187 // Set all parameters in the underlying solver.
188 virtual void SetParameters(MPSolverParameters const& param);
189 // Set each parameter in the underlying solver.
190 virtual void SetRelativeMipGap(double value);
191 virtual void SetPrimalTolerance(double value);
192 virtual void SetDualTolerance(double value);
193 virtual void SetPresolveMode(int value);
194 virtual void SetScalingMode(int value);
195 virtual void SetLpAlgorithm(int value);
196
197 virtual bool ReadParameterFile(std::string const& filename);
198 virtual std::string ValidFileExtensionForParameterFile() const;
199
200 private:
201 // Mark modeling object "out of sync". This implicitly invalidates
202 // solution information as well. It is the counterpart of
203 // MPSolverInterface::InvalidateSolutionSynchronization
InvalidateModelSynchronization()204 void InvalidateModelSynchronization() {
205 mCstat = 0;
206 mRstat = 0;
207 sync_status_ = MUST_RELOAD;
208 }
209
210 // Transform XPRESS basis status to MPSolver basis status.
211 static MPSolver::BasisStatus xformBasisStatus(int xpress_basis_status);
212
213 private:
214 XPRSprob mLp;
215 bool const mMip;
216 // Incremental extraction.
217 // Without incremental extraction we have to re-extract the model every
218 // time we perform a solve. Due to the way the Reset() function is
219 // implemented, this will lose MIP start or basis information from a
220 // previous solve. On the other hand, if there is a significant changes
221 // to the model then just re-extracting everything is usually faster than
222 // keeping the low-level modeling object in sync with the high-level
223 // variables/constraints.
224 // Note that incremental extraction is particularly expensive in function
225 // ExtractNewVariables() since there we must scan _all_ old constraints
226 // and update them with respect to the new variables.
227 bool const supportIncrementalExtraction;
228
229 // Use slow and immediate updates or try to do bulk updates.
230 // For many updates to the model we have the option to either perform
231 // the update immediately with a potentially slow operation or to
232 // just mark the low-level modeling object out of sync and re-extract
233 // the model later.
234 enum SlowUpdates {
235 SlowSetCoefficient = 0x0001,
236 SlowClearConstraint = 0x0002,
237 SlowSetObjectiveCoefficient = 0x0004,
238 SlowClearObjective = 0x0008,
239 SlowSetConstraintBounds = 0x0010,
240 SlowSetVariableInteger = 0x0020,
241 SlowSetVariableBounds = 0x0040,
242 SlowUpdatesAll = 0xffff
243 } const slowUpdates;
244 // XPRESS has no method to query the basis status of a single variable.
245 // Hence we query the status only once and cache the array. This is
246 // much faster in case the basis status of more than one row/column
247 // is required.
248 unique_ptr<int[]> mutable mCstat;
249 unique_ptr<int[]> mutable mRstat;
250
251 // Setup the right-hand side of a constraint from its lower and upper bound.
252 static void MakeRhs(double lb, double ub, double& rhs, char& sense,
253 double& range);
254 };
255
256 /** init XPRESS environment */
init_xpress_env(int xpress_oem_license_key=0)257 int init_xpress_env(int xpress_oem_license_key = 0) {
258 int code;
259
260 const char* xpress_from_env = getenv("XPRESS");
261 std::string xpresspath;
262
263 if (xpress_from_env == nullptr) {
264 #if defined(XPRESS_PATH)
265 std::string path(STRINGIFY(XPRESS_PATH));
266 LOG(WARNING)
267 << "Environment variable XPRESS undefined. Trying compile path "
268 << "'" << path << "'";
269 #if defined(_MSC_VER)
270 // need to remove the enclosing '\"' from the string itself.
271 path.erase(std::remove(path.begin(), path.end(), '\"'), path.end());
272 xpresspath = path + "\\bin";
273 #else // _MSC_VER
274 xpresspath = path + "/bin";
275 #endif // _MSC_VER
276 #else
277 LOG(WARNING)
278 << "XpressInterface Error : Environment variable XPRESS undefined.\n";
279 return -1;
280 #endif
281 } else {
282 xpresspath = xpress_from_env;
283 }
284
285 /** if not an OEM key */
286 if (xpress_oem_license_key == 0) {
287 LOG(WARNING) << "XpressInterface : Initialising xpress-MP with parameter "
288 << xpresspath << std::endl;
289
290 code = XPRSinit(xpresspath.c_str());
291
292 if (!code) {
293 /** XPRSbanner informs about Xpress version, options and error messages */
294 char banner[1000];
295 XPRSgetbanner(banner);
296
297 LOG(WARNING) << "XpressInterface : Xpress banner :\n"
298 << banner << std::endl;
299 return 0;
300 } else {
301 char errmsg[256];
302 XPRSgetlicerrmsg(errmsg, 256);
303
304 VLOG(0) << "XpressInterface : License error : " << errmsg << std::endl;
305 VLOG(0) << "XpressInterface : XPRSinit returned code : " << code << "\n";
306
307 char banner[1000];
308 XPRSgetbanner(banner);
309
310 LOG(ERROR) << "XpressInterface : Xpress banner :\n" << banner << "\n";
311 return -1;
312 }
313 } else {
314 /** if OEM key */
315 LOG(WARNING) << "XpressInterface : Initialising xpress-MP with OEM key "
316 << xpress_oem_license_key << "\n";
317
318 int nvalue = 0;
319 int ierr;
320 char slicmsg[256] = "";
321 char errmsg[256];
322
323 XPRSlicense(&nvalue, slicmsg);
324 VLOG(0) << "XpressInterface : First message from XPRSLicense : " << slicmsg
325 << "\n";
326
327 nvalue = xpress_oem_license_key - ((nvalue * nvalue) / 19);
328 ierr = XPRSlicense(&nvalue, slicmsg);
329
330 VLOG(0) << "XpressInterface : Second message from XPRSLicense : " << slicmsg
331 << "\n";
332 if (ierr == 16) {
333 VLOG(0) << "XpressInterface : Optimizer development software detected\n";
334 } else if (ierr != 0) {
335 /** get the license error message */
336 XPRSgetlicerrmsg(errmsg, 256);
337
338 LOG(ERROR) << "XpressInterface : " << errmsg << "\n";
339 return -1;
340 }
341
342 code = XPRSinit(NULL);
343
344 if (!code) {
345 return 0;
346 } else {
347 LOG(ERROR) << "XPRSinit returned code : " << code << "\n";
348 return -1;
349 }
350 }
351 }
352
353 // Creates a LP/MIP instance.
XpressInterface(MPSolver * const solver,bool mip)354 XpressInterface::XpressInterface(MPSolver* const solver, bool mip)
355 : MPSolverInterface(solver),
356 mLp(0),
357 mMip(mip),
358 supportIncrementalExtraction(false),
359 slowUpdates(static_cast<SlowUpdates>(SlowSetObjectiveCoefficient |
360 SlowClearObjective)),
361 mCstat(),
362 mRstat() {
363 int status = init_xpress_env();
364 CHECK_STATUS(status);
365 status = XPRScreateprob(&mLp);
366 CHECK_STATUS(status);
367 DCHECK(mLp != nullptr); // should not be NULL if status=0
368 CHECK_STATUS(XPRSloadlp(mLp, "newProb", 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
369
370 CHECK_STATUS(
371 XPRSchgobjsense(mLp, maximize_ ? XPRS_OBJ_MAXIMIZE : XPRS_OBJ_MINIMIZE));
372 }
373
~XpressInterface()374 XpressInterface::~XpressInterface() {
375 CHECK_STATUS(XPRSdestroyprob(mLp));
376 CHECK_STATUS(XPRSfree());
377 }
378
SolverVersion() const379 std::string XpressInterface::SolverVersion() const {
380 // We prefer XPRSversionnumber() over XPRSversion() since the
381 // former will never pose any encoding issues.
382 int version = 0;
383 CHECK_STATUS(XPRSgetintcontrol(mLp, XPRS_VERSION, &version));
384
385 int const major = version / 1000000;
386 version -= major * 1000000;
387 int const release = version / 10000;
388 version -= release * 10000;
389 int const mod = version / 100;
390 version -= mod * 100;
391 int const fix = version;
392
393 return absl::StrFormat("XPRESS library version %d.%02d.%02d.%02d", major,
394 release, mod, fix);
395 }
396
397 // ------ Model modifications and extraction -----
398
Reset()399 void XpressInterface::Reset() {
400 // Instead of explicitly clearing all modeling objects we
401 // just delete the problem object and allocate a new one.
402 CHECK_STATUS(XPRSdestroyprob(mLp));
403
404 int status;
405 status = XPRScreateprob(&mLp);
406 CHECK_STATUS(status);
407 DCHECK(mLp != nullptr); // should not be NULL if status=0
408 CHECK_STATUS(XPRSloadlp(mLp, "newProb", 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
409
410 CHECK_STATUS(
411 XPRSchgobjsense(mLp, maximize_ ? XPRS_OBJ_MAXIMIZE : XPRS_OBJ_MINIMIZE));
412
413 ResetExtractionInformation();
414 mCstat = 0;
415 mRstat = 0;
416 }
417
SetOptimizationDirection(bool maximize)418 void XpressInterface::SetOptimizationDirection(bool maximize) {
419 InvalidateSolutionSynchronization();
420 XPRSchgobjsense(mLp, maximize ? XPRS_OBJ_MAXIMIZE : XPRS_OBJ_MINIMIZE);
421 }
422
SetVariableBounds(int var_index,double lb,double ub)423 void XpressInterface::SetVariableBounds(int var_index, double lb, double ub) {
424 InvalidateSolutionSynchronization();
425
426 // Changing the bounds of a variable is fast. However, doing this for
427 // many variables may still be slow. So we don't perform the update by
428 // default. However, if we support incremental extraction
429 // (supportIncrementalExtraction is true) then we MUST perform the
430 // update here or we will lose it.
431
432 if (!supportIncrementalExtraction && !(slowUpdates & SlowSetVariableBounds)) {
433 InvalidateModelSynchronization();
434 } else {
435 if (variable_is_extracted(var_index)) {
436 // Variable has already been extracted, so we must modify the
437 // modeling object.
438 DCHECK_LT(var_index, last_variable_index_);
439 char const lu[2] = {'L', 'U'};
440 double const bd[2] = {lb, ub};
441 int const idx[2] = {var_index, var_index};
442 CHECK_STATUS(XPRSchgbounds(mLp, 2, idx, lu, bd));
443 } else {
444 // Variable is not yet extracted. It is sufficient to just mark
445 // the modeling object "out of sync"
446 InvalidateModelSynchronization();
447 }
448 }
449 }
450
451 // Modifies integrality of an extracted variable.
SetVariableInteger(int var_index,bool integer)452 void XpressInterface::SetVariableInteger(int var_index, bool integer) {
453 InvalidateSolutionSynchronization();
454
455 // NOTE: The type of the model (continuous or mixed integer) is
456 // defined once and for all in the constructor. There are no
457 // dynamic changes to the model type.
458
459 // Changing the type of a variable should be fast. Still, doing all
460 // updates in one big chunk right before solve() is usually faster.
461 // However, if we support incremental extraction
462 // (supportIncrementalExtraction is true) then we MUST change the
463 // type of extracted variables here.
464
465 if (!supportIncrementalExtraction && !slowUpdates &&
466 !SlowSetVariableInteger) {
467 InvalidateModelSynchronization();
468 } else {
469 if (mMip) {
470 if (variable_is_extracted(var_index)) {
471 // Variable is extracted. Change the type immediately.
472 // TODO: Should we check the current type and don't do anything
473 // in case the type does not change?
474 DCHECK_LE(var_index, XPRSgetnumcols(mLp));
475 char const type = integer ? XPRS_INTEGER : XPRS_CONTINUOUS;
476 CHECK_STATUS(XPRSchgcoltype(mLp, 1, &var_index, &type));
477 } else {
478 InvalidateModelSynchronization();
479 }
480 } else {
481 LOG(DFATAL)
482 << "Attempt to change variable to integer in non-MIP problem!";
483 }
484 }
485 }
486
487 // Setup the right-hand side of a constraint.
MakeRhs(double lb,double ub,double & rhs,char & sense,double & range)488 void XpressInterface::MakeRhs(double lb, double ub, double& rhs, char& sense,
489 double& range) {
490 if (lb == ub) {
491 // Both bounds are equal -> this is an equality constraint
492 rhs = lb;
493 range = 0.0;
494 sense = 'E';
495 } else if (lb > XPRS_MINUSINFINITY && ub < XPRS_PLUSINFINITY) {
496 // Both bounds are finite -> this is a ranged constraint
497 // The value of a ranged constraint is allowed to be in
498 // [ rhs[i], rhs[i]+rngval[i] ]
499 // see also the reference documentation for XPRSnewrows()
500 if (ub < lb) {
501 // The bounds for the constraint are contradictory. XPRESS models
502 // a range constraint l <= ax <= u as
503 // ax = l + v
504 // where v is an auxiliary variable the range of which is controlled
505 // by l and u: if l < u then v in [0, u-l]
506 // else v in [u-l, 0]
507 // (the range is specified as the rngval[] argument to XPRSnewrows).
508 // Thus XPRESS cannot represent range constraints with contradictory
509 // bounds and we must error out here.
510 CHECK_STATUS(-1);
511 }
512 rhs = ub;
513 range = ub - lb;
514 sense = 'R';
515 } else if (ub < XPRS_PLUSINFINITY || (std::abs(ub) == XPRS_PLUSINFINITY &&
516 std::abs(lb) > XPRS_PLUSINFINITY)) {
517 // Finite upper, infinite lower bound -> this is a <= constraint
518 rhs = ub;
519 range = 0.0;
520 sense = 'L';
521 } else if (lb > XPRS_MINUSINFINITY || (std::abs(lb) == XPRS_PLUSINFINITY &&
522 std::abs(ub) > XPRS_PLUSINFINITY)) {
523 // Finite lower, infinite upper bound -> this is a >= constraint
524 rhs = lb;
525 range = 0.0;
526 sense = 'G';
527 } else {
528 // Lower and upper bound are both infinite.
529 // This is used for example in .mps files to specify alternate
530 // objective functions.
531 // Note that the case lb==ub was already handled above, so we just
532 // pick the bound with larger magnitude and create a constraint for it.
533 // Note that we replace the infinite bound by XPRS_PLUSINFINITY since
534 // bounds with larger magnitude may cause other XPRESS functions to
535 // fail (for example the export to LP files).
536 DCHECK_GT(std::abs(lb), XPRS_PLUSINFINITY);
537 DCHECK_GT(std::abs(ub), XPRS_PLUSINFINITY);
538 if (std::abs(lb) > std::abs(ub)) {
539 rhs = (lb < 0) ? -XPRS_PLUSINFINITY : XPRS_PLUSINFINITY;
540 sense = 'G';
541 } else {
542 rhs = (ub < 0) ? -XPRS_PLUSINFINITY : XPRS_PLUSINFINITY;
543 sense = 'L';
544 }
545 range = 0.0;
546 }
547 }
548
SetConstraintBounds(int index,double lb,double ub)549 void XpressInterface::SetConstraintBounds(int index, double lb, double ub) {
550 InvalidateSolutionSynchronization();
551
552 // Changing rhs, sense, or range of a constraint is not too slow.
553 // Still, doing all the updates in one large operation is faster.
554 // Note however that if we do not want to re-extract the full model
555 // for each solve (supportIncrementalExtraction is true) then we MUST
556 // update the constraint here, otherwise we lose this update information.
557
558 if (!supportIncrementalExtraction &&
559 !(slowUpdates & SlowSetConstraintBounds)) {
560 InvalidateModelSynchronization();
561 } else {
562 if (constraint_is_extracted(index)) {
563 // Constraint is already extracted, so we must update its bounds
564 // and its type.
565 DCHECK(mLp != NULL);
566 char sense;
567 double range, rhs;
568 MakeRhs(lb, ub, rhs, sense, range);
569 CHECK_STATUS(XPRSchgrhs(mLp, 1, &index, &lb));
570 CHECK_STATUS(XPRSchgrowtype(mLp, 1, &index, &sense));
571 CHECK_STATUS(XPRSchgrhsrange(mLp, 1, &index, &range));
572 } else {
573 // Constraint is not yet extracted. It is sufficient to mark the
574 // modeling object as "out of sync"
575 InvalidateModelSynchronization();
576 }
577 }
578 }
579
AddRowConstraint(MPConstraint * const ct)580 void XpressInterface::AddRowConstraint(MPConstraint* const ct) {
581 // This is currently only invoked when a new constraint is created,
582 // see MPSolver::MakeRowConstraint().
583 // At this point we only have the lower and upper bounds of the
584 // constraint. We could immediately call XPRSaddrows() here but it is
585 // usually much faster to handle the fully populated constraint in
586 // ExtractNewConstraints() right before the solve.
587 InvalidateModelSynchronization();
588 }
589
AddVariable(MPVariable * const ct)590 void XpressInterface::AddVariable(MPVariable* const ct) {
591 // This is currently only invoked when a new variable is created,
592 // see MPSolver::MakeVar().
593 // At this point the variable does not appear in any constraints or
594 // the objective function. We could invoke XPRSaddcols() to immediately
595 // create the variable here but it is usually much faster to handle the
596 // fully setup variable in ExtractNewVariables() right before the solve.
597 InvalidateModelSynchronization();
598 }
599
SetCoefficient(MPConstraint * const constraint,MPVariable const * const variable,double new_value,double)600 void XpressInterface::SetCoefficient(MPConstraint* const constraint,
601 MPVariable const* const variable,
602 double new_value, double) {
603 InvalidateSolutionSynchronization();
604
605 // Changing a single coefficient in the matrix is potentially pretty
606 // slow since that coefficient has to be found in the sparse matrix
607 // representation. So by default we don't perform this update immediately
608 // but instead mark the low-level modeling object "out of sync".
609 // If we want to support incremental extraction then we MUST perform
610 // the modification immediately or we will lose it.
611
612 if (!supportIncrementalExtraction && !(slowUpdates & SlowSetCoefficient)) {
613 InvalidateModelSynchronization();
614 } else {
615 int const row = constraint->index();
616 int const col = variable->index();
617 if (constraint_is_extracted(row) && variable_is_extracted(col)) {
618 // If row and column are both extracted then we can directly
619 // update the modeling object
620 DCHECK_LE(row, last_constraint_index_);
621 DCHECK_LE(col, last_variable_index_);
622 CHECK_STATUS(XPRSchgcoef(mLp, row, col, new_value));
623 } else {
624 // If either row or column is not yet extracted then we can
625 // defer the update to ExtractModel()
626 InvalidateModelSynchronization();
627 }
628 }
629 }
630
ClearConstraint(MPConstraint * const constraint)631 void XpressInterface::ClearConstraint(MPConstraint* const constraint) {
632 int const row = constraint->index();
633 if (!constraint_is_extracted(row))
634 // There is nothing to do if the constraint was not even extracted.
635 return;
636
637 // Clearing a constraint means setting all coefficients in the corresponding
638 // row to 0 (we cannot just delete the row since that would renumber all
639 // the constraints/rows after it).
640 // Modifying coefficients in the matrix is potentially pretty expensive
641 // since they must be found in the sparse matrix representation. That is
642 // why by default we do not modify the coefficients here but only mark
643 // the low-level modeling object "out of sync".
644
645 if (!(slowUpdates & SlowClearConstraint)) {
646 InvalidateModelSynchronization();
647 } else {
648 InvalidateSolutionSynchronization();
649
650 int const len = constraint->coefficients_.size();
651 unique_ptr<int[]> rowind(new int[len]);
652 unique_ptr<int[]> colind(new int[len]);
653 unique_ptr<double[]> val(new double[len]);
654 int j = 0;
655 const auto& coeffs = constraint->coefficients_;
656 for (auto it(coeffs.begin()); it != coeffs.end(); ++it) {
657 int const col = it->first->index();
658 if (variable_is_extracted(col)) {
659 rowind[j] = row;
660 colind[j] = col;
661 val[j] = 0.0;
662 ++j;
663 }
664 }
665 if (j)
666 CHECK_STATUS(XPRSchgmcoef(mLp, j, rowind.get(), colind.get(), val.get()));
667 }
668 }
669
SetObjectiveCoefficient(MPVariable const * const variable,double coefficient)670 void XpressInterface::SetObjectiveCoefficient(MPVariable const* const variable,
671 double coefficient) {
672 int const col = variable->index();
673 if (!variable_is_extracted(col))
674 // Nothing to do if variable was not even extracted
675 return;
676
677 InvalidateSolutionSynchronization();
678
679 // The objective function is stored as a dense vector, so updating a
680 // single coefficient is O(1). So by default we update the low-level
681 // modeling object here.
682 // If we support incremental extraction then we have no choice but to
683 // perform the update immediately.
684
685 if (supportIncrementalExtraction ||
686 (slowUpdates & SlowSetObjectiveCoefficient)) {
687 CHECK_STATUS(XPRSchgobj(mLp, 1, &col, &coefficient));
688 } else {
689 InvalidateModelSynchronization();
690 }
691 }
692
SetObjectiveOffset(double value)693 void XpressInterface::SetObjectiveOffset(double value) {
694 // Changing the objective offset is O(1), so we always do it immediately.
695 InvalidateSolutionSynchronization();
696 CHECK_STATUS(XPRSsetobjoffset(mLp, value));
697 }
698
ClearObjective()699 void XpressInterface::ClearObjective() {
700 InvalidateSolutionSynchronization();
701
702 // Since the objective function is stored as a dense vector updating
703 // it is O(n), so we usually perform the update immediately.
704 // If we want to support incremental extraction then we have no choice
705 // but to perform the update immediately.
706
707 if (supportIncrementalExtraction || (slowUpdates & SlowClearObjective)) {
708 int const cols = XPRSgetnumcols(mLp);
709 unique_ptr<int[]> ind(new int[cols]);
710 unique_ptr<double[]> zero(new double[cols]);
711 int j = 0;
712 const auto& coeffs = solver_->objective_->coefficients_;
713 for (auto it(coeffs.begin()); it != coeffs.end(); ++it) {
714 int const idx = it->first->index();
715 // We only need to reset variables that have been extracted.
716 if (variable_is_extracted(idx)) {
717 DCHECK_LT(idx, cols);
718 ind[j] = idx;
719 zero[j] = 0.0;
720 ++j;
721 }
722 }
723 if (j > 0) CHECK_STATUS(XPRSchgobj(mLp, j, ind.get(), zero.get()));
724 CHECK_STATUS(XPRSsetobjoffset(mLp, 0.0));
725 } else {
726 InvalidateModelSynchronization();
727 }
728 }
729
730 // ------ Query statistics on the solution and the solve ------
731
iterations() const732 int64_t XpressInterface::iterations() const {
733 if (!CheckSolutionIsSynchronized()) return kUnknownNumberOfIterations;
734 return static_cast<int64_t>(XPRSgetitcnt(mLp));
735 }
736
nodes() const737 int64_t XpressInterface::nodes() const {
738 if (mMip) {
739 if (!CheckSolutionIsSynchronized()) return kUnknownNumberOfNodes;
740 return static_cast<int64_t>(XPRSgetnodecnt(mLp));
741 } else {
742 LOG(DFATAL) << "Number of nodes only available for discrete problems";
743 return kUnknownNumberOfNodes;
744 }
745 }
746
747 // Transform a XPRESS basis status to an MPSolver basis status.
xformBasisStatus(int xpress_basis_status)748 MPSolver::BasisStatus XpressInterface::xformBasisStatus(
749 int xpress_basis_status) {
750 switch (xpress_basis_status) {
751 case XPRS_AT_LOWER:
752 return MPSolver::AT_LOWER_BOUND;
753 case XPRS_BASIC:
754 return MPSolver::BASIC;
755 case XPRS_AT_UPPER:
756 return MPSolver::AT_UPPER_BOUND;
757 case XPRS_FREE_SUPER:
758 return MPSolver::FREE;
759 default:
760 LOG(DFATAL) << "Unknown XPRESS basis status";
761 return MPSolver::FREE;
762 }
763 }
764
765 // Returns the basis status of a row.
row_status(int constraint_index) const766 MPSolver::BasisStatus XpressInterface::row_status(int constraint_index) const {
767 if (mMip) {
768 LOG(FATAL) << "Basis status only available for continuous problems";
769 return MPSolver::FREE;
770 }
771
772 if (CheckSolutionIsSynchronized()) {
773 if (!mRstat) {
774 int const rows = XPRSgetnumrows(mLp);
775 unique_ptr<int[]> data(new int[rows]);
776 mRstat.swap(data);
777 CHECK_STATUS(XPRSgetbasis(mLp, 0, mRstat.get()));
778 }
779 } else {
780 mRstat = 0;
781 }
782
783 if (mRstat) {
784 return xformBasisStatus(mRstat[constraint_index]);
785 } else {
786 LOG(FATAL) << "Row basis status not available";
787 return MPSolver::FREE;
788 }
789 }
790
791 // Returns the basis status of a column.
column_status(int variable_index) const792 MPSolver::BasisStatus XpressInterface::column_status(int variable_index) const {
793 if (mMip) {
794 LOG(FATAL) << "Basis status only available for continuous problems";
795 return MPSolver::FREE;
796 }
797
798 if (CheckSolutionIsSynchronized()) {
799 if (!mCstat) {
800 int const cols = XPRSgetnumcols(mLp);
801 unique_ptr<int[]> data(new int[cols]);
802 mCstat.swap(data);
803 CHECK_STATUS(XPRSgetbasis(mLp, mCstat.get(), 0));
804 }
805 } else {
806 mCstat = 0;
807 }
808
809 if (mCstat) {
810 return xformBasisStatus(mCstat[variable_index]);
811 } else {
812 LOG(FATAL) << "Column basis status not available";
813 return MPSolver::FREE;
814 }
815 }
816
817 // Extract all variables that have not yet been extracted.
ExtractNewVariables()818 void XpressInterface::ExtractNewVariables() {
819 // NOTE: The code assumes that a linear expression can never contain
820 // non-zero duplicates.
821
822 InvalidateSolutionSynchronization();
823
824 if (!supportIncrementalExtraction) {
825 // Without incremental extraction ExtractModel() is always called
826 // to extract the full model.
827 CHECK(last_variable_index_ == 0 ||
828 last_variable_index_ == solver_->variables_.size());
829 CHECK(last_constraint_index_ == 0 ||
830 last_constraint_index_ == solver_->constraints_.size());
831 }
832
833 int const last_extracted = last_variable_index_;
834 int const var_count = solver_->variables_.size();
835 int newcols = var_count - last_extracted;
836 if (newcols > 0) {
837 // There are non-extracted variables. Extract them now.
838
839 unique_ptr<double[]> obj(new double[newcols]);
840 unique_ptr<double[]> lb(new double[newcols]);
841 unique_ptr<double[]> ub(new double[newcols]);
842 unique_ptr<char[]> ctype(new char[newcols]);
843 unique_ptr<const char*[]> colname(new const char*[newcols]);
844
845 bool have_names = false;
846 for (int j = 0, varidx = last_extracted; j < newcols; ++j, ++varidx) {
847 MPVariable const* const var = solver_->variables_[varidx];
848 lb[j] = var->lb();
849 ub[j] = var->ub();
850 ctype[j] = var->integer() ? XPRS_INTEGER : XPRS_CONTINUOUS;
851 colname[j] = var->name().empty() ? 0 : var->name().c_str();
852 have_names = have_names || var->name().empty();
853 obj[j] = solver_->objective_->GetCoefficient(var);
854 }
855
856 // Arrays for modifying the problem are setup. Update the index
857 // of variables that will get extracted now. Updating indices
858 // _before_ the actual extraction makes things much simpler in
859 // case we support incremental extraction.
860 // In case of error we just reset the indices.
861 std::vector<MPVariable*> const& variables = solver_->variables();
862 for (int j = last_extracted; j < var_count; ++j) {
863 CHECK(!variable_is_extracted(variables[j]->index()));
864 set_variable_as_extracted(variables[j]->index(), true);
865 }
866
867 try {
868 bool use_newcols = true;
869
870 if (supportIncrementalExtraction) {
871 // If we support incremental extraction then we must
872 // update existing constraints with the new variables.
873 // To do that we use XPRSaddcols() to actually create the
874 // variables. This is supposed to be faster than combining
875 // XPRSnewcols() and XPRSchgcoeflist().
876
877 // For each column count the size of the intersection with
878 // existing constraints.
879 unique_ptr<int[]> collen(new int[newcols]);
880 for (int j = 0; j < newcols; ++j) collen[j] = 0;
881 int nonzeros = 0;
882 // TODO: Use a bitarray to flag the constraints that actually
883 // intersect new variables?
884 for (int i = 0; i < last_constraint_index_; ++i) {
885 MPConstraint const* const ct = solver_->constraints_[i];
886 CHECK(constraint_is_extracted(ct->index()));
887 const auto& coeffs = ct->coefficients_;
888 for (auto it(coeffs.begin()); it != coeffs.end(); ++it) {
889 int const idx = it->first->index();
890 if (variable_is_extracted(idx) && idx > last_variable_index_) {
891 collen[idx - last_variable_index_]++;
892 ++nonzeros;
893 }
894 }
895 }
896
897 if (nonzeros > 0) {
898 // At least one of the new variables did intersect with an
899 // old constraint. We have to create the new columns via
900 // XPRSaddcols().
901 use_newcols = false;
902 unique_ptr<int[]> begin(new int[newcols + 2]);
903 unique_ptr<int[]> cmatind(new int[nonzeros]);
904 unique_ptr<double[]> cmatval(new double[nonzeros]);
905
906 // Here is how cmatbeg[] is setup:
907 // - it is initialized as
908 // [ 0, 0, collen[0], collen[0]+collen[1], ... ]
909 // so that cmatbeg[j+1] tells us where in cmatind[] and
910 // cmatval[] we need to put the next nonzero for column
911 // j
912 // - after nonzeros have been setup the array looks like
913 // [ 0, collen[0], collen[0]+collen[1], ... ]
914 // so that it is the correct input argument for XPRSaddcols
915 int* cmatbeg = begin.get();
916 cmatbeg[0] = 0;
917 cmatbeg[1] = 0;
918 ++cmatbeg;
919 for (int j = 0; j < newcols; ++j)
920 cmatbeg[j + 1] = cmatbeg[j] + collen[j];
921
922 for (int i = 0; i < last_constraint_index_; ++i) {
923 MPConstraint const* const ct = solver_->constraints_[i];
924 int const row = ct->index();
925 const auto& coeffs = ct->coefficients_;
926 for (auto it(coeffs.begin()); it != coeffs.end(); ++it) {
927 int const idx = it->first->index();
928 if (variable_is_extracted(idx) && idx > last_variable_index_) {
929 int const nz = cmatbeg[idx]++;
930 cmatind[nz] = row;
931 cmatval[nz] = it->second;
932 }
933 }
934 }
935 --cmatbeg;
936 CHECK_STATUS(XPRSaddcols(mLp, newcols, nonzeros, obj.get(), cmatbeg,
937 cmatind.get(), cmatval.get(), lb.get(),
938 ub.get()));
939 }
940 }
941
942 if (use_newcols) {
943 // Either incremental extraction is not supported or none of
944 // the new variables did intersect an existing constraint.
945 // We can just use XPRSnewcols() to create the new variables.
946 std::vector<int> collen(newcols, 0);
947 std::vector<int> cmatbeg(newcols, 0);
948 unique_ptr<int[]> cmatind(new int[1]);
949 unique_ptr<double[]> cmatval(new double[1]);
950 cmatind[0] = 0;
951 cmatval[0] = 1.0;
952
953 CHECK_STATUS(XPRSaddcols(mLp, newcols, 0, obj.get(), cmatbeg.data(),
954 cmatind.get(), cmatval.get(), lb.get(),
955 ub.get()));
956 int const cols = XPRSgetnumcols(mLp);
957 unique_ptr<int[]> ind(new int[newcols]);
958 for (int j = 0; j < cols; ++j) ind[j] = j;
959 CHECK_STATUS(
960 XPRSchgcoltype(mLp, cols - last_extracted, ind.get(), ctype.get()));
961 } else {
962 // Incremental extraction: we must update the ctype of the
963 // newly created variables (XPRSaddcols() does not allow
964 // specifying the ctype)
965 if (mMip && XPRSgetnumcols(mLp) > 0) {
966 // Query the actual number of columns in case we did not
967 // manage to extract all columns.
968 int const cols = XPRSgetnumcols(mLp);
969 unique_ptr<int[]> ind(new int[newcols]);
970 for (int j = last_extracted; j < cols; ++j)
971 ind[j - last_extracted] = j;
972 CHECK_STATUS(XPRSchgcoltype(mLp, cols - last_extracted, ind.get(),
973 ctype.get()));
974 }
975 }
976 } catch (...) {
977 // Undo all changes in case of error.
978 int const cols = XPRSgetnumcols(mLp);
979 if (cols > last_extracted) {
980 std::vector<int> colsToDelete;
981 for (int i = last_extracted; i < cols; ++i) colsToDelete.push_back(i);
982 (void)XPRSdelcols(mLp, colsToDelete.size(), colsToDelete.data());
983 }
984 std::vector<MPVariable*> const& variables = solver_->variables();
985 int const size = variables.size();
986 for (int j = last_extracted; j < size; ++j)
987 set_variable_as_extracted(j, false);
988 throw;
989 }
990 }
991 }
992
993 // Extract constraints that have not yet been extracted.
ExtractNewConstraints()994 void XpressInterface::ExtractNewConstraints() {
995 // NOTE: The code assumes that a linear expression can never contain
996 // non-zero duplicates.
997 if (!supportIncrementalExtraction) {
998 // Without incremental extraction ExtractModel() is always called
999 // to extract the full model.
1000 CHECK(last_variable_index_ == 0 ||
1001 last_variable_index_ == solver_->variables_.size());
1002 CHECK(last_constraint_index_ == 0 ||
1003 last_constraint_index_ == solver_->constraints_.size());
1004 }
1005
1006 int const offset = last_constraint_index_;
1007 int const total = solver_->constraints_.size();
1008
1009 if (total > offset) {
1010 // There are constraints that are not yet extracted.
1011
1012 InvalidateSolutionSynchronization();
1013
1014 int newCons = total - offset;
1015 int const cols = XPRSgetnumcols(mLp);
1016 DCHECK_EQ(last_variable_index_, cols);
1017 int const chunk = newCons; // 10; // max number of rows to add in one shot
1018
1019 // Update indices of new constraints _before_ actually extracting
1020 // them. In case of error we will just reset the indices.
1021 for (int c = offset; c < total; ++c) set_constraint_as_extracted(c, true);
1022
1023 try {
1024 unique_ptr<int[]> rmatind(new int[cols]);
1025 unique_ptr<double[]> rmatval(new double[cols]);
1026 unique_ptr<int[]> rmatbeg(new int[chunk]);
1027 unique_ptr<char[]> sense(new char[chunk]);
1028 unique_ptr<double[]> rhs(new double[chunk]);
1029 unique_ptr<char const*[]> name(new char const*[chunk]);
1030 unique_ptr<double[]> rngval(new double[chunk]);
1031 unique_ptr<int[]> rngind(new int[chunk]);
1032 bool haveRanges = false;
1033
1034 // Loop over the new constraints, collecting rows for up to
1035 // CHUNK constraints into the arrays so that adding constraints
1036 // is faster.
1037 for (int c = 0; c < newCons; /* nothing */) {
1038 // Collect up to CHUNK constraints into the arrays.
1039 int nextRow = 0;
1040 int nextNz = 0;
1041 for (/* nothing */; c < newCons && nextRow < chunk; ++c, ++nextRow) {
1042 MPConstraint const* const ct = solver_->constraints_[offset + c];
1043
1044 // Stop if there is not enough room in the arrays
1045 // to add the current constraint.
1046 if (nextNz + ct->coefficients_.size() > cols) {
1047 DCHECK_GT(nextRow, 0);
1048 break;
1049 }
1050
1051 // Setup right-hand side of constraint.
1052 MakeRhs(ct->lb(), ct->ub(), rhs[nextRow], sense[nextRow],
1053 rngval[nextRow]);
1054 haveRanges = haveRanges || (rngval[nextRow] != 0.0);
1055 rngind[nextRow] = offset + c;
1056
1057 // Setup left-hand side of constraint.
1058 rmatbeg[nextRow] = nextNz;
1059 const auto& coeffs = ct->coefficients_;
1060 for (auto it(coeffs.begin()); it != coeffs.end(); ++it) {
1061 int const idx = it->first->index();
1062 if (variable_is_extracted(idx)) {
1063 DCHECK_LT(nextNz, cols);
1064 DCHECK_LT(idx, cols);
1065 rmatind[nextNz] = idx;
1066 rmatval[nextNz] = it->second;
1067 ++nextNz;
1068 }
1069 }
1070
1071 // Finally the name of the constraint.
1072 name[nextRow] = ct->name().empty() ? 0 : ct->name().c_str();
1073 }
1074 if (nextRow > 0) {
1075 CHECK_STATUS(XPRSaddrows(mLp, nextRow, nextNz, sense.get(), rhs.get(),
1076 rngval.get(), rmatbeg.get(), rmatind.get(),
1077 rmatval.get()));
1078 if (haveRanges) {
1079 CHECK_STATUS(
1080 XPRSchgrhsrange(mLp, nextRow, rngind.get(), rngval.get()));
1081 }
1082 }
1083 }
1084 } catch (...) {
1085 // Undo all changes in case of error.
1086 int const rows = XPRSgetnumrows(mLp);
1087 std::vector<int> rowsToDelete;
1088 for (int i = offset; i < rows; ++i) rowsToDelete.push_back(i);
1089 if (rows > offset)
1090 (void)XPRSdelrows(mLp, rowsToDelete.size(), rowsToDelete.data());
1091 std::vector<MPConstraint*> const& constraints = solver_->constraints();
1092 int const size = constraints.size();
1093 for (int i = offset; i < size; ++i) set_constraint_as_extracted(i, false);
1094 throw;
1095 }
1096 }
1097 }
1098
1099 // Extract the objective function.
ExtractObjective()1100 void XpressInterface::ExtractObjective() {
1101 // NOTE: The code assumes that the objective expression does not contain
1102 // any non-zero duplicates.
1103
1104 int const cols = XPRSgetnumcols(mLp);
1105 DCHECK_EQ(last_variable_index_, cols);
1106
1107 unique_ptr<int[]> ind(new int[cols]);
1108 unique_ptr<double[]> val(new double[cols]);
1109 for (int j = 0; j < cols; ++j) {
1110 ind[j] = j;
1111 val[j] = 0.0;
1112 }
1113
1114 const auto& coeffs = solver_->objective_->coefficients_;
1115 for (auto it = coeffs.begin(); it != coeffs.end(); ++it) {
1116 int const idx = it->first->index();
1117 if (variable_is_extracted(idx)) {
1118 DCHECK_LT(idx, cols);
1119 val[idx] = it->second;
1120 }
1121 }
1122
1123 CHECK_STATUS(XPRSchgobj(mLp, cols, ind.get(), val.get()));
1124 CHECK_STATUS(XPRSsetobjoffset(mLp, solver_->Objective().offset()));
1125 }
1126
1127 // ------ Parameters -----
1128
SetParameters(const MPSolverParameters & param)1129 void XpressInterface::SetParameters(const MPSolverParameters& param) {
1130 SetCommonParameters(param);
1131 if (mMip) SetMIPParameters(param);
1132 }
1133
SetRelativeMipGap(double value)1134 void XpressInterface::SetRelativeMipGap(double value) {
1135 if (mMip) {
1136 CHECK_STATUS(XPRSsetdblcontrol(mLp, XPRS_MIPRELSTOP, value));
1137 } else {
1138 LOG(WARNING) << "The relative MIP gap is only available "
1139 << "for discrete problems.";
1140 }
1141 }
1142
SetPrimalTolerance(double value)1143 void XpressInterface::SetPrimalTolerance(double value) {
1144 CHECK_STATUS(XPRSsetdblcontrol(mLp, XPRS_FEASTOL, value));
1145 }
1146
SetDualTolerance(double value)1147 void XpressInterface::SetDualTolerance(double value) {
1148 CHECK_STATUS(XPRSsetdblcontrol(mLp, XPRS_OPTIMALITYTOL, value));
1149 }
1150
SetPresolveMode(int value)1151 void XpressInterface::SetPresolveMode(int value) {
1152 MPSolverParameters::PresolveValues const presolve =
1153 static_cast<MPSolverParameters::PresolveValues>(value);
1154
1155 switch (presolve) {
1156 case MPSolverParameters::PRESOLVE_OFF:
1157 CHECK_STATUS(XPRSsetintcontrol(mLp, XPRS_PRESOLVE, 0));
1158 return;
1159 case MPSolverParameters::PRESOLVE_ON:
1160 CHECK_STATUS(XPRSsetintcontrol(mLp, XPRS_PRESOLVE, 1));
1161 return;
1162 }
1163 SetIntegerParamToUnsupportedValue(MPSolverParameters::PRESOLVE, value);
1164 }
1165
1166 // Sets the scaling mode.
SetScalingMode(int value)1167 void XpressInterface::SetScalingMode(int value) {
1168 MPSolverParameters::ScalingValues const scaling =
1169 static_cast<MPSolverParameters::ScalingValues>(value);
1170
1171 switch (scaling) {
1172 case MPSolverParameters::SCALING_OFF:
1173 CHECK_STATUS(XPRSsetintcontrol(mLp, XPRS_SCALING, 0));
1174 break;
1175 case MPSolverParameters::SCALING_ON:
1176 CHECK_STATUS(XPRSsetdefaultcontrol(mLp, XPRS_SCALING));
1177 // In Xpress, scaling is not a binary on/off control, but a bit vector
1178 // control setting it to 1 would only enable bit 1. Instead we reset it to
1179 // its default (163 for the current version 8.6) Alternatively, we could
1180 // call CHECK_STATUS(XPRSsetintcontrol(mLp, XPRS_SCALING, 163));
1181 break;
1182 }
1183 }
1184
1185 // Sets the LP algorithm : primal, dual or barrier. Note that XPRESS offers
1186 // other LP algorithm (e.g. network) and automatic selection
SetLpAlgorithm(int value)1187 void XpressInterface::SetLpAlgorithm(int value) {
1188 MPSolverParameters::LpAlgorithmValues const algorithm =
1189 static_cast<MPSolverParameters::LpAlgorithmValues>(value);
1190
1191 int alg = 1;
1192
1193 switch (algorithm) {
1194 case MPSolverParameters::DUAL:
1195 alg = 2;
1196 break;
1197 case MPSolverParameters::PRIMAL:
1198 alg = 3;
1199 break;
1200 case MPSolverParameters::BARRIER:
1201 alg = 4;
1202 break;
1203 }
1204
1205 if (alg == XPRS_DEFAULTALG) {
1206 SetIntegerParamToUnsupportedValue(MPSolverParameters::LP_ALGORITHM, value);
1207 } else {
1208 CHECK_STATUS(XPRSsetintcontrol(mLp, XPRS_DEFAULTALG, alg));
1209 }
1210 }
1211
ReadParameterFile(std::string const & filename)1212 bool XpressInterface::ReadParameterFile(std::string const& filename) {
1213 // Return true on success and false on error.
1214 LOG(DFATAL) << "ReadParameterFile not implemented for XPRESS interface";
1215 return false;
1216 }
1217
ValidFileExtensionForParameterFile() const1218 std::string XpressInterface::ValidFileExtensionForParameterFile() const {
1219 return ".prm";
1220 }
1221
Solve(MPSolverParameters const & param)1222 MPSolver::ResultStatus XpressInterface::Solve(MPSolverParameters const& param) {
1223 int status;
1224
1225 // Delete chached information
1226 mCstat = 0;
1227 mRstat = 0;
1228
1229 WallTimer timer;
1230 timer.Start();
1231
1232 // Set incrementality
1233 MPSolverParameters::IncrementalityValues const inc =
1234 static_cast<MPSolverParameters::IncrementalityValues>(
1235 param.GetIntegerParam(MPSolverParameters::INCREMENTALITY));
1236 switch (inc) {
1237 case MPSolverParameters::INCREMENTALITY_OFF: {
1238 Reset(); // This should not be required but re-extracting everything
1239 // may be faster, so we do it.
1240 break;
1241 }
1242 case MPSolverParameters::INCREMENTALITY_ON: {
1243 XPRSsetintcontrol(mLp, XPRS_CRASH, 0);
1244 break;
1245 }
1246 }
1247
1248 // Extract the model to be solved.
1249 // If we don't support incremental extraction and the low-level modeling
1250 // is out of sync then we have to re-extract everything. Note that this
1251 // will lose MIP starts or advanced basis information from a previous
1252 // solve.
1253 if (!supportIncrementalExtraction && sync_status_ == MUST_RELOAD) Reset();
1254 ExtractModel();
1255 VLOG(1) << absl::StrFormat("Model build in %.3f seconds.", timer.Get());
1256
1257 // Set log level.
1258 XPRSsetintcontrol(mLp, XPRS_OUTPUTLOG, quiet() ? 0 : 1);
1259 // Set parameters.
1260 // NOTE: We must invoke SetSolverSpecificParametersAsString() _first_.
1261 // Its current implementation invokes ReadParameterFile() which in
1262 // turn invokes XPRSreadcopyparam(). The latter will _overwrite_
1263 // all current parameter settings in the environment.
1264 solver_->SetSolverSpecificParametersAsString(
1265 solver_->solver_specific_parameter_string_);
1266 SetParameters(param);
1267 if (solver_->time_limit()) {
1268 VLOG(1) << "Setting time limit = " << solver_->time_limit() << " ms.";
1269 // In Xpress, a time limit should usually have a negative sign. With a
1270 // positive sign, the solver will only stop when a solution has been found.
1271 CHECK_STATUS(XPRSsetdblcontrol(mLp, XPRS_MAXTIME,
1272 -1.0 * solver_->time_limit_in_secs()));
1273 }
1274
1275 // Solve.
1276 // Do not CHECK_STATUS here since some errors (for example CPXERR_NO_MEMORY)
1277 // still allow us to query useful information.
1278 timer.Restart();
1279
1280 int xpressstat = 0;
1281 if (mMip) {
1282 if (this->maximize_)
1283 status = XPRSmaxim(mLp, "g");
1284 else
1285 status = XPRSminim(mLp, "g");
1286 XPRSgetintattrib(mLp, XPRS_MIPSTATUS, &xpressstat);
1287 } else {
1288 if (this->maximize_)
1289 status = XPRSmaxim(mLp, "");
1290 else
1291 status = XPRSminim(mLp, "");
1292 XPRSgetintattrib(mLp, XPRS_LPSTATUS, &xpressstat);
1293 }
1294
1295 // Disable screen output right after solve
1296 XPRSsetintcontrol(mLp, XPRS_OUTPUTLOG, 0);
1297
1298 if (status) {
1299 VLOG(1) << absl::StrFormat("Failed to optimize MIP. Error %d", status);
1300 // NOTE: We do not return immediately since there may be information
1301 // to grab (for example an incumbent)
1302 } else {
1303 VLOG(1) << absl::StrFormat("Solved in %.3f seconds.", timer.Get());
1304 }
1305
1306 VLOG(1) << absl::StrFormat("XPRESS solution status %d.", xpressstat);
1307
1308 // Figure out what solution we have.
1309 bool const feasible = (mMip && (xpressstat == XPRS_MIP_OPTIMAL ||
1310 xpressstat == XPRS_MIP_SOLUTION)) ||
1311 (!mMip && xpressstat == XPRS_LP_OPTIMAL);
1312
1313 // Get problem dimensions for solution queries below.
1314 int const rows = XPRSgetnumrows(mLp);
1315 int const cols = XPRSgetnumcols(mLp);
1316 DCHECK_EQ(rows, solver_->constraints_.size());
1317 DCHECK_EQ(cols, solver_->variables_.size());
1318
1319 // Capture objective function value.
1320 objective_value_ = XPRS_NAN;
1321 best_objective_bound_ = XPRS_NAN;
1322 if (feasible) {
1323 if (mMip) {
1324 CHECK_STATUS(XPRSgetdblattrib(mLp, XPRS_MIPOBJVAL, &objective_value_));
1325 CHECK_STATUS(XPRSgetdblattrib(mLp, XPRS_BESTBOUND, &best_objective_bound_));
1326 } else {
1327 CHECK_STATUS(XPRSgetdblattrib(mLp, XPRS_LPOBJVAL, &objective_value_));
1328 }
1329 }
1330 VLOG(1) << "objective=" << objective_value_
1331 << ", bound=" << best_objective_bound_;
1332
1333 // Capture primal and dual solutions
1334 if (mMip) {
1335 // If there is a primal feasible solution then capture it.
1336 if (feasible) {
1337 if (cols > 0) {
1338 unique_ptr<double[]> x(new double[cols]);
1339 CHECK_STATUS(XPRSgetmipsol(mLp, x.get(), 0));
1340 for (int i = 0; i < solver_->variables_.size(); ++i) {
1341 MPVariable* const var = solver_->variables_[i];
1342 var->set_solution_value(x[i]);
1343 VLOG(3) << var->name() << ": value =" << x[i];
1344 }
1345 }
1346 } else {
1347 for (int i = 0; i < solver_->variables_.size(); ++i)
1348 solver_->variables_[i]->set_solution_value(XPRS_NAN);
1349 }
1350
1351 // MIP does not have duals
1352 for (int i = 0; i < solver_->variables_.size(); ++i)
1353 solver_->variables_[i]->set_reduced_cost(XPRS_NAN);
1354 for (int i = 0; i < solver_->constraints_.size(); ++i)
1355 solver_->constraints_[i]->set_dual_value(XPRS_NAN);
1356 } else {
1357 // Continuous problem.
1358 if (cols > 0) {
1359 unique_ptr<double[]> x(new double[cols]);
1360 unique_ptr<double[]> dj(new double[cols]);
1361 if (feasible) CHECK_STATUS(XPRSgetlpsol(mLp, x.get(), 0, 0, dj.get()));
1362 for (int i = 0; i < solver_->variables_.size(); ++i) {
1363 MPVariable* const var = solver_->variables_[i];
1364 var->set_solution_value(x[i]);
1365 bool value = false, dual = false;
1366
1367 if (feasible) {
1368 var->set_solution_value(x[i]);
1369 value = true;
1370 } else {
1371 var->set_solution_value(XPRS_NAN);
1372 }
1373 if (feasible) {
1374 var->set_reduced_cost(dj[i]);
1375 dual = true;
1376 } else {
1377 var->set_reduced_cost(XPRS_NAN);
1378 }
1379 VLOG(3) << var->name() << ":"
1380 << (value ? absl::StrFormat(" value = %f", x[i]) : "")
1381 << (dual ? absl::StrFormat(" reduced cost = %f", dj[i]) : "");
1382 }
1383 }
1384
1385 if (rows > 0) {
1386 unique_ptr<double[]> pi(new double[rows]);
1387 if (feasible) {
1388 CHECK_STATUS(XPRSgetlpsol(mLp, 0, 0, pi.get(), 0));
1389 }
1390 for (int i = 0; i < solver_->constraints_.size(); ++i) {
1391 MPConstraint* const ct = solver_->constraints_[i];
1392 bool dual = false;
1393 if (feasible) {
1394 ct->set_dual_value(pi[i]);
1395 dual = true;
1396 } else {
1397 ct->set_dual_value(XPRS_NAN);
1398 }
1399 VLOG(4) << "row " << ct->index() << ":"
1400 << (dual ? absl::StrFormat(" dual = %f", pi[i]) : "");
1401 }
1402 }
1403 }
1404
1405 // Map XPRESS status to more generic solution status in MPSolver
1406 if (mMip) {
1407 switch (xpressstat) {
1408 case XPRS_MIP_OPTIMAL:
1409 result_status_ = MPSolver::OPTIMAL;
1410 break;
1411 case XPRS_MIP_INFEAS:
1412 result_status_ = MPSolver::INFEASIBLE;
1413 break;
1414 case XPRS_MIP_UNBOUNDED:
1415 result_status_ = MPSolver::UNBOUNDED;
1416 break;
1417 default:
1418 result_status_ = feasible ? MPSolver::FEASIBLE : MPSolver::ABNORMAL;
1419 break;
1420 }
1421 } else {
1422 switch (xpressstat) {
1423 case XPRS_LP_OPTIMAL:
1424 result_status_ = MPSolver::OPTIMAL;
1425 break;
1426 case XPRS_LP_INFEAS:
1427 result_status_ = MPSolver::INFEASIBLE;
1428 break;
1429 case XPRS_LP_UNBOUNDED:
1430 result_status_ = MPSolver::UNBOUNDED;
1431 break;
1432 default:
1433 result_status_ = feasible ? MPSolver::FEASIBLE : MPSolver::ABNORMAL;
1434 break;
1435 }
1436 }
1437
1438 sync_status_ = SOLUTION_SYNCHRONIZED;
1439 return result_status_;
1440 }
1441
BuildXpressInterface(bool mip,MPSolver * const solver)1442 MPSolverInterface* BuildXpressInterface(bool mip, MPSolver* const solver) {
1443 return new XpressInterface(solver, mip);
1444 }
1445
1446 } // namespace operations_research
1447 #endif // #if defined(USE_XPRESS)
1448