From 2075a167f4db86e8be4990b26e6b584a819d5b6e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Thu, 15 Oct 2020 16:14:21 +0200
Subject: [PATCH] Plug PETSc (and various classic solvers and preconditioners)
The 'LinearSolver' class proposes mechanism to choose the used
- library,
- method, and
- preconditioner
Also one can adjust the required convergence epsilon, the maximum
number of iterations and the verbosity.
All of these options can be listed through the 'pgs' script using the
'printLSAvailable()' builtin function.
One drive the linear system solver by changing values in the
script. For instance:
setLSVerbosity(true);
setLSEpsilon(1E-9);
setLSMaxIter(1000);
setLSLibrary("PETSc");
setLSMethod("LU");
setLSPrecond("none");
Finally one can print the current options by using the
'printLSOptions()' function.
---
CMakeLists.txt | 24 +-
src/algebra/BiCGStab.hpp | 26 +-
src/algebra/{PCG.hpp => CG.hpp} | 40 +--
src/algebra/CMakeLists.txt | 3 +-
src/algebra/CRSMatrix.hpp | 15 +-
src/algebra/LinearSolver.cpp | 320 +++++++++++++++++++
src/algebra/LinearSolver.hpp | 36 +++
src/algebra/LinearSolverOptions.cpp | 15 +
src/algebra/LinearSolverOptions.hpp | 233 ++++++++++++++
src/algebra/LinearSystemSolver.cpp | 331 --------------------
src/algebra/LinearSystemSolver.hpp | 79 -----
src/algebra/PETScWrapper.cpp | 1 +
src/algebra/SparseMatrixDescriptor.hpp | 7 +-
src/language/ast/ASTNodeDataType.cpp | 12 +-
src/language/modules/CMakeLists.txt | 1 +
src/language/modules/LinearSolverModule.cpp | 85 +++++
src/language/modules/LinearSolverModule.hpp | 19 ++
src/language/modules/ModuleRepository.cpp | 2 +
tests/CMakeLists.txt | 6 +-
tests/test_BiCGStab.cpp | 4 +-
tests/{test_PCG.cpp => test_CG.cpp} | 10 +-
tests/test_SparseMatrixDescriptor.cpp | 14 +-
22 files changed, 815 insertions(+), 468 deletions(-)
rename src/algebra/{PCG.hpp => CG.hpp} (60%)
create mode 100644 src/algebra/LinearSolver.cpp
create mode 100644 src/algebra/LinearSolver.hpp
create mode 100644 src/algebra/LinearSolverOptions.cpp
create mode 100644 src/algebra/LinearSolverOptions.hpp
delete mode 100644 src/algebra/LinearSystemSolver.cpp
delete mode 100644 src/algebra/LinearSystemSolver.hpp
create mode 100644 src/language/modules/LinearSolverModule.cpp
create mode 100644 src/language/modules/LinearSolverModule.hpp
rename tests/{test_PCG.cpp => test_CG.cpp} (91%)
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 6b4f0a54d..dc0f69398 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -139,8 +139,16 @@ endif()
#------------------------------------------------------
# Check for PETSc
-find_package(PETSc)
-set(PUGS_HAS_PETSC ${PETSC_FOUND})
+# defaults use PETSc
+set(PUGS_ENABLE_PETSC AUTO CACHE STRING
+ "Choose one of: AUTO ON OFF")
+
+if (PUGS_ENABLE_PETSC MATCHES "^(AUTO|ON)$")
+ find_package(PETSc)
+ set(PUGS_HAS_PETSC ${PETSC_FOUND})
+else()
+ unset(PUGS_HAS_PETSC)
+endif()
if (${PETSC_FOUND})
include_directories(SYSTEM ${PETSC_INCLUDES})
@@ -474,14 +482,22 @@ else()
if(NOT PARMETIS_LIBRARIES)
message(" MPI: deactivated: ParMETIS cannot be found!")
else()
- message(" MPI: not found!")
+ if (PUGS_ENABLE_MPI MATCHES "^(AUTO|ON)$")
+ message(" MPI: not found!")
+ else()
+ message(" MPI: explicitly deactivated!")
+ endif()
endif()
endif()
if (PETSC_FOUND)
message(" PETSc: ${PETSC_VERSION}")
else()
- message(" PETSc: not found!")
+ if (PUGS_ENABLE_PETSC MATCHES "^(AUTO|ON)$")
+ message(" PETSc: not found!")
+ else()
+ message(" PETSc: explicitly deactivated!")
+ endif()
endif()
if(CLANG_FORMAT)
diff --git a/src/algebra/BiCGStab.hpp b/src/algebra/BiCGStab.hpp
index bc822668a..ff34b89e1 100644
--- a/src/algebra/BiCGStab.hpp
+++ b/src/algebra/BiCGStab.hpp
@@ -7,16 +7,20 @@
#include <rang.hpp>
-template <bool verbose = true>
struct BiCGStab
{
template <typename VectorType, typename MatrixType>
- BiCGStab(const VectorType& b, const MatrixType& A, VectorType& x, const size_t max_iter, const double epsilon = 1e-6)
+ BiCGStab(const MatrixType& A,
+ VectorType& x,
+ const VectorType& b,
+ const double epsilon,
+ const size_t maximum_iteration,
+ const bool verbose)
{
- if constexpr (verbose) {
+ if (verbose) {
std::cout << "- bi-conjugate gradient stabilized\n";
std::cout << " epsilon = " << epsilon << '\n';
- std::cout << " maximum number of iterations: " << max_iter << '\n';
+ std::cout << " maximum number of iterations: " << maximum_iteration << '\n';
}
VectorType r_k_1{b.size()};
@@ -38,13 +42,13 @@ struct BiCGStab
VectorType r_k{x.size()};
- if constexpr (verbose) {
+ if (verbose) {
std::cout << " initial residu: " << resid0 << '\n';
}
- for (size_t i = 1; i <= max_iter; ++i) {
- if constexpr (verbose) {
- std::cout << " - iteration: " << std::setw(6) << i << "\tresidu: " << residu / resid0
- << "\tabsolute: " << residu << '\n';
+ for (size_t i = 1; i <= maximum_iteration; ++i) {
+ if (verbose) {
+ std::cout << " - iteration: " << std::setw(6) << i << " residu: " << std::scientific << residu / resid0
+ << " absolute: " << std::scientific << residu << '\n';
}
Ap_k = A * p_k;
@@ -77,8 +81,8 @@ struct BiCGStab
<< '\n';
;
std::cout << " - epsilon: " << epsilon << '\n';
- std::cout << " - relative residu : " << residu / resid0 << '\n';
- std::cout << " - absolute residu : " << residu << '\n';
+ std::cout << " - relative residu : " << std::scientific << residu / resid0 << '\n';
+ std::cout << " - absolute residu : " << std::scientific << residu << '\n';
}
}
}
diff --git a/src/algebra/PCG.hpp b/src/algebra/CG.hpp
similarity index 60%
rename from src/algebra/PCG.hpp
rename to src/algebra/CG.hpp
index 36a4635b3..538c77571 100644
--- a/src/algebra/PCG.hpp
+++ b/src/algebra/CG.hpp
@@ -6,30 +6,30 @@
#include <rang.hpp>
-template <bool verbose = true>
-struct PCG
+struct CG
{
- template <typename VectorType, typename MatrixType, typename PreconditionerType>
- PCG(const VectorType& f,
- const MatrixType& A,
- [[maybe_unused]] const PreconditionerType& C,
- VectorType& x,
- const size_t maxiter,
- const double epsilon = 1e-6)
+ template <typename VectorType, typename MatrixType>
+ CG(const MatrixType& A,
+ VectorType& x,
+ const VectorType& f,
+ const double epsilon,
+ const size_t maximum_iteration,
+ const bool verbose = false)
{
- if constexpr (verbose) {
+ if (verbose) {
std::cout << "- conjugate gradient\n";
std::cout << " epsilon = " << epsilon << '\n';
- std::cout << " maximum number of iterations: " << maxiter << '\n';
+ std::cout << " maximum number of iterations: " << maximum_iteration << '\n';
}
VectorType h{f.size()};
VectorType b = copy(f);
- if constexpr (verbose) {
+ if (verbose) {
h = A * x;
h -= f;
- std::cout << "- initial *real* residu : " << (h, h) << '\n';
+ std::cout << "- initial " << rang::style::bold << "real" << rang::style::reset << " residu : " << (h, h)
+ << '\n';
}
VectorType g{b.size()};
@@ -40,7 +40,7 @@ struct PCG
double relativeEpsilon = epsilon;
- for (size_t i = 1; i <= maxiter; ++i) {
+ for (size_t i = 1; i <= maximum_iteration; ++i) {
if (i == 1) {
h = A * x;
@@ -74,13 +74,13 @@ struct PCG
if ((i == 1) && (gcg != 0)) {
relativeEpsilon = epsilon * gcg;
gcg0 = gcg;
- if constexpr (verbose) {
+ if (verbose) {
std::cout << " initial residu: " << gcg << '\n';
}
}
- if constexpr (verbose) {
- std::cout << " - iteration " << std::setw(6) << i << "\tresidu: " << gcg / gcg0;
- std::cout << "\tabsolute: " << gcg << '\n';
+ if (verbose) {
+ std::cout << " - iteration " << std::setw(6) << i << std::scientific << " residu: " << gcg / gcg0;
+ std::cout << " absolute: " << std::scientific << gcg << '\n';
}
if (gcg < relativeEpsilon) {
@@ -96,8 +96,8 @@ struct PCG
if (gcg > relativeEpsilon) {
std::cout << " conjugate gradient: " << rang::fgB::red << "*NOT CONVERGED*" << rang::style::reset << '\n';
std::cout << " - epsilon: " << epsilon << '\n';
- std::cout << " - relative residu : " << gcg / gcg0 << '\n';
- std::cout << " - absolute residu : " << gcg << '\n';
+ std::cout << " - relative residu : " << std::scientific << gcg / gcg0 << '\n';
+ std::cout << " - absolute residu : " << std::scientific << gcg << '\n';
}
}
};
diff --git a/src/algebra/CMakeLists.txt b/src/algebra/CMakeLists.txt
index bbd2d01bf..d1f01ec31 100644
--- a/src/algebra/CMakeLists.txt
+++ b/src/algebra/CMakeLists.txt
@@ -2,5 +2,6 @@
add_library(
PugsAlgebra
- LinearSystemSolver.cpp
+ LinearSolver.cpp
+ LinearSolverOptions.cpp
PETScWrapper.cpp)
diff --git a/src/algebra/CRSMatrix.hpp b/src/algebra/CRSMatrix.hpp
index 97a0a8a8e..d9b03ac93 100644
--- a/src/algebra/CRSMatrix.hpp
+++ b/src/algebra/CRSMatrix.hpp
@@ -30,8 +30,21 @@ class CRSMatrix
return m_host_matrix.numRows();
}
+ auto
+ values() const
+ {
+ return m_values;
+ }
+
+ auto
+ hostMatrix() const
+ {
+ return m_host_matrix;
+ }
+
template <typename DataType2>
- Vector<MutableDataType> operator*(const Vector<DataType2>& x) const
+ Vector<MutableDataType>
+ operator*(const Vector<DataType2>& x) const
{
static_assert(std::is_same<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>(),
"Cannot multiply matrix and vector of different type");
diff --git a/src/algebra/LinearSolver.cpp b/src/algebra/LinearSolver.cpp
new file mode 100644
index 000000000..a275cd581
--- /dev/null
+++ b/src/algebra/LinearSolver.cpp
@@ -0,0 +1,320 @@
+#include <algebra/LinearSolver.hpp>
+#include <utils/pugs_config.hpp>
+
+#include <algebra/BiCGStab.hpp>
+#include <algebra/CG.hpp>
+
+#ifdef PUGS_HAS_PETSC
+#include <petsc.h>
+#endif // PUGS_HAS_PETSC
+
+struct LinearSolver::Internals
+{
+ static bool
+ hasLibrary(const LSLibrary library)
+ {
+ switch (library) {
+ case LSLibrary::builtin: {
+ return true;
+ }
+ case LSLibrary::petsc: {
+#ifdef PUGS_HAS_PETSC
+ return true;
+#else
+ return false;
+#endif
+ }
+ default: {
+ throw UnexpectedError("Linear system library (" + ::name(library) + ") was not set!");
+ }
+ }
+ }
+
+ static void
+ checkHasLibrary(const LSLibrary library)
+ {
+ if (not hasLibrary(library)) {
+ throw NormalError(::name(library) + " is not linked to pugs. Cannot use it!");
+ }
+ }
+
+ static void
+ checkBuiltinMethod(const LSMethod method)
+ {
+ switch (method) {
+ case LSMethod::cg:
+ case LSMethod::bicgstab: {
+ break;
+ }
+ default: {
+ throw NormalError(name(method) + " is not a builtin linear solver!");
+ }
+ }
+ }
+
+ static void
+ checkPETScMethod(const LSMethod method)
+ {
+ switch (method) {
+ case LSMethod::cg:
+ case LSMethod::bicgstab:
+ case LSMethod::bicgstab2:
+ case LSMethod::gmres:
+ case LSMethod::lu:
+ case LSMethod::choleski: {
+ break;
+ }
+ default: {
+ throw NormalError(name(method) + " is not a builtin linear solver!");
+ }
+ }
+ }
+
+ static void
+ checkBuiltinPrecond(const LSPrecond precond)
+ {
+ switch (precond) {
+ case LSPrecond::none: {
+ break;
+ }
+ default: {
+ throw NormalError(name(precond) + " is not a builtin preconditioner!");
+ }
+ }
+ }
+
+ static void
+ checkPETScPrecond(const LSPrecond precond)
+ {
+ switch (precond) {
+ case LSPrecond::none:
+ case LSPrecond::diagonal:
+ case LSPrecond::incomplete_choleski:
+ case LSPrecond::incomplete_LU: {
+ break;
+ }
+ default: {
+ throw NormalError(name(precond) + " is not a PETSc preconditioner!");
+ }
+ }
+ }
+
+ static void
+ checkOptions(const LinearSolverOptions& options)
+ {
+ switch (options.library()) {
+ case LSLibrary::builtin: {
+ checkBuiltinMethod(options.method());
+ checkBuiltinPrecond(options.precond());
+ break;
+ }
+ case LSLibrary::petsc: {
+ checkPETScMethod(options.method());
+ checkPETScPrecond(options.precond());
+ break;
+ }
+ default: {
+ throw UnexpectedError("undefined options compatibility for this library (" + ::name(options.library()) + ")!");
+ }
+ }
+ }
+
+ static void
+ builtinSolveLocalSystem(const CRSMatrix<double, size_t>& A,
+ Vector<double>& x,
+ const Vector<double>& b,
+ const LinearSolverOptions& options)
+ {
+ if (options.precond() != LSPrecond::none) {
+ throw UnexpectedError("builtin linear solver do not allow any preconditioner!");
+ }
+ switch (options.method()) {
+ case LSMethod::cg: {
+ CG{A, x, b, options.epsilon(), options.maximumIteration(), options.verbose()};
+ break;
+ }
+ case LSMethod::bicgstab: {
+ BiCGStab{A, x, b, options.epsilon(), options.maximumIteration(), options.verbose()};
+ break;
+ }
+ default: {
+ throw NotImplementedError("undefined builtin method: " + name(options.method()));
+ }
+ }
+ }
+
+#ifdef PUGS_HAS_PETSC
+
+ static int
+ petscMonitor(KSP, int i, double residu, void*)
+ {
+ std::cout << " - iteration: " << std::setw(6) << i << " residu: " << std::scientific << residu << '\n';
+ return 0;
+ }
+
+ static void
+ petscSolveLocalSystem(const CRSMatrix<double, size_t>& A,
+ Vector<double>& x,
+ const Vector<double>& b,
+ const LinearSolverOptions& options)
+ {
+ Assert(x.size() == b.size() and x.size() == A.numberOfRows());
+
+ Vec petscB;
+ VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, b.size(), b.size(), &b[0], &petscB);
+ Vec petscX;
+ VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, x.size(), x.size(), &x[0], &petscX);
+
+ Mat petscMat;
+ MatCreate(PETSC_COMM_WORLD, &petscMat);
+ MatSetSizes(petscMat, PETSC_DECIDE, PETSC_DECIDE, x.size(), x.size());
+
+ MatSetType(petscMat, MATAIJ);
+
+ Array<PetscScalar> values = copy(A.values());
+ auto hm = A.hostMatrix();
+
+ Array<PetscInt> row_indices{hm.row_map.size()};
+ for (size_t i = 0; i < hm.row_map.size(); ++i) {
+ row_indices[i] = hm.row_map[i];
+ }
+
+ Array<PetscInt> column_indices{values.size()};
+ size_t l = 0;
+ for (size_t i = 0; i < A.numberOfRows(); ++i) {
+ const auto row_i = hm.rowConst(i);
+ for (size_t j = 0; j < row_i.length; ++j) {
+ column_indices[l++] = row_i.colidx(j);
+ }
+ }
+
+ MatCreateSeqAIJWithArrays(PETSC_COMM_WORLD, x.size(), x.size(), &row_indices[0], &column_indices[0], &values[0],
+ &petscMat);
+
+ MatAssemblyBegin(petscMat, MAT_FINAL_ASSEMBLY);
+ MatAssemblyEnd(petscMat, MAT_FINAL_ASSEMBLY);
+
+ KSP ksp;
+ KSPCreate(PETSC_COMM_WORLD, &ksp);
+ KSPSetTolerances(ksp, options.epsilon(), 1E-100, 1E5, options.maximumIteration());
+
+ KSPSetOperators(ksp, petscMat, petscMat);
+
+ PC pc;
+ KSPGetPC(ksp, &pc);
+
+ bool direct_solver = false;
+
+ switch (options.method()) {
+ case LSMethod::bicgstab: {
+ KSPSetType(ksp, KSPBCGS);
+ break;
+ }
+ case LSMethod::bicgstab2: {
+ KSPSetType(ksp, KSPBCGSL);
+ KSPBCGSLSetEll(ksp, 2);
+ break;
+ }
+ case LSMethod::cg: {
+ KSPSetType(ksp, KSPCG);
+ break;
+ }
+ case LSMethod::gmres: {
+ KSPSetType(ksp, KSPGMRES);
+
+ break;
+ }
+ case LSMethod::lu: {
+ KSPSetType(ksp, KSPPREONLY);
+ PCSetType(pc, PCLU);
+ PCFactorSetShiftType(pc, MAT_SHIFT_NONZERO);
+ direct_solver = true;
+ break;
+ }
+ case LSMethod::choleski: {
+ KSPSetType(ksp, KSPPREONLY);
+ PCSetType(pc, PCCHOLESKY);
+ direct_solver = true;
+ break;
+ }
+ default: {
+ throw UnexpectedError("unexpected method: " + name(options.method()));
+ }
+ }
+
+ if (not direct_solver) {
+ switch (options.precond()) {
+ case LSPrecond::diagonal: {
+ PCSetType(pc, PCJACOBI);
+ break;
+ }
+ case LSPrecond::incomplete_LU: {
+ PCSetType(pc, PCILU);
+ break;
+ }
+ case LSPrecond::incomplete_choleski: {
+ PCSetType(pc, PCICC);
+ break;
+ }
+ case LSPrecond::none: {
+ PCSetType(pc, PCNONE);
+ break;
+ }
+ default: {
+ throw UnexpectedError("unexpected preconditioner: " + name(options.precond()));
+ }
+ }
+ }
+ if (options.verbose()) {
+ KSPMonitorSet(ksp, petscMonitor, 0, 0);
+ }
+
+ KSPSolve(ksp, petscB, petscX);
+ MatDestroy(&petscMat);
+ }
+
+#else
+ static void
+ petscSolveLocalSystem(const CRSMatrix<double, size_t>&,
+ Vector<double>&,
+ const Vector<double>&,
+ const LinearSolverOptions&)
+ {
+ checkHasLibrary(LSLibrary::petsc);
+ throw UnexpectedError("unexpected situation should not reach this point!");
+ }
+
+#endif // PUGS_HAS_PETSC
+};
+
+bool
+LinearSolver::hasLibrary(LSLibrary library) const
+{
+ return Internals::hasLibrary(library);
+}
+
+void
+LinearSolver::solveLocalSystem(const CRSMatrix<double, size_t>& A, Vector<double>& x, const Vector<double>& b)
+{
+ switch (m_options.library()) {
+ case LSLibrary::builtin: {
+ Internals::builtinSolveLocalSystem(A, x, b, m_options);
+ break;
+ }
+ case LSLibrary::petsc: {
+ Internals::petscSolveLocalSystem(A, x, b, m_options);
+ break;
+ }
+ default: {
+ throw UnexpectedError(::name(m_options.library()) + " cannot solve local systems for sparse matrices");
+ }
+ }
+}
+
+LinearSolver::LinearSolver(const LinearSolverOptions& options) : m_options{options}
+{
+ Internals::checkHasLibrary(m_options.library());
+ Internals::checkOptions(options);
+}
+
+LinearSolver::LinearSolver() : LinearSolver{LinearSolverOptions::default_options} {}
diff --git a/src/algebra/LinearSolver.hpp b/src/algebra/LinearSolver.hpp
new file mode 100644
index 000000000..b72b09b0f
--- /dev/null
+++ b/src/algebra/LinearSolver.hpp
@@ -0,0 +1,36 @@
+#ifndef LINEAR_SOLVER_HPP
+#define LINEAR_SOLVER_HPP
+
+#include <algebra/CRSMatrix.hpp>
+#include <algebra/LinearSolverOptions.hpp>
+#include <algebra/TinyMatrix.hpp>
+#include <algebra/TinyVector.hpp>
+#include <algebra/Vector.hpp>
+#include <utils/Exceptions.hpp>
+
+class LinearSolver
+{
+ private:
+ struct Internals;
+
+ const LinearSolverOptions m_options;
+
+ void _solveLocalDense(size_t N, const double* A, double* x, const double* b);
+
+ public:
+ bool hasLibrary(LSLibrary library) const;
+
+ void solveLocalSystem(const CRSMatrix<double, size_t>& A, Vector<double>& x, const Vector<double>& b);
+
+ template <size_t N>
+ void
+ solveLocalSystem(const TinyMatrix<N>& A, TinyVector<N>& x, const TinyVector<N>& b)
+ {
+ this->_solveLocalDense(N, &A(0, 0), &x[0], &b[0]);
+ }
+
+ LinearSolver();
+ LinearSolver(const LinearSolverOptions& options);
+};
+
+#endif // LINEAR_SOLVER_HPP
diff --git a/src/algebra/LinearSolverOptions.cpp b/src/algebra/LinearSolverOptions.cpp
new file mode 100644
index 000000000..b40c5fa92
--- /dev/null
+++ b/src/algebra/LinearSolverOptions.cpp
@@ -0,0 +1,15 @@
+#include <algebra/LinearSolverOptions.hpp>
+
+#include <rang.hpp>
+
+std::ostream&
+operator<<(std::ostream& os, const LinearSolverOptions& options)
+{
+ os << " library: " << rang::style::bold << name(options.library()) << rang::style::reset << '\n';
+ os << " method : " << rang::style::bold << name(options.method()) << rang::style::reset << '\n';
+ os << " precond: " << rang::style::bold << name(options.precond()) << rang::style::reset << '\n';
+ os << " epsilon: " << rang::style::bold << options.epsilon() << rang::style::reset << '\n';
+ os << " maxiter: " << rang::style::bold << options.maximumIteration() << rang::style::reset << '\n';
+ os << " verbose: " << rang::style::bold << std::boolalpha << options.verbose() << rang::style::reset << '\n';
+ return os;
+}
diff --git a/src/algebra/LinearSolverOptions.hpp b/src/algebra/LinearSolverOptions.hpp
new file mode 100644
index 000000000..bd7538b53
--- /dev/null
+++ b/src/algebra/LinearSolverOptions.hpp
@@ -0,0 +1,233 @@
+#ifndef LINEAR_SOLVER_OPTIONS_HPP
+#define LINEAR_SOLVER_OPTIONS_HPP
+
+#include <utils/Exceptions.hpp>
+
+#include <iostream>
+
+enum class LSLibrary : int8_t
+{
+ LS__begin = 0,
+ //
+ builtin = LS__begin,
+ petsc,
+ //
+ LS__end
+};
+
+enum class LSMethod : int8_t
+{
+ LS__begin = 0,
+ //
+ cg = LS__begin,
+ bicgstab,
+ bicgstab2,
+ gmres,
+ lu,
+ choleski,
+ //
+ LS__end
+};
+
+enum class LSPrecond : int8_t
+{
+ LS__begin = 0,
+ //
+ none = LS__begin,
+ diagonal,
+ incomplete_choleski,
+ incomplete_LU,
+ //
+ LS__end
+};
+
+inline std::string
+name(const LSLibrary library)
+{
+ switch (library) {
+ case LSLibrary::builtin: {
+ return "builtin";
+ }
+ case LSLibrary::petsc: {
+ return "PETSc";
+ }
+ case LSLibrary::LS__end: {
+ }
+ }
+ throw UnexpectedError("Linear system library name is not defined!");
+}
+
+inline std::string
+name(const LSMethod method)
+{
+ switch (method) {
+ case LSMethod::cg: {
+ return "CG";
+ }
+ case LSMethod::bicgstab: {
+ return "BICGStab";
+ }
+ case LSMethod::bicgstab2: {
+ return "BICGStab2";
+ }
+ case LSMethod::gmres: {
+ return "GMRES";
+ }
+ case LSMethod::lu: {
+ return "LU";
+ }
+ case LSMethod::choleski: {
+ return "Choleski";
+ }
+ case LSMethod::LS__end: {
+ }
+ }
+ throw UnexpectedError("Linear system method name is not defined!");
+}
+
+inline std::string
+name(const LSPrecond precond)
+{
+ switch (precond) {
+ case LSPrecond::none: {
+ return "none";
+ }
+ case LSPrecond::diagonal: {
+ return "diagonal";
+ }
+ case LSPrecond::incomplete_choleski: {
+ return "ICholeski";
+ }
+ case LSPrecond::incomplete_LU: {
+ return "ILU";
+ }
+ case LSPrecond::LS__end: {
+ }
+ }
+ throw UnexpectedError("Linear system preconditioner name is not defined!");
+}
+
+template <typename LSEnumType>
+inline LSEnumType
+getLSEnumFromName(const std::string& enum_name)
+{
+ using BaseT = std::underlying_type_t<LSEnumType>;
+ for (BaseT enum_value = static_cast<BaseT>(LSEnumType::LS__begin);
+ enum_value < static_cast<BaseT>(LSEnumType::LS__end); ++enum_value) {
+ if (name(LSEnumType{enum_value}) == enum_name) {
+ return LSEnumType{enum_value};
+ }
+ }
+ throw NormalError(std::string{"could not find '"} + enum_name + "' associate type!");
+}
+
+template <typename LSEnumType>
+inline void
+printLSEnumListNames(std::ostream& os)
+{
+ using BaseT = std::underlying_type_t<LSEnumType>;
+ for (BaseT enum_value = static_cast<BaseT>(LSEnumType::LS__begin);
+ enum_value < static_cast<BaseT>(LSEnumType::LS__end); ++enum_value) {
+ os << " - " << name(LSEnumType{enum_value}) << '\n';
+ }
+}
+
+class LinearSolverOptions
+{
+ private:
+ LSLibrary m_library = LSLibrary::builtin;
+ LSMethod m_method = LSMethod::bicgstab;
+ LSPrecond m_precond = LSPrecond::none;
+
+ double m_epsilon = 1E-6;
+ size_t m_maximum_iteration = 200;
+
+ bool m_verbose = false;
+
+ public:
+ static LinearSolverOptions default_options;
+
+ friend std::ostream& operator<<(std::ostream& os, const LinearSolverOptions& options);
+
+ LSLibrary&
+ library()
+ {
+ return m_library;
+ }
+
+ LSLibrary
+ library() const
+ {
+ return m_library;
+ }
+
+ LSMethod
+ method() const
+ {
+ return m_method;
+ }
+
+ LSMethod&
+ method()
+ {
+ return m_method;
+ }
+
+ LSPrecond
+ precond() const
+ {
+ return m_precond;
+ }
+
+ LSPrecond&
+ precond()
+ {
+ return m_precond;
+ }
+
+ double
+ epsilon() const
+ {
+ return m_epsilon;
+ }
+
+ double&
+ epsilon()
+ {
+ return m_epsilon;
+ }
+
+ size_t&
+ maximumIteration()
+ {
+ return m_maximum_iteration;
+ }
+
+ size_t
+ maximumIteration() const
+ {
+ return m_maximum_iteration;
+ }
+
+ bool&
+ verbose()
+ {
+ return m_verbose;
+ };
+
+ bool
+ verbose() const
+ {
+ return m_verbose;
+ };
+
+ LinearSolverOptions(const LinearSolverOptions&) = default;
+ LinearSolverOptions(LinearSolverOptions&&) = default;
+
+ LinearSolverOptions() = default;
+ ~LinearSolverOptions() = default;
+};
+
+inline LinearSolverOptions LinearSolverOptions::default_options;
+
+#endif // LINEAR_SOLVER_OPTIONS_HPP
diff --git a/src/algebra/LinearSystemSolver.cpp b/src/algebra/LinearSystemSolver.cpp
deleted file mode 100644
index 488d6763d..000000000
--- a/src/algebra/LinearSystemSolver.cpp
+++ /dev/null
@@ -1,331 +0,0 @@
-#include <algebra/LinearSystemSolver.hpp>
-#include <utils/pugs_config.hpp>
-
-#ifdef PUGS_HAS_PETSC
-#include <petsc.h>
-#endif // PUGS_HAS_PETSC
-
-struct LinearSystemSolver::Internals
-{
- static bool
- hasLibrary(LSLibrary library)
- {
- switch (library) {
- case LSLibrary::builtin: {
- return true;
- }
- case LSLibrary::petsc: {
-#ifdef PUGS_HAS_PETSC
- return true;
-#else
- return false;
-#endif
- }
- default: {
- throw UnexpectedError("Linear system library (" + name(library) + ") was not set!");
- }
- }
- }
-
- static void
- checkHasLibrary(LSLibrary library)
- {
- if (not hasLibrary(library)) {
- throw NormalError(name(library) + " is not linked to pugs. Cannot use it!");
- }
- }
-
- static std::string
- name(LSLibrary library)
- {
- switch (library) {
- case LSLibrary::builtin: {
- return "builtin";
- }
- case LSLibrary::petsc: {
- return "PETSc";
- }
- default: {
- throw UnexpectedError("Linear system library name is not defined!");
- }
- }
- }
-
- static std::string
- name(LSMethod method)
- {
- switch (method) {
- case LSMethod::cg: {
- return "conjugate gradient";
- }
- case LSMethod::bicgstab: {
- return "bicg-stab";
- }
- case LSMethod::gmres: {
- return "gmres";
- }
- default: {
- throw UnexpectedError("Linear system method name is not defined!");
- }
- }
- }
-
- static std::string
- name(LSPrecond precond)
- {
- switch (precond) {
- case LSPrecond::none: {
- return "none";
- }
- case LSPrecond::diagonal: {
- return "Jacobi (diagonal)";
- }
- case LSPrecond::incomplete_choleski: {
- return "incomplete Choleski";
- }
- case LSPrecond::incomplete_LU: {
- return "incomplete LU";
- }
- default: {
- throw UnexpectedError("Linear system preconditioner name is not defined!");
- }
- }
- }
-
- static void
- checkBuiltinMethod(LSMethod method)
- {
- switch (method) {
- case LSMethod::cg:
- case LSMethod::bicgstab: {
- break;
- }
- default: {
- throw NormalError(name(method) + " is not a builtin linear solver!");
- }
- }
- }
-
- static void
- checkPETScMethod(LSMethod method)
- {
- switch (method) {
- case LSMethod::cg:
- case LSMethod::bicgstab:
- case LSMethod::gmres: {
- break;
- }
- default: {
- throw NormalError(name(method) + " is not a builtin linear solver!");
- }
- }
- }
-
- static void
- checkBuiltinPrecond(LSPrecond precond)
- {
- switch (precond) {
- case LSPrecond::none: {
- break;
- }
- default: {
- throw NormalError(name(precond) + " is not a builtin preconditioner!");
- }
- }
- }
-
- static void
- checkPETScPrecond(LSPrecond precond)
- {
- switch (precond) {
- case LSPrecond::none:
- case LSPrecond::diagonal:
- case LSPrecond::incomplete_choleski:
- case LSPrecond::incomplete_LU: {
- break;
- }
- default: {
- throw NormalError(name(precond) + " is not a PETSc preconditioner!");
- }
- }
- }
-
- static void
- checkOptions(LSLibrary library, LSMethod method, LSPrecond precond)
- {
- switch (library) {
- case LSLibrary::builtin: {
- checkBuiltinMethod(method);
- checkBuiltinPrecond(precond);
- break;
- }
- case LSLibrary::petsc: {
- checkPETScMethod(method);
- checkPETScPrecond(precond);
- break;
- }
- default: {
- throw UnexpectedError("Undefined options compatibility for this library (" + name(library) + ")!");
- }
- }
- }
-
-#ifdef PUGS_HAS_PETSC
-
- static int
- petscMonitor(KSP, int it, double norm, void*)
- {
- std::cout << "iteration: " << it << " residual norm = " << norm << '\n';
- return 0;
- }
-
- static void
- petscSolveLocalSystem(const SparseMatrixDescriptor<double, size_t>& A,
- Vector<double>& x,
- const Vector<double>& b,
- LSMethod method,
- LSPrecond precond,
- double epsilon,
- size_t maximum_iteration)
- {
- Assert(x.size() == b.size() and x.size() == A.numberOfRows());
-
- Vec petscB;
- VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, b.size(), b.size(), &b[0], &petscB);
- Vec petscX;
- VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, x.size(), x.size(), &x[0], &petscX);
-
- Mat petscMat;
- MatCreate(PETSC_COMM_WORLD, &petscMat);
- MatSetSizes(petscMat, PETSC_DECIDE, PETSC_DECIDE, x.size(), x.size());
-
- MatSetType(petscMat, MATAIJ);
- MatSetUp(petscMat);
-
- for (size_t i = 0; i < A.numberOfRows(); ++i) {
- const auto& row = A.row(i);
- PetscInt I = i;
- for (auto [j, value] : row.columnIdToValueMap()) {
- PetscInt J = j;
- MatSetValues(petscMat, 1, &I, 1, &J, &value, INSERT_VALUES);
- }
- }
-
- MatAssemblyBegin(petscMat, MAT_FINAL_ASSEMBLY);
- MatAssemblyEnd(petscMat, MAT_FINAL_ASSEMBLY);
-
- PetscPrintf(PETSC_COMM_WORLD, "Matrix _before_ MatScaleUserImpl() operation\n");
- MatView(petscMat, PETSC_VIEWER_STDOUT_SELF);
-
- KSP ksp;
-
- KSPCreate(PETSC_COMM_WORLD, &ksp);
- KSPSetTolerances(ksp, epsilon, 1E-100, 1E5, maximum_iteration);
-
- KSPSetOperators(ksp, petscMat, petscMat);
-
- switch (method) {
- case LSMethod::bicgstab: {
- KSPSetType(ksp, KSPBCGS);
- break;
- }
- case LSMethod::cg: {
- KSPSetType(ksp, KSPCG);
- break;
- }
- case LSMethod::gmres: {
- KSPSetType(ksp, KSPGMRES);
- break;
- }
- default: {
- throw UnexpectedError("unexpected method: " + name(method));
- }
- }
-
- PC pc;
- KSPGetPC(ksp, &pc);
-
- switch (precond) {
- case LSPrecond::diagonal: {
- PCSetType(pc, PCJACOBI);
- break;
- }
- case LSPrecond::incomplete_LU: {
- PCSetType(pc, PCILU);
- break;
- }
- case LSPrecond::incomplete_choleski: {
- PCSetType(pc, PCICC);
- break;
- }
- case LSPrecond::none: {
- PCSetType(pc, PCNONE);
- break;
- }
- default: {
- throw UnexpectedError("unexpected preconditioner: " + name(precond));
- }
- }
-
- KSPMonitorSet(ksp, petscMonitor, 0, 0);
-
- KSPSolve(ksp, petscB, petscX);
-
- VecView(petscX, PETSC_VIEWER_STDOUT_SELF);
-
- MatDestroy(&petscMat);
- }
-
-#else
-
- static void
- petscSolveLocalSystem(const SparseMatrixDescriptor<double, size_t>&,
- Vector<double>&,
- const Vector<double>&,
- LSMethod,
- LSPrecond,
- double,
- size_t)
- {
- checkHasLibrary(LSLibrary::petsc);
- }
-
-#endif // PUGS_HAS_PETSC
-};
-
-bool
-LinearSystemSolver::hasLibrary(LSLibrary library) const
-{
- return Internals::hasLibrary(library);
-}
-
-LinearSystemSolver::LinearSystemSolver(LSLibrary library,
- LSMethod method,
- LSPrecond precond,
- double epsilon,
- size_t maximum_iteration)
- : m_library{library}, m_method{method}, m_precond{precond}, m_epsilon{epsilon}, m_maximum_iteration{maximum_iteration}
-{
- Internals::checkHasLibrary(library);
- Internals::checkOptions(library, method, precond);
-}
-
-void
-LinearSystemSolver::solveLocalSystem(const SparseMatrixDescriptor<double, size_t>& A,
- Vector<double>& x,
- const Vector<double>& b)
-{
- switch (m_library) {
- case LSLibrary::builtin: {
- throw NotImplementedError("builtin solvers are not plugged yet!");
- break;
- }
- case LSLibrary::petsc: {
- Internals::petscSolveLocalSystem(A, x, b, m_method, m_precond, m_epsilon, m_maximum_iteration);
- break;
- }
- default: {
- throw UnexpectedError(Internals::name(m_library) + " cannot solve local systems for sparse matrices");
- }
- }
-}
diff --git a/src/algebra/LinearSystemSolver.hpp b/src/algebra/LinearSystemSolver.hpp
deleted file mode 100644
index f689f96dd..000000000
--- a/src/algebra/LinearSystemSolver.hpp
+++ /dev/null
@@ -1,79 +0,0 @@
-#ifndef LINEAR_SYSTEM_SOLVER_HPP
-#define LINEAR_SYSTEM_SOLVER_HPP
-
-#include <algebra/SparseMatrixDescriptor.hpp>
-#include <algebra/TinyMatrix.hpp>
-#include <algebra/TinyVector.hpp>
-#include <algebra/Vector.hpp>
-
-#include <utils/Exceptions.hpp>
-
-enum class LSLibrary
-{
- builtin,
- petsc
-};
-
-enum class LSMethod
-{
- cg,
- bicgstab,
- gmres
-};
-
-enum class LSPrecond
-{
- none,
- diagonal,
- incomplete_choleski,
- incomplete_LU
-};
-
-class LinearSystemSolver
-{
- private:
- struct Internals;
-
- const LSLibrary m_library;
- const LSMethod m_method;
- const LSPrecond m_precond;
-
- const double m_epsilon;
- const size_t m_maximum_iteration;
-
- bool m_is_verbose = false;
-
- void _solveLocalDense(size_t N, const double* A, double* x, const double* b);
-
- public:
- bool
- isVerbose() const
- {
- return m_is_verbose;
- }
-
- void
- setVerbose(bool verbose)
- {
- m_is_verbose = verbose;
- }
-
- bool hasLibrary(LSLibrary library) const;
-
- void solveLocalSystem(const SparseMatrixDescriptor<double, size_t>& A, Vector<double>& x, const Vector<double>& b);
-
- template <size_t N>
- void
- solveLocalSystem(const TinyMatrix<N>& A, TinyVector<N>& x, const TinyVector<N>& b)
- {
- this->_solveLocalDense(N, &A(0, 0), &x[0], &b[0]);
- }
-
- LinearSystemSolver(LSLibrary library,
- LSMethod method,
- LSPrecond precond,
- double epsilon = 1e-6,
- size_t maximum_iteration = 100);
-};
-
-#endif // LINEAR_SYSTEM_SOLVER_HPP
diff --git a/src/algebra/PETScWrapper.cpp b/src/algebra/PETScWrapper.cpp
index 76e56692a..dd11dc977 100644
--- a/src/algebra/PETScWrapper.cpp
+++ b/src/algebra/PETScWrapper.cpp
@@ -12,6 +12,7 @@ void
initialize([[maybe_unused]] int& argc, [[maybe_unused]] char* argv[])
{
#ifdef PUGS_HAS_PETSC
+ PetscOptionsSetValue(NULL, "-no_signal_handler", "true");
PetscInitialize(&argc, &argv, 0, 0);
#endif // PUGS_HAS_PETSC
}
diff --git a/src/algebra/SparseMatrixDescriptor.hpp b/src/algebra/SparseMatrixDescriptor.hpp
index 09f13085b..d5c14713b 100644
--- a/src/algebra/SparseMatrixDescriptor.hpp
+++ b/src/algebra/SparseMatrixDescriptor.hpp
@@ -151,7 +151,12 @@ class SparseMatrixDescriptor
return values;
}
- SparseMatrixDescriptor(size_t nb_row) : m_row_array{nb_row} {}
+ SparseMatrixDescriptor(size_t nb_row) : m_row_array{nb_row}
+ {
+ for (size_t i = 0; i < nb_row; ++i) {
+ m_row_array[i][i] = 0;
+ }
+ }
~SparseMatrixDescriptor() = default;
};
diff --git a/src/language/ast/ASTNodeDataType.cpp b/src/language/ast/ASTNodeDataType.cpp
index 10ec58b45..ce2bbcc60 100644
--- a/src/language/ast/ASTNodeDataType.cpp
+++ b/src/language/ast/ASTNodeDataType.cpp
@@ -48,11 +48,15 @@ dataTypeName(const ASTNodeDataType& data_type)
case ASTNodeDataType::list_t: {
std::ostringstream data_type_name_list;
const auto& data_type_list = data_type.contentTypeList();
- data_type_name_list << dataTypeName(*data_type_list[0]);
- for (size_t i = 1; i < data_type_list.size(); ++i) {
- data_type_name_list << '*' << dataTypeName(*data_type_list[i]);
+ if (data_type_list.size() > 0) {
+ data_type_name_list << dataTypeName(*data_type_list[0]);
+ for (size_t i = 1; i < data_type_list.size(); ++i) {
+ data_type_name_list << '*' << dataTypeName(*data_type_list[i]);
+ }
+ name = "list(" + data_type_name_list.str() + ")";
+ } else {
+ name = "list(void)";
}
- name = "list(" + data_type_name_list.str() + ")";
break;
}
case ASTNodeDataType::string_t:
diff --git a/src/language/modules/CMakeLists.txt b/src/language/modules/CMakeLists.txt
index f8cf2f0b4..0a996fd67 100644
--- a/src/language/modules/CMakeLists.txt
+++ b/src/language/modules/CMakeLists.txt
@@ -2,6 +2,7 @@
add_library(PugsLanguageModules
BuiltinModule.cpp
+ LinearSolverModule.cpp
MathModule.cpp
MeshModule.cpp
ModuleRepository.cpp
diff --git a/src/language/modules/LinearSolverModule.cpp b/src/language/modules/LinearSolverModule.cpp
new file mode 100644
index 000000000..9e626ba2b
--- /dev/null
+++ b/src/language/modules/LinearSolverModule.cpp
@@ -0,0 +1,85 @@
+#include <language/modules/LinearSolverModule.hpp>
+
+#include <algebra/LinearSolver.hpp>
+#include <language/utils/BuiltinFunctionEmbedder.hpp>
+#include <language/utils/TypeDescriptor.hpp>
+
+LinearSolverModule::LinearSolverModule()
+{
+ this->_addBuiltinFunction("setLSVerbosity", std::make_shared<BuiltinFunctionEmbedder<void(const bool&)>>(
+
+ [](const bool& verbose) -> void {
+ LinearSolverOptions::default_options.verbose() = verbose;
+ }
+
+ ));
+
+ this->_addBuiltinFunction("setLSEpsilon", std::make_shared<BuiltinFunctionEmbedder<void(const double&)>>(
+
+ [](const double& epsilon) -> void {
+ LinearSolverOptions::default_options.epsilon() = epsilon;
+ }
+
+ ));
+
+ this->_addBuiltinFunction("setLSMaxIter", std::make_shared<BuiltinFunctionEmbedder<void(const uint64_t&)>>(
+
+ [](const uint64_t& max_iter) -> void {
+ LinearSolverOptions::default_options.maximumIteration() = max_iter;
+ }
+
+ ));
+
+ this->_addBuiltinFunction("setLSLibrary", std::make_shared<BuiltinFunctionEmbedder<void(const std::string&)>>(
+
+ [](const std::string& library_name) -> void {
+ LinearSolverOptions::default_options.library() =
+ getLSEnumFromName<LSLibrary>(library_name);
+ }
+
+ ));
+
+ this->_addBuiltinFunction("setLSMethod", std::make_shared<BuiltinFunctionEmbedder<void(const std::string&)>>(
+
+ [](const std::string& method_name) -> void {
+ LinearSolverOptions::default_options.method() =
+ getLSEnumFromName<LSMethod>(method_name);
+ }
+
+ ));
+
+ this->_addBuiltinFunction("setLSPrecond", std::make_shared<BuiltinFunctionEmbedder<void(const std::string&)>>(
+
+ [](const std::string& precond_name) -> void {
+ LinearSolverOptions::default_options.precond() =
+ getLSEnumFromName<LSPrecond>(precond_name);
+ }
+
+ ));
+
+ this->_addBuiltinFunction("printLSOptions", std::make_shared<BuiltinFunctionEmbedder<void()>>(
+
+ []() -> void {
+ std::cout << rang::fgB::yellow << "Linear solver options"
+ << rang::style::reset << '\n';
+ std::cout << LinearSolverOptions::default_options;
+ }
+
+ ));
+
+ this->_addBuiltinFunction("printLSAvailable",
+ std::make_shared<BuiltinFunctionEmbedder<void()>>(
+
+ []() -> void {
+ std::cout << rang::fgB::yellow << "Available linear solver options"
+ << rang::style::reset << '\n';
+ std::cout << rang::fgB::blue << " libraries" << rang::style::reset << '\n';
+ printLSEnumListNames<LSLibrary>(std::cout);
+ std::cout << rang::fgB::blue << " methods" << rang::style::reset << '\n';
+ printLSEnumListNames<LSMethod>(std::cout);
+ std::cout << rang::fgB::blue << " preconditioners" << rang::style::reset << '\n';
+ printLSEnumListNames<LSPrecond>(std::cout);
+ }
+
+ ));
+}
diff --git a/src/language/modules/LinearSolverModule.hpp b/src/language/modules/LinearSolverModule.hpp
new file mode 100644
index 000000000..5dc7ae64e
--- /dev/null
+++ b/src/language/modules/LinearSolverModule.hpp
@@ -0,0 +1,19 @@
+#ifndef LINEAR_SOLVER_MODULE_HPP
+#define LINEAR_SOLVER_MODULE_HPP
+
+#include <language/modules/BuiltinModule.hpp>
+
+class LinearSolverModule : public BuiltinModule
+{
+ public:
+ std::string_view
+ name() const final
+ {
+ return "linear_solver";
+ }
+
+ LinearSolverModule();
+ ~LinearSolverModule() = default;
+};
+
+#endif // LINEAR_SOLVER_MODULE_HPP
diff --git a/src/language/modules/ModuleRepository.cpp b/src/language/modules/ModuleRepository.cpp
index 497ba083a..df403fa87 100644
--- a/src/language/modules/ModuleRepository.cpp
+++ b/src/language/modules/ModuleRepository.cpp
@@ -1,6 +1,7 @@
#include <language/modules/ModuleRepository.hpp>
#include <language/ast/ASTNode.hpp>
+#include <language/modules/LinearSolverModule.hpp>
#include <language/modules/MathModule.hpp>
#include <language/modules/MeshModule.hpp>
#include <language/modules/SchemeModule.hpp>
@@ -23,6 +24,7 @@ ModuleRepository::ModuleRepository()
this->_subscribe(std::make_unique<MeshModule>());
this->_subscribe(std::make_unique<VTKModule>());
this->_subscribe(std::make_unique<SchemeModule>());
+ this->_subscribe(std::make_unique<LinearSolverModule>());
}
template <typename NameEmbedderMapT, typename EmbedderTableT>
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 8b2a7c551..8d0344a03 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -47,7 +47,7 @@ add_executable (unit_tests
test_BuiltinFunctionEmbedder.cpp
test_BuiltinFunctionEmbedderTable.cpp
test_BuiltinFunctionProcessor.cpp
- test_MathModule.cpp
+ test_CG.cpp
test_ContinueProcessor.cpp
test_ConcatExpressionProcessor.cpp
test_CRSMatrix.cpp
@@ -65,9 +65,9 @@ add_executable (unit_tests
test_INodeProcessor.cpp
test_ItemType.cpp
test_ListAffectationProcessor.cpp
+ test_MathModule.cpp
test_NameProcessor.cpp
test_OStreamProcessor.cpp
- test_PCG.cpp
test_PugsFunctionAdapter.cpp
test_PugsAssert.cpp
test_RevisionInfo.cpp
@@ -98,6 +98,7 @@ target_link_libraries (unit_tests
kokkos
${PARMETIS_LIBRARIES}
${MPI_CXX_LINK_FLAGS} ${MPI_CXX_LIBRARIES}
+ ${PETSC_LIBRARIES}
Catch2
)
@@ -108,6 +109,7 @@ target_link_libraries (mpi_unit_tests
kokkos
${PARMETIS_LIBRARIES}
${MPI_CXX_LINK_FLAGS} ${MPI_CXX_LIBRARIES}
+ ${PETSC_LIBRARIES}
Catch2
)
diff --git a/tests/test_BiCGStab.cpp b/tests/test_BiCGStab.cpp
index 99b94515b..49a6955be 100644
--- a/tests/test_BiCGStab.cpp
+++ b/tests/test_BiCGStab.cpp
@@ -45,7 +45,7 @@ TEST_CASE("BiCGStab", "[algebra]")
Vector<double> x{5};
x = 0;
- BiCGStab<false>{b, A, x, 10, 1e-12};
+ BiCGStab{A, x, b, 1e-12, 10, false};
Vector error = x - x_exact;
REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x, x)));
}
@@ -88,7 +88,7 @@ TEST_CASE("BiCGStab", "[algebra]")
Vector<double> x{5};
x = 0;
- BiCGStab{b, A, x, 1, 1e-12};
+ BiCGStab{A, x, b, 1e-12, 1, true};
Vector error = x - x_exact;
REQUIRE(std::sqrt((error, error)) > 1E-5 * std::sqrt((x, x)));
}
diff --git a/tests/test_PCG.cpp b/tests/test_CG.cpp
similarity index 91%
rename from tests/test_PCG.cpp
rename to tests/test_CG.cpp
index d4564761d..2ce9febb0 100644
--- a/tests/test_PCG.cpp
+++ b/tests/test_CG.cpp
@@ -1,11 +1,11 @@
#include <catch2/catch.hpp>
+#include <algebra/CG.hpp>
#include <algebra/CRSMatrix.hpp>
-#include <algebra/PCG.hpp>
// clazy:excludeall=non-pod-global-static
-TEST_CASE("PCG", "[algebra]")
+TEST_CASE("CG", "[algebra]")
{
SECTION("no preconditionner")
{
@@ -45,7 +45,7 @@ TEST_CASE("PCG", "[algebra]")
Vector<double> x{5};
x = 0;
- PCG{b, A, A, x, 10, 1e-12};
+ CG{A, x, b, 1e-12, 10};
Vector error = x - x_exact;
REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x, x)));
}
@@ -62,7 +62,7 @@ TEST_CASE("PCG", "[algebra]")
Vector<double> x{5};
x = 0;
- PCG<false>{b, A, A, x, 10, 1e-12};
+ CG{A, x, b, 1e-12, 10, false};
REQUIRE(std::sqrt((x, x)) == 0);
}
@@ -104,7 +104,7 @@ TEST_CASE("PCG", "[algebra]")
Vector<double> x{5};
x = 0;
- PCG<false>{b, A, A, x, 1, 1e-12};
+ CG{A, x, b, 1e-12, 1, false};
Vector error = x - x_exact;
REQUIRE(std::sqrt((error, error)) > 1E-10 * std::sqrt((x, x)));
}
diff --git a/tests/test_SparseMatrixDescriptor.cpp b/tests/test_SparseMatrixDescriptor.cpp
index cd66380b1..f90cc83b3 100644
--- a/tests/test_SparseMatrixDescriptor.cpp
+++ b/tests/test_SparseMatrixDescriptor.cpp
@@ -68,18 +68,18 @@ TEST_CASE("SparseMatrixDescriptor", "[algebra]")
S(4, 1) = 1;
S(4, 4) = -2;
- REQUIRE(S.row(0).numberOfValues() == 1);
+ REQUIRE(S.row(0).numberOfValues() == 2);
REQUIRE(S.row(1).numberOfValues() == 2);
REQUIRE(S.row(2).numberOfValues() == 1);
- REQUIRE(S.row(3).numberOfValues() == 1);
+ REQUIRE(S.row(3).numberOfValues() == 2);
REQUIRE(S.row(4).numberOfValues() == 2);
const auto& const_S = S;
- REQUIRE(const_S.row(0).numberOfValues() == 1);
+ REQUIRE(const_S.row(0).numberOfValues() == 2);
REQUIRE(const_S.row(1).numberOfValues() == 2);
REQUIRE(const_S.row(2).numberOfValues() == 1);
- REQUIRE(const_S.row(3).numberOfValues() == 1);
+ REQUIRE(const_S.row(3).numberOfValues() == 2);
REQUIRE(const_S.row(4).numberOfValues() == 2);
#ifndef NDEBUG
@@ -126,7 +126,7 @@ TEST_CASE("SparseMatrixDescriptor", "[algebra]")
const auto graph = S.graphVector();
REQUIRE(graph.size() == S.numberOfRows());
- REQUIRE(graph[0].size() == 1);
+ REQUIRE(graph[0].size() == 2);
REQUIRE(graph[1].size() == 2);
REQUIRE(graph[2].size() == 1);
REQUIRE(graph[3].size() == 2);
@@ -157,7 +157,7 @@ TEST_CASE("SparseMatrixDescriptor", "[algebra]")
const auto value_array = S.valueArray();
- REQUIRE(value_array.size() == 8);
+ REQUIRE(value_array.size() == 9);
REQUIRE(value_array[0] == 5);
REQUIRE(value_array[1] == 1);
@@ -186,7 +186,7 @@ TEST_CASE("SparseMatrixDescriptor", "[algebra]")
output << '\n' << S;
std::string expected_output = R"(
-0 | 2:5
+0 | 0:0 2:5
1 | 1:1 2:11
2 | 2:4
3 | 1:-3 3:5
--
GitLab