OR-Tools  8.2
glpk_interface.cc
Go to the documentation of this file.
1// Copyright 2010-2018 Google LLC
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//
15
16#if defined(USE_GLPK)
17
18#include <cmath>
19#include <cstddef>
20#include <limits>
21#include <memory>
22#include <string>
23#include <utility>
24#include <vector>
25
26#include "absl/memory/memory.h"
27#include "absl/strings/str_format.h"
29#include "ortools/base/hash.h"
32#include "ortools/base/timer.h"
34
35extern "C" {
36#include "glpk.h"
37}
38
39namespace operations_research {
40// Class to store information gathered in the callback
41class GLPKInformation {
42 public:
43 explicit GLPKInformation(bool maximize) : num_all_nodes_(0) {
44 ResetBestObjectiveBound(maximize);
45 }
46 void Reset(bool maximize) {
47 num_all_nodes_ = 0;
48 ResetBestObjectiveBound(maximize);
49 }
50 void ResetBestObjectiveBound(bool maximize) {
51 if (maximize) {
52 best_objective_bound_ = std::numeric_limits<double>::infinity();
53 } else {
54 best_objective_bound_ = -std::numeric_limits<double>::infinity();
55 }
56 }
57 int num_all_nodes_;
58 double best_objective_bound_;
59};
60
61// Function to be called in the GLPK callback
62void GLPKGatherInformationCallback(glp_tree* tree, void* info) {
63 CHECK(tree != nullptr);
64 CHECK(info != nullptr);
65 GLPKInformation* glpk_info = reinterpret_cast<GLPKInformation*>(info);
66 switch (glp_ios_reason(tree)) {
67 // The best bound and the number of nodes change only when GLPK
68 // branches, generates cuts or finds an integer solution.
69 case GLP_ISELECT:
70 case GLP_IROWGEN:
71 case GLP_IBINGO: {
72 // Get total number of nodes
73 glp_ios_tree_size(tree, nullptr, nullptr, &glpk_info->num_all_nodes_);
74 // Get best bound
75 int node_id = glp_ios_best_node(tree);
76 if (node_id > 0) {
77 glpk_info->best_objective_bound_ = glp_ios_node_bound(tree, node_id);
78 }
79 break;
80 }
81 default:
82 break;
83 }
84}
85
86// ----- GLPK Solver -----
87
88namespace {
89// GLPK indexes its variables and constraints starting at 1.
90int MPSolverIndexToGlpkIndex(int index) { return index + 1; }
91} // namespace
92
93class GLPKInterface : public MPSolverInterface {
94 public:
95 // Constructor that takes a name for the underlying glpk solver.
96 GLPKInterface(MPSolver* const solver, bool mip);
97 ~GLPKInterface() override;
98
99 // Sets the optimization direction (min/max).
100 void SetOptimizationDirection(bool maximize) override;
101
102 // ----- Solve -----
103 // Solve the problem using the parameter values specified.
104 MPSolver::ResultStatus Solve(const MPSolverParameters& param) override;
105
106 // ----- Model modifications and extraction -----
107 // Resets extracted model
108 void Reset() override;
109
110 // Modify bounds.
111 void SetVariableBounds(int mpsolver_var_index, double lb, double ub) override;
112 void SetVariableInteger(int mpsolver_var_index, bool integer) override;
113 void SetConstraintBounds(int mpsolver_constraint_index, double lb,
114 double ub) override;
115
116 // Add Constraint incrementally.
117 void AddRowConstraint(MPConstraint* const ct) override;
118 // Add variable incrementally.
119 void AddVariable(MPVariable* const var) override;
120 // Change a coefficient in a constraint.
121 void SetCoefficient(MPConstraint* const constraint,
122 const MPVariable* const variable, double new_value,
123 double old_value) override;
124 // Clear a constraint from all its terms.
125 void ClearConstraint(MPConstraint* const constraint) override;
126 // Change a coefficient in the linear objective
127 void SetObjectiveCoefficient(const MPVariable* const variable,
128 double coefficient) override;
129 // Change the constant term in the linear objective.
130 void SetObjectiveOffset(double value) override;
131 // Clear the objective from all its terms.
132 void ClearObjective() override;
133
134 // ------ Query statistics on the solution and the solve ------
135 // Number of simplex iterations
136 int64 iterations() const override;
137 // Number of branch-and-bound nodes. Only available for discrete problems.
138 int64 nodes() const override;
139
140 // Returns the basis status of a row.
141 MPSolver::BasisStatus row_status(int constraint_index) const override;
142 // Returns the basis status of a column.
143 MPSolver::BasisStatus column_status(int variable_index) const override;
144
145 // Checks whether a feasible solution exists.
146 bool CheckSolutionExists() const override;
147
148 // ----- Misc -----
149 // Query problem type.
150 bool IsContinuous() const override { return IsLP(); }
151 bool IsLP() const override { return !mip_; }
152 bool IsMIP() const override { return mip_; }
153
154 void ExtractNewVariables() override;
155 void ExtractNewConstraints() override;
156 void ExtractObjective() override;
157
158 std::string SolverVersion() const override {
159 return absl::StrFormat("GLPK %s", glp_version());
160 }
161
162 void* underlying_solver() override { return reinterpret_cast<void*>(lp_); }
163
164 double ComputeExactConditionNumber() const override;
165
166 private:
167 // Configure the solver's parameters.
168 void ConfigureGLPKParameters(const MPSolverParameters& param);
169
170 // Set all parameters in the underlying solver.
171 void SetParameters(const MPSolverParameters& param) override;
172 // Set each parameter in the underlying solver.
173 void SetRelativeMipGap(double value) override;
174 void SetPrimalTolerance(double value) override;
175 void SetDualTolerance(double value) override;
176 void SetPresolveMode(int value) override;
177 void SetScalingMode(int value) override;
178 void SetLpAlgorithm(int value) override;
179
180 void ExtractOldConstraints();
181 void ExtractOneConstraint(MPConstraint* const constraint, int* const indices,
182 double* const coefs);
183 // Transforms basis status from GLPK integer code to MPSolver::BasisStatus.
184 MPSolver::BasisStatus TransformGLPKBasisStatus(int glpk_basis_status) const;
185
186 // Computes the L1-norm of the current scaled basis.
187 // The L1-norm |A| is defined as max_j sum_i |a_ij|
188 // This method is available only for continuous problems.
189 double ComputeScaledBasisL1Norm(int num_rows, int num_cols,
190 double* row_scaling_factor,
191 double* column_scaling_factor) const;
192
193 // Computes the L1-norm of the inverse of the current scaled
194 // basis.
195 // This method is available only for continuous problems.
196 double ComputeInverseScaledBasisL1Norm(int num_rows, int num_cols,
197 double* row_scaling_factor,
198 double* column_scaling_factor) const;
199
200 glp_prob* lp_;
201 bool mip_;
202
203 // Parameters
204 glp_smcp lp_param_;
205 glp_iocp mip_param_;
206 // For the callback
207 std::unique_ptr<GLPKInformation> mip_callback_info_;
208};
209
210// Creates a LP/MIP instance with the specified name and minimization objective.
211GLPKInterface::GLPKInterface(MPSolver* const solver, bool mip)
212 : MPSolverInterface(solver), lp_(nullptr), mip_(mip) {
213 lp_ = glp_create_prob();
214 glp_set_prob_name(lp_, solver_->name_.c_str());
215 glp_set_obj_dir(lp_, GLP_MIN);
216 mip_callback_info_ = absl::make_unique<GLPKInformation>(maximize_);
217}
218
219// Frees the LP memory allocations.
220GLPKInterface::~GLPKInterface() {
221 CHECK(lp_ != nullptr);
222 glp_delete_prob(lp_);
223 lp_ = nullptr;
224}
225
226void GLPKInterface::Reset() {
227 CHECK(lp_ != nullptr);
228 glp_delete_prob(lp_);
229 lp_ = glp_create_prob();
230 glp_set_prob_name(lp_, solver_->name_.c_str());
231 glp_set_obj_dir(lp_, maximize_ ? GLP_MAX : GLP_MIN);
232 ResetExtractionInformation();
233}
234
235// ------ Model modifications and extraction -----
236
237// Not cached
238void GLPKInterface::SetOptimizationDirection(bool maximize) {
239 InvalidateSolutionSynchronization();
240 glp_set_obj_dir(lp_, maximize ? GLP_MAX : GLP_MIN);
241}
242
243void GLPKInterface::SetVariableBounds(int mpsolver_var_index, double lb,
244 double ub) {
245 InvalidateSolutionSynchronization();
246 if (!variable_is_extracted(mpsolver_var_index)) {
247 sync_status_ = MUST_RELOAD;
248 return;
249 }
250 // Not cached if the variable has been extracted.
251 DCHECK(lp_ != nullptr);
252 const double infinity = solver_->infinity();
253 const int glpk_var_index = MPSolverIndexToGlpkIndex(mpsolver_var_index);
254 if (lb != -infinity) {
255 if (ub != infinity) {
256 if (lb == ub) {
257 glp_set_col_bnds(lp_, glpk_var_index, GLP_FX, lb, ub);
258 } else {
259 glp_set_col_bnds(lp_, glpk_var_index, GLP_DB, lb, ub);
260 }
261 } else {
262 glp_set_col_bnds(lp_, glpk_var_index, GLP_LO, lb, 0.0);
263 }
264 } else if (ub != infinity) {
265 glp_set_col_bnds(lp_, glpk_var_index, GLP_UP, 0.0, ub);
266 } else {
267 glp_set_col_bnds(lp_, glpk_var_index, GLP_FR, 0.0, 0.0);
268 }
269}
270
271void GLPKInterface::SetVariableInteger(int mpsolver_var_index, bool integer) {
272 InvalidateSolutionSynchronization();
273 if (mip_) {
274 if (variable_is_extracted(mpsolver_var_index)) {
275 // Not cached if the variable has been extracted.
276 glp_set_col_kind(lp_, MPSolverIndexToGlpkIndex(mpsolver_var_index),
277 integer ? GLP_IV : GLP_CV);
278 } else {
279 sync_status_ = MUST_RELOAD;
280 }
281 }
282}
283
284void GLPKInterface::SetConstraintBounds(int mpsolver_constraint_index,
285 double lb, double ub) {
286 InvalidateSolutionSynchronization();
287 if (!constraint_is_extracted(mpsolver_constraint_index)) {
288 sync_status_ = MUST_RELOAD;
289 return;
290 }
291 // Not cached if the row has been extracted
292 const int glpk_constraint_index =
293 MPSolverIndexToGlpkIndex(mpsolver_constraint_index);
294 DCHECK(lp_ != nullptr);
295 const double infinity = solver_->infinity();
296 if (lb != -infinity) {
297 if (ub != infinity) {
298 if (lb == ub) {
299 glp_set_row_bnds(lp_, glpk_constraint_index, GLP_FX, lb, ub);
300 } else {
301 glp_set_row_bnds(lp_, glpk_constraint_index, GLP_DB, lb, ub);
302 }
303 } else {
304 glp_set_row_bnds(lp_, glpk_constraint_index, GLP_LO, lb, 0.0);
305 }
306 } else if (ub != infinity) {
307 glp_set_row_bnds(lp_, glpk_constraint_index, GLP_UP, 0.0, ub);
308 } else {
309 glp_set_row_bnds(lp_, glpk_constraint_index, GLP_FR, 0.0, 0.0);
310 }
311}
312
313void GLPKInterface::SetCoefficient(MPConstraint* const constraint,
314 const MPVariable* const variable,
315 double new_value, double old_value) {
316 InvalidateSolutionSynchronization();
317 // GLPK does not allow to modify one coefficient at a time, so we
318 // extract the whole constraint again, if it has been extracted
319 // already and if it does not contain new variables. Otherwise, we
320 // cache the modification.
321 if (constraint_is_extracted(constraint->index()) &&
322 (sync_status_ == MODEL_SYNCHRONIZED ||
323 !constraint->ContainsNewVariables())) {
324 const int size = constraint->coefficients_.size();
325 std::unique_ptr<int[]> indices(new int[size + 1]);
326 std::unique_ptr<double[]> coefs(new double[size + 1]);
327 ExtractOneConstraint(constraint, indices.get(), coefs.get());
328 }
329}
330
331// Not cached
332void GLPKInterface::ClearConstraint(MPConstraint* const constraint) {
333 InvalidateSolutionSynchronization();
334 // Constraint may have not been extracted yet.
335 if (constraint_is_extracted(constraint->index())) {
336 glp_set_mat_row(lp_, MPSolverIndexToGlpkIndex(constraint->index()), 0,
337 nullptr, nullptr);
338 }
339}
340
341// Cached
342void GLPKInterface::SetObjectiveCoefficient(const MPVariable* const variable,
343 double coefficient) {
344 sync_status_ = MUST_RELOAD;
345}
346
347// Cached
348void GLPKInterface::SetObjectiveOffset(double value) {
349 sync_status_ = MUST_RELOAD;
350}
351
352// Clear objective of all its terms (linear)
353void GLPKInterface::ClearObjective() {
354 InvalidateSolutionSynchronization();
355 for (const auto& entry : solver_->objective_->coefficients_) {
356 const int mpsolver_var_index = entry.first->index();
357 // Variable may have not been extracted yet.
358 if (!variable_is_extracted(mpsolver_var_index)) {
359 DCHECK_NE(MODEL_SYNCHRONIZED, sync_status_);
360 } else {
361 glp_set_obj_coef(lp_, MPSolverIndexToGlpkIndex(mpsolver_var_index), 0.0);
362 }
363 }
364 // Constant term.
365 glp_set_obj_coef(lp_, 0, 0.0);
366}
367
368void GLPKInterface::AddRowConstraint(MPConstraint* const ct) {
369 sync_status_ = MUST_RELOAD;
370}
371
372void GLPKInterface::AddVariable(MPVariable* const var) {
373 sync_status_ = MUST_RELOAD;
374}
375
376// Define new variables and add them to existing constraints.
377void GLPKInterface::ExtractNewVariables() {
378 int total_num_vars = solver_->variables_.size();
379 if (total_num_vars > last_variable_index_) {
380 glp_add_cols(lp_, total_num_vars - last_variable_index_);
381 for (int j = last_variable_index_; j < solver_->variables_.size(); ++j) {
382 MPVariable* const var = solver_->variables_[j];
383 set_variable_as_extracted(j, true);
384 if (!var->name().empty()) {
385 glp_set_col_name(lp_, MPSolverIndexToGlpkIndex(j), var->name().c_str());
386 }
387 SetVariableBounds(/*mpsolver_var_index=*/j, var->lb(), var->ub());
388 SetVariableInteger(/*mpsolver_var_index=*/j, var->integer());
389
390 // The true objective coefficient will be set later in ExtractObjective.
391 double tmp_obj_coef = 0.0;
392 glp_set_obj_coef(lp_, MPSolverIndexToGlpkIndex(j), tmp_obj_coef);
393 }
394 // Add new variables to the existing constraints.
395 ExtractOldConstraints();
396 }
397}
398
399// Extract again existing constraints if they contain new variables.
400void GLPKInterface::ExtractOldConstraints() {
401 const int max_constraint_size =
402 solver_->ComputeMaxConstraintSize(0, last_constraint_index_);
403 // The first entry in the following arrays is dummy, to be
404 // consistent with glpk API.
405 std::unique_ptr<int[]> indices(new int[max_constraint_size + 1]);
406 std::unique_ptr<double[]> coefs(new double[max_constraint_size + 1]);
407
408 for (int i = 0; i < last_constraint_index_; ++i) {
409 MPConstraint* const ct = solver_->constraints_[i];
410 DCHECK(constraint_is_extracted(i));
411 const int size = ct->coefficients_.size();
412 if (size == 0) {
413 continue;
414 }
415 // Update the constraint's coefficients if it contains new variables.
416 if (ct->ContainsNewVariables()) {
417 ExtractOneConstraint(ct, indices.get(), coefs.get());
418 }
419 }
420}
421
422// Extract one constraint. Arrays indices and coefs must be
423// preallocated to have enough space to contain the constraint's
424// coefficients.
425void GLPKInterface::ExtractOneConstraint(MPConstraint* const constraint,
426 int* const indices,
427 double* const coefs) {
428 // GLPK convention is to start indexing at 1.
429 int k = 1;
430 for (const auto& entry : constraint->coefficients_) {
431 DCHECK(variable_is_extracted(entry.first->index()));
432 indices[k] = MPSolverIndexToGlpkIndex(entry.first->index());
433 coefs[k] = entry.second;
434 ++k;
435 }
436 glp_set_mat_row(lp_, MPSolverIndexToGlpkIndex(constraint->index()), k - 1,
437 indices, coefs);
438}
439
440// Define new constraints on old and new variables.
441void GLPKInterface::ExtractNewConstraints() {
442 int total_num_rows = solver_->constraints_.size();
443 if (last_constraint_index_ < total_num_rows) {
444 // Define new constraints
445 glp_add_rows(lp_, total_num_rows - last_constraint_index_);
446 int num_coefs = 0;
447 for (int i = last_constraint_index_; i < total_num_rows; ++i) {
448 MPConstraint* ct = solver_->constraints_[i];
449 set_constraint_as_extracted(i, true);
450 if (ct->name().empty()) {
451 glp_set_row_name(lp_, MPSolverIndexToGlpkIndex(i),
452 absl::StrFormat("ct_%i", i).c_str());
453 } else {
454 glp_set_row_name(lp_, MPSolverIndexToGlpkIndex(i), ct->name().c_str());
455 }
456 // All constraints are set to be of the type <= limit_ .
457 SetConstraintBounds(/*mpsolver_constraint_index=*/i, ct->lb(), ct->ub());
458 num_coefs += ct->coefficients_.size();
459 }
460
461 // Fill new constraints with coefficients
462 if (last_variable_index_ == 0 && last_constraint_index_ == 0) {
463 // Faster extraction when nothing has been extracted yet: build
464 // and load whole matrix at once instead of constructing rows
465 // separately.
466
467 // The first entry in the following arrays is dummy, to be
468 // consistent with glpk API.
469 std::unique_ptr<int[]> variable_indices(new int[num_coefs + 1]);
470 std::unique_ptr<int[]> constraint_indices(new int[num_coefs + 1]);
471 std::unique_ptr<double[]> coefs(new double[num_coefs + 1]);
472 int k = 1;
473 for (int i = 0; i < solver_->constraints_.size(); ++i) {
474 MPConstraint* ct = solver_->constraints_[i];
475 for (const auto& entry : ct->coefficients_) {
476 DCHECK(variable_is_extracted(entry.first->index()));
477 constraint_indices[k] = MPSolverIndexToGlpkIndex(ct->index());
478 variable_indices[k] = MPSolverIndexToGlpkIndex(entry.first->index());
479 coefs[k] = entry.second;
480 ++k;
481 }
482 }
483 CHECK_EQ(num_coefs + 1, k);
484 glp_load_matrix(lp_, num_coefs, constraint_indices.get(),
485 variable_indices.get(), coefs.get());
486 } else {
487 // Build each new row separately.
488 int max_constraint_size = solver_->ComputeMaxConstraintSize(
489 last_constraint_index_, total_num_rows);
490 // The first entry in the following arrays is dummy, to be
491 // consistent with glpk API.
492 std::unique_ptr<int[]> indices(new int[max_constraint_size + 1]);
493 std::unique_ptr<double[]> coefs(new double[max_constraint_size + 1]);
494 for (int i = last_constraint_index_; i < total_num_rows; i++) {
495 ExtractOneConstraint(solver_->constraints_[i], indices.get(),
496 coefs.get());
497 }
498 }
499 }
500}
501
502void GLPKInterface::ExtractObjective() {
503 // Linear objective: set objective coefficients for all variables
504 // (some might have been modified).
505 for (const auto& entry : solver_->objective_->coefficients_) {
506 glp_set_obj_coef(lp_, MPSolverIndexToGlpkIndex(entry.first->index()),
507 entry.second);
508 }
509 // Constant term.
510 glp_set_obj_coef(lp_, 0, solver_->Objective().offset());
511}
512
513// Solve the problem using the parameter values specified.
514MPSolver::ResultStatus GLPKInterface::Solve(const MPSolverParameters& param) {
515 WallTimer timer;
516 timer.Start();
517
518 // Note that GLPK provides incrementality for LP but not for MIP.
519 if (param.GetIntegerParam(MPSolverParameters::INCREMENTALITY) ==
520 MPSolverParameters::INCREMENTALITY_OFF) {
521 Reset();
522 }
523
524 // Set log level.
525 if (quiet_) {
526 glp_term_out(GLP_OFF);
527 } else {
528 glp_term_out(GLP_ON);
529 }
530
531 ExtractModel();
532 VLOG(1) << absl::StrFormat("Model built in %.3f seconds.", timer.Get());
533
534 // Configure parameters at every solve, even when the model has not
535 // been changed, in case some of the parameters such as the time
536 // limit have been changed since the last solve.
537 ConfigureGLPKParameters(param);
538
539 // Solve
540 timer.Restart();
541 int solver_status = glp_simplex(lp_, &lp_param_);
542 if (mip_) {
543 // glp_intopt requires to solve the root LP separately.
544 // If the root LP was solved successfully, solve the MIP.
545 if (solver_status == 0) {
546 solver_status = glp_intopt(lp_, &mip_param_);
547 } else {
548 // Something abnormal occurred during the root LP solve. It is
549 // highly unlikely that an integer feasible solution is
550 // available at this point, so we don't put any effort in trying
551 // to recover it.
552 result_status_ = MPSolver::ABNORMAL;
553 if (solver_status == GLP_ETMLIM) {
554 result_status_ = MPSolver::NOT_SOLVED;
555 }
556 sync_status_ = SOLUTION_SYNCHRONIZED;
557 return result_status_;
558 }
559 }
560 VLOG(1) << absl::StrFormat("GLPK Status: %i (time spent: %.3f seconds).",
561 solver_status, timer.Get());
562
563 // Get the results.
564 if (mip_) {
565 objective_value_ = glp_mip_obj_val(lp_);
566 best_objective_bound_ = mip_callback_info_->best_objective_bound_;
567 } else {
568 objective_value_ = glp_get_obj_val(lp_);
569 }
570 VLOG(1) << "objective=" << objective_value_
571 << ", bound=" << best_objective_bound_;
572 for (int i = 0; i < solver_->variables_.size(); ++i) {
573 MPVariable* const var = solver_->variables_[i];
574 double val;
575 if (mip_) {
576 val = glp_mip_col_val(lp_, MPSolverIndexToGlpkIndex(i));
577 } else {
578 val = glp_get_col_prim(lp_, MPSolverIndexToGlpkIndex(i));
579 }
580 var->set_solution_value(val);
581 VLOG(3) << var->name() << ": value =" << val;
582 if (!mip_) {
583 double reduced_cost;
584 reduced_cost = glp_get_col_dual(lp_, MPSolverIndexToGlpkIndex(i));
585 var->set_reduced_cost(reduced_cost);
586 VLOG(4) << var->name() << ": reduced cost = " << reduced_cost;
587 }
588 }
589 for (int i = 0; i < solver_->constraints_.size(); ++i) {
590 MPConstraint* const ct = solver_->constraints_[i];
591 if (!mip_) {
592 const double dual_value =
593 glp_get_row_dual(lp_, MPSolverIndexToGlpkIndex(i));
594 ct->set_dual_value(dual_value);
595 VLOG(4) << "row " << MPSolverIndexToGlpkIndex(i)
596 << ": dual value = " << dual_value;
597 }
598 }
599
600 // Check the status: optimal, infeasible, etc.
601 if (mip_) {
602 int tmp_status = glp_mip_status(lp_);
603 VLOG(1) << "GLPK result status: " << tmp_status;
604 if (tmp_status == GLP_OPT) {
605 result_status_ = MPSolver::OPTIMAL;
606 } else if (tmp_status == GLP_FEAS) {
607 result_status_ = MPSolver::FEASIBLE;
608 } else if (tmp_status == GLP_NOFEAS) {
609 // For infeasible problems, GLPK actually seems to return
610 // GLP_UNDEF. So this is never (?) reached. Return infeasible
611 // in case GLPK returns a correct status in future versions.
612 result_status_ = MPSolver::INFEASIBLE;
613 } else if (solver_status == GLP_ETMLIM) {
614 result_status_ = MPSolver::NOT_SOLVED;
615 } else {
616 result_status_ = MPSolver::ABNORMAL;
617 // GLPK does not have a status code for unbounded MIP models, so
618 // we return an abnormal status instead.
619 }
620 } else {
621 int tmp_status = glp_get_status(lp_);
622 VLOG(1) << "GLPK result status: " << tmp_status;
623 if (tmp_status == GLP_OPT) {
624 result_status_ = MPSolver::OPTIMAL;
625 } else if (tmp_status == GLP_FEAS) {
626 result_status_ = MPSolver::FEASIBLE;
627 } else if (tmp_status == GLP_NOFEAS || tmp_status == GLP_INFEAS) {
628 // For infeasible problems, GLPK actually seems to return
629 // GLP_UNDEF. So this is never (?) reached. Return infeasible
630 // in case GLPK returns a correct status in future versions.
631 result_status_ = MPSolver::INFEASIBLE;
632 } else if (tmp_status == GLP_UNBND) {
633 // For unbounded problems, GLPK actually seems to return
634 // GLP_UNDEF. So this is never (?) reached. Return unbounded
635 // in case GLPK returns a correct status in future versions.
636 result_status_ = MPSolver::UNBOUNDED;
637 } else if (solver_status == GLP_ETMLIM) {
638 result_status_ = MPSolver::NOT_SOLVED;
639 } else {
640 result_status_ = MPSolver::ABNORMAL;
641 }
642 }
643
644 sync_status_ = SOLUTION_SYNCHRONIZED;
645
646 return result_status_;
647}
648
649MPSolver::BasisStatus GLPKInterface::TransformGLPKBasisStatus(
650 int glpk_basis_status) const {
651 switch (glpk_basis_status) {
652 case GLP_BS:
653 return MPSolver::BASIC;
654 case GLP_NL:
655 return MPSolver::AT_LOWER_BOUND;
656 case GLP_NU:
657 return MPSolver::AT_UPPER_BOUND;
658 case GLP_NF:
659 return MPSolver::FREE;
660 case GLP_NS:
661 return MPSolver::FIXED_VALUE;
662 default:
663 LOG(FATAL) << "Unknown GLPK basis status";
664 return MPSolver::FREE;
665 }
666}
667
668// ------ Query statistics on the solution and the solve ------
669
670int64 GLPKInterface::iterations() const {
671#if GLP_MAJOR_VERSION == 4 && GLP_MINOR_VERSION < 49
672 if (!mip_ && CheckSolutionIsSynchronized()) {
673 return lpx_get_int_parm(lp_, LPX_K_ITCNT);
674 }
675#elif GLP_MAJOR_VERSION == 4 && GLP_MINOR_VERSION >= 53
676 if (!mip_ && CheckSolutionIsSynchronized()) {
677 return glp_get_it_cnt(lp_);
678 }
679#endif
680 LOG(WARNING) << "Total number of iterations is not available";
681 return kUnknownNumberOfIterations;
682}
683
684int64 GLPKInterface::nodes() const {
685 if (mip_) {
686 if (!CheckSolutionIsSynchronized()) return kUnknownNumberOfNodes;
687 return mip_callback_info_->num_all_nodes_;
688 } else {
689 LOG(DFATAL) << "Number of nodes only available for discrete problems";
690 return kUnknownNumberOfNodes;
691 }
692}
693
694MPSolver::BasisStatus GLPKInterface::row_status(int constraint_index) const {
695 DCHECK_GE(constraint_index, 0);
696 DCHECK_LT(constraint_index, last_constraint_index_);
697 const int glpk_basis_status =
698 glp_get_row_stat(lp_, MPSolverIndexToGlpkIndex(constraint_index));
699 return TransformGLPKBasisStatus(glpk_basis_status);
700}
701
702MPSolver::BasisStatus GLPKInterface::column_status(int variable_index) const {
703 DCHECK_GE(variable_index, 0);
704 DCHECK_LT(variable_index, last_variable_index_);
705 const int glpk_basis_status =
706 glp_get_col_stat(lp_, MPSolverIndexToGlpkIndex(variable_index));
707 return TransformGLPKBasisStatus(glpk_basis_status);
708}
709
710bool GLPKInterface::CheckSolutionExists() const {
711 if (result_status_ == MPSolver::ABNORMAL) {
712 LOG(WARNING) << "Ignoring ABNORMAL status from GLPK: This status may or may"
713 << " not indicate that a solution exists.";
714 return true;
715 } else {
716 // Call default implementation
717 return MPSolverInterface::CheckSolutionExists();
718 }
719}
720
721double GLPKInterface::ComputeExactConditionNumber() const {
722 if (!IsContinuous()) {
723 // TODO(user): support MIP.
724 LOG(DFATAL) << "ComputeExactConditionNumber not implemented for"
725 << " GLPK_MIXED_INTEGER_PROGRAMMING";
726 return 0.0;
727 }
728 if (!CheckSolutionIsSynchronized()) return 0.0;
729 // Simplex is the only LP algorithm supported in the wrapper for
730 // GLPK, so when a solution exists, a basis exists.
731 CheckSolutionExists();
732 const int num_rows = glp_get_num_rows(lp_);
733 const int num_cols = glp_get_num_cols(lp_);
734 // GLPK indexes everything starting from 1 instead of 0.
735 std::unique_ptr<double[]> row_scaling_factor(new double[num_rows + 1]);
736 std::unique_ptr<double[]> column_scaling_factor(new double[num_cols + 1]);
737 for (int row = 1; row <= num_rows; ++row) {
738 row_scaling_factor[row] = glp_get_rii(lp_, row);
739 }
740 for (int col = 1; col <= num_cols; ++col) {
741 column_scaling_factor[col] = glp_get_sjj(lp_, col);
742 }
743 return ComputeInverseScaledBasisL1Norm(num_rows, num_cols,
744 row_scaling_factor.get(),
745 column_scaling_factor.get()) *
746 ComputeScaledBasisL1Norm(num_rows, num_cols, row_scaling_factor.get(),
747 column_scaling_factor.get());
748}
749
750double GLPKInterface::ComputeScaledBasisL1Norm(
751 int num_rows, int num_cols, double* row_scaling_factor,
752 double* column_scaling_factor) const {
753 double norm = 0.0;
754 std::unique_ptr<double[]> values(new double[num_rows + 1]);
755 std::unique_ptr<int[]> indices(new int[num_rows + 1]);
756 for (int col = 1; col <= num_cols; ++col) {
757 const int glpk_basis_status = glp_get_col_stat(lp_, col);
758 // Take into account only basic columns.
759 if (glpk_basis_status == GLP_BS) {
760 // Compute L1-norm of column 'col': sum_row |a_row,col|.
761 const int num_nz = glp_get_mat_col(lp_, col, indices.get(), values.get());
762 double column_norm = 0.0;
763 for (int k = 1; k <= num_nz; k++) {
764 column_norm += fabs(values[k] * row_scaling_factor[indices[k]]);
765 }
766 column_norm *= fabs(column_scaling_factor[col]);
767 // Compute max_col column_norm
768 norm = std::max(norm, column_norm);
769 }
770 }
771 // Slack variables.
772 for (int row = 1; row <= num_rows; ++row) {
773 const int glpk_basis_status = glp_get_row_stat(lp_, row);
774 // Take into account only basic slack variables.
775 if (glpk_basis_status == GLP_BS) {
776 // Only one non-zero coefficient: +/- 1.0 in the corresponding
777 // row. The row has a scaling coefficient but the slack variable
778 // is never scaled on top of that.
779 const double column_norm = fabs(row_scaling_factor[row]);
780 // Compute max_col column_norm
781 norm = std::max(norm, column_norm);
782 }
783 }
784 return norm;
785}
786
787double GLPKInterface::ComputeInverseScaledBasisL1Norm(
788 int num_rows, int num_cols, double* row_scaling_factor,
789 double* column_scaling_factor) const {
790 // Compute the LU factorization if it doesn't exist yet.
791 if (!glp_bf_exists(lp_)) {
792 const int factorize_status = glp_factorize(lp_);
793 switch (factorize_status) {
794 case GLP_EBADB: {
795 LOG(FATAL) << "Not able to factorize: error GLP_EBADB.";
796 break;
797 }
798 case GLP_ESING: {
799 LOG(WARNING)
800 << "Not able to factorize: "
801 << "the basis matrix is singular within the working precision.";
802 return MPSolver::infinity();
803 }
804 case GLP_ECOND: {
805 LOG(WARNING)
806 << "Not able to factorize: the basis matrix is ill-conditioned.";
807 return MPSolver::infinity();
808 }
809 default:
810 break;
811 }
812 }
813 std::unique_ptr<double[]> right_hand_side(new double[num_rows + 1]);
814 double norm = 0.0;
815 // Iteratively solve B x = e_k, where e_k is the kth unit vector.
816 // The result of this computation is the kth column of B^-1.
817 // glp_ftran works on original matrix. Scale input and result to
818 // obtain the norm of the kth column in the inverse scaled
819 // matrix. See glp_ftran documentation in glpapi12.c for how the
820 // scaling is done: inv(B'') = inv(SB) * inv(B) * inv(R) where:
821 // o B'' is the scaled basis
822 // o B is the original basis
823 // o R is the diagonal row scaling matrix
824 // o SB consists of the basic columns of the augmented column
825 // scaling matrix (auxiliary variables then structural variables):
826 // S~ = diag(inv(R) | S).
827 for (int k = 1; k <= num_rows; ++k) {
828 for (int row = 1; row <= num_rows; ++row) {
829 right_hand_side[row] = 0.0;
830 }
831 right_hand_side[k] = 1.0;
832 // Multiply input by inv(R).
833 for (int row = 1; row <= num_rows; ++row) {
834 right_hand_side[row] /= row_scaling_factor[row];
835 }
836 glp_ftran(lp_, right_hand_side.get());
837 // glp_ftran stores the result in the same vector where the right
838 // hand side was provided.
839 // Multiply result by inv(SB).
840 for (int row = 1; row <= num_rows; ++row) {
841 const int k = glp_get_bhead(lp_, row);
842 if (k <= num_rows) {
843 // Auxiliary variable.
844 right_hand_side[row] *= row_scaling_factor[k];
845 } else {
846 // Structural variable.
847 right_hand_side[row] /= column_scaling_factor[k - num_rows];
848 }
849 }
850 // Compute sum_row |vector_row|.
851 double column_norm = 0.0;
852 for (int row = 1; row <= num_rows; ++row) {
853 column_norm += fabs(right_hand_side[row]);
854 }
855 // Compute max_col column_norm
856 norm = std::max(norm, column_norm);
857 }
858 return norm;
859}
860
861// ------ Parameters ------
862
863void GLPKInterface::ConfigureGLPKParameters(const MPSolverParameters& param) {
864 if (mip_) {
865 glp_init_iocp(&mip_param_);
866 // Time limit
867 if (solver_->time_limit()) {
868 VLOG(1) << "Setting time limit = " << solver_->time_limit() << " ms.";
869 mip_param_.tm_lim = solver_->time_limit();
870 }
871 // Initialize structures related to the callback.
872 mip_param_.cb_func = GLPKGatherInformationCallback;
873 mip_callback_info_->Reset(maximize_);
874 mip_param_.cb_info = mip_callback_info_.get();
875 // TODO(user): switch some cuts on? All cuts are off by default!?
876 }
877
878 // Configure LP parameters in all cases since they will be used to
879 // solve the root LP in the MIP case.
880 glp_init_smcp(&lp_param_);
881 // Time limit
882 if (solver_->time_limit()) {
883 VLOG(1) << "Setting time limit = " << solver_->time_limit() << " ms.";
884 lp_param_.tm_lim = solver_->time_limit();
885 }
886
887 // Should give a numerically better representation of the problem.
888 glp_scale_prob(lp_, GLP_SF_AUTO);
889
890 // Use advanced initial basis (options: standard / advanced / Bixby's).
891 glp_adv_basis(lp_, 0);
892
893 // Set parameters specified by the user.
894 SetParameters(param);
895}
896
897void GLPKInterface::SetParameters(const MPSolverParameters& param) {
898 SetCommonParameters(param);
899 if (mip_) {
900 SetMIPParameters(param);
901 }
902}
903
904void GLPKInterface::SetRelativeMipGap(double value) {
905 if (mip_) {
906 mip_param_.mip_gap = value;
907 } else {
908 LOG(WARNING) << "The relative MIP gap is only available "
909 << "for discrete problems.";
910 }
911}
912
913void GLPKInterface::SetPrimalTolerance(double value) {
914 lp_param_.tol_bnd = value;
915}
916
917void GLPKInterface::SetDualTolerance(double value) { lp_param_.tol_dj = value; }
918
919void GLPKInterface::SetPresolveMode(int value) {
920 switch (value) {
921 case MPSolverParameters::PRESOLVE_OFF: {
922 mip_param_.presolve = GLP_OFF;
923 lp_param_.presolve = GLP_OFF;
924 break;
925 }
926 case MPSolverParameters::PRESOLVE_ON: {
927 mip_param_.presolve = GLP_ON;
928 lp_param_.presolve = GLP_ON;
929 break;
930 }
931 default: {
932 SetIntegerParamToUnsupportedValue(MPSolverParameters::PRESOLVE, value);
933 }
934 }
935}
936
937void GLPKInterface::SetScalingMode(int value) {
938 SetUnsupportedIntegerParam(MPSolverParameters::SCALING);
939}
940
941void GLPKInterface::SetLpAlgorithm(int value) {
942 switch (value) {
943 case MPSolverParameters::DUAL: {
944 // Use dual, and if it fails, switch to primal.
945 lp_param_.meth = GLP_DUALP;
946 break;
947 }
948 case MPSolverParameters::PRIMAL: {
949 lp_param_.meth = GLP_PRIMAL;
950 break;
951 }
952 case MPSolverParameters::BARRIER:
953 default: {
954 SetIntegerParamToUnsupportedValue(MPSolverParameters::LP_ALGORITHM,
955 value);
956 }
957 }
958}
959
960MPSolverInterface* BuildGLPKInterface(bool mip, MPSolver* const solver) {
961 return new GLPKInterface(solver, mip);
962}
963
964} // namespace operations_research
965#endif // #if defined(USE_GLPK)
int64 max
Definition: alldiff_cst.cc:139
#define CHECK(condition)
Definition: base/logging.h:495
#define DCHECK_NE(val1, val2)
Definition: base/logging.h:886
#define CHECK_EQ(val1, val2)
Definition: base/logging.h:697
#define DCHECK_GE(val1, val2)
Definition: base/logging.h:889
#define DCHECK_LT(val1, val2)
Definition: base/logging.h:888
#define LOG(severity)
Definition: base/logging.h:420
#define DCHECK(condition)
Definition: base/logging.h:884
#define VLOG(verboselevel)
Definition: base/logging.h:978
void Start()
Definition: timer.h:31
void Restart()
Definition: timer.h:35
double Get() const
Definition: timer.h:45
ResultStatus
The status of solving the problem.
BasisStatus
Advanced usage: possible basis status values for a variable and the slack variable of a linear constr...
const Constraint * ct
int64 value
IntVar * var
Definition: expr_array.cc:1858
int64_t int64
A C++ wrapper that provides a simple and unified interface to several linear programming and mixed in...
const int WARNING
Definition: log_severity.h:31
const int FATAL
Definition: log_severity.h:32
ColIndex col
Definition: markowitz.cc:176
RowIndex row
Definition: markowitz.cc:175
CpSolverResponse Solve(const CpModelProto &model_proto)
Solves the given CpModelProto and returns an instance of CpSolverResponse.
The vehicle routing library lets one model and solve generic vehicle routing problems ranging from th...
int index
Definition: pack.cc:508
int64 coefficient
const bool maximize_
Definition: search.cc:2499