Skip to content

Commit

Permalink
Cleaned up MA27 factorization function
Browse files Browse the repository at this point in the history
  • Loading branch information
@cvanaret
cvanaret committed Dec 1, 2024
1 parent 7ff8c68 commit d3a62b3
Show file tree
Hide file tree
Showing 3 changed files with 58 additions and 68 deletions.
2 changes: 1 addition & 1 deletion bindings/AMPL/AMPLModel.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -260,7 +260,7 @@ namespace uno {
const Collection<size_t>& AMPLModel::get_upper_bounded_variables() const {
return this->upper_bounded_variables_collection;
}

void AMPLModel::evaluate_lagrangian_hessian(const Vector<double>& x, double objective_multiplier, const Vector<double>& multipliers,
SymmetricMatrix<size_t, double>& hessian) const {
assert(hessian.capacity() >= this->number_asl_hessian_nonzeros);
Expand Down
122 changes: 56 additions & 66 deletions uno/solvers/MA27/MA27Solver.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -104,7 +104,7 @@ 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.
};


Expand Down Expand Up @@ -165,58 +165,41 @@ namespace uno {
std::copy(matrix.data_pointer(), matrix.data_pointer() + matrix.number_nonzeros(), factor.begin());

// 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);
// 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));

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;
}
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
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");
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 @@ -237,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 @@ -250,34 +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::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::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
2 changes: 1 addition & 1 deletion uno/solvers/MA27/MA27Solver.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -58,7 +58,7 @@ class MA27Solver

// 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 d3a62b3

Please sign in to comment.