diff --git a/FEATURES.md b/FEATURES.md index bc7bf6d82c..eb17f2ed97 100644 --- a/FEATURES.md +++ b/FEATURES.md @@ -16,3 +16,4 @@ Fixed bug in presolve when pointers stored in HighsMatrixSlice get invalidated w Primal and dual residual tolerances - applied following IPM or PDLP solution - now documented as options +Highs::getCols (Highs::getRows) now runs in linear time if the internal constraint matrix is stored column-wise (row-wise). Added ensureColwise/Rowwise to the Highs class, the C API and highspy so that users can set the internal constraint matrix storage orientation diff --git a/check/Avgas.cpp b/check/Avgas.cpp index fb64f348fd..96db2cd0ba 100644 --- a/check/Avgas.cpp +++ b/check/Avgas.cpp @@ -15,324 +15,231 @@ #include #include // For printf +#include "lp_data/HConst.h" + const bool dev_run = false; -void Avgas::row(HighsInt row, HighsInt& num_row, HighsInt& num_row_nz, - std::vector& rowLower, std::vector& rowUpper, - std::vector& ARstart, std::vector& ARindex, - std::vector& ARvalue) { - rowLower.resize(num_row + 1); - rowUpper.resize(num_row + 1); - ARstart.resize(num_row + 1); - ARstart[num_row] = num_row_nz; +void Avgas::addRow(const HighsInt row, HighsInt& num_row, HighsInt& num_row_nz, + std::vector& rowLower, std::vector& rowUpper, + std::vector& ARstart, + std::vector& ARindex, + std::vector& ARvalue) { + double lower; + double upper; + std::vector index; + std::vector value; + getRow(row, lower, upper, index, value); + rowLower.push_back(lower); + rowUpper.push_back(upper); + HighsInt num_nz = index.size(); + ARstart.push_back(num_row_nz); + assert(HighsInt(ARstart.size()) == num_row + 1); + for (HighsInt iEl = 0; iEl < num_nz; iEl++) { + ARvalue.push_back(value[iEl]); + ARindex.push_back(index[iEl]); + } + num_row++; + num_row_nz += num_nz; +} + +void Avgas::getRow(const HighsInt row, double& lower, double& upper, + std::vector& index, std::vector& value) { + index.clear(); + value.clear(); + upper = kHighsInf; if (row == 0) { - rowLower[num_row] = -1; - rowUpper[num_row] = 1e31; - HighsInt num_new_nz = 2; - ARindex.resize(num_row_nz + num_new_nz); - ARvalue.resize(num_row_nz + num_new_nz); - ARindex[num_row_nz] = 0; - ARvalue[num_row_nz] = -1; - num_row_nz++; - ARindex[num_row_nz] = 1; - ARvalue[num_row_nz] = -1; - num_row_nz++; + lower = -1; + index.push_back(0); + value.push_back(-1); + index.push_back(1); + value.push_back(-1); } else if (row == 1) { - rowLower[num_row] = -1; - rowUpper[num_row] = 1e31; - HighsInt num_new_nz = 2; - ARindex.resize(num_row_nz + num_new_nz); - ARvalue.resize(num_row_nz + num_new_nz); - ARindex[num_row_nz] = 2; - ARvalue[num_row_nz] = -1; - num_row_nz++; - ARindex[num_row_nz] = 3; - ARvalue[num_row_nz] = -1; - num_row_nz++; + lower = -1; + index.push_back(2); + value.push_back(-1); + index.push_back(3); + value.push_back(-1); } else if (row == 2) { - rowLower[num_row] = -1; - rowUpper[num_row] = 1e31; - HighsInt num_new_nz = 2; - ARindex.resize(num_row_nz + num_new_nz); - ARvalue.resize(num_row_nz + num_new_nz); - ARindex[num_row_nz] = 4; - ARvalue[num_row_nz] = -1; - num_row_nz++; - ARindex[num_row_nz] = 5; - ARvalue[num_row_nz] = -1; - num_row_nz++; + lower = -1; + index.push_back(4); + value.push_back(-1); + index.push_back(5); + value.push_back(-1); } else if (row == 3) { - rowLower[num_row] = -1; - rowUpper[num_row] = 1e31; - HighsInt num_new_nz = 2; - ARindex.resize(num_row_nz + num_new_nz); - ARvalue.resize(num_row_nz + num_new_nz); - ARindex[num_row_nz] = 6; - ARvalue[num_row_nz] = -1; - num_row_nz++; - ARindex[num_row_nz] = 7; - ARvalue[num_row_nz] = -1; - num_row_nz++; + lower = -1; + index.push_back(6); + value.push_back(-1); + index.push_back(7); + value.push_back(-1); } else if (row == 4) { - rowLower[num_row] = -2; - rowUpper[num_row] = 1e31; - HighsInt num_new_nz = 4; - ARindex.resize(num_row_nz + num_new_nz); - ARvalue.resize(num_row_nz + num_new_nz); - ARindex[num_row_nz] = 0; - ARvalue[num_row_nz] = -1; - num_row_nz++; - ARindex[num_row_nz] = 2; - ARvalue[num_row_nz] = -1; - num_row_nz++; - ARindex[num_row_nz] = 4; - ARvalue[num_row_nz] = -1; - num_row_nz++; - ARindex[num_row_nz] = 6; - ARvalue[num_row_nz] = -1; - num_row_nz++; + lower = -2; + index.push_back(0); + value.push_back(-1); + index.push_back(2); + value.push_back(-1); + index.push_back(4); + value.push_back(-1); + index.push_back(6); + value.push_back(-1); } else if (row == 5) { - rowLower[num_row] = -2; - rowUpper[num_row] = 1e31; - HighsInt num_new_nz = 4; - ARindex.resize(num_row_nz + num_new_nz); - ARvalue.resize(num_row_nz + num_new_nz); - ARindex[num_row_nz] = 1; - ARvalue[num_row_nz] = -1; - num_row_nz++; - ARindex[num_row_nz] = 3; - ARvalue[num_row_nz] = -1; - num_row_nz++; - ARindex[num_row_nz] = 5; - ARvalue[num_row_nz] = -1; - num_row_nz++; - ARindex[num_row_nz] = 7; - ARvalue[num_row_nz] = -1; - num_row_nz++; + lower = -2; + index.push_back(1); + value.push_back(-1); + index.push_back(3); + value.push_back(-1); + index.push_back(5); + value.push_back(-1); + index.push_back(7); + value.push_back(-1); } else if (row == 6) { - rowLower[num_row] = 0; - rowUpper[num_row] = 1e31; - HighsInt num_new_nz = 3; - ARindex.resize(num_row_nz + num_new_nz); - ARvalue.resize(num_row_nz + num_new_nz); - ARindex[num_row_nz] = 0; - ARvalue[num_row_nz] = 2; - num_row_nz++; - ARindex[num_row_nz] = 2; - ARvalue[num_row_nz] = 1; - num_row_nz++; - ARindex[num_row_nz] = 6; - ARvalue[num_row_nz] = -1; - num_row_nz++; + lower = 0; + index.push_back(0); + value.push_back(2); + index.push_back(2); + value.push_back(1); + index.push_back(6); + value.push_back(-1); } else if (row == 7) { - rowLower[num_row] = 0; - rowUpper[num_row] = 1e31; - HighsInt num_new_nz = 4; - ARindex.resize(num_row_nz + num_new_nz); - ARvalue.resize(num_row_nz + num_new_nz); - ARindex[num_row_nz] = 0; - ARvalue[num_row_nz] = 5; - num_row_nz++; - ARindex[num_row_nz] = 2; - ARvalue[num_row_nz] = 3; - num_row_nz++; - ARindex[num_row_nz] = 4; - ARvalue[num_row_nz] = -3; - num_row_nz++; - ARindex[num_row_nz] = 6; - ARvalue[num_row_nz] = -1; - num_row_nz++; + lower = 0; + index.push_back(0); + value.push_back(5); + index.push_back(2); + value.push_back(3); + index.push_back(4); + value.push_back(-3); + index.push_back(6); + value.push_back(-1); } else if (row == 8) { - rowLower[num_row] = 0; - rowUpper[num_row] = 1e31; - HighsInt num_new_nz = 4; - ARindex.resize(num_row_nz + num_new_nz); - ARvalue.resize(num_row_nz + num_new_nz); - ARindex[num_row_nz] = 1; - ARvalue[num_row_nz] = 1; - num_row_nz++; - ARindex[num_row_nz] = 3; - ARvalue[num_row_nz] = -1; - num_row_nz++; - ARindex[num_row_nz] = 5; - ARvalue[num_row_nz] = -3; - num_row_nz++; - ARindex[num_row_nz] = 7; - ARvalue[num_row_nz] = -5; - num_row_nz++; + lower = 0; + index.push_back(1); + value.push_back(1); + index.push_back(3); + value.push_back(-1); + index.push_back(5); + value.push_back(-3); + index.push_back(7); + value.push_back(-5); } else if (row == 9) { - rowLower[num_row] = 0; - rowUpper[num_row] = 1e31; - HighsInt num_new_nz = 3; - ARindex.resize(num_row_nz + num_new_nz); - ARvalue.resize(num_row_nz + num_new_nz); - ARindex[num_row_nz] = 1; - ARvalue[num_row_nz] = 1; - num_row_nz++; - ARindex[num_row_nz] = 5; - ARvalue[num_row_nz] = -3; - num_row_nz++; - ARindex[num_row_nz] = 7; - ARvalue[num_row_nz] = -2; - num_row_nz++; + lower = 0; + index.push_back(1); + value.push_back(1); + index.push_back(5); + value.push_back(-3); + index.push_back(7); + value.push_back(-2); } else { - if (dev_run) printf("Avgas: row %" HIGHSINT_FORMAT " out of range\n", row); + if (dev_run) printf("Avgas: row %d out of range\n", HighsInt(row)); } - num_row++; } -void Avgas::col(HighsInt col, HighsInt& num_col, HighsInt& num_col_nz, - std::vector& colCost, std::vector& colLower, - std::vector& colUpper, std::vector& Astart, - std::vector& Aindex, std::vector& Avalue) { - colCost.resize(num_col + 1); - colLower.resize(num_col + 1); - colUpper.resize(num_col + 1); - Astart.resize(num_col + 1); - Astart[num_col] = num_col_nz; - HighsInt num_new_nz = 4; +void Avgas::addCol(HighsInt col, HighsInt& num_col, HighsInt& num_col_nz, + std::vector& colCost, std::vector& colLower, + std::vector& colUpper, std::vector& Astart, + std::vector& Aindex, std::vector& Avalue) { + double cost; + double lower; + double upper; + std::vector index; + std::vector value; + getCol(col, cost, lower, upper, index, value); + colCost.push_back(cost); + colLower.push_back(lower); + colUpper.push_back(upper); + HighsInt num_nz = index.size(); + Astart.push_back(num_col_nz); + assert(HighsInt(Astart.size()) == num_col + 1); + for (HighsInt iEl = 0; iEl < num_nz; iEl++) { + Avalue.push_back(value[iEl]); + Aindex.push_back(index[iEl]); + } + num_col++; + num_col_nz += num_nz; +} + +void Avgas::getCol(const HighsInt col, double& cost, double& lower, + double& upper, std::vector& index, + std::vector& value) { + index.clear(); + value.clear(); + lower = 0; + upper = 1; if (col == 0) { - colCost[num_col] = 0; - colLower[num_col] = 0; - colUpper[num_col] = 1; - Aindex.resize(num_col_nz + num_new_nz); - Avalue.resize(num_col_nz + num_new_nz); - Aindex[num_col_nz] = 0; - Avalue[num_col_nz] = -1; - num_col_nz++; - Aindex[num_col_nz] = 4; - Avalue[num_col_nz] = -1; - num_col_nz++; - Aindex[num_col_nz] = 6; - Avalue[num_col_nz] = 2; - num_col_nz++; - Aindex[num_col_nz] = 7; - Avalue[num_col_nz] = 5; - num_col_nz++; + cost = 0; + index.push_back(0); + value.push_back(-1); + index.push_back(4); + value.push_back(-1); + index.push_back(6); + value.push_back(2); + index.push_back(7); + value.push_back(5); } else if (col == 1) { - colCost[num_col] = -2; - colLower[num_col] = 0; - colUpper[num_col] = 1; - Aindex.resize(num_col_nz + num_new_nz); - Avalue.resize(num_col_nz + num_new_nz); - Aindex[num_col_nz] = 0; - Avalue[num_col_nz] = -1; - num_col_nz++; - Aindex[num_col_nz] = 5; - Avalue[num_col_nz] = -1; - num_col_nz++; - Aindex[num_col_nz] = 8; - Avalue[num_col_nz] = 1; - num_col_nz++; - Aindex[num_col_nz] = 9; - Avalue[num_col_nz] = 1; - num_col_nz++; + cost = -2; + index.push_back(0); + value.push_back(-1); + index.push_back(5); + value.push_back(-1); + index.push_back(8); + value.push_back(1); + index.push_back(9); + value.push_back(1); } else if (col == 2) { - colCost[num_col] = -1; - colLower[num_col] = 0; - colUpper[num_col] = 1; - Aindex.resize(num_col_nz + num_new_nz); - Avalue.resize(num_col_nz + num_new_nz); - Aindex[num_col_nz] = 1; - Avalue[num_col_nz] = -1; - num_col_nz++; - Aindex[num_col_nz] = 4; - Avalue[num_col_nz] = -1; - num_col_nz++; - Aindex[num_col_nz] = 6; - Avalue[num_col_nz] = 1; - num_col_nz++; - Aindex[num_col_nz] = 7; - Avalue[num_col_nz] = 3; - num_col_nz++; + cost = -1; + index.push_back(1); + value.push_back(-1); + index.push_back(4); + value.push_back(-1); + index.push_back(6); + value.push_back(1); + index.push_back(7); + value.push_back(3); } else if (col == 3) { - num_new_nz = 3; - colCost[num_col] = -3; - colLower[num_col] = 0; - colUpper[num_col] = 1; - Aindex.resize(num_col_nz + num_new_nz); - Avalue.resize(num_col_nz + num_new_nz); - Aindex[num_col_nz] = 1; - Avalue[num_col_nz] = -1; - num_col_nz++; - Aindex[num_col_nz] = 5; - Avalue[num_col_nz] = -1; - num_col_nz++; - Aindex[num_col_nz] = 8; - Avalue[num_col_nz] = -1; - num_col_nz++; + cost = -3; + index.push_back(1); + value.push_back(-1); + index.push_back(5); + value.push_back(-1); + index.push_back(8); + value.push_back(-1); } else if (col == 4) { - num_new_nz = 3; - colCost[num_col] = -2; - colLower[num_col] = 0; - colUpper[num_col] = 1; - Aindex.resize(num_col_nz + num_new_nz); - Avalue.resize(num_col_nz + num_new_nz); - Aindex[num_col_nz] = 2; - Avalue[num_col_nz] = -1; - num_col_nz++; - Aindex[num_col_nz] = 4; - Avalue[num_col_nz] = -1; - num_col_nz++; - Aindex[num_col_nz] = 7; - Avalue[num_col_nz] = -3; - num_col_nz++; + cost = -2; + index.push_back(2); + value.push_back(-1); + index.push_back(4); + value.push_back(-1); + index.push_back(7); + value.push_back(-3); } else if (col == 5) { - colCost[num_col] = -4; - colLower[num_col] = 0; - colUpper[num_col] = 1; - Aindex.resize(num_col_nz + num_new_nz); - Avalue.resize(num_col_nz + num_new_nz); - Aindex[num_col_nz] = 2; - Avalue[num_col_nz] = -1; - num_col_nz++; - Aindex[num_col_nz] = 5; - Avalue[num_col_nz] = -1; - num_col_nz++; - Aindex[num_col_nz] = 8; - Avalue[num_col_nz] = -3; - num_col_nz++; - Aindex[num_col_nz] = 9; - Avalue[num_col_nz] = -3; - num_col_nz++; + cost = -4; + index.push_back(2); + value.push_back(-1); + index.push_back(5); + value.push_back(-1); + index.push_back(8); + value.push_back(-3); + index.push_back(9); + value.push_back(-3); } else if (col == 6) { - colCost[num_col] = -3; - colLower[num_col] = 0; - colUpper[num_col] = 1; - Aindex.resize(num_col_nz + num_new_nz); - Avalue.resize(num_col_nz + num_new_nz); - Aindex[num_col_nz] = 3; - Avalue[num_col_nz] = -1; - num_col_nz++; - Aindex[num_col_nz] = 4; - Avalue[num_col_nz] = -1; - num_col_nz++; - Aindex[num_col_nz] = 6; - Avalue[num_col_nz] = -1; - num_col_nz++; - Aindex[num_col_nz] = 7; - Avalue[num_col_nz] = -1; - num_col_nz++; + cost = -3; + index.push_back(3); + value.push_back(-1); + index.push_back(4); + value.push_back(-1); + index.push_back(6); + value.push_back(-1); + index.push_back(7); + value.push_back(-1); } else if (col == 7) { - colCost[num_col] = -5; - colLower[num_col] = 0; - colUpper[num_col] = 1; - Aindex.resize(num_col_nz + num_new_nz); - Avalue.resize(num_col_nz + num_new_nz); - Aindex[num_col_nz] = 3; - Avalue[num_col_nz] = -1; - num_col_nz++; - Aindex[num_col_nz] = 5; - Avalue[num_col_nz] = -1; - num_col_nz++; - Aindex[num_col_nz] = 8; - Avalue[num_col_nz] = -5; - num_col_nz++; - Aindex[num_col_nz] = 9; - Avalue[num_col_nz] = -2; - num_col_nz++; + cost = -5; + index.push_back(3); + value.push_back(-1); + index.push_back(5); + value.push_back(-1); + index.push_back(8); + value.push_back(-5); + index.push_back(9); + value.push_back(-2); } else { - if (dev_run) printf("Avgas: col %" HIGHSINT_FORMAT " out of range\n", col); + if (dev_run) printf("Avgas: col %d out of range\n", HighsInt(col)); } - num_col++; } diff --git a/check/Avgas.h b/check/Avgas.h index e34bbb1b61..6e7154cfc2 100644 --- a/check/Avgas.h +++ b/check/Avgas.h @@ -17,19 +17,28 @@ #include "util/HighsInt.h" +const HighsInt avgas_num_col = 8; +const HighsInt avgas_num_row = 10; + /** * @brief Utilities for tests with AVGAS */ class Avgas { public: - void row(HighsInt row, HighsInt& num_row, HighsInt& num_row_nz, - std::vector& rowLower, std::vector& rowUpper, - std::vector& ARstart, std::vector& ARindex, - std::vector& ARvalue); + void addRow(const HighsInt row, HighsInt& num_row, HighsInt& num_row_nz, + std::vector& rowLower, std::vector& rowUpper, + std::vector& ARstart, std::vector& ARindex, + std::vector& ARvalue); + + void getRow(const HighsInt row, double& lower, double& upper, + std::vector& index, std::vector& value); + + void addCol(const HighsInt col, HighsInt& num_col, HighsInt& num_col_nz, + std::vector& colCost, std::vector& colLower, + std::vector& colUpper, std::vector& Astart, + std::vector& Aindex, std::vector& Avalue); - void col(HighsInt col, HighsInt& num_col, HighsInt& num_col_nz, - std::vector& colCost, std::vector& colLower, - std::vector& colUpper, std::vector& Astart, - std::vector& Aindex, std::vector& Avalue); + void getCol(const HighsInt col, double& cost, double& lower, double& upper, + std::vector& index, std::vector& value); }; #endif /* SIMPLEX_AVGAS_H_ */ diff --git a/check/TestLpModification.cpp b/check/TestLpModification.cpp index ea55d16264..30bdac0f06 100644 --- a/check/TestLpModification.cpp +++ b/check/TestLpModification.cpp @@ -5,6 +5,7 @@ #include "catch.hpp" #include "lp_data/HighsLpUtils.h" #include "util/HighsRandom.h" +#include "util/HighsTimer.h" #include "util/HighsUtils.h" const bool dev_run = false; @@ -38,6 +39,11 @@ bool areLpRowEqual(const HighsInt num_row0, const double* rowLower0, bool areLpEqual(const HighsLp lp0, const HighsLp lp1, const double infinite_bound); +bool equalSparseVectors(const HighsInt dim, const HighsInt num_nz0, + const HighsInt* index0, const double* value0, + const HighsInt num_nz1, const HighsInt* index1, + const double* value1); + void testDeleteKeep(const HighsIndexCollection& index_collection); bool testAllDeleteKeep(HighsInt num_row); @@ -425,8 +431,6 @@ TEST_CASE("LP-modification", "[highs_data]") { // options.log_dev_level = kHighsLogDevLevelVerbose; Avgas avgas; - const HighsInt avgas_num_col = 8; - const HighsInt avgas_num_row = 10; HighsInt num_row = 0; HighsInt num_row_nz = 0; std::vector rowLower; @@ -436,8 +440,8 @@ TEST_CASE("LP-modification", "[highs_data]") { std::vector ARvalue; for (HighsInt row = 0; row < avgas_num_row; row++) { - avgas.row(row, num_row, num_row_nz, rowLower, rowUpper, ARstart, ARindex, - ARvalue); + avgas.addRow(row, num_row, num_row_nz, rowLower, rowUpper, ARstart, ARindex, + ARvalue); } HighsInt num_col = 0; @@ -449,8 +453,8 @@ TEST_CASE("LP-modification", "[highs_data]") { std::vector Aindex; std::vector Avalue; for (HighsInt col = 0; col < avgas_num_col; col++) { - avgas.col(col, num_col, num_col_nz, colCost, colLower, colUpper, Astart, - Aindex, Avalue); + avgas.addCol(col, num_col, num_col_nz, colCost, colLower, colUpper, Astart, + Aindex, Avalue); } HighsStatus return_status; @@ -1876,6 +1880,39 @@ TEST_CASE("mod-duplicate-indices", "[highs_data]") { REQUIRE(objective0 == -7.75); } +bool equalSparseVectors(const HighsInt dim, const HighsInt num_nz0, + const HighsInt* index0, const double* value0, + const HighsInt num_nz1, const HighsInt* index1, + const double* value1) { + if (num_nz0 != num_nz1) { + if (dev_run) printf("num_nz0 != num_nz1\n"); + return false; + } + std::vector full_vector; + full_vector.assign(dim, 0); + for (HighsInt iEl = 0; iEl < num_nz0; iEl++) + full_vector[index0[iEl]] = value0[iEl]; + for (HighsInt iEl = 0; iEl < num_nz1; iEl++) { + HighsInt iRow = index1[iEl]; + if (full_vector[iRow] != value1[iEl]) { + if (dev_run) + printf("vector0[%d] = %g <> %g = vector1[%d]\n", int(iRow), + full_vector[iRow], value1[iEl], int(iRow)); + return false; + } + + full_vector[iRow] = 0; + } + for (HighsInt iRow = 0; iRow < dim; iRow++) + if (full_vector[iRow]) { + if (dev_run) + printf("Full vector[%d] = %g, not zero\n", int(iRow), + full_vector[iRow]); + return false; + } + return true; +} + TEST_CASE("resize-integrality", "[highs_data]") { Highs highs; highs.setOptionValue("output_flag", dev_run); @@ -1954,3 +1991,286 @@ TEST_CASE("zero-matrix-entries", "[highs_data]") { lp.a_matrix_.value_ = {1, 0, 0, 1}; REQUIRE(highs.passModel(lp) == HighsStatus::kOk); } + +void testAvgasGetRow(Highs& h) { + Avgas avgas; + double cost; + double lower; + double upper; + std::vector index; + std::vector value; + HighsInt get_num; + HighsInt lp_nnz; + std::vector lp_cost(1); + std::vector lp_lower(1); + std::vector lp_upper(1); + std::vector lp_start(1); + std::vector lp_index(avgas_num_col); + std::vector lp_value(avgas_num_col); + std::vector set(1); + std::vector mask(avgas_num_row); + for (HighsInt row = 0; row < avgas_num_row; row++) { + avgas.getRow(row, lower, upper, index, value); + HighsInt avgas_nnz = index.size(); + h.getRows(row, row, get_num, lp_lower.data(), lp_upper.data(), lp_nnz, + lp_start.data(), lp_index.data(), lp_value.data()); + REQUIRE(lp_lower[0] == lower); + REQUIRE(lp_upper[0] == upper); + REQUIRE(equalSparseVectors(avgas_num_col, avgas_nnz, index.data(), + value.data(), lp_nnz, lp_index.data(), + lp_value.data())); + } + HighsInt from_row = 2; + HighsInt to_row = 5; + HighsInt num_row = to_row - from_row + 1; + lp_lower.resize(num_row); + lp_upper.resize(num_row); + lp_start.resize(num_row); + lp_index.resize(num_row * avgas_num_col); + lp_value.resize(num_row * avgas_num_col); + h.getRows(from_row, to_row, get_num, lp_lower.data(), lp_upper.data(), lp_nnz, + lp_start.data(), lp_index.data(), lp_value.data()); + REQUIRE(get_num == num_row); + for (HighsInt row = 0; row < num_row; row++) { + HighsInt avgas_row = from_row + row; + avgas.getRow(avgas_row, lower, upper, index, value); + HighsInt avgas_nnz = index.size(); + REQUIRE(lp_lower[row] == lower); + REQUIRE(lp_upper[row] == upper); + HighsInt from_el = lp_start[row]; + HighsInt lp_col_nnz = + row < num_row - 1 ? lp_start[row + 1] - from_el : lp_nnz - from_el; + REQUIRE(equalSparseVectors(avgas_num_col, avgas_nnz, index.data(), + value.data(), lp_col_nnz, &lp_index[from_el], + &lp_value[from_el])); + } + set = {1, 2, 3, 6, 7}; + num_row = set.size(); + lp_lower.resize(num_row); + lp_upper.resize(num_row); + lp_start.resize(num_row); + lp_index.resize(num_row * avgas_num_col); + lp_value.resize(num_row * avgas_num_col); + h.getRows(num_row, set.data(), get_num, lp_lower.data(), lp_upper.data(), + lp_nnz, lp_start.data(), lp_index.data(), lp_value.data()); + REQUIRE(get_num == num_row); + for (HighsInt row = 0; row < num_row; row++) { + HighsInt avgas_row = set[row]; + avgas.getRow(avgas_row, lower, upper, index, value); + HighsInt avgas_nnz = index.size(); + REQUIRE(lp_lower[row] == lower); + REQUIRE(lp_upper[row] == upper); + HighsInt from_el = lp_start[row]; + HighsInt lp_col_nnz = + row < num_row - 1 ? lp_start[row + 1] - from_el : lp_nnz - from_el; + REQUIRE(equalSparseVectors(avgas_num_col, avgas_nnz, index.data(), + value.data(), lp_col_nnz, &lp_index[from_el], + &lp_value[from_el])); + } + mask[0] = 1; + mask[1] = 1; + mask[4] = 1; + mask[6] = 1; + mask[7] = 1; + num_row = 5; + lp_lower.resize(num_row); + lp_upper.resize(num_row); + lp_start.resize(num_row); + lp_index.resize(num_row * avgas_num_col); + lp_value.resize(num_row * avgas_num_col); + h.getRows(mask.data(), get_num, lp_lower.data(), lp_upper.data(), lp_nnz, + lp_start.data(), lp_index.data(), lp_value.data()); + REQUIRE(get_num == num_row); + HighsInt row = 0; + for (HighsInt iRow = 0; iRow < avgas_num_row; iRow++) { + if (!mask[iRow]) continue; + HighsInt avgas_row = iRow; + avgas.getRow(avgas_row, lower, upper, index, value); + HighsInt avgas_nnz = index.size(); + REQUIRE(lp_lower[row] == lower); + REQUIRE(lp_upper[row] == upper); + HighsInt from_el = lp_start[row]; + HighsInt lp_col_nnz = + row < num_row - 1 ? lp_start[row + 1] - from_el : lp_nnz - from_el; + REQUIRE(equalSparseVectors(avgas_num_col, avgas_nnz, index.data(), + value.data(), lp_col_nnz, &lp_index[from_el], + &lp_value[from_el])); + row++; + } +} + +void testAvgasGetCol(Highs& h) { + Avgas avgas; + double cost; + double lower; + double upper; + std::vector index; + std::vector value; + HighsInt get_num; + HighsInt lp_nnz; + std::vector lp_cost(1); + std::vector lp_lower(1); + std::vector lp_upper(1); + std::vector lp_start(1); + std::vector lp_index(avgas_num_row); + std::vector lp_value(avgas_num_row); + std::vector set(1); + std::vector mask(avgas_num_col); + mask.assign(avgas_num_col, 0); + for (HighsInt col = 0; col < avgas_num_col; col++) { + avgas.getCol(col, cost, lower, upper, index, value); + HighsInt avgas_nnz = index.size(); + h.getCols(col, col, get_num, lp_cost.data(), lp_lower.data(), + lp_upper.data(), lp_nnz, lp_start.data(), lp_index.data(), + lp_value.data()); + REQUIRE(lp_cost[0] == cost); + REQUIRE(lp_lower[0] == lower); + REQUIRE(lp_upper[0] == upper); + REQUIRE(equalSparseVectors(avgas_num_row, avgas_nnz, index.data(), + value.data(), lp_nnz, lp_index.data(), + lp_value.data())); + set[0] = col; + h.getCols(1, set.data(), get_num, lp_cost.data(), lp_lower.data(), + lp_upper.data(), lp_nnz, lp_start.data(), lp_index.data(), + lp_value.data()); + REQUIRE(lp_cost[0] == cost); + REQUIRE(lp_lower[0] == lower); + REQUIRE(lp_upper[0] == upper); + REQUIRE(equalSparseVectors(avgas_num_row, avgas_nnz, index.data(), + value.data(), lp_nnz, lp_index.data(), + lp_value.data())); + mask[col] = 1; + h.getCols(mask.data(), get_num, lp_cost.data(), lp_lower.data(), + lp_upper.data(), lp_nnz, lp_start.data(), lp_index.data(), + lp_value.data()); + REQUIRE(lp_cost[0] == cost); + REQUIRE(lp_lower[0] == lower); + REQUIRE(lp_upper[0] == upper); + REQUIRE(equalSparseVectors(avgas_num_row, avgas_nnz, index.data(), + value.data(), lp_nnz, lp_index.data(), + lp_value.data())); + mask[col] = 0; + } + HighsInt from_col = 2; + HighsInt to_col = 5; + HighsInt num_col = to_col - from_col + 1; + lp_cost.resize(num_col); + lp_lower.resize(num_col); + lp_upper.resize(num_col); + lp_start.resize(num_col); + lp_index.resize(num_col * avgas_num_row); + lp_value.resize(num_col * avgas_num_row); + h.getCols(from_col, to_col, get_num, lp_cost.data(), lp_lower.data(), + lp_upper.data(), lp_nnz, lp_start.data(), lp_index.data(), + lp_value.data()); + REQUIRE(get_num == num_col); + for (HighsInt col = 0; col < num_col; col++) { + HighsInt avgas_col = from_col + col; + avgas.getCol(avgas_col, cost, lower, upper, index, value); + HighsInt avgas_nnz = index.size(); + REQUIRE(lp_cost[col] == cost); + REQUIRE(lp_lower[col] == lower); + REQUIRE(lp_upper[col] == upper); + HighsInt from_el = lp_start[col]; + HighsInt lp_row_nnz = + col < num_col - 1 ? lp_start[col + 1] - from_el : lp_nnz - from_el; + REQUIRE(equalSparseVectors(avgas_num_row, avgas_nnz, index.data(), + value.data(), lp_row_nnz, &lp_index[from_el], + &lp_value[from_el])); + } + set = {1, 2, 3, 6, 7}; + num_col = set.size(); + lp_cost.resize(num_col); + lp_lower.resize(num_col); + lp_upper.resize(num_col); + lp_start.resize(num_col); + lp_index.resize(num_col * avgas_num_row); + lp_value.resize(num_col * avgas_num_row); + h.getCols(num_col, set.data(), get_num, lp_cost.data(), lp_lower.data(), + lp_upper.data(), lp_nnz, lp_start.data(), lp_index.data(), + lp_value.data()); + REQUIRE(get_num == num_col); + for (HighsInt col = 0; col < num_col; col++) { + HighsInt avgas_col = set[col]; + avgas.getCol(avgas_col, cost, lower, upper, index, value); + HighsInt avgas_nnz = index.size(); + REQUIRE(lp_cost[col] == cost); + REQUIRE(lp_lower[col] == lower); + REQUIRE(lp_upper[col] == upper); + HighsInt from_el = lp_start[col]; + HighsInt lp_row_nnz = + col < num_col - 1 ? lp_start[col + 1] - from_el : lp_nnz - from_el; + REQUIRE(equalSparseVectors(avgas_num_row, avgas_nnz, index.data(), + value.data(), lp_row_nnz, &lp_index[from_el], + &lp_value[from_el])); + } + mask[0] = 1; + mask[1] = 1; + mask[4] = 1; + mask[6] = 1; + mask[7] = 1; + num_col = 5; + lp_cost.resize(num_col); + lp_lower.resize(num_col); + lp_upper.resize(num_col); + lp_start.resize(num_col); + lp_index.resize(num_col * avgas_num_row); + lp_value.resize(num_col * avgas_num_row); + h.getCols(mask.data(), get_num, lp_cost.data(), lp_lower.data(), + lp_upper.data(), lp_nnz, lp_start.data(), lp_index.data(), + lp_value.data()); + REQUIRE(get_num == num_col); + HighsInt col = 0; + for (HighsInt iCol = 0; iCol < avgas_num_col; iCol++) { + if (!mask[iCol]) continue; + HighsInt avgas_col = iCol; + avgas.getCol(avgas_col, cost, lower, upper, index, value); + HighsInt avgas_nnz = index.size(); + REQUIRE(lp_cost[col] == cost); + REQUIRE(lp_lower[col] == lower); + REQUIRE(lp_upper[col] == upper); + HighsInt from_el = lp_start[col]; + HighsInt lp_row_nnz = + col < num_col - 1 ? lp_start[col + 1] - from_el : lp_nnz - from_el; + REQUIRE(equalSparseVectors(avgas_num_row, avgas_nnz, index.data(), + value.data(), lp_row_nnz, &lp_index[from_el], + &lp_value[from_el])); + col++; + } +} + +TEST_CASE("row-wise-get-row-avgas", "[highs_data]") { + Avgas avgas; + const HighsInt avgas_num_col = 8; + const HighsInt avgas_num_row = 10; + + Highs h; + h.setOptionValue("output_flag", dev_run); + double cost; + double lower; + double upper; + std::vector index; + std::vector value; + + for (HighsInt col = 0; col < avgas_num_col; col++) { + avgas.getCol(col, cost, lower, upper, index, value); + REQUIRE(h.addCol(cost, lower, upper, 0, nullptr, nullptr) == + HighsStatus::kOk); + } + for (HighsInt row = 0; row < avgas_num_row; row++) { + avgas.getRow(row, lower, upper, index, value); + HighsInt avgas_nnz = index.size(); + REQUIRE(h.addRow(lower, upper, avgas_nnz, index.data(), value.data()) == + HighsStatus::kOk); + } + + // Test extraction of rows and columns, with rowwise and colwise + // internal storage + h.ensureRowwise(); + testAvgasGetRow(h); + testAvgasGetCol(h); + + h.ensureColwise(); + testAvgasGetRow(h); + testAvgasGetCol(h); +} diff --git a/check/TestLpOrientation.cpp b/check/TestLpOrientation.cpp index e94297b9bc..a43a6d5b20 100644 --- a/check/TestLpOrientation.cpp +++ b/check/TestLpOrientation.cpp @@ -10,8 +10,6 @@ const bool dev_run = false; // No commas in test case name. TEST_CASE("LP-orientation", "[lp_orientation]") { Avgas avgas; - const HighsInt avgas_num_col = 8; - const HighsInt avgas_num_row = 10; HighsInt num_row = 0; HighsInt num_row_nz = 0; vector rowLower; @@ -21,11 +19,10 @@ TEST_CASE("LP-orientation", "[lp_orientation]") { vector ARvalue; for (HighsInt row = 0; row < avgas_num_row; row++) - avgas.row(row, num_row, num_row_nz, rowLower, rowUpper, ARstart, ARindex, - ARvalue); + avgas.addRow(row, num_row, num_row_nz, rowLower, rowUpper, ARstart, ARindex, + ARvalue); - ARstart.resize(num_row + 1); - ARstart[num_row] = num_row_nz; + ARstart.push_back(num_row_nz); HighsInt num_col = 0; HighsInt num_col_nz = 0; @@ -36,10 +33,10 @@ TEST_CASE("LP-orientation", "[lp_orientation]") { vector Aindex; vector Avalue; for (HighsInt col = 0; col < avgas_num_col; col++) - avgas.col(col, num_col, num_col_nz, colCost, colLower, colUpper, Astart, - Aindex, Avalue); - Astart.resize(num_col + 1); - Astart[num_col] = num_col_nz; + avgas.addCol(col, num_col, num_col_nz, colCost, colLower, colUpper, Astart, + Aindex, Avalue); + Astart.push_back(num_col_nz); + assert(num_col_nz == num_row_nz); double optimal_objective_function_value = -7.75; @@ -102,20 +99,18 @@ TEST_CASE("LP-orientation", "[lp_orientation]") { highs.clearModel(); highs.addCols(num_col, colCost.data(), colLower.data(), colUpper.data(), 0, NULL, NULL, NULL); - vector one_row_Lower; - vector one_row_Upper; - vector one_row_start; + double one_row_lower; + double one_row_upper; vector one_row_index; vector one_row_value; + HighsInt one_row_num_nz; for (HighsInt row = 0; row < avgas_num_row; row++) { - HighsInt one_row_numnz = 0; - HighsInt one_row_numrow = 0; - avgas.row(row, one_row_numrow, one_row_numnz, one_row_Lower, one_row_Upper, - one_row_start, one_row_index, one_row_value); - REQUIRE(highs.addRows(1, one_row_Lower.data(), one_row_Upper.data(), - one_row_numnz, one_row_start.data(), - one_row_index.data(), - one_row_value.data()) == HighsStatus::kOk); + avgas.getRow(row, one_row_lower, one_row_upper, one_row_index, + one_row_value); + one_row_num_nz = one_row_index.size(); + REQUIRE(highs.addRow(one_row_lower, one_row_upper, one_row_num_nz, + one_row_index.data(), + one_row_value.data()) == HighsStatus::kOk); } highs.run(); REQUIRE(info.objective_function_value == optimal_objective_function_value); diff --git a/check/TestLpValidation.cpp b/check/TestLpValidation.cpp index 473fdeb681..f048d505ed 100644 --- a/check/TestLpValidation.cpp +++ b/check/TestLpValidation.cpp @@ -171,10 +171,9 @@ TEST_CASE("LP-validation", "[highs_data]") { vector ARstart; vector ARindex; vector ARvalue; - for (HighsInt row = 0; row < avgas_num_row; row++) { - avgas.row(row, num_row, num_row_nz, rowLower, rowUpper, ARstart, ARindex, - ARvalue); + avgas.addRow(row, num_row, num_row_nz, rowLower, rowUpper, ARstart, ARindex, + ARvalue); } HighsInt num_col = 0; @@ -186,8 +185,8 @@ TEST_CASE("LP-validation", "[highs_data]") { vector Aindex; vector Avalue; for (HighsInt col = 0; col < avgas_num_col; col++) { - avgas.col(col, num_col, num_col_nz, colCost, colLower, colUpper, Astart, - Aindex, Avalue); + avgas.addCol(col, num_col, num_col_nz, colCost, colLower, colUpper, Astart, + Aindex, Avalue); } return_status = assessLp(lp, options); diff --git a/docs/src/interfaces/python/example-py.md b/docs/src/interfaces/python/example-py.md index 7f026fb903..ddb72b24e8 100644 --- a/docs/src/interfaces/python/example-py.md +++ b/docs/src/interfaces/python/example-py.md @@ -208,6 +208,8 @@ print('Basis validity = ', h.basisValidityToString(info.basis_validity)) ## Modify model data + * `EnsureColwise` + * `EnsureRowwise` * `changeObjectiveSense` * `changeColCost` * `changeColBounds` diff --git a/src/Highs.h b/src/Highs.h index dfd6eab44d..137fc86f78 100644 --- a/src/Highs.h +++ b/src/Highs.h @@ -1048,6 +1048,16 @@ class Highs { const HighsInt* starts, const HighsInt* indices, const double* values); + HighsStatus ensureColwise() { + this->model_.lp_.ensureColwise(); + return HighsStatus::kOk; + } + + HighsStatus ensureRowwise() { + this->model_.lp_.ensureRowwise(); + return HighsStatus::kOk; + } + /** * @brief Delete multiple columns from the incumbent model given by an * interval [from_col, to_col] @@ -1306,6 +1316,14 @@ class Highs { return ekk_instance_.primal_phase1_dual_; } + /** + * @brief Development methods + */ + HighsInt defineClock(const char* name) { + return this->timer_.clock_def(name); + } + void writeAllClocks() { this->timer_.writeAllClocks(); } + // Start of deprecated methods std::string compilationDate() const { return "deprecated"; } @@ -1550,15 +1568,14 @@ class Highs { void deleteRowsInterface(HighsIndexCollection& index_collection); void getColsInterface(const HighsIndexCollection& index_collection, - HighsInt& num_col, double* col_cost, double* col_lower, - double* col_upper, HighsInt& num_nz, - HighsInt* col_matrix_start, HighsInt* col_matrix_index, - double* col_matrix_value); + HighsInt& num_col, double* cost, double* lower, + double* upper, HighsInt& num_nz, HighsInt* start, + HighsInt* index, double* value); void getRowsInterface(const HighsIndexCollection& index_collection, - HighsInt& num_row, double* row_lower, double* row_upper, - HighsInt& num_nz, HighsInt* row_matrix_start, - HighsInt* row_matrix_index, double* row_matrix_value); + HighsInt& num_row, double* lower, double* upper, + HighsInt& num_nz, HighsInt* start, HighsInt* index, + double* value); void getCoefficientInterface(const HighsInt ext_row, const HighsInt ext_col, double& value); diff --git a/src/highs_bindings.cpp b/src/highs_bindings.cpp index d4c7ad0ebe..2ea1c35817 100644 --- a/src/highs_bindings.cpp +++ b/src/highs_bindings.cpp @@ -1330,6 +1330,8 @@ PYBIND11_MODULE(_core, m, py::mod_gil_not_used()) { .def("addCols", &highs_addCols) .def("addVar", &highs_addVar) .def("addVars", &highs_addVars) + .def("ensureColwise", &Highs::ensureColwise) + .def("ensureRowwise", &Highs::ensureRowwise) .def("changeColsCost", &highs_changeColsCost) .def("changeColsBounds", &highs_changeColsBounds) .def("changeColsIntegrality", &highs_changeColsIntegrality) diff --git a/src/interfaces/highs_c_api.cpp b/src/interfaces/highs_c_api.cpp index af354daf02..e5fefce606 100644 --- a/src/interfaces/highs_c_api.cpp +++ b/src/interfaces/highs_c_api.cpp @@ -802,6 +802,14 @@ HighsInt Highs_addRows(void* highs, const HighsInt num_new_row, ->addRows(num_new_row, lower, upper, num_new_nz, starts, index, value); } +HighsInt Highs_ensureColwise(void* highs) { + return (HighsInt)((Highs*)highs)->ensureColwise(); +} + +HighsInt Highs_ensureRowwise(void* highs) { + return (HighsInt)((Highs*)highs)->ensureRowwise(); +} + HighsInt Highs_changeObjectiveSense(void* highs, const HighsInt sense) { ObjSense pass_sense = ObjSense::kMinimize; if (sense == (HighsInt)ObjSense::kMaximize) pass_sense = ObjSense::kMaximize; diff --git a/src/interfaces/highs_c_api.h b/src/interfaces/highs_c_api.h index 89f6fa4f79..02fe1930b4 100644 --- a/src/interfaces/highs_c_api.h +++ b/src/interfaces/highs_c_api.h @@ -1426,6 +1426,25 @@ HighsInt Highs_addRows(void* highs, const HighsInt num_new_row, const HighsInt num_new_nz, const HighsInt* starts, const HighsInt* index, const double* value); +/** + * Ensure that the constraint matrix of the incumbent model is stored + * column-wise. + * + * @param highs A pointer to the Highs instance. + * + * @returns A `kHighsStatus` constant indicating whether the call succeeded. + */ +HighsInt Highs_ensureColwise(void* highs); + +/** + * Ensure that the constraint matrix of the incumbent model is stored row-wise. + * + * @param highs A pointer to the Highs instance. + * + * @returns A `kHighsStatus` constant indicating whether the call succeeded. + */ +HighsInt Highs_ensureRowwise(void* highs); + /** * Change the objective sense of the model. * diff --git a/src/lp_data/HighsInterface.cpp b/src/lp_data/HighsInterface.cpp index c36dcf89c3..e522e46db9 100644 --- a/src/lp_data/HighsInterface.cpp +++ b/src/lp_data/HighsInterface.cpp @@ -666,181 +666,36 @@ void Highs::deleteRowsInterface(HighsIndexCollection& index_collection) { } void Highs::getColsInterface(const HighsIndexCollection& index_collection, - HighsInt& get_num_col, double* col_cost, - double* col_lower, double* col_upper, - HighsInt& get_num_nz, HighsInt* col_matrix_start, - HighsInt* col_matrix_index, - double* col_matrix_value) { - HighsLp& lp = model_.lp_; - // Ensure that the LP is column-wise - lp.ensureColwise(); - assert(ok(index_collection)); - HighsInt from_k; - HighsInt to_k; - limits(index_collection, from_k, to_k); - // Surely this is checked elsewhere - assert(0 <= from_k && to_k < lp.num_col_); - assert(from_k <= to_k); - HighsInt out_from_col; - HighsInt out_to_col; - HighsInt in_from_col; - HighsInt in_to_col = -1; - HighsInt current_set_entry = 0; - HighsInt col_dim = lp.num_col_; - get_num_col = 0; - get_num_nz = 0; - for (HighsInt k = from_k; k <= to_k; k++) { - updateOutInIndex(index_collection, out_from_col, out_to_col, in_from_col, - in_to_col, current_set_entry); - assert(out_to_col < col_dim); - assert(in_to_col < col_dim); - for (HighsInt iCol = out_from_col; iCol <= out_to_col; iCol++) { - if (col_cost != NULL) col_cost[get_num_col] = lp.col_cost_[iCol]; - if (col_lower != NULL) col_lower[get_num_col] = lp.col_lower_[iCol]; - if (col_upper != NULL) col_upper[get_num_col] = lp.col_upper_[iCol]; - if (col_matrix_start != NULL) - col_matrix_start[get_num_col] = get_num_nz + lp.a_matrix_.start_[iCol] - - lp.a_matrix_.start_[out_from_col]; - get_num_col++; - } - for (HighsInt iEl = lp.a_matrix_.start_[out_from_col]; - iEl < lp.a_matrix_.start_[out_to_col + 1]; iEl++) { - if (col_matrix_index != NULL) - col_matrix_index[get_num_nz] = lp.a_matrix_.index_[iEl]; - if (col_matrix_value != NULL) - col_matrix_value[get_num_nz] = lp.a_matrix_.value_[iEl]; - get_num_nz++; - } - if (out_to_col == col_dim - 1 || in_to_col == col_dim - 1) break; + HighsInt& num_col, double* cost, double* lower, + double* upper, HighsInt& num_nz, HighsInt* start, + HighsInt* index, double* value) { + const HighsLp& lp = model_.lp_; + if (lp.a_matrix_.isColwise()) { + getSubVectors(index_collection, lp.num_col_, lp.col_cost_.data(), + lp.col_lower_.data(), lp.col_upper_.data(), lp.a_matrix_, + num_col, cost, lower, upper, num_nz, start, index, value); + } else { + getSubVectorsTranspose(index_collection, lp.num_col_, lp.col_cost_.data(), + lp.col_lower_.data(), lp.col_upper_.data(), + lp.a_matrix_, num_col, cost, lower, upper, num_nz, + start, index, value); } } void Highs::getRowsInterface(const HighsIndexCollection& index_collection, - HighsInt& get_num_row, double* row_lower, - double* row_upper, HighsInt& get_num_nz, - HighsInt* row_matrix_start, - HighsInt* row_matrix_index, - double* row_matrix_value) { - HighsLp& lp = model_.lp_; - // Ensure that the LP is column-wise - lp.ensureColwise(); - assert(ok(index_collection)); - HighsInt from_k; - HighsInt to_k; - limits(index_collection, from_k, to_k); - // Surely this is checked elsewhere - assert(0 <= from_k && to_k < lp.num_row_); - assert(from_k <= to_k); - // "Out" means not in the set to be extracted - // "In" means in the set to be extracted - HighsInt out_from_row; - HighsInt out_to_row; - HighsInt in_from_row; - HighsInt in_to_row = -1; - HighsInt current_set_entry = 0; - HighsInt row_dim = lp.num_row_; - // Ensure that the LP is column-wise - lp.ensureColwise(); - // Set up a row mask so that entries to be got from the column-wise - // matrix can be identified and have their correct row index. - vector new_index; - new_index.resize(lp.num_row_); - - get_num_row = 0; - get_num_nz = 0; - if (!index_collection.is_mask_) { - out_to_row = -1; - current_set_entry = 0; - for (HighsInt k = from_k; k <= to_k; k++) { - updateOutInIndex(index_collection, in_from_row, in_to_row, out_from_row, - out_to_row, current_set_entry); - if (k == from_k) { - // Account for any initial rows not being extracted - for (HighsInt iRow = 0; iRow < in_from_row; iRow++) { - new_index[iRow] = -1; - } - } - for (HighsInt iRow = in_from_row; iRow <= in_to_row; iRow++) { - new_index[iRow] = get_num_row; - get_num_row++; - } - for (HighsInt iRow = out_from_row; iRow <= out_to_row; iRow++) { - new_index[iRow] = -1; - } - if (out_to_row >= row_dim - 1) break; - } + HighsInt& num_row, double* lower, double* upper, + HighsInt& num_nz, HighsInt* start, HighsInt* index, + double* value) { + const HighsLp& lp = model_.lp_; + if (lp.a_matrix_.isColwise()) { + getSubVectorsTranspose(index_collection, lp.num_row_, nullptr, + lp.row_lower_.data(), lp.row_upper_.data(), + lp.a_matrix_, num_row, nullptr, lower, upper, num_nz, + start, index, value); } else { - for (HighsInt iRow = 0; iRow < lp.num_row_; iRow++) { - if (index_collection.mask_[iRow]) { - new_index[iRow] = get_num_row; - get_num_row++; - } else { - new_index[iRow] = -1; - } - } - } - - // Bail out if no rows are to be extracted - if (get_num_row == 0) return; - - for (HighsInt iRow = 0; iRow < lp.num_row_; iRow++) { - HighsInt new_iRow = new_index[iRow]; - if (new_iRow >= 0) { - assert(new_iRow < get_num_row); - if (row_lower != NULL) row_lower[new_iRow] = lp.row_lower_[iRow]; - if (row_upper != NULL) row_upper[new_iRow] = lp.row_upper_[iRow]; - } - } - const bool extract_start = row_matrix_start != NULL; - const bool extract_index = row_matrix_index != NULL; - const bool extract_value = row_matrix_value != NULL; - const bool extract_matrix = extract_index || extract_value; - // Allocate an array of lengths for the row-wise matrix to be - // extracted: necessary even if just the number of nonzeros is - // required - vector row_matrix_length; - row_matrix_length.assign(get_num_row, 0); - // Identify the lengths of the rows in the row-wise matrix to be extracted - for (HighsInt col = 0; col < lp.num_col_; col++) { - for (HighsInt iEl = lp.a_matrix_.start_[col]; - iEl < lp.a_matrix_.start_[col + 1]; iEl++) { - HighsInt iRow = lp.a_matrix_.index_[iEl]; - HighsInt new_iRow = new_index[iRow]; - if (new_iRow >= 0) row_matrix_length[new_iRow]++; - } - } - if (!extract_start) { - // bail out if no matrix starts are to be extracted, but only after - // computing the number of nonzeros - for (HighsInt iRow = 0; iRow < get_num_row; iRow++) - get_num_nz += row_matrix_length[iRow]; - return; - } - // Allocate an array of lengths for the row-wise matrix to be extracted - row_matrix_start[0] = 0; - for (HighsInt iRow = 0; iRow < get_num_row - 1; iRow++) { - row_matrix_start[iRow + 1] = - row_matrix_start[iRow] + row_matrix_length[iRow]; - row_matrix_length[iRow] = row_matrix_start[iRow]; - } - HighsInt iRow = get_num_row - 1; - get_num_nz = row_matrix_start[iRow] + row_matrix_length[iRow]; - // Bail out if matrix indices and values are not required - if (!extract_matrix) return; - row_matrix_length[iRow] = row_matrix_start[iRow]; - // Fill the row-wise matrix with indices and values - for (HighsInt col = 0; col < lp.num_col_; col++) { - for (HighsInt iEl = lp.a_matrix_.start_[col]; - iEl < lp.a_matrix_.start_[col + 1]; iEl++) { - HighsInt iRow = lp.a_matrix_.index_[iEl]; - HighsInt new_iRow = new_index[iRow]; - if (new_iRow >= 0) { - HighsInt row_iEl = row_matrix_length[new_iRow]; - if (extract_index) row_matrix_index[row_iEl] = col; - if (extract_value) row_matrix_value[row_iEl] = lp.a_matrix_.value_[iEl]; - row_matrix_length[new_iRow]++; - } - } + getSubVectors(index_collection, lp.num_row_, nullptr, lp.row_lower_.data(), + lp.row_upper_.data(), lp.a_matrix_, num_row, nullptr, lower, + upper, num_nz, start, index, value); } } diff --git a/src/lp_data/HighsLpUtils.cpp b/src/lp_data/HighsLpUtils.cpp index 253d8a26f6..6b6042fc03 100644 --- a/src/lp_data/HighsLpUtils.cpp +++ b/src/lp_data/HighsLpUtils.cpp @@ -23,7 +23,6 @@ #include "util/HighsCDouble.h" #include "util/HighsMatrixUtils.h" #include "util/HighsSort.h" -#include "util/HighsTimer.h" using std::fabs; using std::max; @@ -3077,3 +3076,196 @@ void removeRowsOfCountOne(const HighsLogOptions& log_options, HighsLp& lp) { highsLogUser(log_options, HighsLogType::kWarning, "Removed %d rows of count 1\n", (int)num_row_count_1); } + +void getSubVectors(const HighsIndexCollection& index_collection, + const HighsInt data_dim, const double* data0, + const double* data1, const double* data2, + const HighsSparseMatrix& matrix, HighsInt& num_sub_vector, + double* sub_vector_data0, double* sub_vector_data1, + double* sub_vector_data2, HighsInt& sub_matrix_num_nz, + HighsInt* sub_matrix_start, HighsInt* sub_matrix_index, + double* sub_matrix_value) { + // Ensure that if there's no data0 then it's not required in the + // sub-vector + if (data0 == nullptr) assert(sub_vector_data0 == nullptr); + assert(ok(index_collection)); + HighsInt from_k; + HighsInt to_k; + limits(index_collection, from_k, to_k); + // Surely this is checked elsewhere + assert(0 <= from_k && to_k < data_dim); + assert(from_k <= to_k); + HighsInt out_from_vector; + HighsInt out_to_vector; + HighsInt in_from_vector; + HighsInt in_to_vector = -1; + HighsInt current_set_entry = 0; + num_sub_vector = 0; + sub_matrix_num_nz = 0; + for (HighsInt k = from_k; k <= to_k; k++) { + updateOutInIndex(index_collection, out_from_vector, out_to_vector, + in_from_vector, in_to_vector, current_set_entry); + assert(out_to_vector < data_dim); + assert(in_to_vector < data_dim); + for (HighsInt iVector = out_from_vector; iVector <= out_to_vector; + iVector++) { + if (sub_vector_data0 != nullptr) + sub_vector_data0[num_sub_vector] = data0[iVector]; + if (sub_vector_data1 != nullptr) + sub_vector_data1[num_sub_vector] = data1[iVector]; + if (sub_vector_data2 != nullptr) + sub_vector_data2[num_sub_vector] = data2[iVector]; + if (sub_matrix_start != nullptr) + sub_matrix_start[num_sub_vector] = sub_matrix_num_nz + + matrix.start_[iVector] - + matrix.start_[out_from_vector]; + num_sub_vector++; + } + for (HighsInt iEl = matrix.start_[out_from_vector]; + iEl < matrix.start_[out_to_vector + 1]; iEl++) { + if (sub_matrix_index != nullptr) + sub_matrix_index[sub_matrix_num_nz] = matrix.index_[iEl]; + if (sub_matrix_value != nullptr) + sub_matrix_value[sub_matrix_num_nz] = matrix.value_[iEl]; + sub_matrix_num_nz++; + } + if (out_to_vector == data_dim - 1 || in_to_vector == data_dim - 1) break; + } +} + +void getSubVectorsTranspose(const HighsIndexCollection& index_collection, + const HighsInt data_dim, const double* data0, + const double* data1, const double* data2, + const HighsSparseMatrix& matrix, + HighsInt& num_sub_vector, double* sub_vector_data0, + double* sub_vector_data1, double* sub_vector_data2, + HighsInt& sub_matrix_num_nz, + HighsInt* sub_matrix_start, + HighsInt* sub_matrix_index, + double* sub_matrix_value) { + // Ensure that if there's no data0 then it's not required in the + // sub-vector + if (data0 == nullptr) assert(sub_vector_data0 == nullptr); + assert(ok(index_collection)); + HighsInt from_k; + HighsInt to_k; + limits(index_collection, from_k, to_k); + // Surely this is checked elsewhere + assert(0 <= from_k && to_k < data_dim); + assert(from_k <= to_k); + // "Out" means not in the set to be extracted + // "In" means in the set to be extracted + HighsInt out_from_vector; + HighsInt out_to_vector; + HighsInt in_from_vector; + HighsInt in_to_vector = -1; + HighsInt current_set_entry = 0; + // Set up a mask so that entries to be got from the matrix can be + // identified and have their correct index. + vector new_index; + new_index.resize(data_dim); + + num_sub_vector = 0; + sub_matrix_num_nz = 0; + if (!index_collection.is_mask_) { + out_to_vector = -1; + current_set_entry = 0; + for (HighsInt k = from_k; k <= to_k; k++) { + updateOutInIndex(index_collection, in_from_vector, in_to_vector, + out_from_vector, out_to_vector, current_set_entry); + if (k == from_k) { + // Account for any initial vectors not being extracted + for (HighsInt iVector = 0; iVector < in_from_vector; iVector++) { + new_index[iVector] = -1; + } + } + for (HighsInt iVector = in_from_vector; iVector <= in_to_vector; + iVector++) { + new_index[iVector] = num_sub_vector; + num_sub_vector++; + } + for (HighsInt iVector = out_from_vector; iVector <= out_to_vector; + iVector++) { + new_index[iVector] = -1; + } + if (out_to_vector >= data_dim - 1) break; + } + } else { + for (HighsInt iVector = 0; iVector < data_dim; iVector++) { + if (index_collection.mask_[iVector]) { + new_index[iVector] = num_sub_vector; + num_sub_vector++; + } else { + new_index[iVector] = -1; + } + } + } + + // Bail out if no vectors are to be extracted + if (num_sub_vector == 0) return; + + for (HighsInt iVector = 0; iVector < data_dim; iVector++) { + HighsInt new_iVector = new_index[iVector]; + if (new_iVector >= 0) { + assert(new_iVector < num_sub_vector); + if (sub_vector_data0 != NULL) + sub_vector_data0[new_iVector] = data0[iVector]; + if (sub_vector_data1 != NULL) + sub_vector_data1[new_iVector] = data1[iVector]; + if (sub_vector_data2 != NULL) + sub_vector_data2[new_iVector] = data2[iVector]; + } + } + const bool extract_start = sub_matrix_start != NULL; + const bool extract_index = sub_matrix_index != NULL; + const bool extract_value = sub_matrix_value != NULL; + const bool extract_matrix = extract_index || extract_value; + // Allocate an array of lengths for the sub-matrix to be + // extracted: necessary even if just the number of nonzeros is + // required + vector sub_matrix_length; + sub_matrix_length.assign(num_sub_vector, 0); + // Identify the lengths of the vectors in the sub-matrix to be extracted + HighsInt num_vector = matrix.start_.size() - 1; + for (HighsInt vector = 0; vector < num_vector; vector++) { + for (HighsInt iEl = matrix.start_[vector]; iEl < matrix.start_[vector + 1]; + iEl++) { + HighsInt iVector = matrix.index_[iEl]; + HighsInt new_iVector = new_index[iVector]; + if (new_iVector >= 0) sub_matrix_length[new_iVector]++; + } + } + if (!extract_start) { + // bail out if no matrix starts are to be extracted, but only after + // computing the number of nonzeros + for (HighsInt iVector = 0; iVector < num_sub_vector; iVector++) + sub_matrix_num_nz += sub_matrix_length[iVector]; + return; + } + // Allocate an array of lengths for the sub-matrix to be extracted + sub_matrix_start[0] = 0; + for (HighsInt iVector = 0; iVector < num_sub_vector - 1; iVector++) { + sub_matrix_start[iVector + 1] = + sub_matrix_start[iVector] + sub_matrix_length[iVector]; + sub_matrix_length[iVector] = sub_matrix_start[iVector]; + } + HighsInt iVector = num_sub_vector - 1; + sub_matrix_num_nz = sub_matrix_start[iVector] + sub_matrix_length[iVector]; + // Bail out if matrix indices and values are not required + if (!extract_matrix) return; + sub_matrix_length[iVector] = sub_matrix_start[iVector]; + // Fill the row-wise matrix with indices and values + for (HighsInt vector = 0; vector < num_vector; vector++) { + for (HighsInt iEl = matrix.start_[vector]; iEl < matrix.start_[vector + 1]; + iEl++) { + HighsInt iVector = matrix.index_[iEl]; + HighsInt new_iVector = new_index[iVector]; + if (new_iVector >= 0) { + HighsInt row_iEl = sub_matrix_length[new_iVector]; + if (extract_index) sub_matrix_index[row_iEl] = vector; + if (extract_value) sub_matrix_value[row_iEl] = matrix.value_[iEl]; + sub_matrix_length[new_iVector]++; + } + } + } +} diff --git a/src/lp_data/HighsLpUtils.h b/src/lp_data/HighsLpUtils.h index c6cff7a8bb..867e92c12c 100644 --- a/src/lp_data/HighsLpUtils.h +++ b/src/lp_data/HighsLpUtils.h @@ -263,4 +263,36 @@ HighsLp withoutSemiVariables(const HighsLp& lp, HighsSolution& solution, void removeRowsOfCountOne(const HighsLogOptions& log_options, HighsLp& lp); +// Get subvectors from data structure of data0, data1, data2 and +// matrix, where the storage of the matrix is compatible with the +// vectors to be extracted from it +// +// Data to be extracted is given by sub_* being nullpointer +// +// * cost, lower and upper bounds for columns, and column-wise LP constraint +// matrix +// +// * lower and upper bounds for rows, and row-wise LP constraint +// * matrix. "cost" is nullptr, and so must be sub_vector_data0 +// +void getSubVectors(const HighsIndexCollection& index_collection, + const HighsInt data_dim, const double* data0, + const double* data1, const double* data2, + const HighsSparseMatrix& matrix, HighsInt& num_sub_vector, + double* sub_vector_data0, double* sub_vector_data1, + double* sub_vector_data2, HighsInt& sub_matrix_num_nz, + HighsInt* sub_matrix_start, HighsInt* sub_matrix_index, + double* sub_matrix_value); + +void getSubVectorsTranspose(const HighsIndexCollection& index_collection, + const HighsInt data_dim, const double* data0, + const double* data1, const double* data2, + const HighsSparseMatrix& matrix, + HighsInt& num_sub_vector, double* sub_vector_data0, + double* sub_vector_data1, double* sub_vector_data2, + HighsInt& sub_matrix_num_nz, + HighsInt* sub_matrix_start, + HighsInt* sub_matrix_index, + double* sub_matrix_value); + #endif // LP_DATA_HIGHSLPUTILS_H_ diff --git a/src/util/HighsTimer.h b/src/util/HighsTimer.h index 71f1c63223..a171866fe5 100644 --- a/src/util/HighsTimer.h +++ b/src/util/HighsTimer.h @@ -126,6 +126,16 @@ class HighsTimer { } } + /** + * @brief write all clocks + */ + void writeAllClocks() { + for (HighsInt i = 0; i < num_clock; i++) + if (clock_num_call[i]) + printf("Time %7.5f for %9d calls of clock %3d: %s\n", clock_time[i], + int(clock_num_call[i]), int(i), clock_names[i].c_str()); + } + /** * @brief Start a clock */