Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Linear solvers: renamed symbolic factorization into symbolic analysis #118

Merged
merged 1 commit into from
Dec 2, 2024
+79 −104
Merged

Linear solvers: renamed symbolic factorization into symbolic analysis #118

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 5 additions & 5 deletions uno/ingredients/hessian_models/ConvexifiedHessian.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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) {
Expand All @@ -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) {
Expand Down
2 changes: 1 addition & 1 deletion uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -92,7 +92,7 @@ namespace uno {
void SymmetricIndefiniteLinearSystem<ElementType>::factorize_matrix(const Model& /*model*/,
DirectSymmetricIndefiniteLinearSolver<size_t, ElementType>& 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++;
Expand Down
3 changes: 2 additions & 1 deletion uno/preprocessing/Preprocessing.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -59,7 +59,8 @@ namespace uno {

/* solve the system */
Vector<double> 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
Expand Down
3 changes: 1 addition & 2 deletions uno/solvers/DirectSymmetricIndefiniteLinearSolver.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,8 +13,7 @@ namespace uno {
explicit DirectSymmetricIndefiniteLinearSolver(size_t dimension) : SymmetricIndefiniteLinearSolver<IndexType, ElementType>(dimension) { };
virtual ~DirectSymmetricIndefiniteLinearSolver() = default;

virtual void factorize(const SymmetricMatrix<IndexType, ElementType>& matrix) = 0;
virtual void do_symbolic_factorization(const SymmetricMatrix<IndexType, ElementType>& matrix) = 0;
virtual void do_symbolic_analysis(const SymmetricMatrix<IndexType, ElementType>& matrix) = 0;
virtual void do_numerical_factorization(const SymmetricMatrix<IndexType, ElementType>& matrix) = 0;

[[nodiscard]] virtual std::tuple<size_t, size_t, size_t> get_inertia() const = 0;
Expand Down
14 changes: 4 additions & 10 deletions uno/solvers/MA27/MA27Solver.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -110,7 +110,7 @@ namespace uno {

MA27Solver::MA27Solver(size_t max_dimension, size_t max_number_nonzeros):
DirectSymmetricIndefiniteLinearSolver<size_t, double>(max_dimension),
nz_max(static_cast<int>(max_number_nonzeros)), n(static_cast<int>(max_dimension)), nnz(static_cast<int>(max_number_nonzeros)),
n(static_cast<int>(max_dimension)), nnz(static_cast<int>(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) {
Expand All @@ -124,13 +124,7 @@ namespace uno {
icntl[eICNTL::LDIAG] = 0;
}

void MA27Solver::factorize(const SymmetricMatrix<size_t, double>& matrix) {
// general factorization method: symbolic factorization and numerical factorization
do_symbolic_factorization(matrix);
do_numerical_factorization(matrix);
}

void MA27Solver::do_symbolic_factorization(const SymmetricMatrix<size_t, double>& matrix) {
void MA27Solver::do_symbolic_analysis(const SymmetricMatrix<size_t, double>& 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");
Expand All @@ -141,7 +135,7 @@ namespace uno {
n = static_cast<int>(matrix.dimension());
nnz = static_cast<int>(matrix.number_nonzeros());

// symbolic factorization
// symbolic analysis
int liw = static_cast<int>(iw.size());
MA27AD(&n, &nnz, /* size info */
irn.data(), icn.data(), /* matrix indices */
Expand All @@ -151,7 +145,7 @@ namespace uno {
// resize the factor by at least INFO(5) (here, 50% more)
factor.resize(static_cast<size_t>(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';
}
Expand Down
103 changes: 49 additions & 54 deletions uno/solvers/MA27/MA27Solver.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,57 +9,52 @@
#include "solvers/DirectSymmetricIndefiniteLinearSolver.hpp"

namespace uno {
// forward declaration
template <typename ElementType>
class Vector;


class MA27Solver
: public DirectSymmetricIndefiniteLinearSolver<size_t, double>
{
public:
explicit MA27Solver(size_t max_dimension, size_t max_number_nonzeros);
~MA27Solver() override = default;

void factorize(const SymmetricMatrix<size_t, double>& matrix) override;
void do_symbolic_factorization(const SymmetricMatrix<size_t, double>& matrix) override;
void do_numerical_factorization(const SymmetricMatrix<size_t, double>& matrix) override;
void solve_indefinite_system(const SymmetricMatrix<size_t, double>& matrix, const Vector<double>& rhs, Vector<double>& result) override;


[[nodiscard]] std::tuple<size_t, size_t, size_t> 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<int, 30> icntl{}; // integer array of length 30; integer control values
std::array<double, 5> cntl{}; // double array of length 5; double control values

std::vector<int> irn{}; // row index of input
std::vector<int> icn{}; // col index of input

std::vector<int> iw{}; // integer workspace of length liw
std::vector<int> ikeep{}; // integer array of 3*n; pivot sequence
std::vector<int> 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<int, 20> info{}; // integer array of length 20
double ops{}; // double, operations count

std::vector<double> factor{}; // data array of length la;
int maxfrt{}; // integer, to be set by ma27
std::vector<double> w{}; // double workspace


// 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 check_factorization_status();
};

} // namespace uno
#endif // UNO_MA27SOLVER_H
// forward declaration
template <typename ElementType>
class Vector;

class MA27Solver: public DirectSymmetricIndefiniteLinearSolver<size_t, double> {
public:
explicit MA27Solver(size_t max_dimension, size_t max_number_nonzeros);
~MA27Solver() override = default;

void do_symbolic_analysis(const SymmetricMatrix<size_t, double>& matrix) override;
void do_numerical_factorization(const SymmetricMatrix<size_t, double>& matrix) override;
void solve_indefinite_system(const SymmetricMatrix<size_t, double>& matrix, const Vector<double>& rhs, Vector<double>& result) override;


[[nodiscard]] std::tuple<size_t, size_t, size_t> 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<int, 30> icntl{}; // integer array of length 30; integer control values
std::array<double, 5> cntl{}; // double array of length 5; double control values

std::vector<int> irn{}; // row index of input
std::vector<int> icn{}; // col index of input

std::vector<int> iw{}; // integer workspace of length liw
std::vector<int> ikeep{}; // integer array of 3*n; pivot sequence
std::vector<int> 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<int, 20> info{}; // integer array of length 20
double ops{}; // double, operations count

std::vector<double> factor{}; // data array of length la;
int maxfrt{}; // integer, to be set by ma27
std::vector<double> w{}; // double workspace


// 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 check_factorization_status();
};
} // namespace

#endif // UNO_MA27SOLVER_H
16 changes: 5 additions & 11 deletions uno/solvers/MA57/MA57Solver.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -50,13 +50,7 @@ namespace uno {
this->icntl[8] = 1;
}

// general factorization method: symbolic factorization and numerical factorization
void MA57Solver::factorize(const SymmetricMatrix<size_t, double>& matrix) {
this->do_symbolic_factorization(matrix);
this->do_numerical_factorization(matrix);
}

void MA57Solver::do_symbolic_factorization(const SymmetricMatrix<size_t, double>& matrix) {
void MA57Solver::do_symbolic_analysis(const SymmetricMatrix<size_t, double>& 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");
Expand All @@ -67,7 +61,7 @@ namespace uno {
const int n = static_cast<int>(matrix.dimension());
const int nnz = static_cast<int>(matrix.number_nonzeros());

// symbolic factorization
// symbolic analysis
MA57AD(/* const */ &n,
/* const */ &nnz,
/* const */ this->row_indices.data(),
Expand All @@ -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';
}
Expand All @@ -90,7 +84,7 @@ namespace uno {
this->fact.resize(static_cast<size_t>(lfact));
this->ifact.resize(static_cast<size_t>(lifact));

// store the sizes of the symbolic factorization
// store the sizes of the symbolic analysis
this->factorization = {n, nnz, lfact, lifact};
}

Expand Down
3 changes: 1 addition & 2 deletions uno/solvers/MA57/MA57Solver.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -33,8 +33,7 @@ namespace uno {
MA57Solver(size_t dimension, size_t number_nonzeros);
~MA57Solver() override = default;

void factorize(const SymmetricMatrix<size_t, double>& matrix) override;
void do_symbolic_factorization(const SymmetricMatrix<size_t, double>& matrix) override;
void do_symbolic_analysis(const SymmetricMatrix<size_t, double>& matrix) override;
void do_numerical_factorization(const SymmetricMatrix<size_t, double>& matrix) override;
void solve_indefinite_system(const SymmetricMatrix<size_t, double>& matrix, const Vector<double>& rhs, Vector<double>& result) override;

Expand Down
8 changes: 1 addition & 7 deletions uno/solvers/MUMPS/MUMPSSolver.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -37,13 +37,7 @@ namespace uno {
dmumps_c(&this->mumps_structure);
}

void MUMPSSolver::factorize(const SymmetricMatrix<size_t, double>& 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<size_t, double>& matrix) {
void MUMPSSolver::do_symbolic_analysis(const SymmetricMatrix<size_t, double>& matrix) {
this->save_matrix_to_local_format(matrix);
this->mumps_structure.n = static_cast<int>(matrix.dimension());
this->mumps_structure.nnz = static_cast<int>(matrix.number_nonzeros());
Expand Down
3 changes: 1 addition & 2 deletions uno/solvers/MUMPS/MUMPSSolver.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,8 +14,7 @@ namespace uno {
explicit MUMPSSolver(size_t dimension, size_t number_nonzeros);
~MUMPSSolver() override;

void factorize(const SymmetricMatrix<size_t, double>& matrix) override;
void do_symbolic_factorization(const SymmetricMatrix<size_t, double>& matrix) override;
void do_symbolic_analysis(const SymmetricMatrix<size_t, double>& matrix) override;
void do_numerical_factorization(const SymmetricMatrix<size_t, double>& matrix) override;
void solve_indefinite_system(const SymmetricMatrix<size_t, double>& matrix, const Vector<double>& rhs,
Vector<double>& result) override;
Expand Down
6 changes: 3 additions & 3 deletions unotest/functional_tests/MA27SolverTests.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,7 @@ TEST(MA27Solver, SystemSize5) {
const std::array<double, n> 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);

Expand All @@ -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();
Expand All @@ -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)
Expand Down
6 changes: 3 additions & 3 deletions unotest/functional_tests/MA57SolverTests.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,7 @@ TEST(MA57Solver, SystemSize5) {
const std::array<double, n> 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);

Expand All @@ -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();
Expand All @@ -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)
Expand Down
6 changes: 3 additions & 3 deletions unotest/functional_tests/MUMPSSolverTests.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,7 @@ TEST(MUMPSSolver, SystemSize5) {
const std::array<double, n> 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);

Expand Down Expand Up @@ -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();
Expand All @@ -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)
Expand Down
Loading