Skip to content

Commit

Permalink
Merge pull request #116 from cvanaret/ma27
Browse files Browse the repository at this point in the history 
(WIP) Review of the MA27 interface
  • Loading branch information
@cvanaret
cvanaret authored Dec 1, 2024
2 parents f3d4bd9 + b70d50b commit 48554eb
Show file tree
Hide file tree
Showing 5 changed files with 85 additions and 90 deletions.
4 changes: 4 additions & 0 deletions bindings/AMPL/AMPLModel.cpp
View file
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
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
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
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
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

0 comments on commit 48554eb

Please sign in to comment.