diff --git a/uno/ingredients/hessian_models/ConvexifiedHessian.cpp b/uno/ingredients/hessian_models/ConvexifiedHessian.cpp index cd9ce44a..fe5a532d 100644 --- a/uno/ingredients/hessian_models/ConvexifiedHessian.cpp +++ b/uno/ingredients/hessian_models/ConvexifiedHessian.cpp @@ -40,7 +40,7 @@ namespace uno { double regularization_factor = (smallest_diagonal_entry > 0.) ? 0. : this->regularization_initial_value - smallest_diagonal_entry; bool good_inertia = false; - bool symbolic_factorization_performed = false; + bool symbolic_analysis_performed = false; while (not good_inertia) { DEBUG << "Testing factorization with regularization factor " << regularization_factor << '\n'; if (0. < regularization_factor) { @@ -50,10 +50,10 @@ namespace uno { } DEBUG << "Current Hessian:\n" << hessian << '\n'; - // perform the symbolic factorization only once - if (not symbolic_factorization_performed) { - this->linear_solver->do_symbolic_factorization(hessian); - symbolic_factorization_performed = true; + // perform the symbolic analysis only once + if (not symbolic_analysis_performed) { + this->linear_solver->do_symbolic_analysis(hessian); + symbolic_analysis_performed = true; } this->linear_solver->do_numerical_factorization(hessian); if (this->linear_solver->rank() == number_original_variables && this->linear_solver->number_negative_eigenvalues() == 0) { diff --git a/uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp b/uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp index a468b74c..69176875 100644 --- a/uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp +++ b/uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp @@ -92,7 +92,7 @@ namespace uno { void SymmetricIndefiniteLinearSystem::factorize_matrix(const Model& /*model*/, DirectSymmetricIndefiniteLinearSolver& linear_solver) { if (true || this->number_factorizations == 0) { - linear_solver.do_symbolic_factorization(this->matrix); + linear_solver.do_symbolic_analysis(this->matrix); } linear_solver.do_numerical_factorization(this->matrix); this->number_factorizations++; diff --git a/uno/preprocessing/Preprocessing.cpp b/uno/preprocessing/Preprocessing.cpp index 1fd658a4..00584c19 100644 --- a/uno/preprocessing/Preprocessing.cpp +++ b/uno/preprocessing/Preprocessing.cpp @@ -59,7 +59,8 @@ namespace uno { /* solve the system */ Vector solution(matrix.dimension()); - linear_solver.factorize(matrix); + linear_solver.do_symbolic_analysis(matrix); + linear_solver.do_numerical_factorization(matrix); linear_solver.solve_indefinite_system(matrix, rhs, solution); // if least-square multipliers too big, discard them. Otherwise, keep them diff --git a/uno/solvers/DirectSymmetricIndefiniteLinearSolver.hpp b/uno/solvers/DirectSymmetricIndefiniteLinearSolver.hpp index 0d7a4059..df5ed7c9 100644 --- a/uno/solvers/DirectSymmetricIndefiniteLinearSolver.hpp +++ b/uno/solvers/DirectSymmetricIndefiniteLinearSolver.hpp @@ -13,8 +13,7 @@ namespace uno { explicit DirectSymmetricIndefiniteLinearSolver(size_t dimension) : SymmetricIndefiniteLinearSolver(dimension) { }; virtual ~DirectSymmetricIndefiniteLinearSolver() = default; - virtual void factorize(const SymmetricMatrix& matrix) = 0; - virtual void do_symbolic_factorization(const SymmetricMatrix& matrix) = 0; + virtual void do_symbolic_analysis(const SymmetricMatrix& matrix) = 0; virtual void do_numerical_factorization(const SymmetricMatrix& matrix) = 0; [[nodiscard]] virtual std::tuple get_inertia() const = 0; diff --git a/uno/solvers/MA27/MA27Solver.cpp b/uno/solvers/MA27/MA27Solver.cpp index cda96f2f..8e2c5a59 100644 --- a/uno/solvers/MA27/MA27Solver.cpp +++ b/uno/solvers/MA27/MA27Solver.cpp @@ -110,7 +110,7 @@ namespace uno { 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)), + 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), // 20% more than 2*nnz + 3*n + 1 ikeep(3 * max_dimension), iw1(2 * max_dimension) { @@ -124,13 +124,7 @@ namespace uno { icntl[eICNTL::LDIAG] = 0; } - void MA27Solver::factorize(const SymmetricMatrix& matrix) { - // general factorization method: symbolic factorization and numerical factorization - do_symbolic_factorization(matrix); - do_numerical_factorization(matrix); - } - - void MA27Solver::do_symbolic_factorization(const SymmetricMatrix& matrix) { + void MA27Solver::do_symbolic_analysis(const SymmetricMatrix& matrix) { assert(matrix.dimension() <= iw1.capacity() && "MA27Solver: the dimension of the matrix is larger than the preallocated size"); assert(matrix.number_nonzeros() <= irn.capacity() && "MA27Solver: the number of nonzeros of the matrix is larger than the preallocated size"); @@ -141,7 +135,7 @@ namespace uno { n = static_cast(matrix.dimension()); nnz = static_cast(matrix.number_nonzeros()); - // symbolic factorization + // symbolic analysis int liw = static_cast(iw.size()); MA27AD(&n, &nnz, /* size info */ irn.data(), icn.data(), /* matrix indices */ @@ -151,7 +145,7 @@ namespace uno { // resize the factor by at least INFO(5) (here, 50% more) factor.resize(static_cast(3 * info[eINFO::NRLNEC] / 2)); - assert(info[eINFO::IFLAG] == eIFLAG::SUCCESS && "MA27: the symbolic factorization failed"); + assert(info[eINFO::IFLAG] == eIFLAG::SUCCESS && "MA27: the symbolic analysis failed"); if (info[eINFO::IFLAG] != eIFLAG::SUCCESS) { WARNING << "MA27 has issued a warning: IFLAG = " << info[eINFO::IFLAG] << " additional info, IERROR = " << info[eINFO::IERROR] << '\n'; } diff --git a/uno/solvers/MA27/MA27Solver.hpp b/uno/solvers/MA27/MA27Solver.hpp index 51d961ca..2ef8f776 100644 --- a/uno/solvers/MA27/MA27Solver.hpp +++ b/uno/solvers/MA27/MA27Solver.hpp @@ -9,57 +9,52 @@ #include "solvers/DirectSymmetricIndefiniteLinearSolver.hpp" namespace uno { -// forward declaration -template -class Vector; - - -class MA27Solver - : public DirectSymmetricIndefiniteLinearSolver -{ -public: - explicit MA27Solver(size_t max_dimension, size_t max_number_nonzeros); - ~MA27Solver() override = default; - - void factorize(const SymmetricMatrix& matrix) override; - void do_symbolic_factorization(const SymmetricMatrix& matrix) override; - void do_numerical_factorization(const SymmetricMatrix& matrix) override; - void solve_indefinite_system(const SymmetricMatrix& matrix, const Vector& rhs, Vector& result) override; - - - [[nodiscard]] std::tuple get_inertia() const override; - [[nodiscard]] size_t number_negative_eigenvalues() const override; - // [[nodiscard]] bool matrix_is_positive_definite() const override; - [[nodiscard]] bool matrix_is_singular() const override; - [[nodiscard]] size_t rank() const override; - -private: - 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::vector irn{}; // row index of input - std::vector icn{}; // col index of input - - std::vector iw{}; // integer workspace of length liw - std::vector ikeep{}; // integer array of 3*n; pivot sequence - 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 - 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 check_factorization_status(); -}; - -} // namespace uno -#endif // UNO_MA27SOLVER_H + // forward declaration + template + class Vector; + + class MA27Solver: public DirectSymmetricIndefiniteLinearSolver { + public: + explicit MA27Solver(size_t max_dimension, size_t max_number_nonzeros); + ~MA27Solver() override = default; + + void do_symbolic_analysis(const SymmetricMatrix& matrix) override; + void do_numerical_factorization(const SymmetricMatrix& matrix) override; + void solve_indefinite_system(const SymmetricMatrix& matrix, const Vector& rhs, Vector& result) override; + + + [[nodiscard]] std::tuple get_inertia() const override; + [[nodiscard]] size_t number_negative_eigenvalues() const override; + // [[nodiscard]] bool matrix_is_positive_definite() const override; + [[nodiscard]] bool matrix_is_singular() const override; + [[nodiscard]] size_t rank() const override; + + private: + 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::vector irn{}; // row index of input + std::vector icn{}; // col index of input + + std::vector iw{}; // integer workspace of length liw + std::vector ikeep{}; // integer array of 3*n; pivot sequence + 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 + 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 check_factorization_status(); + }; +} // namespace + +#endif // UNO_MA27SOLVER_H \ No newline at end of file diff --git a/uno/solvers/MA57/MA57Solver.cpp b/uno/solvers/MA57/MA57Solver.cpp index 6e482d48..0688f66a 100644 --- a/uno/solvers/MA57/MA57Solver.cpp +++ b/uno/solvers/MA57/MA57Solver.cpp @@ -19,7 +19,7 @@ namespace uno { // MA57 // default values of controlling parameters void MA57ID(double cntl[], int icntl[]); - // symbolic factorization + // symbolic analysis void MA57AD(const int* n, const int* ne, const int irn[], const int jcn[], const int* lkeep, int keep[], int iwork[], int icntl[], int info[], double rinfo[]); // numerical factorization @@ -50,13 +50,7 @@ namespace uno { this->icntl[8] = 1; } - // general factorization method: symbolic factorization and numerical factorization - void MA57Solver::factorize(const SymmetricMatrix& matrix) { - this->do_symbolic_factorization(matrix); - this->do_numerical_factorization(matrix); - } - - void MA57Solver::do_symbolic_factorization(const SymmetricMatrix& matrix) { + void MA57Solver::do_symbolic_analysis(const SymmetricMatrix& matrix) { assert(matrix.dimension() <= this->dimension && "MA57Solver: the dimension of the matrix is larger than the preallocated size"); assert(matrix.number_nonzeros() <= this->row_indices.capacity() && "MA57Solver: the number of nonzeros of the matrix is larger than the preallocated size"); @@ -67,7 +61,7 @@ namespace uno { const int n = static_cast(matrix.dimension()); const int nnz = static_cast(matrix.number_nonzeros()); - // symbolic factorization + // symbolic analysis MA57AD(/* const */ &n, /* const */ &nnz, /* const */ this->row_indices.data(), @@ -79,7 +73,7 @@ namespace uno { /* out */ this->info.data(), /* out */ this->rinfo.data()); - assert(0 <= info[0] && "MA57: the symbolic factorization failed"); + assert(0 <= info[0] && "MA57: the symbolic analysis failed"); if (0 < info[0]) { WARNING << "MA57 has issued a warning: info(1) = " << info[0] << '\n'; } @@ -90,7 +84,7 @@ namespace uno { this->fact.resize(static_cast(lfact)); this->ifact.resize(static_cast(lifact)); - // store the sizes of the symbolic factorization + // store the sizes of the symbolic analysis this->factorization = {n, nnz, lfact, lifact}; } diff --git a/uno/solvers/MA57/MA57Solver.hpp b/uno/solvers/MA57/MA57Solver.hpp index 9ccded90..33f9c524 100644 --- a/uno/solvers/MA57/MA57Solver.hpp +++ b/uno/solvers/MA57/MA57Solver.hpp @@ -33,8 +33,7 @@ namespace uno { MA57Solver(size_t dimension, size_t number_nonzeros); ~MA57Solver() override = default; - void factorize(const SymmetricMatrix& matrix) override; - void do_symbolic_factorization(const SymmetricMatrix& matrix) override; + void do_symbolic_analysis(const SymmetricMatrix& matrix) override; void do_numerical_factorization(const SymmetricMatrix& matrix) override; void solve_indefinite_system(const SymmetricMatrix& matrix, const Vector& rhs, Vector& result) override; diff --git a/uno/solvers/MUMPS/MUMPSSolver.cpp b/uno/solvers/MUMPS/MUMPSSolver.cpp index 8a754c3d..1218cbbd 100644 --- a/uno/solvers/MUMPS/MUMPSSolver.cpp +++ b/uno/solvers/MUMPS/MUMPSSolver.cpp @@ -37,13 +37,7 @@ namespace uno { dmumps_c(&this->mumps_structure); } - void MUMPSSolver::factorize(const SymmetricMatrix& matrix) { - // general factorization method: symbolic factorization and numerical factorization - this->do_symbolic_factorization(matrix); - this->do_numerical_factorization(matrix); - } - - void MUMPSSolver::do_symbolic_factorization(const SymmetricMatrix& matrix) { + void MUMPSSolver::do_symbolic_analysis(const SymmetricMatrix& matrix) { this->save_matrix_to_local_format(matrix); this->mumps_structure.n = static_cast(matrix.dimension()); this->mumps_structure.nnz = static_cast(matrix.number_nonzeros()); diff --git a/uno/solvers/MUMPS/MUMPSSolver.hpp b/uno/solvers/MUMPS/MUMPSSolver.hpp index 83962b36..0c7e6ca4 100644 --- a/uno/solvers/MUMPS/MUMPSSolver.hpp +++ b/uno/solvers/MUMPS/MUMPSSolver.hpp @@ -14,8 +14,7 @@ namespace uno { explicit MUMPSSolver(size_t dimension, size_t number_nonzeros); ~MUMPSSolver() override; - void factorize(const SymmetricMatrix& matrix) override; - void do_symbolic_factorization(const SymmetricMatrix& matrix) override; + void do_symbolic_analysis(const SymmetricMatrix& matrix) override; void do_numerical_factorization(const SymmetricMatrix& matrix) override; void solve_indefinite_system(const SymmetricMatrix& matrix, const Vector& rhs, Vector& result) override; diff --git a/unotest/functional_tests/MA27SolverTests.cpp b/unotest/functional_tests/MA27SolverTests.cpp index 4c1f6386..a123a817 100644 --- a/unotest/functional_tests/MA27SolverTests.cpp +++ b/unotest/functional_tests/MA27SolverTests.cpp @@ -24,7 +24,7 @@ TEST(MA27Solver, SystemSize5) { const std::array reference{1., 2., 3., 4., 5.}; MA27Solver solver(n, nnz); - solver.do_symbolic_factorization(matrix); + solver.do_symbolic_analysis(matrix); solver.do_numerical_factorization(matrix); solver.solve_indefinite_system(matrix, rhs, result); @@ -46,7 +46,7 @@ TEST(MA27Solver, Inertia) { matrix.insert(1., 4, 4); MA27Solver solver(n, nnz); - solver.do_symbolic_factorization(matrix); + solver.do_symbolic_analysis(matrix); solver.do_numerical_factorization(matrix); const auto [number_positive, number_negative, number_zero] = solver.get_inertia(); @@ -68,7 +68,7 @@ TEST(MA27Solver, SingularMatrix) { matrix.insert(0., 2, 2); matrix.insert(0., 3, 3); MA27Solver solver(n, nnz); - solver.do_symbolic_factorization(matrix); + solver.do_symbolic_analysis(matrix); solver.do_numerical_factorization(matrix); // expected inertia (1, 1, 2) diff --git a/unotest/functional_tests/MA57SolverTests.cpp b/unotest/functional_tests/MA57SolverTests.cpp index 839dbd23..a75eb355 100644 --- a/unotest/functional_tests/MA57SolverTests.cpp +++ b/unotest/functional_tests/MA57SolverTests.cpp @@ -24,7 +24,7 @@ TEST(MA57Solver, SystemSize5) { const std::array reference{1., 2., 3., 4., 5.}; MA57Solver solver(n, nnz); - solver.do_symbolic_factorization(matrix); + solver.do_symbolic_analysis(matrix); solver.do_numerical_factorization(matrix); solver.solve_indefinite_system(matrix, rhs, result); @@ -46,7 +46,7 @@ TEST(MA57Solver, Inertia) { matrix.insert(1., 4, 4); MA57Solver solver(n, nnz); - solver.do_symbolic_factorization(matrix); + solver.do_symbolic_analysis(matrix); solver.do_numerical_factorization(matrix); const auto [number_positive, number_negative, number_zero] = solver.get_inertia(); @@ -68,7 +68,7 @@ TEST(MA57Solver, SingularMatrix) { matrix.insert(0., 2, 2); matrix.insert(0., 3, 3); MA57Solver solver(n, nnz); - solver.do_symbolic_factorization(matrix); + solver.do_symbolic_analysis(matrix); solver.do_numerical_factorization(matrix); // expected inertia (1, 1, 2) diff --git a/unotest/functional_tests/MUMPSSolverTests.cpp b/unotest/functional_tests/MUMPSSolverTests.cpp index 44b3c297..a886c9a0 100644 --- a/unotest/functional_tests/MUMPSSolverTests.cpp +++ b/unotest/functional_tests/MUMPSSolverTests.cpp @@ -26,7 +26,7 @@ TEST(MUMPSSolver, SystemSize5) { const std::array reference{1., 2., 3., 4., 5.}; MUMPSSolver solver(n, nnz); - solver.do_symbolic_factorization(matrix); + solver.do_symbolic_analysis(matrix); solver.do_numerical_factorization(matrix); solver.solve_indefinite_system(matrix, rhs, result); @@ -65,7 +65,7 @@ TEST(MUMPSSolver, Inertia) { matrix.insert(1., 4, 4); MUMPSSolver solver(n, nnz); - solver.do_symbolic_factorization(matrix); + solver.do_symbolic_analysis(matrix); solver.do_numerical_factorization(matrix); const auto [number_positive, number_negative, number_zero] = solver.get_inertia(); @@ -87,7 +87,7 @@ TEST(MUMPSSolver, SingularMatrix) { matrix.insert(0., 2, 2); matrix.insert(0., 3, 3); MUMPSSolver solver(n, nnz); - solver.do_symbolic_factorization(matrix); + solver.do_symbolic_analysis(matrix); solver.do_numerical_factorization(matrix); // expected inertia (1, 1, 2)