1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the HiGHS linear optimization suite */
4 /* */
5 /* Written and engineered 2008-2021 at the University of Edinburgh */
6 /* */
7 /* Available as open-source under the MIT License */
8 /* */
9 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
10 /**@file lp_data/HighsUtils.cpp
11 * @brief Class-independent utilities for HiGHS
12 * @author Julian Hall, Ivet Galabova, Qi Huangfu and Michael Feldmeier
13 */
14 #include "lp_data/HighsLpUtils.h"
15
16 #include <algorithm>
17 #include <cassert>
18
19 #include "HConfig.h"
20 #include "io/Filereader.h"
21 #include "io/HMPSIO.h"
22 #include "io/HighsIO.h"
23 #include "lp_data/HighsModelUtils.h"
24 #include "lp_data/HighsStatus.h"
25 #include "util/HighsSort.h"
26 #include "util/HighsTimer.h"
27
assessLp(HighsLp & lp,const HighsOptions & options)28 HighsStatus assessLp(HighsLp& lp, const HighsOptions& options) {
29 HighsStatus return_status = HighsStatus::OK;
30 HighsStatus call_status;
31 // Assess the LP dimensions and vector sizes, returning on error
32 call_status = assessLpDimensions(options, lp);
33 return_status =
34 interpretCallStatus(call_status, return_status, "assessLpDimensions");
35 if (return_status == HighsStatus::Error) return return_status;
36
37 // If the LP has no columns there is nothing left to test
38 // NB assessLpDimensions returns HighsStatus::Error if lp.numCol_ < 0
39 if (lp.numCol_ == 0) return HighsStatus::OK;
40
41 // From here, any LP has lp.numCol_ > 0 and lp.Astart_[lp.numCol_] exists (as
42 // the number of nonzeros)
43 assert(lp.numCol_ > 0);
44
45 // Assess the LP column costs
46 HighsIndexCollection index_collection;
47 index_collection.dimension_ = lp.numCol_;
48 index_collection.is_interval_ = true;
49 index_collection.from_ = 0;
50 index_collection.to_ = lp.numCol_ - 1;
51 call_status = assessCosts(options, 0, index_collection, lp.colCost_,
52 options.infinite_cost);
53 return_status =
54 interpretCallStatus(call_status, return_status, "assessCosts");
55 if (return_status == HighsStatus::Error) return return_status;
56 // Assess the LP column bounds
57 call_status = assessBounds(options, "Col", 0, index_collection, lp.colLower_,
58 lp.colUpper_, options.infinite_bound);
59 return_status =
60 interpretCallStatus(call_status, return_status, "assessBounds");
61 if (return_status == HighsStatus::Error) return return_status;
62 if (lp.numRow_) {
63 // Assess the LP row bounds
64 index_collection.dimension_ = lp.numRow_;
65 index_collection.is_interval_ = true;
66 index_collection.from_ = 0;
67 index_collection.to_ = lp.numRow_ - 1;
68 call_status =
69 assessBounds(options, "Row", 0, index_collection, lp.rowLower_,
70 lp.rowUpper_, options.infinite_bound);
71 return_status =
72 interpretCallStatus(call_status, return_status, "assessBounds");
73 if (return_status == HighsStatus::Error) return return_status;
74 // Assess the LP matrix
75 //
76 // The start of column 0 must be zero. It's possible to make
77 // everything consistent with the start being positive but, in
78 // particular, this means that the number of nonzeros is no longer
79 // lp.Astart_[lp.numCol_], a real banana skin to avoid just for
80 // the few wierdos who want lp.Astart_[0] to be positive. Hence
81 // it's not permitted!
82 if (lp.Astart_[0]) {
83 HighsLogMessage(options.logfile, HighsMessageType::ERROR,
84 "LP has nonzero value (%d) for the start of column 0\n",
85 lp.Astart_[0]);
86 return HighsStatus::Error;
87 }
88 call_status = assessMatrix(
89 options, lp.numRow_, lp.numCol_, lp.Astart_, lp.Aindex_, lp.Avalue_,
90 options.small_matrix_value, options.large_matrix_value);
91 return_status =
92 interpretCallStatus(call_status, return_status, "assessMatrix");
93 if (return_status == HighsStatus::Error) return return_status;
94 int lp_num_nz = lp.Astart_[lp.numCol_];
95 // If entries have been removed from the matrix, resize the index
96 // and value vectors to prevent bug in presolve
97 if ((int)lp.Aindex_.size() > lp_num_nz) lp.Aindex_.resize(lp_num_nz);
98 if ((int)lp.Avalue_.size() > lp_num_nz) lp.Avalue_.resize(lp_num_nz);
99 }
100 if (return_status == HighsStatus::Error)
101 return_status = HighsStatus::Error;
102 else
103 return_status = HighsStatus::OK;
104 #ifdef HiGHSDEV
105 HighsLogMessage(options.logfile, HighsMessageType::INFO,
106 "assess_lp returns HighsStatus = %s",
107 HighsStatusToString(return_status).c_str());
108 #endif
109 return return_status;
110 }
111
assessLpDimensions(const HighsOptions & options,const HighsLp & lp)112 HighsStatus assessLpDimensions(const HighsOptions& options, const HighsLp& lp) {
113 HighsStatus return_status = HighsStatus::OK;
114
115 // Use error_found to track whether an error has been found in multiple tests
116 bool error_found = false;
117
118 // Don't expect the matrix_start_size to be legal if there are no columns
119 bool check_matrix_start_size = lp.numCol_ > 0;
120
121 // Assess column-related dimensions
122 bool legal_num_col = lp.numCol_ >= 0;
123 if (!legal_num_col) {
124 HighsLogMessage(options.logfile, HighsMessageType::ERROR,
125 "LP has illegal number of cols = %d\n", lp.numCol_);
126 error_found = true;
127 } else {
128 // Check the size of the column vectors
129 int col_cost_size = lp.colCost_.size();
130 int col_lower_size = lp.colLower_.size();
131 int col_upper_size = lp.colUpper_.size();
132 int matrix_start_size = lp.Astart_.size();
133 bool legal_col_cost_size = col_cost_size >= lp.numCol_;
134 bool legal_col_lower_size = col_lower_size >= lp.numCol_;
135 bool legal_col_upper_size = col_lower_size >= lp.numCol_;
136
137 if (!legal_col_cost_size) {
138 HighsLogMessage(options.logfile, HighsMessageType::ERROR,
139 "LP has illegal colCost size = %d < %d\n", col_cost_size,
140 lp.numCol_);
141 error_found = true;
142 }
143 if (!legal_col_lower_size) {
144 HighsLogMessage(options.logfile, HighsMessageType::ERROR,
145 "LP has illegal colLower size = %d < %d\n",
146 col_lower_size, lp.numCol_);
147 error_found = true;
148 }
149 if (!legal_col_upper_size) {
150 HighsLogMessage(options.logfile, HighsMessageType::ERROR,
151 "LP has illegal colUpper size = %d < %d\n",
152 col_upper_size, lp.numCol_);
153 error_found = true;
154 }
155 if (check_matrix_start_size) {
156 bool legal_matrix_start_size = matrix_start_size >= lp.numCol_ + 1;
157 if (!legal_matrix_start_size) {
158 HighsLogMessage(options.logfile, HighsMessageType::ERROR,
159 "LP has illegal Astart size = %d < %d\n",
160 matrix_start_size, lp.numCol_ + 1);
161 error_found = true;
162 }
163 }
164 }
165
166 // Assess row-related dimensions
167 bool legal_num_row = lp.numRow_ >= 0;
168 if (!legal_num_row) {
169 HighsLogMessage(options.logfile, HighsMessageType::ERROR,
170 "LP has illegal number of rows = %d\n", lp.numRow_);
171 error_found = true;
172 } else {
173 int row_lower_size = lp.rowLower_.size();
174 int row_upper_size = lp.rowUpper_.size();
175 bool legal_row_lower_size = row_lower_size >= lp.numRow_;
176 bool legal_row_upper_size = row_lower_size >= lp.numRow_;
177 if (!legal_row_lower_size) {
178 HighsLogMessage(options.logfile, HighsMessageType::ERROR,
179 "LP has illegal rowLower size = %d < %d\n",
180 row_lower_size, lp.numRow_);
181 error_found = true;
182 }
183 if (!legal_row_upper_size) {
184 HighsLogMessage(options.logfile, HighsMessageType::ERROR,
185 "LP has illegal rowUpper size = %d < %d\n",
186 row_upper_size, lp.numRow_);
187 error_found = true;
188 }
189 }
190
191 // Assess matrix-related dimensions
192 if (check_matrix_start_size) {
193 int lp_num_nz = lp.Astart_[lp.numCol_];
194 bool legal_num_nz = lp_num_nz >= 0;
195 if (!legal_num_nz) {
196 HighsLogMessage(options.logfile, HighsMessageType::ERROR,
197 "LP has illegal number of nonzeros = %d\n", lp_num_nz);
198 error_found = true;
199 } else {
200 int matrix_index_size = lp.Aindex_.size();
201 int matrix_value_size = lp.Avalue_.size();
202 bool legal_matrix_index_size = matrix_index_size >= lp_num_nz;
203 bool legal_matrix_value_size = matrix_value_size >= lp_num_nz;
204 if (!legal_matrix_index_size) {
205 HighsLogMessage(options.logfile, HighsMessageType::ERROR,
206 "LP has illegal Aindex size = %d < %d\n",
207 matrix_index_size, lp_num_nz);
208 error_found = true;
209 }
210 if (!legal_matrix_value_size) {
211 HighsLogMessage(options.logfile, HighsMessageType::ERROR,
212 "LP has illegal Avalue size = %d < %d\n",
213 matrix_value_size, lp_num_nz);
214 error_found = true;
215 }
216 }
217 }
218 if (error_found)
219 return_status = HighsStatus::Error;
220 else
221 return_status = HighsStatus::OK;
222
223 return return_status;
224 }
225
assessCosts(const HighsOptions & options,const int ml_col_os,const HighsIndexCollection & index_collection,vector<double> & cost,const double infinite_cost)226 HighsStatus assessCosts(const HighsOptions& options, const int ml_col_os,
227 const HighsIndexCollection& index_collection,
228 vector<double>& cost, const double infinite_cost) {
229 HighsStatus return_status = HighsStatus::OK;
230 // Check parameters for technique and, if OK set the loop limits
231 if (!assessIndexCollection(options, index_collection))
232 return interpretCallStatus(HighsStatus::Error, return_status,
233 "assessIndexCollection");
234 int from_k;
235 int to_k;
236 if (!limitsForIndexCollection(options, index_collection, from_k, to_k))
237 return interpretCallStatus(HighsStatus::Error, return_status,
238 "limitsForIndexCollection");
239 if (from_k > to_k) return return_status;
240
241 return_status = HighsStatus::OK;
242 bool error_found = false;
243 // Work through the data to be assessed.
244 //
245 // Loop is k \in [from_k...to_k) covering the entries in the
246 // interval, set or mask to be considered.
247 //
248 // For an interval or mask, these values of k are the columns to be
249 // considered in a local sense, as well as the entries in the
250 // cost data to be assessed
251 //
252 // For a set, these values of k are the indices in the set, from
253 // which the columns to be considered in a local sense are
254 // drawn. The entries in the cost data to be assessed correspond
255 // to the values of k
256 //
257 // Adding the value of ml_col_os to local_col yields the value of
258 // ml_col, being the column in a global (whole-model) sense. This is
259 // necessary when assessing the costs of columns being added to a
260 // model, since they are specified using an interval
261 // [0...num_new_col) which must be offset by the current number of
262 // columns in the model.
263 //
264 int local_col;
265 int data_col;
266 int ml_col;
267 for (int k = from_k; k < to_k + 1; k++) {
268 if (index_collection.is_interval_ || index_collection.is_mask_) {
269 local_col = k;
270 data_col = k;
271 } else {
272 local_col = index_collection.set_[k];
273 data_col = k;
274 }
275 ml_col = ml_col_os + local_col;
276 if (index_collection.is_mask_ && !index_collection.mask_[local_col])
277 continue;
278 double abs_cost = fabs(cost[data_col]);
279 bool legal_cost = abs_cost < infinite_cost;
280 if (!legal_cost) {
281 HighsLogMessage(options.logfile, HighsMessageType::ERROR,
282 "Col %12d has |cost| of %12g >= %12g", ml_col, abs_cost,
283 infinite_cost);
284 error_found = !allow_infinite_costs;
285 }
286 }
287 if (error_found)
288 return_status = HighsStatus::Error;
289 else
290 return_status = HighsStatus::OK;
291
292 return return_status;
293 }
294
assessBounds(const HighsOptions & options,const char * type,const int ml_ix_os,const HighsIndexCollection & index_collection,vector<double> & lower,vector<double> & upper,const double infinite_bound)295 HighsStatus assessBounds(const HighsOptions& options, const char* type,
296 const int ml_ix_os,
297 const HighsIndexCollection& index_collection,
298 vector<double>& lower, vector<double>& upper,
299 const double infinite_bound) {
300 HighsStatus return_status = HighsStatus::OK;
301 // Check parameters for technique and, if OK set the loop limits
302 if (!assessIndexCollection(options, index_collection))
303 return interpretCallStatus(HighsStatus::Error, return_status,
304 "assessIndexCollection");
305 int from_k;
306 int to_k;
307 if (!limitsForIndexCollection(options, index_collection, from_k, to_k))
308 return interpretCallStatus(HighsStatus::Error, return_status,
309 "limitsForIndexCollection");
310 if (from_k > to_k) return HighsStatus::OK;
311
312 return_status = HighsStatus::OK;
313 bool error_found = false;
314 bool warning_found = false;
315 // Work through the data to be assessed.
316 //
317 // Loop is k \in [from_k...to_k) covering the entries in the
318 // interval, set or mask to be considered.
319 //
320 // For an interval or mask, these values of k are the row/column
321 // indices to be considered in a local sense, as well as the entries
322 // in the lower and upper bound data to be assessed
323 //
324 // For a set, these values of k are the indices in the set, from
325 // which the indices to be considered in a local sense are
326 // drawn. The entries in the lower and
327 // upper bound data to be assessed correspond to the values of
328 // k.
329 //
330 // Adding the value of ml_ix_os to local_ix yields the value of
331 // ml_ix, being the index in a global (whole-model) sense. This is
332 // necessary when assessing the bounds of rows/columns being added
333 // to a model, since they are specified using an interval
334 // [0...num_new_row/col) which must be offset by the current number
335 // of rows/columns (generically indices) in the model.
336 //
337 int num_infinite_lower_bound = 0;
338 int num_infinite_upper_bound = 0;
339 int local_ix;
340 int data_ix;
341 int ml_ix;
342 for (int k = from_k; k < to_k + 1; k++) {
343 if (index_collection.is_interval_ || index_collection.is_mask_) {
344 local_ix = k;
345 data_ix = k;
346 } else {
347 local_ix = index_collection.set_[k];
348 data_ix = k;
349 }
350 ml_ix = ml_ix_os + local_ix;
351 if (index_collection.is_mask_ && !index_collection.mask_[local_ix])
352 continue;
353
354 if (!highs_isInfinity(-lower[data_ix])) {
355 // Check whether a finite lower bound will be treated as -Infinity
356 bool infinite_lower_bound = lower[data_ix] <= -infinite_bound;
357 if (infinite_lower_bound) {
358 lower[data_ix] = -HIGHS_CONST_INF;
359 num_infinite_lower_bound++;
360 }
361 }
362 if (!highs_isInfinity(upper[data_ix])) {
363 // Check whether a finite upper bound will be treated as Infinity
364 bool infinite_upper_bound = upper[data_ix] >= infinite_bound;
365 if (infinite_upper_bound) {
366 upper[data_ix] = HIGHS_CONST_INF;
367 num_infinite_upper_bound++;
368 }
369 }
370 // Check that the lower bound does not exceed the upper bound
371 bool legalLowerUpperBound = lower[data_ix] <= upper[data_ix];
372 if (!legalLowerUpperBound) {
373 // Leave inconsistent bounds to be used to deduce infeasibility
374 HighsLogMessage(options.logfile, HighsMessageType::WARNING,
375 "%3s %12d has inconsistent bounds [%12g, %12g]", type,
376 ml_ix, lower[data_ix], upper[data_ix]);
377 warning_found = true;
378 }
379 // Check that the lower bound is not as much as +Infinity
380 bool legalLowerBound = lower[data_ix] < infinite_bound;
381 if (!legalLowerBound) {
382 HighsLogMessage(options.logfile, HighsMessageType::ERROR,
383 "%3s %12d has lower bound of %12g >= %12g", type, ml_ix,
384 lower[data_ix], infinite_bound);
385 error_found = true;
386 }
387 // Check that the upper bound is not as little as -Infinity
388 bool legalUpperBound = upper[data_ix] > -infinite_bound;
389 if (!legalUpperBound) {
390 HighsLogMessage(options.logfile, HighsMessageType::ERROR,
391 "%3s %12d has upper bound of %12g <= %12g", type, ml_ix,
392 upper[data_ix], -infinite_bound);
393 error_found = true;
394 }
395 }
396 if (num_infinite_lower_bound) {
397 HighsLogMessage(
398 options.logfile, HighsMessageType::INFO,
399 "%3ss:%12d lower bounds exceeding %12g are treated as -Infinity", type,
400 num_infinite_lower_bound, -infinite_bound);
401 }
402 if (num_infinite_upper_bound) {
403 HighsLogMessage(
404 options.logfile, HighsMessageType::INFO,
405 "%3ss:%12d upper bounds exceeding %12g are treated as +Infinity", type,
406 num_infinite_upper_bound, infinite_bound);
407 }
408
409 if (error_found)
410 return_status = HighsStatus::Error;
411 else if (warning_found)
412 return_status = HighsStatus::Warning;
413 else
414 return_status = HighsStatus::OK;
415
416 return return_status;
417 }
418
assessMatrix(const HighsOptions & options,const int vec_dim,const int num_vec,vector<int> & Astart,vector<int> & Aindex,vector<double> & Avalue,const double small_matrix_value,const double large_matrix_value)419 HighsStatus assessMatrix(const HighsOptions& options, const int vec_dim,
420 const int num_vec, vector<int>& Astart,
421 vector<int>& Aindex, vector<double>& Avalue,
422 const double small_matrix_value,
423 const double large_matrix_value) {
424 int num_nz = Astart[num_vec];
425 if (num_nz > 0 && vec_dim <= 0) return HighsStatus::Error;
426 if (num_nz <= 0) return HighsStatus::OK;
427
428 HighsStatus return_status = HighsStatus::OK;
429 bool error_found = false;
430 bool warning_found = false;
431
432 // Return a error if the first start is not zero
433 if (Astart[0]) {
434 HighsLogMessage(options.logfile, HighsMessageType::WARNING,
435 "Matrix starts do not begin with 0");
436 return HighsStatus::Error;
437 }
438 // Assess the starts
439 // Set up previous_start for a fictitious previous empty packed vector
440 int previous_start = Astart[0];
441 for (int ix = 0; ix < num_vec; ix++) {
442 int this_start = Astart[ix];
443 bool this_start_too_small = this_start < previous_start;
444 if (this_start_too_small) {
445 HighsLogMessage(options.logfile, HighsMessageType::ERROR,
446 "Matrix packed vector %d has illegal start of %d < %d = "
447 "previous start",
448 ix, this_start, previous_start);
449 return HighsStatus::Error;
450 }
451 bool this_start_too_big = this_start > num_nz;
452 if (this_start_too_big) {
453 HighsLogMessage(options.logfile, HighsMessageType::ERROR,
454 "Matrix packed vector %d has illegal start of %d > %d = "
455 "number of nonzeros",
456 ix, this_start, num_nz);
457 return HighsStatus::Error;
458 }
459 }
460
461 // Assess the indices and values
462 // Count the number of acceptable indices/values
463 int num_new_nz = 0;
464 int num_small_values = 0;
465 double max_small_value = 0;
466 double min_small_value = HIGHS_CONST_INF;
467 // Set up a zeroed vector to detect duplicate indices
468 vector<int> check_vector;
469 if (vec_dim > 0) check_vector.assign(vec_dim, 0);
470 for (int ix = 0; ix < num_vec; ix++) {
471 int from_el = Astart[ix];
472 int to_el = Astart[ix + 1];
473 // Account for any index-value pairs removed so far
474 Astart[ix] = num_new_nz;
475 for (int el = from_el; el < to_el; el++) {
476 // Check the index
477 int component = Aindex[el];
478 // Check that the index is non-negative
479 bool legal_component = component >= 0;
480 if (!legal_component) {
481 HighsLogMessage(
482 options.logfile, HighsMessageType::ERROR,
483 "Matrix packed vector %d, entry %d, is illegal index %d", ix, el,
484 component);
485 return HighsStatus::Error;
486 }
487 // Check that the index does not exceed the vector dimension
488 legal_component = component < vec_dim;
489 if (!legal_component) {
490 HighsLogMessage(options.logfile, HighsMessageType::ERROR,
491 "Matrix packed vector %d, entry %d, is illegal index "
492 "%12d >= %d = vector dimension",
493 ix, el, component, vec_dim);
494 return HighsStatus::Error;
495 }
496 // Check that the index has not already ocurred
497 legal_component = check_vector[component] == 0;
498 if (!legal_component) {
499 HighsLogMessage(
500 options.logfile, HighsMessageType::ERROR,
501 "Matrix packed vector %d, entry %d, is duplicate index %d", ix, el,
502 component);
503 return HighsStatus::Error;
504 }
505 // Indicate that the index has occurred
506 check_vector[component] = 1;
507 // Check the value
508 double abs_value = fabs(Avalue[el]);
509 /*
510 // Check that the value is not zero
511 bool zero_value = abs_value == 0;
512 if (zero_value) {
513 HighsLogMessage(options.logfile, HighsMessageType::ERROR,
514 "Matrix packed vector %d, entry %d, is zero", ix, el);
515 return HighsStatus::Error;
516 }
517 */
518 // Check that the value is not too large
519 bool large_value = abs_value > large_matrix_value;
520 if (large_value) {
521 HighsLogMessage(
522 options.logfile, HighsMessageType::ERROR,
523 "Matrix packed vector %d, entry %d, is large value |%g| >= %g", ix,
524 el, abs_value, large_matrix_value);
525 return HighsStatus::Error;
526 }
527 bool ok_value = abs_value > small_matrix_value;
528 if (!ok_value) {
529 if (max_small_value < abs_value) max_small_value = abs_value;
530 if (min_small_value > abs_value) min_small_value = abs_value;
531 num_small_values++;
532 }
533 if (ok_value) {
534 // Shift the index and value of the OK entry to the new
535 // position in the index and value vectors, and increment
536 // the new number of nonzeros
537 Aindex[num_new_nz] = Aindex[el];
538 Avalue[num_new_nz] = Avalue[el];
539 num_new_nz++;
540 } else {
541 // Zero the check_vector entry since the small value
542 // _hasn't_ occurred
543 check_vector[component] = 0;
544 }
545 }
546 // Zero check_vector
547 for (int el = Astart[ix]; el < num_new_nz; el++)
548 check_vector[Aindex[el]] = 0;
549 #ifdef HiGHSDEV
550 // NB This is very expensive so shouldn't be true
551 const bool check_check_vector = false;
552 if (check_check_vector) {
553 // Check zeroing of check vector
554 for (int component = 0; component < vec_dim; component++) {
555 if (check_vector[component]) error_found = true;
556 }
557 if (error_found)
558 HighsLogMessage(options.logfile, HighsMessageType::ERROR,
559 "assessMatrix: check_vector not zeroed");
560 }
561 #endif
562 }
563 if (num_small_values) {
564 HighsLogMessage(options.logfile, HighsMessageType::WARNING,
565 "Matrix packed vector contains %d |values| in [%g, %g] "
566 "less than %g: ignored",
567 num_small_values, min_small_value, max_small_value,
568 small_matrix_value);
569 warning_found = true;
570 }
571 Astart[num_vec] = num_new_nz;
572 if (error_found)
573 return_status = HighsStatus::Error;
574 else if (warning_found)
575 return_status = HighsStatus::Warning;
576 else
577 return_status = HighsStatus::OK;
578
579 return return_status;
580 }
581
cleanBounds(const HighsOptions & options,HighsLp & lp)582 HighsStatus cleanBounds(const HighsOptions& options, HighsLp& lp) {
583 double max_residual = 0;
584 int num_change = 0;
585 for (int iCol = 0; iCol < lp.numCol_; iCol++) {
586 double residual = lp.colLower_[iCol] - lp.colUpper_[iCol];
587 if (residual > options.primal_feasibility_tolerance) {
588 HighsLogMessage(options.logfile, HighsMessageType::ERROR,
589 "Column %d has inconsistent bounds [%g, %g] (residual = "
590 "%g) after presolve ",
591 iCol, lp.colLower_[iCol], lp.colUpper_[iCol], residual);
592 return HighsStatus::Error;
593 } else if (residual > 0) {
594 num_change++;
595 max_residual = std::max(residual, max_residual);
596 double mid = 0.5 * (lp.colLower_[iCol] + lp.colUpper_[iCol]);
597 lp.colLower_[iCol] = mid;
598 lp.colUpper_[iCol] = mid;
599 }
600 }
601 for (int iRow = 0; iRow < lp.numRow_; iRow++) {
602 double residual = lp.rowLower_[iRow] - lp.rowUpper_[iRow];
603 if (residual > options.primal_feasibility_tolerance) {
604 HighsLogMessage(options.logfile, HighsMessageType::ERROR,
605 "Row %d has inconsistent bounds [%g, %g] (residual = %g) "
606 "after presolve ",
607 iRow, lp.rowLower_[iRow], lp.rowUpper_[iRow], residual);
608 return HighsStatus::Error;
609 } else if (residual > 0) {
610 num_change++;
611 max_residual = std::max(residual, max_residual);
612 double mid = 0.5 * (lp.rowLower_[iRow] + lp.rowUpper_[iRow]);
613 lp.rowLower_[iRow] = mid;
614 lp.rowUpper_[iRow] = mid;
615 }
616 }
617 if (num_change) {
618 HighsLogMessage(options.logfile, HighsMessageType::WARNING,
619 "Resolved %d inconsistent bounds (maximum residual = "
620 "%9.4g) after presolve ",
621 num_change, max_residual);
622 return HighsStatus::Warning;
623 }
624 return HighsStatus::OK;
625 }
626
applyScalingToLp(const HighsOptions & options,HighsLp & lp,const HighsScale & scale)627 HighsStatus applyScalingToLp(const HighsOptions& options, HighsLp& lp,
628 const HighsScale& scale) {
629 if (!scale.is_scaled_) return HighsStatus::OK;
630 if ((int)scale.col_.size() < lp.numCol_) return HighsStatus::Error;
631 if ((int)scale.row_.size() < lp.numRow_) return HighsStatus::Error;
632 bool scale_error = false;
633 // Set up column and row index collections for scaling
634 HighsIndexCollection all_cols;
635 all_cols.is_interval_ = true;
636 all_cols.dimension_ = lp.numCol_;
637 all_cols.from_ = 0;
638 all_cols.to_ = lp.numCol_ - 1;
639 HighsIndexCollection all_rows;
640 all_rows.is_interval_ = true;
641 all_rows.dimension_ = lp.numRow_;
642 all_rows.from_ = 0;
643 all_rows.to_ = lp.numRow_ - 1;
644
645 scale_error = applyScalingToLpColCost(options, lp, scale.col_, all_cols) !=
646 HighsStatus::OK ||
647 scale_error;
648 scale_error = applyScalingToLpColBounds(options, lp, scale.col_, all_cols) !=
649 HighsStatus::OK ||
650 scale_error;
651 scale_error = applyScalingToLpRowBounds(options, lp, scale.row_, all_rows) !=
652 HighsStatus::OK ||
653 scale_error;
654 scale_error = applyScalingToLpMatrix(options, lp, &scale.col_[0],
655 &scale.row_[0], 0, lp.numCol_ - 1, 0,
656 lp.numRow_ - 1) != HighsStatus::OK ||
657 scale_error;
658 if (scale_error) return HighsStatus::Error;
659 return HighsStatus::OK;
660 }
661
applyScalingToLpColCost(const HighsOptions & options,HighsLp & lp,const vector<double> & colScale,const HighsIndexCollection & index_collection)662 HighsStatus applyScalingToLpColCost(
663 const HighsOptions& options, HighsLp& lp, const vector<double>& colScale,
664 const HighsIndexCollection& index_collection) {
665 HighsStatus return_status = HighsStatus::OK;
666 // Check parameters for technique and, if OK set the loop limits
667 if (!assessIndexCollection(options, index_collection))
668 return interpretCallStatus(HighsStatus::Error, return_status,
669 "assessIndexCollection");
670
671 int from_k;
672 int to_k;
673 if (!limitsForIndexCollection(options, index_collection, from_k, to_k))
674 return interpretCallStatus(HighsStatus::Error, return_status,
675 "limitsForIndexCollection");
676 if (from_k > to_k) return HighsStatus::OK;
677
678 const bool& interval = index_collection.is_interval_;
679 const bool& mask = index_collection.is_mask_;
680 const int* col_set = index_collection.set_;
681 const int* col_mask = index_collection.mask_;
682
683 int local_col;
684 int ml_col;
685 const int ml_col_os = 0;
686 for (int k = from_k; k < to_k + 1; k++) {
687 if (interval || mask) {
688 local_col = k;
689 } else {
690 local_col = col_set[k];
691 }
692 ml_col = ml_col_os + local_col;
693 if (mask && !col_mask[local_col]) continue;
694 lp.colCost_[ml_col] *= colScale[ml_col];
695 }
696
697 return HighsStatus::OK;
698 }
699
applyScalingToLpColBounds(const HighsOptions & options,HighsLp & lp,const vector<double> & colScale,const HighsIndexCollection & index_collection)700 HighsStatus applyScalingToLpColBounds(
701 const HighsOptions& options, HighsLp& lp, const vector<double>& colScale,
702 const HighsIndexCollection& index_collection) {
703 HighsStatus return_status = HighsStatus::OK;
704 // Check parameters for technique and, if OK set the loop limits
705 if (!assessIndexCollection(options, index_collection))
706 return interpretCallStatus(HighsStatus::Error, return_status,
707 "assessIndexCollection");
708
709 int from_k;
710 int to_k;
711 if (!limitsForIndexCollection(options, index_collection, from_k, to_k))
712 return interpretCallStatus(HighsStatus::Error, return_status,
713 "limitsForIndexCollection");
714 if (from_k > to_k) return HighsStatus::OK;
715
716 const bool& interval = index_collection.is_interval_;
717 const bool& mask = index_collection.is_mask_;
718 const int* col_set = index_collection.set_;
719 const int* col_mask = index_collection.mask_;
720
721 int local_col;
722 int ml_col;
723 const int ml_col_os = 0;
724 for (int k = from_k; k < to_k + 1; k++) {
725 if (interval || mask) {
726 local_col = k;
727 } else {
728 local_col = col_set[k];
729 }
730 ml_col = ml_col_os + local_col;
731 if (mask && !col_mask[local_col]) continue;
732 if (!highs_isInfinity(-lp.colLower_[ml_col]))
733 lp.colLower_[ml_col] /= colScale[ml_col];
734 if (!highs_isInfinity(lp.colUpper_[ml_col]))
735 lp.colUpper_[ml_col] /= colScale[ml_col];
736 }
737
738 return HighsStatus::OK;
739 }
740
applyScalingToLpRowBounds(const HighsOptions & options,HighsLp & lp,const vector<double> & rowScale,const HighsIndexCollection & index_collection)741 HighsStatus applyScalingToLpRowBounds(
742 const HighsOptions& options, HighsLp& lp, const vector<double>& rowScale,
743 const HighsIndexCollection& index_collection) {
744 HighsStatus return_status = HighsStatus::OK;
745 // Check parameters for technique and, if OK set the loop limits
746 if (!assessIndexCollection(options, index_collection))
747 return interpretCallStatus(HighsStatus::Error, return_status,
748 "assessIndexCollection");
749
750 int from_k;
751 int to_k;
752 if (!limitsForIndexCollection(options, index_collection, from_k, to_k))
753 return interpretCallStatus(HighsStatus::Error, return_status,
754 "limitsForIndexCollection");
755 if (from_k > to_k) return HighsStatus::OK;
756
757 const bool& interval = index_collection.is_interval_;
758 const bool& mask = index_collection.is_mask_;
759 const int* row_set = index_collection.set_;
760 const int* row_mask = index_collection.mask_;
761
762 int local_row;
763 int ml_row;
764 const int ml_row_os = 0;
765 for (int k = from_k; k < to_k + 1; k++) {
766 if (interval || mask) {
767 local_row = k;
768 } else {
769 local_row = row_set[k];
770 }
771 ml_row = ml_row_os + local_row;
772 if (mask && !row_mask[local_row]) continue;
773 if (!highs_isInfinity(-lp.rowLower_[ml_row]))
774 lp.rowLower_[ml_row] *= rowScale[ml_row];
775 if (!highs_isInfinity(lp.rowUpper_[ml_row]))
776 lp.rowUpper_[ml_row] *= rowScale[ml_row];
777 }
778
779 return HighsStatus::OK;
780 }
781
applyScalingToLpMatrix(const HighsOptions & options,HighsLp & lp,const double * colScale,const double * rowScale,const int from_col,const int to_col,const int from_row,const int to_row)782 HighsStatus applyScalingToLpMatrix(const HighsOptions& options, HighsLp& lp,
783 const double* colScale,
784 const double* rowScale, const int from_col,
785 const int to_col, const int from_row,
786 const int to_row) {
787 if (from_col < 0) return HighsStatus::Error;
788 if (to_col >= lp.numCol_) return HighsStatus::Error;
789 if (from_row < 0) return HighsStatus::Error;
790 if (to_row >= lp.numRow_) return HighsStatus::Error;
791 if (colScale != NULL) {
792 if (rowScale != NULL) {
793 for (int iCol = from_col; iCol <= to_col; iCol++) {
794 for (int iEl = lp.Astart_[iCol]; iEl < lp.Astart_[iCol + 1]; iEl++) {
795 int iRow = lp.Aindex_[iEl];
796 if (iRow < from_row || iRow > to_row) continue;
797 lp.Avalue_[iEl] *= (colScale[iCol] * rowScale[iRow]);
798 }
799 }
800 } else {
801 // No row scaling
802 for (int iCol = from_col; iCol <= to_col; iCol++) {
803 for (int iEl = lp.Astart_[iCol]; iEl < lp.Astart_[iCol + 1]; iEl++) {
804 int iRow = lp.Aindex_[iEl];
805 if (iRow < from_row || iRow > to_row) continue;
806 lp.Avalue_[iEl] *= colScale[iCol];
807 }
808 }
809 }
810 } else {
811 // No column scaling
812 if (rowScale != NULL) {
813 for (int iCol = from_col; iCol <= to_col; iCol++) {
814 for (int iEl = lp.Astart_[iCol]; iEl < lp.Astart_[iCol + 1]; iEl++) {
815 int iRow = lp.Aindex_[iEl];
816 if (iRow < from_row || iRow > to_row) continue;
817 lp.Avalue_[iEl] *= rowScale[iRow];
818 }
819 }
820 }
821 }
822 return HighsStatus::OK;
823 }
824
applyRowScalingToMatrix(const vector<double> & rowScale,const int numCol,const vector<int> & Astart,const vector<int> & Aindex,vector<double> & Avalue)825 void applyRowScalingToMatrix(const vector<double>& rowScale, const int numCol,
826 const vector<int>& Astart,
827 const vector<int>& Aindex,
828 vector<double>& Avalue) {
829 for (int iCol = 0; iCol < numCol; iCol++) {
830 for (int el = Astart[iCol]; el < Astart[iCol + 1]; el++) {
831 Avalue[el] *= rowScale[Aindex[el]];
832 }
833 }
834 }
835
colScaleMatrix(const int max_scale_factor_exponent,double * colScale,const int numCol,const vector<int> & Astart,const vector<int> & Aindex,vector<double> & Avalue)836 void colScaleMatrix(const int max_scale_factor_exponent, double* colScale,
837 const int numCol, const vector<int>& Astart,
838 const vector<int>& Aindex, vector<double>& Avalue) {
839 const double log2 = log(2.0);
840 const double max_allow_scale = pow(2.0, max_scale_factor_exponent);
841 const double min_allow_scale = 1 / max_allow_scale;
842
843 const double min_allow_col_scale = min_allow_scale;
844 const double max_allow_col_scale = max_allow_scale;
845
846 for (int iCol = 0; iCol < numCol; iCol++) {
847 double col_max_value = 0;
848 for (int k = Astart[iCol]; k < Astart[iCol + 1]; k++)
849 col_max_value = max(fabs(Avalue[k]), col_max_value);
850 if (col_max_value) {
851 double col_scale_value = 1 / col_max_value;
852 // Convert the col scale factor to the nearest power of two, and
853 // ensure that it is not excessively large or small
854 col_scale_value = pow(2.0, floor(log(col_scale_value) / log2 + 0.5));
855 col_scale_value =
856 min(max(min_allow_col_scale, col_scale_value), max_allow_col_scale);
857 colScale[iCol] = col_scale_value;
858 // Scale the column
859 for (int k = Astart[iCol]; k < Astart[iCol + 1]; k++)
860 Avalue[k] *= colScale[iCol];
861 } else {
862 // Empty column
863 colScale[iCol] = 1;
864 }
865 }
866 }
867
applyScalingToLpCol(const HighsOptions & options,HighsLp & lp,const int col,const double colScale)868 HighsStatus applyScalingToLpCol(const HighsOptions& options, HighsLp& lp,
869 const int col, const double colScale) {
870 if (col < 0) return HighsStatus::Error;
871 if (col >= lp.numCol_) return HighsStatus::Error;
872 if (!colScale) return HighsStatus::Error;
873
874 for (int el = lp.Astart_[col]; el < lp.Astart_[col + 1]; el++)
875 lp.Avalue_[el] *= colScale;
876 lp.colCost_[col] *= colScale;
877 if (colScale > 0) {
878 lp.colLower_[col] /= colScale;
879 lp.colUpper_[col] /= colScale;
880 } else {
881 const double new_upper = lp.colLower_[col] / colScale;
882 lp.colLower_[col] = lp.colUpper_[col] / colScale;
883 lp.colUpper_[col] = new_upper;
884 }
885 return HighsStatus::OK;
886 }
887
applyScalingToLpRow(const HighsOptions & options,HighsLp & lp,const int row,const double rowScale)888 HighsStatus applyScalingToLpRow(const HighsOptions& options, HighsLp& lp,
889 const int row, const double rowScale) {
890 if (row < 0) return HighsStatus::Error;
891 if (row >= lp.numRow_) return HighsStatus::Error;
892 if (!rowScale) return HighsStatus::Error;
893
894 for (int col = 0; col < lp.numCol_; col++) {
895 for (int el = lp.Astart_[col]; el < lp.Astart_[col + 1]; el++) {
896 if (lp.Aindex_[el] == row) lp.Avalue_[el] *= rowScale;
897 }
898 }
899 if (rowScale > 0) {
900 lp.rowLower_[row] /= rowScale;
901 lp.rowUpper_[row] /= rowScale;
902 } else {
903 const double new_upper = lp.rowLower_[row] / rowScale;
904 lp.rowLower_[row] = lp.rowUpper_[row] / rowScale;
905 lp.rowUpper_[row] = new_upper;
906 }
907 return HighsStatus::OK;
908 }
909
appendColsToLpVectors(HighsLp & lp,const int num_new_col,const vector<double> & colCost,const vector<double> & colLower,const vector<double> & colUpper)910 HighsStatus appendColsToLpVectors(HighsLp& lp, const int num_new_col,
911 const vector<double>& colCost,
912 const vector<double>& colLower,
913 const vector<double>& colUpper) {
914 if (num_new_col < 0) return HighsStatus::Error;
915 if (num_new_col == 0) return HighsStatus::OK;
916 int new_num_col = lp.numCol_ + num_new_col;
917 lp.colCost_.resize(new_num_col);
918 lp.colLower_.resize(new_num_col);
919 lp.colUpper_.resize(new_num_col);
920 bool have_names = lp.col_names_.size();
921 if (have_names) lp.col_names_.resize(new_num_col);
922 for (int new_col = 0; new_col < num_new_col; new_col++) {
923 int iCol = lp.numCol_ + new_col;
924 lp.colCost_[iCol] = colCost[new_col];
925 lp.colLower_[iCol] = colLower[new_col];
926 lp.colUpper_[iCol] = colUpper[new_col];
927 // Cannot guarantee to create unique names, so name is blank
928 if (have_names) lp.col_names_[iCol] = "";
929 }
930 return HighsStatus::OK;
931 }
932
appendRowsToLpVectors(HighsLp & lp,const int num_new_row,const vector<double> & rowLower,const vector<double> & rowUpper)933 HighsStatus appendRowsToLpVectors(HighsLp& lp, const int num_new_row,
934 const vector<double>& rowLower,
935 const vector<double>& rowUpper) {
936 if (num_new_row < 0) return HighsStatus::Error;
937 if (num_new_row == 0) return HighsStatus::OK;
938 int new_num_row = lp.numRow_ + num_new_row;
939 lp.rowLower_.resize(new_num_row);
940 lp.rowUpper_.resize(new_num_row);
941 bool have_names = lp.row_names_.size();
942 if (have_names) lp.row_names_.resize(new_num_row);
943
944 for (int new_row = 0; new_row < num_new_row; new_row++) {
945 int iRow = lp.numRow_ + new_row;
946 lp.rowLower_[iRow] = rowLower[new_row];
947 lp.rowUpper_[iRow] = rowUpper[new_row];
948 // Cannot guarantee to create unique names, so name is blank
949 if (have_names) lp.row_names_[iRow] = "";
950 }
951 return HighsStatus::OK;
952 }
953
appendColsToLpMatrix(HighsLp & lp,const int num_new_col,const int num_new_nz,const int * XAstart,const int * XAindex,const double * XAvalue)954 HighsStatus appendColsToLpMatrix(HighsLp& lp, const int num_new_col,
955 const int num_new_nz, const int* XAstart,
956 const int* XAindex, const double* XAvalue) {
957 if (num_new_col < 0) return HighsStatus::Error;
958 if (num_new_col == 0) return HighsStatus::OK;
959 // Check that nonzeros aren't being appended to a matrix with no rows
960 if (num_new_nz > 0 && lp.numRow_ <= 0) return HighsStatus::Error;
961 // Determine the new number of columns in the matrix and resize the
962 // starts accordingly.
963 int new_num_col = lp.numCol_ + num_new_col;
964 lp.Astart_.resize(new_num_col + 1);
965 // If adding columns to an empty LP then introduce the start for the
966 // fictitious column 0
967 if (lp.numCol_ == 0) lp.Astart_[0] = 0;
968
969 // Determine the current number of nonzeros and the new number of nonzeros
970 int current_num_nz = lp.Astart_[lp.numCol_];
971 int new_num_nz = current_num_nz + num_new_nz;
972
973 // Append the starts of the new columns
974 if (num_new_nz) {
975 // Nontrivial number of nonzeros being added, so use XAstart
976 assert(XAstart != NULL);
977 for (int col = 0; col < num_new_col; col++)
978 lp.Astart_[lp.numCol_ + col] = current_num_nz + XAstart[col];
979 } else {
980 // No nonzeros being added, so XAstart may be null, but entries of
981 // zero are implied.
982 for (int col = 0; col < num_new_col; col++)
983 lp.Astart_[lp.numCol_ + col] = current_num_nz;
984 }
985 lp.Astart_[lp.numCol_ + num_new_col] = new_num_nz;
986
987 // If no nonzeros are being added then there's nothing else to do
988 if (num_new_nz <= 0) return HighsStatus::OK;
989
990 // Adding a non-trivial matrix: resize the column-wise matrix arrays
991 // accordingly
992 lp.Aindex_.resize(new_num_nz);
993 lp.Avalue_.resize(new_num_nz);
994 // Copy in the new indices and values
995 for (int el = 0; el < num_new_nz; el++) {
996 lp.Aindex_[current_num_nz + el] = XAindex[el];
997 lp.Avalue_[current_num_nz + el] = XAvalue[el];
998 }
999 return HighsStatus::OK;
1000 }
1001
appendRowsToLpMatrix(HighsLp & lp,const int num_new_row,const int num_new_nz,const int * XARstart,const int * XARindex,const double * XARvalue)1002 HighsStatus appendRowsToLpMatrix(HighsLp& lp, const int num_new_row,
1003 const int num_new_nz, const int* XARstart,
1004 const int* XARindex, const double* XARvalue) {
1005 if (num_new_row < 0) return HighsStatus::Error;
1006 if (num_new_row == 0) return HighsStatus::OK;
1007 // Check that nonzeros aren't being appended to a matrix with no columns
1008 if (num_new_nz > 0 && lp.numCol_ <= 0) return HighsStatus::Error;
1009 if (num_new_nz == 0) return HighsStatus::OK;
1010 int current_num_nz = lp.Astart_[lp.numCol_];
1011 vector<int> Alength;
1012 Alength.assign(lp.numCol_, 0);
1013 for (int el = 0; el < num_new_nz; el++) Alength[XARindex[el]]++;
1014 // Determine the new number of nonzeros and resize the column-wise matrix
1015 // arrays
1016 int new_num_nz = current_num_nz + num_new_nz;
1017 lp.Aindex_.resize(new_num_nz);
1018 lp.Avalue_.resize(new_num_nz);
1019
1020 // Append the new rows
1021 // Shift the existing columns to make space for the new entries
1022 int new_el = new_num_nz;
1023 for (int col = lp.numCol_ - 1; col >= 0; col--) {
1024 int start_col_plus_1 = new_el;
1025 new_el -= Alength[col];
1026 for (int el = lp.Astart_[col + 1] - 1; el >= lp.Astart_[col]; el--) {
1027 new_el--;
1028 lp.Aindex_[new_el] = lp.Aindex_[el];
1029 lp.Avalue_[new_el] = lp.Avalue_[el];
1030 }
1031 lp.Astart_[col + 1] = start_col_plus_1;
1032 }
1033 assert(new_el == 0);
1034
1035 // Insert the new entries
1036 for (int row = 0; row < num_new_row; row++) {
1037 int first_el = XARstart[row];
1038 int last_el = (row < num_new_row - 1 ? XARstart[row + 1] : num_new_nz);
1039 for (int el = first_el; el < last_el; el++) {
1040 int col = XARindex[el];
1041 new_el = lp.Astart_[col + 1] - Alength[col];
1042 Alength[col]--;
1043 lp.Aindex_[new_el] = lp.numRow_ + row;
1044 lp.Avalue_[new_el] = XARvalue[el];
1045 }
1046 }
1047 return HighsStatus::OK;
1048 }
1049
deleteLpCols(const HighsOptions & options,HighsLp & lp,const HighsIndexCollection & index_collection)1050 HighsStatus deleteLpCols(const HighsOptions& options, HighsLp& lp,
1051 const HighsIndexCollection& index_collection) {
1052 int new_num_col;
1053 HighsStatus call_status;
1054 call_status =
1055 deleteColsFromLpVectors(options, lp, new_num_col, index_collection);
1056 if (call_status != HighsStatus::OK) return call_status;
1057 call_status = deleteColsFromLpMatrix(options, lp, index_collection);
1058 if (call_status != HighsStatus::OK) return call_status;
1059 lp.numCol_ = new_num_col;
1060 return HighsStatus::OK;
1061 }
1062
deleteColsFromLpVectors(const HighsOptions & options,HighsLp & lp,int & new_num_col,const HighsIndexCollection & index_collection)1063 HighsStatus deleteColsFromLpVectors(
1064 const HighsOptions& options, HighsLp& lp, int& new_num_col,
1065 const HighsIndexCollection& index_collection) {
1066 HighsStatus return_status = HighsStatus::OK;
1067 if (!assessIndexCollection(options, index_collection))
1068 return interpretCallStatus(HighsStatus::Error, return_status,
1069 "assessIndexCollection");
1070 int from_k;
1071 int to_k;
1072 if (!limitsForIndexCollection(options, index_collection, from_k, to_k))
1073 return interpretCallStatus(HighsStatus::Error, return_status,
1074 "limitsForIndexCollection");
1075 if (index_collection.is_set_) {
1076 // For deletion by set it must be increasing
1077 if (!increasingSetOk(index_collection.set_,
1078 index_collection.set_num_entries_, 0, lp.numCol_ - 1,
1079 true))
1080 return HighsStatus::Error;
1081 }
1082 // Initialise new_num_col in case none is removed due to from_k > to_k
1083 new_num_col = lp.numCol_;
1084 if (from_k > to_k) return HighsStatus::OK;
1085
1086 int delete_from_col;
1087 int delete_to_col;
1088 int keep_from_col;
1089 int keep_to_col = -1;
1090 int current_set_entry = 0;
1091
1092 int col_dim = lp.numCol_;
1093 new_num_col = 0;
1094 bool have_names = lp.col_names_.size();
1095 for (int k = from_k; k <= to_k; k++) {
1096 updateIndexCollectionOutInIndex(index_collection, delete_from_col,
1097 delete_to_col, keep_from_col, keep_to_col,
1098 current_set_entry);
1099 // Account for the initial columns being kept
1100 if (k == from_k) new_num_col = delete_from_col;
1101 if (delete_to_col >= col_dim - 1) break;
1102 assert(delete_to_col < col_dim);
1103 for (int col = keep_from_col; col <= keep_to_col; col++) {
1104 lp.colCost_[new_num_col] = lp.colCost_[col];
1105 lp.colLower_[new_num_col] = lp.colLower_[col];
1106 lp.colUpper_[new_num_col] = lp.colUpper_[col];
1107 if (have_names) lp.col_names_[new_num_col] = lp.col_names_[col];
1108 new_num_col++;
1109 }
1110 if (keep_to_col >= col_dim - 1) break;
1111 }
1112 lp.colCost_.resize(new_num_col);
1113 lp.colLower_.resize(new_num_col);
1114 lp.colUpper_.resize(new_num_col);
1115 if (have_names) lp.col_names_.resize(new_num_col);
1116 return HighsStatus::OK;
1117 }
1118
deleteColsFromLpMatrix(const HighsOptions & options,HighsLp & lp,const HighsIndexCollection & index_collection)1119 HighsStatus deleteColsFromLpMatrix(
1120 const HighsOptions& options, HighsLp& lp,
1121 const HighsIndexCollection& index_collection) {
1122 HighsStatus return_status = HighsStatus::OK;
1123 if (!assessIndexCollection(options, index_collection))
1124 return interpretCallStatus(HighsStatus::Error, return_status,
1125 "assessIndexCollection");
1126 int from_k;
1127 int to_k;
1128 if (!limitsForIndexCollection(options, index_collection, from_k, to_k))
1129 return interpretCallStatus(HighsStatus::Error, return_status,
1130 "limitsForIndexCollection");
1131 if (index_collection.is_set_) {
1132 // For deletion by set it must be increasing
1133 if (!increasingSetOk(index_collection.set_,
1134 index_collection.set_num_entries_, 0, lp.numCol_ - 1,
1135 true))
1136 return HighsStatus::Error;
1137 }
1138 if (from_k > to_k) return HighsStatus::OK;
1139
1140 int delete_from_col;
1141 int delete_to_col;
1142 int keep_from_col;
1143 int keep_to_col = -1;
1144 int current_set_entry = 0;
1145
1146 int col_dim = lp.numCol_;
1147 int new_num_col = 0;
1148 int new_num_nz = 0;
1149 for (int k = from_k; k <= to_k; k++) {
1150 updateIndexCollectionOutInIndex(index_collection, delete_from_col,
1151 delete_to_col, keep_from_col, keep_to_col,
1152 current_set_entry);
1153 if (k == from_k) {
1154 // Account for the initial columns being kept
1155 new_num_col = delete_from_col;
1156 new_num_nz = lp.Astart_[delete_from_col];
1157 }
1158 // Ensure that the starts of the deleted columns are zeroed to
1159 // avoid redundant start information for columns whose indices
1160 // are't used after the deletion takes place. In particular, if
1161 // all columns are deleted then something must be done to ensure
1162 // that the matrix isn't magially recreated by increasing the
1163 // number of columns from zero when there are no rows in the LP.
1164 for (int col = delete_from_col; col <= delete_to_col; col++)
1165 lp.Astart_[col] = 0;
1166 // Shift the starts - both in place and value - to account for the
1167 // columns and nonzeros removed
1168 const int keep_from_el = lp.Astart_[keep_from_col];
1169 for (int col = keep_from_col; col <= keep_to_col; col++) {
1170 lp.Astart_[new_num_col] = new_num_nz + lp.Astart_[col] - keep_from_el;
1171 new_num_col++;
1172 }
1173 for (int el = keep_from_el; el < lp.Astart_[keep_to_col + 1]; el++) {
1174 lp.Aindex_[new_num_nz] = lp.Aindex_[el];
1175 lp.Avalue_[new_num_nz] = lp.Avalue_[el];
1176 new_num_nz++;
1177 }
1178 if (keep_to_col >= col_dim - 1) break;
1179 }
1180 // Ensure that the start of the spurious last column is zeroed so
1181 // that it doesn't give a positive number of matrix entries if the
1182 // number of columns in the LP is increased when there are no rows
1183 // in the LP.
1184 lp.Astart_[lp.numCol_] = 0;
1185 lp.Astart_[new_num_col] = new_num_nz;
1186 lp.Astart_.resize(new_num_col + 1);
1187 lp.Aindex_.resize(new_num_nz);
1188 lp.Avalue_.resize(new_num_nz);
1189 return HighsStatus::OK;
1190 }
1191
deleteLpRows(const HighsOptions & options,HighsLp & lp,const HighsIndexCollection & index_collection)1192 HighsStatus deleteLpRows(const HighsOptions& options, HighsLp& lp,
1193 const HighsIndexCollection& index_collection) {
1194 HighsStatus return_status = HighsStatus::OK;
1195 HighsStatus call_status;
1196 int new_num_row;
1197 call_status =
1198 deleteRowsFromLpVectors(options, lp, new_num_row, index_collection);
1199 return_status = interpretCallStatus(call_status, return_status,
1200 "deleteRowsFromLpVectors");
1201 if (return_status == HighsStatus::Error) return return_status;
1202 call_status = deleteRowsFromLpMatrix(options, lp, index_collection);
1203 return_status =
1204 interpretCallStatus(call_status, return_status, "deleteRowsFromLpMatrix");
1205 if (return_status == HighsStatus::Error) return return_status;
1206 lp.numRow_ = new_num_row;
1207 return HighsStatus::OK;
1208 }
1209
deleteRowsFromLpVectors(const HighsOptions & options,HighsLp & lp,int & new_num_row,const HighsIndexCollection & index_collection)1210 HighsStatus deleteRowsFromLpVectors(
1211 const HighsOptions& options, HighsLp& lp, int& new_num_row,
1212 const HighsIndexCollection& index_collection) {
1213 HighsStatus return_status = HighsStatus::OK;
1214 if (!assessIndexCollection(options, index_collection))
1215 return interpretCallStatus(HighsStatus::Error, return_status,
1216 "assessIndexCollection");
1217 int from_k;
1218 int to_k;
1219 if (!limitsForIndexCollection(options, index_collection, from_k, to_k))
1220 return interpretCallStatus(HighsStatus::Error, return_status,
1221 "limitsForIndexCollection");
1222 if (index_collection.is_set_) {
1223 // For deletion by set it must be increasing
1224 if (!increasingSetOk(index_collection.set_,
1225 index_collection.set_num_entries_, 0, lp.numRow_ - 1,
1226 true))
1227 return HighsStatus::Error;
1228 }
1229 // Initialise new_num_row in case none is removed due to from_k > to_k
1230 new_num_row = lp.numRow_;
1231 if (from_k > to_k) return HighsStatus::OK;
1232
1233 int delete_from_row;
1234 int delete_to_row;
1235 int keep_from_row;
1236 int keep_to_row = -1;
1237 int current_set_entry = 0;
1238
1239 int row_dim = lp.numRow_;
1240 new_num_row = 0;
1241 bool have_names = (int)lp.row_names_.size() > 0;
1242 for (int k = from_k; k <= to_k; k++) {
1243 updateIndexCollectionOutInIndex(index_collection, delete_from_row,
1244 delete_to_row, keep_from_row, keep_to_row,
1245 current_set_entry);
1246 if (k == from_k) {
1247 // Account for the initial rows being kept
1248 new_num_row = delete_from_row;
1249 }
1250 if (delete_to_row >= row_dim - 1) break;
1251 assert(delete_to_row < row_dim);
1252 for (int row = keep_from_row; row <= keep_to_row; row++) {
1253 lp.rowLower_[new_num_row] = lp.rowLower_[row];
1254 lp.rowUpper_[new_num_row] = lp.rowUpper_[row];
1255 if (have_names) lp.row_names_[new_num_row] = lp.row_names_[row];
1256 new_num_row++;
1257 }
1258 if (keep_to_row >= row_dim - 1) break;
1259 }
1260 lp.rowLower_.resize(new_num_row);
1261 lp.rowUpper_.resize(new_num_row);
1262 if (have_names) lp.row_names_.resize(new_num_row);
1263 return HighsStatus::OK;
1264 }
1265
deleteRowsFromLpMatrix(const HighsOptions & options,HighsLp & lp,const HighsIndexCollection & index_collection)1266 HighsStatus deleteRowsFromLpMatrix(
1267 const HighsOptions& options, HighsLp& lp,
1268 const HighsIndexCollection& index_collection) {
1269 HighsStatus return_status = HighsStatus::OK;
1270 if (!assessIndexCollection(options, index_collection))
1271 return interpretCallStatus(HighsStatus::Error, return_status,
1272 "assessIndexCollection");
1273 int from_k;
1274 int to_k;
1275 if (!limitsForIndexCollection(options, index_collection, from_k, to_k))
1276 return interpretCallStatus(HighsStatus::Error, return_status,
1277 "limitsForIndexCollection");
1278 if (index_collection.is_set_) {
1279 // For deletion by set it must be increasing
1280 if (!increasingSetOk(index_collection.set_,
1281 index_collection.set_num_entries_, 0, lp.numRow_ - 1,
1282 true))
1283 return HighsStatus::Error;
1284 }
1285 if (from_k > to_k) return HighsStatus::OK;
1286
1287 int delete_from_row;
1288 int delete_to_row;
1289 int keep_from_row;
1290 int row_dim = lp.numRow_;
1291 int keep_to_row = -1;
1292 int current_set_entry = 0;
1293
1294 // Set up a row mask to indicate the new row index of kept rows and
1295 // -1 for deleted rows so that the kept entries in the column-wise
1296 // matrix can be identified and have their correct row index.
1297 vector<int> new_index;
1298 new_index.resize(lp.numRow_);
1299 int new_num_row = 0;
1300 bool mask = index_collection.is_mask_;
1301 const int* row_mask = index_collection.mask_;
1302 if (!mask) {
1303 keep_to_row = -1;
1304 current_set_entry = 0;
1305 for (int k = from_k; k <= to_k; k++) {
1306 updateIndexCollectionOutInIndex(index_collection, delete_from_row,
1307 delete_to_row, keep_from_row, keep_to_row,
1308 current_set_entry);
1309 if (k == from_k) {
1310 // Account for any initial rows being kept
1311 for (int row = 0; row < delete_from_row; row++) {
1312 new_index[row] = new_num_row;
1313 new_num_row++;
1314 }
1315 }
1316 for (int row = delete_from_row; row <= delete_to_row; row++) {
1317 new_index[row] = -1;
1318 }
1319 for (int row = keep_from_row; row <= keep_to_row; row++) {
1320 new_index[row] = new_num_row;
1321 new_num_row++;
1322 }
1323 if (keep_to_row >= row_dim - 1) break;
1324 }
1325 } else {
1326 for (int row = 0; row < lp.numRow_; row++) {
1327 if (row_mask[row]) {
1328 new_index[row] = -1;
1329 } else {
1330 new_index[row] = new_num_row;
1331 new_num_row++;
1332 }
1333 }
1334 }
1335 int new_num_nz = 0;
1336 for (int col = 0; col < lp.numCol_; col++) {
1337 int from_el = lp.Astart_[col];
1338 lp.Astart_[col] = new_num_nz;
1339 for (int el = from_el; el < lp.Astart_[col + 1]; el++) {
1340 int row = lp.Aindex_[el];
1341 int new_row = new_index[row];
1342 if (new_row >= 0) {
1343 lp.Aindex_[new_num_nz] = new_row;
1344 lp.Avalue_[new_num_nz] = lp.Avalue_[el];
1345 new_num_nz++;
1346 }
1347 }
1348 }
1349 lp.Astart_[lp.numCol_] = new_num_nz;
1350 lp.Astart_.resize(lp.numCol_ + 1);
1351 lp.Aindex_.resize(new_num_nz);
1352 lp.Avalue_.resize(new_num_nz);
1353 return HighsStatus::OK;
1354 }
1355
changeLpMatrixCoefficient(HighsLp & lp,const int row,const int col,const double new_value)1356 HighsStatus changeLpMatrixCoefficient(HighsLp& lp, const int row, const int col,
1357 const double new_value) {
1358 if (row < 0 || row > lp.numRow_) return HighsStatus::Error;
1359 if (col < 0 || col > lp.numCol_) return HighsStatus::Error;
1360 int changeElement = -1;
1361 for (int el = lp.Astart_[col]; el < lp.Astart_[col + 1]; el++) {
1362 if (lp.Aindex_[el] == row) {
1363 changeElement = el;
1364 break;
1365 }
1366 }
1367 if (changeElement < 0) {
1368 changeElement = lp.Astart_[col + 1];
1369 int new_num_nz = lp.Astart_[lp.numCol_] + 1;
1370 lp.Aindex_.resize(new_num_nz);
1371 lp.Avalue_.resize(new_num_nz);
1372 for (int i = col + 1; i <= lp.numCol_; i++) lp.Astart_[i]++;
1373 for (int el = new_num_nz - 1; el > changeElement; el--) {
1374 lp.Aindex_[el] = lp.Aindex_[el - 1];
1375 lp.Avalue_[el] = lp.Avalue_[el - 1];
1376 }
1377 }
1378 lp.Aindex_[changeElement] = row;
1379 lp.Avalue_[changeElement] = new_value;
1380
1381 return HighsStatus::OK;
1382 }
1383
changeLpCosts(const HighsOptions & options,HighsLp & lp,const HighsIndexCollection & index_collection,const vector<double> & new_col_cost)1384 HighsStatus changeLpCosts(const HighsOptions& options, HighsLp& lp,
1385 const HighsIndexCollection& index_collection,
1386 const vector<double>& new_col_cost) {
1387 HighsStatus return_status = HighsStatus::OK;
1388 // Check parameters for technique and, if OK set the loop limits
1389 if (!assessIndexCollection(options, index_collection))
1390 return interpretCallStatus(HighsStatus::Error, return_status,
1391 "assessIndexCollection");
1392 int from_k;
1393 int to_k;
1394 if (!limitsForIndexCollection(options, index_collection, from_k, to_k))
1395 return interpretCallStatus(HighsStatus::Error, return_status,
1396 "limitsForIndexCollection");
1397 if (from_k > to_k) return HighsStatus::OK;
1398
1399 const bool& interval = index_collection.is_interval_;
1400 const bool& mask = index_collection.is_mask_;
1401 const int* col_set = index_collection.set_;
1402 const int* col_mask = index_collection.mask_;
1403
1404 // Change the costs to the user-supplied costs, according to the technique
1405 int usr_col;
1406 for (int k = from_k; k < to_k + 1; k++) {
1407 if (interval || mask) {
1408 usr_col = k;
1409 } else {
1410 usr_col = col_set[k];
1411 }
1412 int col = usr_col;
1413 if (mask && !col_mask[col]) continue;
1414 lp.colCost_[col] = new_col_cost[k];
1415 }
1416 return HighsStatus::OK;
1417 }
1418
changeLpColBounds(const HighsOptions & options,HighsLp & lp,const HighsIndexCollection & index_collection,const vector<double> & new_col_lower,const vector<double> & new_col_upper)1419 HighsStatus changeLpColBounds(const HighsOptions& options, HighsLp& lp,
1420 const HighsIndexCollection& index_collection,
1421 const vector<double>& new_col_lower,
1422 const vector<double>& new_col_upper) {
1423 return changeBounds(options, lp.colLower_, lp.colUpper_, index_collection,
1424 new_col_lower, new_col_upper);
1425 }
1426
changeLpRowBounds(const HighsOptions & options,HighsLp & lp,const HighsIndexCollection & index_collection,const vector<double> & new_row_lower,const vector<double> & new_row_upper)1427 HighsStatus changeLpRowBounds(const HighsOptions& options, HighsLp& lp,
1428 const HighsIndexCollection& index_collection,
1429 const vector<double>& new_row_lower,
1430 const vector<double>& new_row_upper) {
1431 return changeBounds(options, lp.rowLower_, lp.rowUpper_, index_collection,
1432 new_row_lower, new_row_upper);
1433 }
1434
changeBounds(const HighsOptions & options,vector<double> & lower,vector<double> & upper,const HighsIndexCollection & index_collection,const vector<double> & new_lower,const vector<double> & new_upper)1435 HighsStatus changeBounds(const HighsOptions& options, vector<double>& lower,
1436 vector<double>& upper,
1437 const HighsIndexCollection& index_collection,
1438 const vector<double>& new_lower,
1439 const vector<double>& new_upper) {
1440 HighsStatus return_status = HighsStatus::OK;
1441 // Check parameters for technique and, if OK set the loop limits
1442 if (!assessIndexCollection(options, index_collection))
1443 return interpretCallStatus(HighsStatus::Error, return_status,
1444 "assessIndexCollection");
1445 int from_k;
1446 int to_k;
1447 if (!limitsForIndexCollection(options, index_collection, from_k, to_k))
1448 return interpretCallStatus(HighsStatus::Error, return_status,
1449 "limitsForIndexCollection");
1450 if (from_k > to_k) return HighsStatus::OK;
1451
1452 const bool& interval = index_collection.is_interval_;
1453 const bool& mask = index_collection.is_mask_;
1454 const int* ix_set = index_collection.set_;
1455 const int* ix_mask = index_collection.mask_;
1456
1457 // Change the bounds to the user-supplied bounds, according to the technique
1458 int usr_ix;
1459 for (int k = from_k; k < to_k + 1; k++) {
1460 if (interval || mask) {
1461 usr_ix = k;
1462 } else {
1463 usr_ix = ix_set[k];
1464 }
1465 int ix = usr_ix;
1466 if (mask && !ix_mask[ix]) continue;
1467 lower[ix] = new_lower[k];
1468 upper[ix] = new_upper[k];
1469 }
1470 return HighsStatus::OK;
1471 }
1472
getNumInt(const HighsLp & lp)1473 int getNumInt(const HighsLp& lp) {
1474 int num_int = 0;
1475 if (lp.integrality_.size()) {
1476 for (int iCol = 0; iCol < lp.numCol_; iCol++)
1477 if (lp.integrality_[iCol] == HighsVarType::INTEGER) num_int++;
1478 }
1479 return num_int;
1480 }
1481
getLpCosts(const HighsLp & lp,const int from_col,const int to_col,double * XcolCost)1482 HighsStatus getLpCosts(const HighsLp& lp, const int from_col, const int to_col,
1483 double* XcolCost) {
1484 if (from_col < 0 || to_col >= lp.numCol_) return HighsStatus::Error;
1485 if (from_col > to_col) return HighsStatus::OK;
1486 for (int col = from_col; col < to_col + 1; col++)
1487 XcolCost[col - from_col] = lp.colCost_[col];
1488 return HighsStatus::OK;
1489 }
1490
getLpColBounds(const HighsLp & lp,const int from_col,const int to_col,double * XcolLower,double * XcolUpper)1491 HighsStatus getLpColBounds(const HighsLp& lp, const int from_col,
1492 const int to_col, double* XcolLower,
1493 double* XcolUpper) {
1494 if (from_col < 0 || to_col >= lp.numCol_) return HighsStatus::Error;
1495 if (from_col > to_col) return HighsStatus::OK;
1496 for (int col = from_col; col < to_col + 1; col++) {
1497 if (XcolLower != NULL) XcolLower[col - from_col] = lp.colLower_[col];
1498 if (XcolUpper != NULL) XcolUpper[col - from_col] = lp.colUpper_[col];
1499 }
1500 return HighsStatus::OK;
1501 }
1502
getLpRowBounds(const HighsLp & lp,const int from_row,const int to_row,double * XrowLower,double * XrowUpper)1503 HighsStatus getLpRowBounds(const HighsLp& lp, const int from_row,
1504 const int to_row, double* XrowLower,
1505 double* XrowUpper) {
1506 if (from_row < 0 || to_row >= lp.numRow_) return HighsStatus::Error;
1507 if (from_row > to_row) return HighsStatus::OK;
1508 for (int row = from_row; row < to_row + 1; row++) {
1509 if (XrowLower != NULL) XrowLower[row - from_row] = lp.rowLower_[row];
1510 if (XrowUpper != NULL) XrowUpper[row - from_row] = lp.rowUpper_[row];
1511 }
1512 return HighsStatus::OK;
1513 }
1514
1515 // Get a single coefficient from the matrix
getLpMatrixCoefficient(const HighsLp & lp,const int Xrow,const int Xcol,double * val)1516 HighsStatus getLpMatrixCoefficient(const HighsLp& lp, const int Xrow,
1517 const int Xcol, double* val) {
1518 if (Xrow < 0 || Xrow >= lp.numRow_) return HighsStatus::Error;
1519 if (Xcol < 0 || Xcol >= lp.numCol_) return HighsStatus::Error;
1520
1521 int get_el = -1;
1522 for (int el = lp.Astart_[Xcol]; el < lp.Astart_[Xcol + 1]; el++) {
1523 if (lp.Aindex_[el] == Xrow) {
1524 get_el = el;
1525 break;
1526 }
1527 }
1528 if (get_el < 0) {
1529 *val = 0;
1530 } else {
1531 *val = lp.Avalue_[get_el];
1532 }
1533 return HighsStatus::OK;
1534 }
1535
1536 // Methods for reporting an LP, including its row and column data and matrix
1537 //
1538 // Report the whole LP
reportLp(const HighsOptions & options,const HighsLp & lp,const int report_level)1539 void reportLp(const HighsOptions& options, const HighsLp& lp,
1540 const int report_level) {
1541 reportLpBrief(options, lp);
1542 if (report_level >= 1) {
1543 reportLpColVectors(options, lp);
1544 reportLpRowVectors(options, lp);
1545 if (report_level >= 2) reportLpColMatrix(options, lp);
1546 }
1547 }
1548
1549 // Report the LP briefly
reportLpBrief(const HighsOptions & options,const HighsLp & lp)1550 void reportLpBrief(const HighsOptions& options, const HighsLp& lp) {
1551 reportLpDimensions(options, lp);
1552 reportLpObjSense(options, lp);
1553 }
1554
1555 // Report the LP dimensions
reportLpDimensions(const HighsOptions & options,const HighsLp & lp)1556 void reportLpDimensions(const HighsOptions& options, const HighsLp& lp) {
1557 int lp_num_nz;
1558 if (lp.numCol_ == 0)
1559 lp_num_nz = 0;
1560 else
1561 lp_num_nz = lp.Astart_[lp.numCol_];
1562 HighsPrintMessage(options.output, options.message_level, ML_MINIMAL,
1563 "LP has %d columns, %d rows", lp.numCol_, lp.numRow_);
1564 int num_int = getNumInt(lp);
1565 if (num_int) {
1566 HighsPrintMessage(options.output, options.message_level, ML_MINIMAL,
1567 ", %d nonzeros and %d integer columns\n", lp_num_nz,
1568 num_int);
1569 } else {
1570 HighsPrintMessage(options.output, options.message_level, ML_MINIMAL,
1571 " and %d nonzeros\n", lp_num_nz, num_int);
1572 }
1573 }
1574
1575 // Report the LP objective sense
reportLpObjSense(const HighsOptions & options,const HighsLp & lp)1576 void reportLpObjSense(const HighsOptions& options, const HighsLp& lp) {
1577 if (lp.sense_ == ObjSense::MINIMIZE)
1578 HighsPrintMessage(options.output, options.message_level, ML_MINIMAL,
1579 "Objective sense is minimize\n");
1580 else if (lp.sense_ == ObjSense::MAXIMIZE)
1581 HighsPrintMessage(options.output, options.message_level, ML_MINIMAL,
1582 "Objective sense is maximize\n");
1583 else
1584 HighsPrintMessage(options.output, options.message_level, ML_MINIMAL,
1585 "Objective sense is ill-defined as %d\n", lp.sense_);
1586 }
1587
getBoundType(const double lower,const double upper)1588 std::string getBoundType(const double lower, const double upper) {
1589 std::string type;
1590 if (highs_isInfinity(-lower)) {
1591 if (highs_isInfinity(upper)) {
1592 type = "FR";
1593 } else {
1594 type = "UB";
1595 }
1596 } else {
1597 if (highs_isInfinity(upper)) {
1598 type = "LB";
1599 } else {
1600 if (lower < upper) {
1601 type = "BX";
1602 } else {
1603 type = "FX";
1604 }
1605 }
1606 }
1607 return type;
1608 }
1609
1610 // Report the vectors of LP column data
reportLpColVectors(const HighsOptions & options,const HighsLp & lp)1611 void reportLpColVectors(const HighsOptions& options, const HighsLp& lp) {
1612 if (lp.numCol_ <= 0) return;
1613 std::string type;
1614 int count;
1615 bool have_integer_columns = getNumInt(lp);
1616 bool have_col_names = lp.col_names_.size();
1617
1618 HighsPrintMessage(options.output, options.message_level, ML_VERBOSE,
1619 " Column Lower Upper Cost "
1620 "Type Count");
1621 if (have_integer_columns)
1622 HighsPrintMessage(options.output, options.message_level, ML_VERBOSE,
1623 " Discrete");
1624 if (have_col_names)
1625 HighsPrintMessage(options.output, options.message_level, ML_VERBOSE,
1626 " Name");
1627 HighsPrintMessage(options.output, options.message_level, ML_VERBOSE, "\n");
1628
1629 for (int iCol = 0; iCol < lp.numCol_; iCol++) {
1630 type = getBoundType(lp.colLower_[iCol], lp.colUpper_[iCol]);
1631 count = lp.Astart_[iCol + 1] - lp.Astart_[iCol];
1632 HighsPrintMessage(options.output, options.message_level, ML_VERBOSE,
1633 "%8d %12g %12g %12g %2s %12d", iCol,
1634 lp.colLower_[iCol], lp.colUpper_[iCol], lp.colCost_[iCol],
1635 type.c_str(), count);
1636 if (have_integer_columns) {
1637 std::string integer_column = "";
1638 if (lp.integrality_[iCol] == HighsVarType::INTEGER) {
1639 if (lp.colLower_[iCol] == 0 && lp.colUpper_[iCol] == 1) {
1640 integer_column = "Binary";
1641 } else {
1642 integer_column = "Integer";
1643 }
1644 }
1645 HighsPrintMessage(options.output, options.message_level, ML_VERBOSE,
1646 " %-8s", integer_column.c_str());
1647 }
1648 if (have_col_names)
1649 HighsPrintMessage(options.output, options.message_level, ML_VERBOSE,
1650 " %-s", lp.col_names_[iCol].c_str());
1651 HighsPrintMessage(options.output, options.message_level, ML_VERBOSE, "\n");
1652 }
1653 }
1654
1655 // Report the vectors of LP row data
reportLpRowVectors(const HighsOptions & options,const HighsLp & lp)1656 void reportLpRowVectors(const HighsOptions& options, const HighsLp& lp) {
1657 if (lp.numRow_ <= 0) return;
1658 std::string type;
1659 vector<int> count;
1660 bool have_row_names = lp.row_names_.size();
1661
1662 count.resize(lp.numRow_, 0);
1663 if (lp.numCol_ > 0) {
1664 for (int el = 0; el < lp.Astart_[lp.numCol_]; el++) count[lp.Aindex_[el]]++;
1665 }
1666
1667 HighsPrintMessage(
1668 options.output, options.message_level, ML_VERBOSE,
1669 " Row Lower Upper Type Count");
1670 if (have_row_names)
1671 HighsPrintMessage(options.output, options.message_level, ML_VERBOSE,
1672 " Name");
1673 HighsPrintMessage(options.output, options.message_level, ML_VERBOSE, "\n");
1674
1675 for (int iRow = 0; iRow < lp.numRow_; iRow++) {
1676 type = getBoundType(lp.rowLower_[iRow], lp.rowUpper_[iRow]);
1677 std::string name = "";
1678 HighsPrintMessage(options.output, options.message_level, ML_VERBOSE,
1679 "%8d %12g %12g %2s %12d", iRow,
1680 lp.rowLower_[iRow], lp.rowUpper_[iRow], type.c_str(),
1681 count[iRow]);
1682 if (have_row_names)
1683 HighsPrintMessage(options.output, options.message_level, ML_VERBOSE,
1684 " %-s", lp.row_names_[iRow].c_str());
1685 HighsPrintMessage(options.output, options.message_level, ML_VERBOSE, "\n");
1686 }
1687 }
1688
1689 // Report the LP column-wise matrix
reportLpColMatrix(const HighsOptions & options,const HighsLp & lp)1690 void reportLpColMatrix(const HighsOptions& options, const HighsLp& lp) {
1691 if (lp.numCol_ <= 0) return;
1692 if (lp.numRow_) {
1693 // With postitive number of rows, can assume that there are index and value
1694 // vectors to pass
1695 reportMatrix(options, "Column", lp.numCol_, lp.Astart_[lp.numCol_],
1696 &lp.Astart_[0], &lp.Aindex_[0], &lp.Avalue_[0]);
1697 } else {
1698 // With no rows, can's assume that there are index and value vectors to pass
1699 reportMatrix(options, "Column", lp.numCol_, lp.Astart_[lp.numCol_],
1700 &lp.Astart_[0], NULL, NULL);
1701 }
1702 }
1703
reportMatrix(const HighsOptions & options,const std::string message,const int num_col,const int num_nz,const int * start,const int * index,const double * value)1704 void reportMatrix(const HighsOptions& options, const std::string message,
1705 const int num_col, const int num_nz, const int* start,
1706 const int* index, const double* value) {
1707 if (num_col <= 0) return;
1708 HighsPrintMessage(options.output, options.message_level, ML_VERBOSE,
1709 "%6s Index Value\n", message.c_str());
1710 for (int col = 0; col < num_col; col++) {
1711 HighsPrintMessage(options.output, options.message_level, ML_VERBOSE,
1712 " %8d Start %10d\n", col, start[col]);
1713 int to_el = (col < num_col - 1 ? start[col + 1] : num_nz);
1714 for (int el = start[col]; el < to_el; el++)
1715 HighsPrintMessage(options.output, options.message_level, ML_VERBOSE,
1716 " %8d %12g\n", index[el], value[el]);
1717 }
1718 HighsPrintMessage(options.output, options.message_level, ML_VERBOSE,
1719 " Start %10d\n", num_nz);
1720 }
1721
1722 #ifdef HiGHSDEV
analyseLp(const HighsLp & lp,const std::string message)1723 void analyseLp(const HighsLp& lp, const std::string message) {
1724 vector<double> min_colBound;
1725 vector<double> min_rowBound;
1726 vector<double> colRange;
1727 vector<double> rowRange;
1728 min_colBound.resize(lp.numCol_);
1729 min_rowBound.resize(lp.numRow_);
1730 colRange.resize(lp.numCol_);
1731 rowRange.resize(lp.numRow_);
1732 for (int col = 0; col < lp.numCol_; col++)
1733 min_colBound[col] = min(fabs(lp.colLower_[col]), fabs(lp.colUpper_[col]));
1734 for (int row = 0; row < lp.numRow_; row++)
1735 min_rowBound[row] = min(fabs(lp.rowLower_[row]), fabs(lp.rowUpper_[row]));
1736 for (int col = 0; col < lp.numCol_; col++)
1737 colRange[col] = lp.colUpper_[col] - lp.colLower_[col];
1738 for (int row = 0; row < lp.numRow_; row++)
1739 rowRange[row] = lp.rowUpper_[row] - lp.rowLower_[row];
1740
1741 printf("\n%s model data: Analysis\n", message.c_str());
1742 analyseVectorValues("Column costs", lp.numCol_, lp.colCost_);
1743 analyseVectorValues("Column lower bounds", lp.numCol_, lp.colLower_);
1744 analyseVectorValues("Column upper bounds", lp.numCol_, lp.colUpper_);
1745 analyseVectorValues("Column min abs bound", lp.numCol_, min_colBound);
1746 analyseVectorValues("Column range", lp.numCol_, colRange);
1747 analyseVectorValues("Row lower bounds", lp.numRow_, lp.rowLower_);
1748 analyseVectorValues("Row upper bounds", lp.numRow_, lp.rowUpper_);
1749 analyseVectorValues("Row min abs bound", lp.numRow_, min_rowBound);
1750 analyseVectorValues("Row range", lp.numRow_, rowRange);
1751 analyseVectorValues("Matrix sparsity", lp.Astart_[lp.numCol_], lp.Avalue_,
1752 true, lp.model_name_);
1753 analyseMatrixSparsity("Constraint matrix", lp.numCol_, lp.numRow_, lp.Astart_,
1754 lp.Aindex_);
1755 analyseModelBounds("Column", lp.numCol_, lp.colLower_, lp.colUpper_);
1756 analyseModelBounds("Row", lp.numRow_, lp.rowLower_, lp.rowUpper_);
1757 }
1758 #endif
1759
writeSolutionToFile(FILE * file,const HighsLp & lp,const HighsBasis & basis,const HighsSolution & solution,const bool pretty)1760 void writeSolutionToFile(FILE* file, const HighsLp& lp, const HighsBasis& basis,
1761 const HighsSolution& solution, const bool pretty) {
1762 if (pretty) {
1763 writeModelBoundSol(file, true, lp.numCol_, lp.colLower_, lp.colUpper_,
1764 lp.col_names_, solution.col_value, solution.col_dual,
1765 basis.col_status);
1766 writeModelBoundSol(file, false, lp.numRow_, lp.rowLower_, lp.rowUpper_,
1767 lp.row_names_, solution.row_value, solution.row_dual,
1768 basis.row_status);
1769 } else {
1770 fprintf(file,
1771 "%d %d : Number of columns and rows for primal and dual solution "
1772 "and basis\n",
1773 lp.numCol_, lp.numRow_);
1774 const bool with_basis = basis.valid_;
1775 if (with_basis) {
1776 fprintf(file, "T\n");
1777 } else {
1778 fprintf(file, "F\n");
1779 }
1780 for (int iCol = 0; iCol < lp.numCol_; iCol++) {
1781 fprintf(file, "%g %g", solution.col_value[iCol], solution.col_dual[iCol]);
1782 if (with_basis) fprintf(file, " %d", (int)basis.col_status[iCol]);
1783 fprintf(file, " \n");
1784 }
1785 for (int iRow = 0; iRow < lp.numRow_; iRow++) {
1786 fprintf(file, "%g %g", solution.row_value[iRow], solution.row_dual[iRow]);
1787 if (with_basis) fprintf(file, " %d", (int)basis.row_status[iRow]);
1788 fprintf(file, " \n");
1789 }
1790 }
1791 }
1792
writeBasisFile(const HighsOptions & options,const HighsBasis & basis,const std::string filename)1793 HighsStatus writeBasisFile(const HighsOptions& options, const HighsBasis& basis,
1794 const std::string filename) {
1795 HighsStatus return_status = HighsStatus::OK;
1796 if (basis.valid_ == false) {
1797 HighsLogMessage(options.logfile, HighsMessageType::ERROR,
1798 "writeBasisFile: Cannot write an invalid basis");
1799 return HighsStatus::Error;
1800 }
1801 std::ofstream outFile(filename);
1802 if (outFile.fail()) {
1803 HighsLogMessage(options.logfile, HighsMessageType::ERROR,
1804 "writeBasisFile: Cannot open writeable file \"%s\"",
1805 filename.c_str());
1806 return HighsStatus::Error;
1807 }
1808 outFile << "HiGHS Version " << HIGHS_VERSION_MAJOR << std::endl;
1809 outFile << basis.col_status.size() << " " << basis.row_status.size()
1810 << std::endl;
1811 for (const auto& status : basis.col_status) {
1812 outFile << (int)status << " ";
1813 }
1814 outFile << std::endl;
1815 for (const auto& status : basis.row_status) {
1816 outFile << (int)status << " ";
1817 }
1818 outFile << std::endl;
1819 outFile << std::endl;
1820 outFile.close();
1821 return return_status;
1822 }
1823
readBasisFile(const HighsOptions & options,HighsBasis & basis,const std::string filename)1824 HighsStatus readBasisFile(const HighsOptions& options, HighsBasis& basis,
1825 const std::string filename) {
1826 // Reads a basis file, returning an error if what's read is
1827 // inconsistent with the sizes of the HighsBasis passed in
1828 HighsStatus return_status = HighsStatus::OK;
1829 std::ifstream inFile(filename);
1830 if (inFile.fail()) {
1831 HighsLogMessage(options.logfile, HighsMessageType::ERROR,
1832 "readBasisFile: Cannot open readable file \"%s\"",
1833 filename.c_str());
1834 return HighsStatus::Error;
1835 }
1836 std::string string_highs, string_version;
1837 int highs_version_number;
1838 inFile >> string_highs >> string_version >> highs_version_number;
1839 if (highs_version_number == 1) {
1840 int numCol, numRow;
1841 inFile >> numCol >> numRow;
1842 int basis_numCol = (int)basis.col_status.size();
1843 int basis_numRow = (int)basis.row_status.size();
1844 if (numCol != basis_numCol) {
1845 HighsLogMessage(options.logfile, HighsMessageType::ERROR,
1846 "readBasisFile: Basis file is for %d columns, not %d",
1847 numCol, basis_numCol);
1848 return HighsStatus::Error;
1849 }
1850 if (numRow != basis_numRow) {
1851 HighsLogMessage(options.logfile, HighsMessageType::ERROR,
1852 "readBasisFile: Basis file is for %d rows, not %d",
1853 numRow, basis_numRow);
1854 return HighsStatus::Error;
1855 }
1856 int int_status;
1857 for (int iCol = 0; iCol < numCol; iCol++) {
1858 inFile >> int_status;
1859 basis.col_status[iCol] = (HighsBasisStatus)int_status;
1860 }
1861 for (int iRow = 0; iRow < numRow; iRow++) {
1862 inFile >> int_status;
1863 basis.row_status[iRow] = (HighsBasisStatus)int_status;
1864 }
1865 if (inFile.eof()) {
1866 HighsLogMessage(
1867 options.logfile, HighsMessageType::ERROR,
1868 "readBasisFile: Reached end of file before reading complete basis");
1869 return_status = HighsStatus::Error;
1870 }
1871 } else {
1872 HighsLogMessage(
1873 options.logfile, HighsMessageType::ERROR,
1874 "readBasisFile: Cannot read basis file for HiGHS version %d",
1875 highs_version_number);
1876 return_status = HighsStatus::Error;
1877 }
1878 inFile.close();
1879 return return_status;
1880 }
1881
calculateColDuals(const HighsLp & lp,HighsSolution & solution)1882 HighsStatus calculateColDuals(const HighsLp& lp, HighsSolution& solution) {
1883 assert(solution.row_dual.size() > 0);
1884 if (!isSolutionRightSize(lp, solution)) return HighsStatus::Error;
1885
1886 solution.col_dual.assign(lp.numCol_, 0);
1887
1888 for (int col = 0; col < lp.numCol_; col++) {
1889 for (int i = lp.Astart_[col]; i < lp.Astart_[col + 1]; i++) {
1890 const int row = lp.Aindex_[i];
1891 assert(row >= 0);
1892 assert(row < lp.numRow_);
1893
1894 solution.col_dual[col] -= solution.row_dual[row] * lp.Avalue_[i];
1895 }
1896 solution.col_dual[col] += lp.colCost_[col];
1897 }
1898
1899 return HighsStatus::OK;
1900 }
1901
calculateRowValues(const HighsLp & lp,HighsSolution & solution)1902 HighsStatus calculateRowValues(const HighsLp& lp, HighsSolution& solution) {
1903 assert(solution.col_value.size() > 0);
1904 if (int(solution.col_value.size()) != lp.numCol_) return HighsStatus::Error;
1905
1906 solution.row_value.clear();
1907 solution.row_value.assign(lp.numRow_, 0);
1908
1909 for (int col = 0; col < lp.numCol_; col++) {
1910 for (int i = lp.Astart_[col]; i < lp.Astart_[col + 1]; i++) {
1911 const int row = lp.Aindex_[i];
1912 assert(row >= 0);
1913 assert(row < lp.numRow_);
1914
1915 solution.row_value[row] += solution.col_value[col] * lp.Avalue_[i];
1916 }
1917 }
1918
1919 return HighsStatus::OK;
1920 }
1921
calculateObjective(const HighsLp & lp,HighsSolution & solution)1922 double calculateObjective(const HighsLp& lp, HighsSolution& solution) {
1923 assert(isSolutionRightSize(lp, solution));
1924 double sum = 0;
1925 for (int col = 0; col < lp.numCol_; col++)
1926 sum += lp.colCost_[col] * solution.col_value[col];
1927
1928 return sum;
1929 }
1930
isColDataNull(const HighsOptions & options,const double * usr_col_cost,const double * usr_col_lower,const double * usr_col_upper)1931 bool isColDataNull(const HighsOptions& options, const double* usr_col_cost,
1932 const double* usr_col_lower, const double* usr_col_upper) {
1933 bool null_data = false;
1934 null_data =
1935 doubleUserDataNotNull(options.logfile, usr_col_cost, "column costs") ||
1936 null_data;
1937 null_data = doubleUserDataNotNull(options.logfile, usr_col_lower,
1938 "column lower bounds") ||
1939 null_data;
1940 null_data = doubleUserDataNotNull(options.logfile, usr_col_upper,
1941 "column upper bounds") ||
1942 null_data;
1943 return null_data;
1944 }
1945
isRowDataNull(const HighsOptions & options,const double * usr_row_lower,const double * usr_row_upper)1946 bool isRowDataNull(const HighsOptions& options, const double* usr_row_lower,
1947 const double* usr_row_upper) {
1948 bool null_data = false;
1949 null_data = doubleUserDataNotNull(options.logfile, usr_row_lower,
1950 "row lower bounds") ||
1951 null_data;
1952 null_data = doubleUserDataNotNull(options.logfile, usr_row_upper,
1953 "row upper bounds") ||
1954 null_data;
1955 return null_data;
1956 }
1957
isMatrixDataNull(const HighsOptions & options,const int * usr_matrix_start,const int * usr_matrix_index,const double * usr_matrix_value)1958 bool isMatrixDataNull(const HighsOptions& options, const int* usr_matrix_start,
1959 const int* usr_matrix_index,
1960 const double* usr_matrix_value) {
1961 bool null_data = false;
1962 null_data =
1963 intUserDataNotNull(options.logfile, usr_matrix_start, "matrix starts") ||
1964 null_data;
1965 null_data =
1966 intUserDataNotNull(options.logfile, usr_matrix_index, "matrix indices") ||
1967 null_data;
1968 null_data = doubleUserDataNotNull(options.logfile, usr_matrix_value,
1969 "matrix values") ||
1970 null_data;
1971 return null_data;
1972 }
1973
transformIntoEqualityProblem(const HighsLp & lp,HighsLp & equality_lp)1974 HighsStatus transformIntoEqualityProblem(const HighsLp& lp,
1975 HighsLp& equality_lp) {
1976 // Copy lp.
1977 equality_lp = lp;
1978
1979 // Add slacks for each row with more than one bound.
1980 std::vector<double> rhs(lp.numRow_, 0);
1981
1982 for (int row = 0; row < lp.numRow_; row++) {
1983 assert(equality_lp.Astart_[equality_lp.numCol_] ==
1984 (int)equality_lp.Avalue_.size());
1985 assert((int)equality_lp.Aindex_.size() == (int)equality_lp.Avalue_.size());
1986 const int nnz = equality_lp.Astart_[equality_lp.numCol_];
1987
1988 if (lp.rowLower_[row] <= -HIGHS_CONST_INF &&
1989 lp.rowUpper_[row] >= HIGHS_CONST_INF) {
1990 // free row
1991 equality_lp.Astart_.push_back(nnz + 1);
1992 equality_lp.Aindex_.push_back(row);
1993 equality_lp.Avalue_.push_back(1.0);
1994
1995 equality_lp.numCol_++;
1996 equality_lp.colLower_.push_back(-HIGHS_CONST_INF);
1997 equality_lp.colUpper_.push_back(HIGHS_CONST_INF);
1998 equality_lp.colCost_.push_back(0);
1999 } else if (lp.rowLower_[row] > -HIGHS_CONST_INF &&
2000 lp.rowUpper_[row] >= HIGHS_CONST_INF) {
2001 // only lower bound
2002 rhs[row] = lp.rowLower_[row];
2003
2004 equality_lp.Astart_.push_back(nnz + 1);
2005 equality_lp.Aindex_.push_back(row);
2006 equality_lp.Avalue_.push_back(-1.0);
2007
2008 equality_lp.numCol_++;
2009 equality_lp.colLower_.push_back(0);
2010 equality_lp.colUpper_.push_back(HIGHS_CONST_INF);
2011 equality_lp.colCost_.push_back(0);
2012 } else if (lp.rowLower_[row] <= -HIGHS_CONST_INF &&
2013 lp.rowUpper_[row] < HIGHS_CONST_INF) {
2014 // only upper bound
2015 rhs[row] = lp.rowUpper_[row];
2016
2017 equality_lp.Astart_.push_back(nnz + 1);
2018 equality_lp.Aindex_.push_back(row);
2019 equality_lp.Avalue_.push_back(1.0);
2020
2021 equality_lp.numCol_++;
2022 equality_lp.colLower_.push_back(0);
2023 equality_lp.colUpper_.push_back(HIGHS_CONST_INF);
2024 equality_lp.colCost_.push_back(0);
2025 } else if (lp.rowLower_[row] > -HIGHS_CONST_INF &&
2026 lp.rowUpper_[row] < HIGHS_CONST_INF &&
2027 lp.rowLower_[row] != lp.rowUpper_[row]) {
2028 // both lower and upper bound that are different
2029 double rhs_value, coefficient;
2030 double difference = lp.rowUpper_[row] - lp.rowLower_[row];
2031 if (fabs(lp.rowLower_[row]) < fabs(lp.rowUpper_[row])) {
2032 rhs_value = lp.rowLower_[row];
2033 coefficient = -1;
2034 } else {
2035 rhs_value = lp.rowUpper_[row];
2036 coefficient = 1;
2037 }
2038 rhs[row] = rhs_value;
2039
2040 equality_lp.Astart_.push_back(nnz + 1);
2041 equality_lp.Aindex_.push_back(row);
2042 equality_lp.Avalue_.push_back(coefficient);
2043
2044 equality_lp.numCol_++;
2045 equality_lp.colLower_.push_back(0);
2046 equality_lp.colUpper_.push_back(difference);
2047 equality_lp.colCost_.push_back(0);
2048 } else if (lp.rowLower_[row] == lp.rowUpper_[row]) {
2049 // equality row
2050 rhs[row] = lp.rowLower_[row];
2051 } else {
2052 #ifdef HiGHSDEV
2053 printf(
2054 "transformIntoEqualityProblem: Unknown row type when adding slacks");
2055 #endif
2056 return HighsStatus::Error;
2057 }
2058 }
2059 equality_lp.rowLower_ = rhs;
2060 equality_lp.rowUpper_ = rhs;
2061
2062 return HighsStatus::OK;
2063 }
2064
2065 // Given (P) returns (D) for the pair
2066 // (P)
2067 // min c'x st Ax=b
2068 // st l <= x <= u
2069 // (D)
2070 // max b'y + l'zl - u'zu
2071 // st A'y + zl - zu = c
2072 // y free, zl >=0, zu >= 0
dualizeEqualityProblem(const HighsLp & lp,HighsLp & dual)2073 HighsStatus dualizeEqualityProblem(const HighsLp& lp, HighsLp& dual) {
2074 std::vector<double> colCost = lp.colCost_;
2075 if (lp.sense_ != ObjSense::MINIMIZE) {
2076 for (int col = 0; col < lp.numCol_; col++) colCost[col] = -colCost[col];
2077 }
2078
2079 assert(lp.rowLower_ == lp.rowUpper_);
2080
2081 const int ncols = lp.numRow_;
2082 const int nrows = lp.numCol_;
2083
2084 dual.numRow_ = nrows;
2085 dual.rowLower_ = colCost;
2086 dual.rowUpper_ = colCost;
2087
2088 // Add columns (y)
2089 dual.numCol_ = ncols;
2090 dual.colLower_.resize(ncols);
2091 dual.colUpper_.resize(ncols);
2092 dual.colCost_.resize(ncols);
2093
2094 for (int col = 0; col < ncols; col++) {
2095 dual.colLower_[col] = -HIGHS_CONST_INF;
2096 dual.colUpper_[col] = HIGHS_CONST_INF;
2097 // cost b'y
2098 dual.colCost_[col] = lp.rowLower_[col];
2099 }
2100
2101 // Get transpose of A
2102 int i, k;
2103 vector<int> iwork(lp.numRow_, 0);
2104 dual.Astart_.resize(lp.numRow_ + 1, 0);
2105 int AcountX = lp.Aindex_.size();
2106 dual.Aindex_.resize(AcountX);
2107 dual.Avalue_.resize(AcountX);
2108 for (int k = 0; k < AcountX; k++) iwork.at(lp.Aindex_.at(k))++;
2109 for (i = 1; i <= lp.numRow_; i++)
2110 dual.Astart_.at(i) = dual.Astart_.at(i - 1) + iwork.at(i - 1);
2111 for (i = 0; i < lp.numRow_; i++) iwork.at(i) = dual.Astart_.at(i);
2112 for (int iCol = 0; iCol < lp.numCol_; iCol++) {
2113 for (k = lp.Astart_.at(iCol); k < lp.Astart_.at(iCol + 1); k++) {
2114 int iRow = lp.Aindex_.at(k);
2115 int iPut = iwork.at(iRow)++;
2116 dual.Aindex_.at(iPut) = iCol;
2117 dual.Avalue_.at(iPut) = lp.Avalue_[k];
2118 }
2119 }
2120
2121 // Add columns (zl)
2122 for (int col = 0; col < lp.numCol_; col++) {
2123 if (lp.colLower_[col] > -HIGHS_CONST_INF) {
2124 const int nnz = dual.Astart_[dual.numCol_];
2125
2126 dual.colLower_.push_back(0);
2127 dual.colUpper_.push_back(HIGHS_CONST_INF);
2128
2129 dual.colCost_.push_back(lp.colLower_[col]);
2130
2131 // Add constaints
2132 dual.Astart_.push_back(nnz + 1);
2133 dual.Aindex_.push_back(col);
2134 dual.Avalue_.push_back(1.0);
2135
2136 dual.numCol_++;
2137 }
2138 }
2139
2140 // Add columns (zu)
2141 for (int col = 0; col < lp.numCol_; col++) {
2142 if (lp.colUpper_[col] < HIGHS_CONST_INF) {
2143 const int nnz = dual.Astart_[dual.numCol_];
2144
2145 dual.colLower_.push_back(0);
2146 dual.colUpper_.push_back(HIGHS_CONST_INF);
2147
2148 dual.colCost_.push_back(-lp.colUpper_[col]);
2149
2150 // Add constaints
2151 dual.Astart_.push_back(nnz + 1);
2152 dual.Aindex_.push_back(col);
2153 dual.Avalue_.push_back(-1.0);
2154
2155 dual.numCol_++;
2156 }
2157 }
2158
2159 dual.sense_ = ObjSense::MINIMIZE;
2160 for (int col = 0; col < dual.numCol_; col++) {
2161 dual.colCost_[col] = -dual.colCost_[col];
2162 }
2163
2164 dual.model_name_ = lp.model_name_ + "_dualized";
2165
2166 #ifdef HiGHSDEV
2167 printf("Dualized equality LP\n");
2168 #endif
2169
2170 return HighsStatus::OK;
2171 }
2172
reportPresolveReductions(const HighsOptions & options,const HighsLp & lp,const HighsLp & presolve_lp)2173 void reportPresolveReductions(const HighsOptions& options, const HighsLp& lp,
2174 const HighsLp& presolve_lp) {
2175 int num_col_from = lp.numCol_;
2176 int num_row_from = lp.numRow_;
2177 int num_els_from = lp.Astart_[num_col_from];
2178 int num_col_to = presolve_lp.numCol_;
2179 int num_row_to = presolve_lp.numRow_;
2180 int num_els_to;
2181 if (num_col_to) {
2182 num_els_to = presolve_lp.Astart_[num_col_to];
2183 } else {
2184 num_els_to = 0;
2185 }
2186 char elemsignchar = '-';
2187 int elemdelta = num_els_from - num_els_to;
2188 if (num_els_from < num_els_to) {
2189 elemdelta = -elemdelta;
2190 elemsignchar = '+';
2191 }
2192 HighsPrintMessage(options.logfile, options.message_level, ML_ALWAYS,
2193 "Presolve : Reductions: rows %d(-%d); columns %d(-%d); "
2194 "elements %d(%c%d)\n",
2195 num_row_to, (num_row_from - num_row_to), num_col_to,
2196 (num_col_from - num_col_to), num_els_to, elemsignchar,
2197 elemdelta);
2198 }
2199
reportPresolveReductions(const HighsOptions & options,const HighsLp & lp,const bool presolve_to_empty)2200 void reportPresolveReductions(const HighsOptions& options, const HighsLp& lp,
2201 const bool presolve_to_empty) {
2202 int num_col_from = lp.numCol_;
2203 int num_row_from = lp.numRow_;
2204 int num_els_from = lp.Astart_[num_col_from];
2205 int num_col_to;
2206 int num_row_to;
2207 int num_els_to;
2208 std::string message;
2209 if (presolve_to_empty) {
2210 num_col_to = 0;
2211 num_row_to = 0;
2212 num_els_to = 0;
2213 message = "- Reduced to empty";
2214 } else {
2215 num_col_to = num_col_from;
2216 num_row_to = num_row_from;
2217 num_els_to = num_els_from;
2218 message = "- Not reduced";
2219 }
2220 HighsPrintMessage(options.logfile, options.message_level, ML_ALWAYS,
2221 "Presolve : Reductions: rows %d(-%d); columns %d(-%d); "
2222 "elements %d(-%d) %s\n",
2223 num_row_to, (num_row_from - num_row_to), num_col_to,
2224 (num_col_from - num_col_to), num_els_to,
2225 (num_els_from - num_els_to), message.c_str());
2226 }
2227
isLessInfeasibleDSECandidate(const HighsOptions & options,const HighsLp & lp)2228 bool isLessInfeasibleDSECandidate(const HighsOptions& options,
2229 const HighsLp& lp) {
2230 int max_col_num_en = -1;
2231 const int max_allowed_col_num_en = 24;
2232 const int max_assess_col_num_en = std::max(9, max_allowed_col_num_en);
2233 const int max_average_col_num_en = 6;
2234 vector<int> col_length_k;
2235 col_length_k.resize(1 + max_assess_col_num_en, 0);
2236 bool LiDSE_candidate = true;
2237 bool all_unit_nonzeros = true;
2238 for (int col = 0; col < lp.numCol_; col++) {
2239 // Check limit on number of entries in the column has not been breached
2240 int col_num_en = lp.Astart_[col + 1] - lp.Astart_[col];
2241 max_col_num_en = std::max(col_num_en, max_col_num_en);
2242 if (col_num_en > max_assess_col_num_en) {
2243 #ifdef HiGHSDEV
2244 if (LiDSE_candidate)
2245 printf("Column %d has %d > %d entries so LP is not LiDSE candidate\n",
2246 col, col_num_en, max_allowed_col_num_en);
2247 LiDSE_candidate = false;
2248 #else
2249 LiDSE_candidate = false;
2250 return LiDSE_candidate;
2251 #endif
2252 } else {
2253 col_length_k[col_num_en]++;
2254 }
2255 for (int en = lp.Astart_[col]; en < lp.Astart_[col + 1]; en++) {
2256 double value = lp.Avalue_[en];
2257 // All nonzeros must be +1 or -1
2258 if (fabs(value) != 1) {
2259 all_unit_nonzeros = false;
2260 #ifdef HiGHSDEV
2261 if (LiDSE_candidate)
2262 printf(
2263 "Column %d has entry %d with value %g so LP is not LiDSE "
2264 "candidate\n",
2265 col, en - lp.Astart_[col], value);
2266 LiDSE_candidate = false;
2267 #else
2268 LiDSE_candidate = false;
2269 return LiDSE_candidate;
2270 #endif
2271 }
2272 }
2273 }
2274 #ifdef HiGHSDEV
2275 /*
2276 printf("LP has\n");
2277 int to_num_en = std::min(max_assess_col_num_en, max_col_num_en);
2278 for (int col_num_en = 0; col_num_en < to_num_en+1; col_num_en++)
2279 printf("%7d columns of count %1d\n", col_length_k[col_num_en], col_num_en);
2280 */
2281 #endif
2282 double average_col_num_en = lp.Astart_[lp.numCol_];
2283 average_col_num_en = average_col_num_en / lp.numCol_;
2284 LiDSE_candidate =
2285 LiDSE_candidate && average_col_num_en <= max_average_col_num_en;
2286 std::string logic0 = "has";
2287 if (!all_unit_nonzeros) logic0 = "does not have";
2288 std::string logic1 = "is not";
2289 if (LiDSE_candidate) logic1 = "is";
2290 HighsLogMessage(
2291 options.logfile, HighsMessageType::INFO,
2292 "LP %s %s all |entries|=1; max column count = %d (limit %d); average "
2293 "column count = %0.2g (limit %d): So %s a candidate for LiDSE",
2294 lp.model_name_.c_str(), logic0.c_str(), max_col_num_en,
2295 max_allowed_col_num_en, average_col_num_en, max_average_col_num_en,
2296 logic1.c_str());
2297 #ifdef HiGHSDEV
2298 int int_average_col_num_en = average_col_num_en;
2299 printf("grep_count_distrib,%s,%d,%d,%d\n", lp.model_name_.c_str(),
2300 max_col_num_en, int_average_col_num_en, LiDSE_candidate);
2301 #endif
2302 return LiDSE_candidate;
2303 }
2304
convertToMinimization(HighsLp & lp)2305 void convertToMinimization(HighsLp& lp) {
2306 if (lp.sense_ != ObjSense::MINIMIZE) {
2307 for (int col = 0; col < lp.numCol_; col++)
2308 lp.colCost_[col] = -lp.colCost_[col];
2309 }
2310 }
2311
isEqualityProblem(const HighsLp & lp)2312 bool isEqualityProblem(const HighsLp& lp) {
2313 for (int row = 0; row < lp.numRow_; row++)
2314 if (lp.rowLower_[row] != lp.rowUpper_[row]) return false;
2315
2316 return true;
2317 }
2318
vectorProduct(const std::vector<double> & v1,const std::vector<double> & v2)2319 double vectorProduct(const std::vector<double>& v1,
2320 const std::vector<double>& v2) {
2321 assert(v1.size() == v2.size());
2322 double sum = 0;
2323 for (int i = 0; i < (int)v1.size(); i++) sum += v1[i] * v2[i];
2324 return sum;
2325 }
2326
calculateResidual(const HighsLp & lp,HighsSolution & solution,std::vector<double> & residual)2327 HighsStatus calculateResidual(const HighsLp& lp, HighsSolution& solution,
2328 std::vector<double>& residual) {
2329 HighsStatus status = calculateRowValues(lp, solution);
2330 if (status != HighsStatus::OK) return status;
2331
2332 residual.clear();
2333 residual.resize(lp.numRow_);
2334
2335 for (int row = 0; row < lp.numRow_; row++) {
2336 if (solution.row_value[row] < lp.rowLower_[row]) {
2337 residual[row] = lp.rowLower_[row] - solution.row_value[row];
2338 } else if (solution.row_value[row] > lp.rowUpper_[row]) {
2339 residual[row] = solution.row_value[row] - lp.rowUpper_[row];
2340 }
2341 }
2342
2343 return status;
2344 }
2345