Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Review of the MA27 interface #116

Merged
merged 6 commits into from
Dec 1, 2024
Merged

Review of the MA27 interface #116

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
d319c89
Reviewing the MA27 interface (in progress): fixed two bugs, which sol…
cvanaret Dec 1, 2024
4717794
MA27Solver: create the w vector (double workspace) once and resize it…
cvanaret Dec 1, 2024
93dab07
Model: removed unused fixed_hessian_sparsity flag + cleaned up
cvanaret Dec 1, 2024
7ff8c68
Merge branch 'main' into ma27
cvanaret Dec 1, 2024
d3a62b3
Cleaned up MA27 factorization function
cvanaret Dec 1, 2024
b70d50b
Fixed small bug that I introduced in assertion
cvanaret Dec 1, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 4 additions & 0 deletions bindings/AMPL/AMPLModel.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -249,6 +249,10 @@ namespace uno {
return this->single_upper_bounded_variables_collection;
}

const Vector<size_t>& AMPLModel::get_fixed_variables() const {
return this->fixed_variables;
}

const Collection<size_t>& AMPLModel::get_lower_bounded_variables() const {
return this->lower_bounded_variables_collection;
}
Expand Down
2 changes: 1 addition & 1 deletion bindings/AMPL/AMPLModel.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,7 @@ namespace uno {
[[nodiscard]] const SparseVector<size_t>& get_slacks() const override;
[[nodiscard]] const Collection<size_t>& get_single_lower_bounded_variables() const override;
[[nodiscard]] const Collection<size_t>& get_single_upper_bounded_variables() const override;
[[nodiscard]] const Vector<size_t>& get_fixed_variables() const override { return this->fixed_variables; }
[[nodiscard]] const Vector<size_t>& get_fixed_variables() const override;

[[nodiscard]] double constraint_lower_bound(size_t constraint_index) const override;
[[nodiscard]] double constraint_upper_bound(size_t constraint_index) const override;
Expand Down
2 changes: 0 additions & 2 deletions uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -91,8 +91,6 @@ namespace uno {
template <typename ElementType>
void SymmetricIndefiniteLinearSystem<ElementType>::factorize_matrix(const Model& /*model*/,
DirectSymmetricIndefiniteLinearSolver<size_t, ElementType>& linear_solver) {
// compute the symbolic factorization only when:
// the problem has a non-constant augmented system (ie is not an LP or a QP) or it is the first factorization
if (true || this->number_factorizations == 0) {
linear_solver.do_symbolic_factorization(this->matrix);
}
Expand Down
158 changes: 75 additions & 83 deletions uno/solvers/MA27/MA27Solver.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -104,26 +104,28 @@ namespace uno {
SUCCESS = 0, // Successful completion.
IDXOUTOFRANGE = 1, // ndex (in IRN or ICN) out of range. Action taken by subroutine is to ignore any such entries and continue (MA27A/AD and MA27B/BD entries). INFO(2) is set to the number of faulty entries. Details of the first ten are printed on unit ICNTL(2).
FALSEDEFINITENESS, // Pivots have different signs when factorizing a supposedly definite matrix (when the value of U in CNTL(1) is zero) (MA27B/BD entry only). INFO(2) is set to the number of sign changes. Note that this warning will overwrite an INFO(1)=1 warning. Details of the first ten are printed on unit ICNTL(2).
RANKDEFECT, // Matrix is rank deficient. In this case, a decomposition will still have been produced which will enable the subsequent solution of consistent equations (MA27B/BD entry only). INFO(2) will be set to the rank of the matrix. Note that this warning will overwrite an INFO(1)=1 or INFO(1)=2 warning.
RANK_DEFICIENT, // Matrix is rank deficient. In this case, a decomposition will still have been produced which will enable the subsequent solution of consistent equations (MA27B/BD entry only). INFO(2) will be set to the rank of the matrix. Note that this warning will overwrite an INFO(1)=1 or INFO(1)=2 warning.
};


MA27Solver::MA27Solver(size_t max_dimension, size_t max_number_nonzeros):
DirectSymmetricIndefiniteLinearSolver<size_t, double>(max_dimension),
nz_max(static_cast<int>(max_number_nonzeros)), n(static_cast<int>(max_dimension)), nnz(static_cast<int>(max_number_nonzeros)),
irn(max_number_nonzeros), icn(max_number_nonzeros), iw((2 * max_number_nonzeros + 3 * max_dimension + 1) * 6 / 5),
ikeep(3 * max_number_nonzeros), iw1(2 * max_dimension) {
iflag = 0;
// set the default values of the controlling parameters
irn(max_number_nonzeros), icn(max_number_nonzeros),
iw((2 * max_number_nonzeros + 3 * max_dimension + 1) * 6 / 5), // 20% more than 2*nnz + 3*n + 1
ikeep(3 * max_dimension), iw1(2 * max_dimension) {
// initialization: set the default values of the controlling parameters
MA27ID(icntl.data(), cntl.data());
// a suitable pivot order is to be chosen automatically
iflag = 0;
// suppress warning messages
icntl[eICNTL::LP] = 0;
icntl[eICNTL::MP] = 0;
icntl[eICNTL::LDIAG] = 0;
}

void MA27Solver::factorize(const SymmetricMatrix<size_t, double>& matrix) {
// // general factorization method: symbolic factorization and numerical factorization
// general factorization method: symbolic factorization and numerical factorization
do_symbolic_factorization(matrix);
do_numerical_factorization(matrix);
}
Expand All @@ -144,102 +146,61 @@ namespace uno {
MA27AD(&n, &nnz, /* size info */
irn.data(), icn.data(), /* matrix indices */
iw.data(), &liw, ikeep.data(), iw1.data(), /* solver workspace */
&nsteps, &iflag,
icntl.data(), cntl.data(),
info.data(), &ops);
&nsteps, &iflag, icntl.data(), cntl.data(), info.data(), &ops);

// resize the factor by at least INFO(5) (here, 50% more)
factor.resize(static_cast<size_t>(3 * info[eINFO::NRLNEC] / 2));

std::copy(matrix.data_pointer(), matrix.data_pointer() + matrix.number_nonzeros(), factor.begin());

assert(eIFLAG::SUCCESS == info[eINFO::IFLAG] && "MA27: the symbolic factorization failed");
if (eIFLAG::SUCCESS != info[eINFO::IFLAG]) {
assert(info[eINFO::IFLAG] == eIFLAG::SUCCESS && "MA27: the symbolic factorization failed");
if (info[eINFO::IFLAG] != eIFLAG::SUCCESS) {
WARNING << "MA27 has issued a warning: IFLAG = " << info[eINFO::IFLAG] << " additional info, IERROR = " << info[eINFO::IERROR] << '\n';
}
}

void MA27Solver::repeat_factorization_after_resizing([[maybe_unused]]const SymmetricMatrix<size_t, double>& matrix) {
if (eIFLAG::INSUFFICIENTINTEGER == info[eINFO::IFLAG]) {
INFO << "MA27: insufficient integer workspace, resizing and retrying. \n";
// increase the size of iw
iw.resize(static_cast<size_t>(info[eINFO::IERROR]));
}

if (eIFLAG::INSUFFICIENTREAL == info[eINFO::IFLAG]) {
INFO << "MA27: insufficient real workspace, resizing and retrying. \n";
// increase the size of factor
factor.resize(static_cast<size_t>(info[eINFO::IERROR]));
}

int la = static_cast<int>(factor.size());
int liw = static_cast<int>(iw.size());

MA27BD(&n, &nnz, irn.data(), icn.data(), factor.data(), &la, iw.data(), &liw,
ikeep.data(), &nsteps, &maxfrt, iw1.data(), icntl.data(), cntl.data(), info.data());

if (eIFLAG::INSUFFICIENTINTEGER == info[eINFO::IFLAG] || eIFLAG::INSUFFICIENTREAL == info[eINFO::IFLAG]) {
repeat_factorization_after_resizing(matrix);
}
}

void MA27Solver::do_numerical_factorization([[maybe_unused]]const SymmetricMatrix<size_t, double>& matrix) {
assert(matrix.dimension() <= iw1.capacity() && "MA27Solver: the dimension of the matrix is larger than the preallocated size");
assert(nnz == static_cast<int>(matrix.number_nonzeros()) && "MA27Solver: the numbers of nonzeros do not match");

// numerical factorization
int la = static_cast<int>(factor.size());
int liw = static_cast<int>(iw.size());
MA27BD(&n, &nnz, irn.data(), icn.data(), factor.data(), &la, iw.data(), &liw,
ikeep.data(), &nsteps, &maxfrt, iw1.data(), icntl.data(), cntl.data(), info.data());

if (eIFLAG::INSUFFICIENTINTEGER == info[eINFO::IFLAG] || eIFLAG::INSUFFICIENTREAL == info[eINFO::IFLAG]) {
repeat_factorization_after_resizing(matrix);
}
// initialize factor with the entries of the matrix. It will be modified by MA27BD
std::copy(matrix.data_pointer(), matrix.data_pointer() + matrix.number_nonzeros(), factor.begin());

switch (info[eINFO::IFLAG]) {
case NSTEPS:
WARNING << "MA27BD: Value of NSTEPS outside the range 1 ≤ NSTEPS ≤ N" << '\n';
break;
case PIVOTSIGN:
WARNING << "MA27BD: A change of sign of pivots has been detected when U was negative. Detected at pivot step " << info[eINFO::IERROR]
<< '\n';
break;
case SINGULAR:
DEBUG << "MA27BD: Matrix is singular. Singularity detected during pivot step " << info[eINFO::IERROR] << '\n';
break;
case NZOUTOFRANGE:
WARNING << "MA27BD: Value of NZ out of range. NZ < 0." << '\n';
break;
case NOUTOFRANGE:
WARNING << "MA27BD: Value of N out of range. N < 1." << '\n';
break;
case IDXOUTOFRANGE:
WARNING << "MA27BD: Index (in IRN or ICN) out of range. " << info[eINFO::IERROR] << " indices affected." << '\n';
break;
case FALSEDEFINITENESS:
WARNING << "MA27BD: Matrix was supposed to be definite, but pivots have different signs when factorizing. Detected "
<< info[eINFO::IERROR] << " sign changes." << '\n';
break;
case RANKDEFECT:
DEBUG << "MA27BD: Matrix is rank deficient. Rank: " << info[eINFO::IERROR] << " whereas dimension " << n << '\n';
break;
// numerical factorization
// may fail because of insufficient space. In this case, more memory is allocated and the factorization tried again
bool factorization_done = false;
while (not factorization_done) {
int la = static_cast<int>(factor.size());
int liw = static_cast<int>(iw.size());
MA27BD(&n, &nnz, irn.data(), icn.data(), factor.data(), &la, iw.data(), &liw, ikeep.data(), &nsteps, &maxfrt, iw1.data(), icntl.data(),
cntl.data(), info.data());
factorization_done = true;

if (info[eINFO::IFLAG] == eIFLAG::INSUFFICIENTINTEGER) {
INFO << "MA27: insufficient integer workspace, resizing and retrying. \n";
// increase the size of iw
iw.resize(static_cast<size_t>(info[eINFO::IERROR]));
factorization_done = false;
}
if (info[eINFO::IFLAG] == eIFLAG::INSUFFICIENTREAL) {
INFO << "MA27: insufficient real workspace, resizing and retrying. \n";
// increase the size of factor
factor.resize(static_cast<size_t>(info[eINFO::IERROR]));
factorization_done = false;
}
}

this->w.resize(static_cast<size_t>(maxfrt));
this->check_factorization_status();
}

void MA27Solver::solve_indefinite_system([[maybe_unused]]const SymmetricMatrix<size_t, double>& matrix, const Vector<double>& rhs,
Vector<double>& result) {
// solve
std::vector<double> w(maxfrt); // double workspace
void MA27Solver::solve_indefinite_system(const SymmetricMatrix<size_t, double>& /*matrix*/, const Vector<double>& rhs, Vector<double>& result) {
int la = static_cast<int>(factor.size());
int liw = static_cast<int>(iw.size());

result = rhs;

MA27CD(&n, factor.data(), &la, iw.data(), &liw, w.data(), &maxfrt, result.data(), iw1.data(), &nsteps, icntl.data(), info.data());

assert(info[eINFO::IFLAG] == eIFLAG::SUCCESS && "MA27: the solution failed");
if (eIFLAG::SUCCESS != info[eINFO::IFLAG]) {
assert(info[eINFO::IFLAG] == eIFLAG::SUCCESS && "MA27: the linear solve failed");
if (info[eINFO::IFLAG] != eIFLAG::SUCCESS) {
WARNING << "MA27 has issued a warning: IFLAG = " << info[eINFO::IFLAG] << " additional info, IERROR = " << info[eINFO::IERROR] << '\n';
}
}
Expand All @@ -259,11 +220,11 @@ namespace uno {
}

bool MA27Solver::matrix_is_singular() const {
return (info[eINFO::IFLAG] == eIFLAG::SINGULAR || info[eINFO::IFLAG] == eIFLAG::RANKDEFECT);
return (info[eINFO::IFLAG] == eIFLAG::SINGULAR || info[eINFO::IFLAG] == eIFLAG::RANK_DEFICIENT);
}

size_t MA27Solver::rank() const {
return info[eINFO::IFLAG] == eIFLAG::RANKDEFECT ? static_cast<size_t>(info[eINFO::IERROR]) : static_cast<size_t>(n);
return (info[eINFO::IFLAG] == eIFLAG::RANK_DEFICIENT) ? static_cast<size_t>(info[eINFO::IERROR]) : static_cast<size_t>(n);
}

void MA27Solver::save_matrix_to_local_format(const SymmetricMatrix<size_t, double>& matrix) {
Expand All @@ -272,10 +233,41 @@ namespace uno {
icn.clear();
factor.clear();
constexpr auto fortran_shift = 1;
for (const auto[row_index, column_index, element]: matrix) {
for (const auto [row_index, column_index, element]: matrix) {
irn.emplace_back(static_cast<int>(row_index + fortran_shift));
icn.emplace_back(static_cast<int>(column_index + fortran_shift));
factor.emplace_back(element);
}
}

void MA27Solver::check_factorization_status() {
switch (info[eINFO::IFLAG]) {
case NSTEPS:
WARNING << "MA27BD: Value of NSTEPS outside the range 1 ≤ NSTEPS ≤ N" << '\n';
break;
case PIVOTSIGN:
WARNING << "MA27BD: A change of sign of pivots has been detected when U was negative. Detected at pivot step " << info[eINFO::IERROR]
<< '\n';
break;
case SINGULAR:
DEBUG << "MA27BD: Matrix is singular. Singularity detected during pivot step " << info[eINFO::IERROR] << '\n';
break;
case NZOUTOFRANGE:
WARNING << "MA27BD: Value of NZ out of range. NZ < 0." << '\n';
break;
case NOUTOFRANGE:
WARNING << "MA27BD: Value of N out of range. N < 1." << '\n';
break;
case IDXOUTOFRANGE:
WARNING << "MA27BD: Index (in IRN or ICN) out of range. " << info[eINFO::IERROR] << " indices affected." << '\n';
break;
case FALSEDEFINITENESS:
WARNING << "MA27BD: Matrix was supposed to be definite, but pivots have different signs when factorizing. Detected "
<< info[eINFO::IERROR] << " sign changes." << '\n';
break;
case RANK_DEFICIENT:
DEBUG << "MA27BD: Matrix is rank deficient. Rank: " << info[eINFO::IERROR] << " whereas dimension " << n << '\n';
break;
}
}
} // namespace
9 changes: 5 additions & 4 deletions uno/solvers/MA27/MA27Solver.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -37,8 +37,8 @@ class MA27Solver
int nz_max{}; // maximal number of nonzeros entries
int n{}; // dimension of current factorisation (maximal value here is <= max_dimension)
int nnz{}; // number of nonzeros of current factorisation
std::array<int,30> icntl{}; // integer array of length 30; integer control values
std::array<double,5> cntl{}; // double array of length 5; double control values
std::array<int, 30> icntl{}; // integer array of length 30; integer control values
std::array<double, 5> cntl{}; // double array of length 5; double control values

std::vector<int> irn{}; // row index of input
std::vector<int> icn{}; // col index of input
Expand All @@ -48,16 +48,17 @@ class MA27Solver
std::vector<int> iw1{}; // integer workspace array of length n
int nsteps{}; // integer, to be set by ma27
int iflag{}; // integer; 0 if pivot order chosen automatically; 1 if the pivot order set by ikeep
std::array<int,20> info{}; // integer array of length 20
std::array<int, 20> info{}; // integer array of length 20
double ops{}; // double, operations count

std::vector<double> factor{}; // data array of length la;
int maxfrt{}; // integer, to be set by ma27
std::vector<double> w{}; // double workspace


// bool use_iterative_refinement{false}; // Not sure how to do this with ma27
void save_matrix_to_local_format(const SymmetricMatrix<size_t, double>& matrix);
void repeat_factorization_after_resizing(const SymmetricMatrix<size_t, double> &matrix);
void check_factorization_status();
};

} // namespace uno
Expand Down
Loading