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