diff --git a/bindings/AMPL/AMPLModel.cpp b/bindings/AMPL/AMPLModel.cpp index f7af1062..5882c4a6 100644 --- a/bindings/AMPL/AMPLModel.cpp +++ b/bindings/AMPL/AMPLModel.cpp @@ -249,6 +249,10 @@ namespace uno { return this->single_upper_bounded_variables_collection; } + const Vector& AMPLModel::get_fixed_variables() const { + return this->fixed_variables; + } + const Collection& AMPLModel::get_lower_bounded_variables() const { return this->lower_bounded_variables_collection; } diff --git a/bindings/AMPL/AMPLModel.hpp b/bindings/AMPL/AMPLModel.hpp index da9aecce..f93fbd34 100644 --- a/bindings/AMPL/AMPLModel.hpp +++ b/bindings/AMPL/AMPLModel.hpp @@ -47,7 +47,7 @@ namespace uno { [[nodiscard]] const SparseVector& get_slacks() const override; [[nodiscard]] const Collection& get_single_lower_bounded_variables() const override; [[nodiscard]] const Collection& get_single_upper_bounded_variables() const override; - [[nodiscard]] const Vector& get_fixed_variables() const override { return this->fixed_variables; } + [[nodiscard]] const Vector& 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; diff --git a/uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp b/uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp index d06c2847..a468b74c 100644 --- a/uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp +++ b/uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp @@ -91,8 +91,6 @@ namespace uno { template void SymmetricIndefiniteLinearSystem::factorize_matrix(const Model& /*model*/, DirectSymmetricIndefiniteLinearSolver& 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); } diff --git a/uno/solvers/MA27/MA27Solver.cpp b/uno/solvers/MA27/MA27Solver.cpp index e0c5e56e..cda96f2f 100644 --- a/uno/solvers/MA27/MA27Solver.cpp +++ b/uno/solvers/MA27/MA27Solver.cpp @@ -104,18 +104,20 @@ 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(max_dimension), nz_max(static_cast(max_number_nonzeros)), n(static_cast(max_dimension)), nnz(static_cast(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; @@ -123,7 +125,7 @@ namespace uno { } void MA27Solver::factorize(const SymmetricMatrix& 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); } @@ -144,93 +146,52 @@ 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(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& 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::do_numerical_factorization([[maybe_unused]]const SymmetricMatrix& matrix) { assert(matrix.dimension() <= iw1.capacity() && "MA27Solver: the dimension of the matrix is larger than the preallocated size"); assert(nnz == static_cast(matrix.number_nonzeros()) && "MA27Solver: the numbers of nonzeros do not match"); - // 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); - } + // 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(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)); + this->check_factorization_status(); } - void MA27Solver::solve_indefinite_system([[maybe_unused]]const SymmetricMatrix& matrix, const Vector& rhs, - Vector& result) { - // solve - std::vector w(maxfrt); // double workspace + 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()); @@ -238,8 +199,8 @@ 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"); - 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'; } } @@ -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(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) { @@ -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(row_index + fortran_shift)); icn.emplace_back(static_cast(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 diff --git a/uno/solvers/MA27/MA27Solver.hpp b/uno/solvers/MA27/MA27Solver.hpp index 149f68cb..51d961ca 100644 --- a/uno/solvers/MA27/MA27Solver.hpp +++ b/uno/solvers/MA27/MA27Solver.hpp @@ -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 icntl{}; // integer array of length 30; integer control values - std::array cntl{}; // double array of length 5; double control values + std::array icntl{}; // integer array of length 30; integer control values + std::array cntl{}; // double array of length 5; double control values std::vector irn{}; // row index of input std::vector icn{}; // col index of input @@ -48,16 +48,17 @@ class MA27Solver std::vector 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 info{}; // integer array of length 20 + std::array info{}; // integer array of length 20 double ops{}; // double, operations count std::vector factor{}; // data array of length la; int maxfrt{}; // integer, to be set by ma27 + std::vector w{}; // double workspace // 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