diff --git a/CMakeLists.txt b/CMakeLists.txt
index 6b4f0a54d493750762038fb9ca4230b4293f5655..dc0f693980d6771d005c961cfae8f023ae81d174 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 bc822668a7f63c4d5d7cc4c45f4d1b1d70084fc0..ff34b89e148d844215cf7fffc21c2ad78bca5ec9 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 36a4635b313d2bba13d97e1b7a4ea1882e2a66a2..538c77571dac4d551d7476ea79bf855e8026fb53 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 bbd2d01bf417baaacb13c6c7362a825337a180be..d1f01ec310858be0e775514f934ce0a7c0e9da28 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 97a0a8a8ed44d660f9524eac269abee0324cd136..d9b03ac93fb2f1c4baeff98c90daac17547ff7cf 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 0000000000000000000000000000000000000000..a275cd5816d6ce04715fd1d160008b4dd8508051
--- /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 0000000000000000000000000000000000000000..b72b09b0f6a3255f914c96ee3df7717923931294
--- /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 0000000000000000000000000000000000000000..b40c5fa92f3af73015afa133eaae99b91feb6fda
--- /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 0000000000000000000000000000000000000000..bd7538b53df8740ee6c151e0d2e72a106fdc3c8c
--- /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 488d6763d1838cc4d5edc596589310dc2ef66038..0000000000000000000000000000000000000000
--- 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 f689f96ddba12ca8f7f4e4c19bbd78621dd1981f..0000000000000000000000000000000000000000
--- 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 76e56692a67e86bb61f6333c09680a30bc9e07ae..dd11dc9772f11ad29ca4970c5b88656c623a3508 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 09f13085b4881d2d602d41a42e0d8d0918cca633..d5c14713bb6b91235b796a56f40b9a2f5a605805 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 10ec58b451ead8402a328744cb977c29557dbbbf..ce2bbcc60182e0fe025d4836ee43fe0cb6d37a59 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 f8cf2f0b445882a672d54e60dbc75858bbb246c7..0a996fd6764b4d2eed57c02e8682283f1e640901 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 0000000000000000000000000000000000000000..9e626ba2b0e4707df0495de2057a00d81ea4b118
--- /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 0000000000000000000000000000000000000000..5dc7ae64efd3ed43f554e821693f761b0a66d412
--- /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 497ba083ac2d9d2ca0359c5d73e85e094b1fed0e..df403fa87f0c0d39778414a92ec7aa991b3d8d5b 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 8b2a7c5516394444debece5e32413dde7804c56d..8d0344a03785fc30e30a20b71b498d7cd2b28d4a 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 99b94515ba7795903434db761d509b9052fcaf71..49a6955be4ae1da7c45e036977a3036be0ccf6d3 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 d4564761d8c564082e4668360e5a2e2a8ecce008..2ce9febb0a31b087d2a7651e4442245b9d6c51ce 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 cd66380b1d1beb6bcb9ba1a9c0d9c495947a3235..f90cc83b3778e39a3319824f53a34e2041dadcf7 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