diff --git a/bindings/AMPL/AMPLModel.cpp b/bindings/AMPL/AMPLModel.cpp index e846e2f3..5882c4a6 100644 --- a/bindings/AMPL/AMPLModel.cpp +++ b/bindings/AMPL/AMPLModel.cpp @@ -260,7 +260,7 @@ namespace uno { const Collection& AMPLModel::get_upper_bounded_variables() const { return this->upper_bounded_variables_collection; } - + void AMPLModel::evaluate_lagrangian_hessian(const Vector& x, double objective_multiplier, const Vector& multipliers, SymmetricMatrix& hessian) const { assert(hessian.capacity() >= this->number_asl_hessian_nonzeros); diff --git a/uno/solvers/MA27/MA27Solver.cpp b/uno/solvers/MA27/MA27Solver.cpp index a976eece..a82890d5 100644 --- a/uno/solvers/MA27/MA27Solver.cpp +++ b/uno/solvers/MA27/MA27Solver.cpp @@ -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. }; @@ -165,50 +165,33 @@ namespace uno { std::copy(matrix.data_pointer(), matrix.data_pointer() + matrix.number_nonzeros(), factor.begin()); // numerical factorization - int la = static_cast(factor.size()); - int liw = static_cast(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(factor.size()); + int liw = static_cast(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(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(info[eINFO::IERROR])); + factorization_done = false; + } } - this->w.resize(static_cast(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& matrix, const Vector& rhs, - Vector& result) { - // solve + void MA27Solver::solve_indefinite_system(const SymmetricMatrix& /*matrix*/, const Vector& rhs, Vector& result) { int la = static_cast(factor.size()); int liw = static_cast(iw.size()); @@ -216,7 +199,7 @@ namespace uno { 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'; } @@ -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(info[eINFO::IERROR]) : static_cast(n); + return (info[eINFO::IFLAG] == eIFLAG::RANK_DEFICIENT) ? static_cast(info[eINFO::IERROR]) : static_cast(n); } void MA27Solver::save_matrix_to_local_format(const SymmetricMatrix& matrix) { @@ -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(row_index + fortran_shift)); icn.emplace_back(static_cast(column_index + fortran_shift)); factor.emplace_back(element); } } - void MA27Solver::repeat_factorization_after_resizing([[maybe_unused]] const SymmetricMatrix& 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(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(info[eINFO::IERROR])); - } - - int la = static_cast(factor.size()); - int liw = static_cast(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 diff --git a/uno/solvers/MA27/MA27Solver.hpp b/uno/solvers/MA27/MA27Solver.hpp index 9561053a..51d961ca 100644 --- a/uno/solvers/MA27/MA27Solver.hpp +++ b/uno/solvers/MA27/MA27Solver.hpp @@ -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& matrix); - void repeat_factorization_after_resizing(const SymmetricMatrix &matrix); + void check_factorization_status(); }; } // namespace uno