From 4717794c4d0899ca1b65e4fb2d82de1167ecf870 Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Sun, 1 Dec 2024 20:09:37 +0100 Subject: [PATCH] MA27Solver: create the w vector (double workspace) once and resize it at every iteration + minor changes --- uno/solvers/MA27/MA27Solver.cpp | 60 ++++++++++++++++----------------- uno/solvers/MA27/MA27Solver.hpp | 1 + 2 files changed, 31 insertions(+), 30 deletions(-) diff --git a/uno/solvers/MA27/MA27Solver.cpp b/uno/solvers/MA27/MA27Solver.cpp index e0505e28..a976eece 100644 --- a/uno/solvers/MA27/MA27Solver.cpp +++ b/uno/solvers/MA27/MA27Solver.cpp @@ -125,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); } @@ -151,36 +151,12 @@ namespace uno { // resize the factor by at least INFO(5) (here, 50% more) factor.resize(static_cast(3 * info[eINFO::NRLNEC] / 2)); - 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"); @@ -198,6 +174,8 @@ namespace uno { repeat_factorization_after_resizing(matrix); } + this->w.resize(static_cast(maxfrt)); + switch (info[eINFO::IFLAG]) { case NSTEPS: WARNING << "MA27BD: Value of NSTEPS outside the range 1 ≤ NSTEPS ≤ N" << '\n'; @@ -226,13 +204,11 @@ namespace uno { DEBUG << "MA27BD: Matrix is rank deficient. Rank: " << info[eINFO::IERROR] << " whereas dimension " << n << '\n'; break; } - } void MA27Solver::solve_indefinite_system([[maybe_unused]]const SymmetricMatrix& matrix, const Vector& rhs, Vector& result) { // solve - std::vector w(maxfrt); // double workspace int la = static_cast(factor.size()); int liw = static_cast(iw.size()); @@ -241,7 +217,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"); - if (eIFLAG::SUCCESS != info[eINFO::IFLAG]) { + if (info[eINFO::IFLAG] != eIFLAG::SUCCESS) { WARNING << "MA27 has issued a warning: IFLAG = " << info[eINFO::IFLAG] << " additional info, IERROR = " << info[eINFO::IERROR] << '\n'; } } @@ -280,4 +256,28 @@ namespace uno { 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); + } + } } // namespace diff --git a/uno/solvers/MA27/MA27Solver.hpp b/uno/solvers/MA27/MA27Solver.hpp index e8e208cd..9561053a 100644 --- a/uno/solvers/MA27/MA27Solver.hpp +++ b/uno/solvers/MA27/MA27Solver.hpp @@ -53,6 +53,7 @@ class MA27Solver 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