diff --git a/src/algebra/DenseMatrixWrapper.hpp b/src/algebra/DenseMatrixWrapper.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..dee57f27eff3d4c906cb256754d86ca6eeafbdfa
--- /dev/null
+++ b/src/algebra/DenseMatrixWrapper.hpp
@@ -0,0 +1,59 @@
+#ifndef DENSE_MATRIX_WRAPPER_HPP
+#define DENSE_MATRIX_WRAPPER_HPP
+
+#include <algebra/SmallMatrix.hpp>
+#include <algebra/TinyMatrix.hpp>
+
+template <typename DataType>
+class DenseMatrixWrapper
+{
+ private:
+ const size_t m_number_of_rows;
+ const size_t m_number_of_columns;
+ const DataType* const m_matrix_ptr;
+
+ public:
+ PUGS_INLINE
+ bool
+ isSquare() const
+ {
+ return (m_number_of_rows == m_number_of_columns);
+ }
+
+ PUGS_INLINE
+ size_t
+ numberOfRows() const
+ {
+ return m_number_of_rows;
+ }
+
+ PUGS_INLINE
+ size_t
+ numberOfColumns() const
+ {
+ return m_number_of_columns;
+ }
+
+ PUGS_INLINE
+ const DataType* const&
+ ptr() const
+ {
+ return m_matrix_ptr;
+ }
+
+ DenseMatrixWrapper(const SmallMatrix<DataType>& A)
+ : m_number_of_rows{A.numberOfRows()},
+ m_number_of_columns{A.numberOfColumns()},
+ m_matrix_ptr{(m_number_of_columns * m_number_of_rows > 0) ? (&A(0, 0)) : nullptr}
+ {}
+
+ template <size_t M, size_t N>
+ DenseMatrixWrapper(const TinyMatrix<M, N, DataType>& A)
+ : m_number_of_rows{M}, m_number_of_columns{N}, m_matrix_ptr{(M * N > 0) ? (&A(0, 0)) : nullptr}
+ {}
+
+ DenseMatrixWrapper(const DenseMatrixWrapper&) = delete;
+ DenseMatrixWrapper(DenseMatrixWrapper&&) = delete;
+};
+
+#endif // DENSE_MATRIX_WRAPPER_HPP
diff --git a/src/algebra/Eigen3Utils.hpp b/src/algebra/Eigen3Utils.hpp
index 5224a4b0bb10dc252b0f299cc9cbbbfe9fa87712..d9d61e02746c77e22002fb1c1ecfb93db6c1a9bb 100644
--- a/src/algebra/Eigen3Utils.hpp
+++ b/src/algebra/Eigen3Utils.hpp
@@ -6,8 +6,7 @@
#ifdef PUGS_HAS_EIGEN3
#include <algebra/CRSMatrix.hpp>
-#include <algebra/SmallMatrix.hpp>
-#include <algebra/TinyMatrix.hpp>
+#include <algebra/DenseMatrixWrapper.hpp>
#include <eigen3/Eigen/Eigen>
@@ -25,6 +24,13 @@ class Eigen3DenseMatrixEmbedder
{}
public:
+ PUGS_INLINE
+ size_t
+ isSquare() const
+ {
+ return (m_eigen_matrix.rows() == m_eigen_matrix.cols());
+ }
+
PUGS_INLINE
size_t
numberOfRows() const
@@ -53,12 +59,8 @@ class Eigen3DenseMatrixEmbedder
return m_eigen_matrix;
}
- template <size_t N>
- Eigen3DenseMatrixEmbedder(const TinyMatrix<N>& A) : Eigen3DenseMatrixEmbedder{N, N, &A(0, 0)}
- {}
-
- Eigen3DenseMatrixEmbedder(const SmallMatrix<double>& A)
- : Eigen3DenseMatrixEmbedder{A.numberOfRows(), A.numberOfColumns(), &A(0, 0)}
+ Eigen3DenseMatrixEmbedder(const DenseMatrixWrapper<double>& A)
+ : Eigen3DenseMatrixEmbedder{A.numberOfRows(), A.numberOfColumns(), A.ptr()}
{}
~Eigen3DenseMatrixEmbedder() = default;
@@ -88,6 +90,13 @@ class Eigen3SparseMatrixEmbedder
return m_eigen_matrix.cols();
}
+ PUGS_INLINE
+ size_t
+ isSquare() const
+ {
+ return (m_eigen_matrix.rows() == m_eigen_matrix.cols());
+ }
+
PUGS_INLINE
Eigen3MapMatrixType&
matrix()
@@ -113,6 +122,24 @@ class Eigen3SparseMatrixEmbedder
~Eigen3SparseMatrixEmbedder() = default;
};
+#else // PUGS_HAS_EIGEN3
+
+class Eigen3DenseMatrixEmbedder
+{
+ public:
+ Eigen3DenseMatrixEmbedder(const DenseMatrixWrapper<double>&) {}
+
+ ~Eigen3DenseMatrixEmbedder() = default;
+};
+
+class Eigen3SparseMatrixEmbedder
+{
+ public:
+ Eigen3SparseMatrixEmbedder(const CRSMatrix<double>&) {}
+
+ ~Eigen3SparseMatrixEmbedder() = default;
+};
+
#endif // PUGS_HAS_EIGEN3
#endif // EIGEN3_UTILS_HPP
diff --git a/src/algebra/LinearSolver.cpp b/src/algebra/LinearSolver.cpp
index 1ebb59663405883aac2405979ff60bdb7ccb3e7e..52bf8fc3969d2001685c5548f31c30ce7af20dd0 100644
--- a/src/algebra/LinearSolver.cpp
+++ b/src/algebra/LinearSolver.cpp
@@ -187,10 +187,10 @@ struct LinearSolver::Internals
template <typename MatrixType, typename SolutionVectorType, typename RHSVectorType>
static void
- builtinSolveLocalSystem(const MatrixType& A,
- SolutionVectorType& x,
- const RHSVectorType& b,
- const LinearSolverOptions& options)
+ _builtinSolveLocalSystem(const MatrixType& A,
+ SolutionVectorType& x,
+ const RHSVectorType& b,
+ const LinearSolverOptions& options)
{
if (options.precond() != LSPrecond::none) {
// LCOV_EXCL_START
@@ -215,18 +215,15 @@ struct LinearSolver::Internals
}
#ifdef PUGS_HAS_EIGEN3
- template <typename PreconditionerType,
- typename Eigen3MatrixEmbedderType,
- typename SolutionVectorType,
- typename RHSVectorType>
+ template <typename PreconditionerType, typename Eigen3MatrixEmbedderType>
static void
_preconditionedEigen3SolveLocalSystem(const Eigen3MatrixEmbedderType& eigen3_A,
- SolutionVectorType& x,
- const RHSVectorType& b,
+ VectorWrapper<double>& x,
+ const VectorWrapper<const double>& b,
const LinearSolverOptions& options)
{
- Eigen::Map<Eigen::VectorX<double>> eigen3_x{(x.size() > 0) ? (&x[0]) : nullptr, static_cast<int>(x.size())};
- Eigen::Map<const Eigen::VectorX<double>> eigen3_b{(b.size() > 0) ? (&b[0]) : nullptr, static_cast<int>(b.size())};
+ Eigen::Map<Eigen::VectorX<double>> eigen3_x{x.ptr(), static_cast<int>(x.size())};
+ Eigen::Map<const Eigen::VectorX<double>> eigen3_b{b.ptr(), static_cast<int>(b.size())};
using Eigen3MatrixType = typename Eigen3MatrixEmbedderType::Eigen3MatrixType;
@@ -287,28 +284,29 @@ struct LinearSolver::Internals
}
}
- template <typename Eigen3MatrixEmbedderType, typename SolutionVectorType, typename RHSVectorType>
+ template <typename Eigen3MatrixEmbedderType>
static void
_eigen3SolveLocalSystem(const Eigen3MatrixEmbedderType& eigen3_A,
- SolutionVectorType& x,
- const RHSVectorType& b,
+ VectorWrapper<double>& x,
+ const VectorWrapper<const double>& b,
const LinearSolverOptions& options)
{
switch (options.precond()) {
case LSPrecond::none: {
- _preconditionedEigen3SolveLocalSystem<Eigen::IdentityPreconditioner, Eigen3MatrixEmbedderType, SolutionVectorType,
- RHSVectorType>(eigen3_A, x, b, options);
+ _preconditionedEigen3SolveLocalSystem<Eigen::IdentityPreconditioner, Eigen3MatrixEmbedderType>(eigen3_A, x, b,
+ options);
break;
}
case LSPrecond::diagonal: {
- _preconditionedEigen3SolveLocalSystem<Eigen::DiagonalPreconditioner<double>, Eigen3MatrixEmbedderType,
- SolutionVectorType, RHSVectorType>(eigen3_A, x, b, options);
+ _preconditionedEigen3SolveLocalSystem<Eigen::DiagonalPreconditioner<double>, Eigen3MatrixEmbedderType>(eigen3_A,
+ x, b,
+ options);
break;
}
case LSPrecond::incomplete_cholesky: {
if constexpr (std::is_same_v<Eigen3SparseMatrixEmbedder, Eigen3MatrixEmbedderType>) {
- _preconditionedEigen3SolveLocalSystem<Eigen::IncompleteCholesky<double>, Eigen3MatrixEmbedderType,
- SolutionVectorType, RHSVectorType>(eigen3_A, x, b, options);
+ _preconditionedEigen3SolveLocalSystem<Eigen::IncompleteCholesky<double>, Eigen3MatrixEmbedderType>(eigen3_A, x,
+ b, options);
} else {
throw NormalError("incomplete cholesky is not available for dense matrices in Eigen3");
}
@@ -316,8 +314,8 @@ struct LinearSolver::Internals
}
case LSPrecond::incomplete_LU: {
if constexpr (std::is_same_v<Eigen3SparseMatrixEmbedder, Eigen3MatrixEmbedderType>) {
- _preconditionedEigen3SolveLocalSystem<Eigen::IncompleteLUT<double>, Eigen3MatrixEmbedderType,
- SolutionVectorType, RHSVectorType>(eigen3_A, x, b, options);
+ _preconditionedEigen3SolveLocalSystem<Eigen::IncompleteLUT<double>, Eigen3MatrixEmbedderType>(eigen3_A, x, b,
+ options);
} else {
throw NormalError("incomplete LU is not available for dense matrices in Eigen3");
}
@@ -330,39 +328,12 @@ struct LinearSolver::Internals
// LCOV_EXCL_STOP
}
}
-
- template <typename DataType, typename SolutionVectorType, typename RHSVectorType>
- static void
- eigen3SolveLocalSystem(const CRSMatrix<DataType>& A,
- SolutionVectorType& x,
- const RHSVectorType& b,
- const LinearSolverOptions& options)
- {
- Assert((x.size() == b.size()) and (x.size() - A.numberOfColumns() == 0) and A.isSquare());
-
- Eigen3SparseMatrixEmbedder eigen3_A{A};
- _eigen3SolveLocalSystem(eigen3_A, x, b, options);
- }
-
- template <typename DataType, typename SolutionVectorType, typename RHSVectorType>
- static void
- eigen3SolveLocalSystem(const SmallMatrix<DataType>& A,
- SolutionVectorType& x,
- const RHSVectorType& b,
- const LinearSolverOptions& options)
- {
- Assert((x.size() == b.size()) and (x.size() - A.numberOfColumns() == 0) and A.isSquare());
-
- Eigen3DenseMatrixEmbedder eigen3_A{A};
- _eigen3SolveLocalSystem(eigen3_A, x, b, options);
- }
-
#else // PUGS_HAS_EIGEN3
// LCOV_EXCL_START
template <typename MatrixType, typename SolutionVectorType, typename RHSVectorType>
static void
- eigen3SolveLocalSystem(const MatrixType&, SolutionVectorType&, const RHSVectorType&, const LinearSolverOptions&)
+ _eigen3SolveLocalSystem(const MatrixType&, SolutionVectorType&, const RHSVectorType&, const LinearSolverOptions&)
{
checkHasLibrary(LSLibrary::eigen3);
throw UnexpectedError("unexpected situation should not reach this point!");
@@ -379,19 +350,17 @@ struct LinearSolver::Internals
return 0;
}
- template <typename MatrixType, typename SolutionVectorType, typename RHSVectorType>
+ template <typename MatrixType>
static void
- petscSolveLocalSystem(const MatrixType& A,
- SolutionVectorType& x,
- const RHSVectorType& b,
- const LinearSolverOptions& options)
+ _petscSolveLocalSystem(const MatrixType& A,
+ VectorWrapper<double>& x,
+ const VectorWrapper<const double>& b,
+ const LinearSolverOptions& options)
{
- Assert((x.size() == b.size()) and (x.size() - A.numberOfColumns() == 0) and A.isSquare());
-
Vec petscB;
- VecCreateMPIWithArray(PETSC_COMM_SELF, 1, b.size(), b.size(), &b[0], &petscB);
+ VecCreateMPIWithArray(PETSC_COMM_SELF, 1, b.size(), b.size(), b.ptr(), &petscB);
Vec petscX;
- VecCreateMPIWithArray(PETSC_COMM_SELF, 1, x.size(), x.size(), &x[0], &petscX);
+ VecCreateMPIWithArray(PETSC_COMM_SELF, 1, x.size(), x.size(), x.ptr(), &petscX);
PETScAijMatrixEmbedder petscA(A);
@@ -491,7 +460,7 @@ struct LinearSolver::Internals
// LCOV_EXCL_START
template <typename MatrixType, typename SolutionVectorType, typename RHSVectorType>
static void
- petscSolveLocalSystem(const MatrixType&, SolutionVectorType&, const RHSVectorType&, const LinearSolverOptions&)
+ _petscSolveLocalSystem(const MatrixType&, SolutionVectorType&, const RHSVectorType&, const LinearSolverOptions&)
{
checkHasLibrary(LSLibrary::petsc);
throw UnexpectedError("unexpected situation should not reach this point!");
@@ -499,40 +468,6 @@ struct LinearSolver::Internals
// LCOV_EXCL_STOP
#endif // PUGS_HAS_PETSC
-
- template <typename MatrixType, typename SolutionVectorType, typename RHSVectorType>
- static void
- solveLocalSystem(const MatrixType& A,
- SolutionVectorType& x,
- const RHSVectorType& b,
- const LinearSolverOptions& options)
- {
- switch (options.library()) {
- case LSLibrary::builtin: {
- builtinSolveLocalSystem(A, x, b, options);
- break;
- }
- // LCOV_EXCL_START
- case LSLibrary::eigen3: {
- // not covered since if Eigen3 is not linked this point is
- // unreachable: LinearSolver throws an exception at construction
- // in this case.
- eigen3SolveLocalSystem(A, x, b, options);
- break;
- }
- case LSLibrary::petsc: {
- // not covered since if PETSc is not linked this point is
- // unreachable: LinearSolver throws an exception at construction
- // in this case.
- petscSolveLocalSystem(A, x, b, options);
- break;
- }
- default: {
- throw UnexpectedError(::name(options.library()) + " cannot solve local systems");
- }
- // LCOV_EXCL_STOP
- }
- }
};
bool
@@ -548,15 +483,52 @@ LinearSolver::checkOptions(const LinearSolverOptions& options) const
}
void
-LinearSolver::solveLocalSystem(const CRSMatrix<double, int>& A, Vector<double>& x, const Vector<const double>& b)
+LinearSolver::_builtinSolveLocalSystem(const CRSMatrix<double, int>& A,
+ Vector<double>& x,
+ const Vector<const double>& b)
+{
+ Internals::_builtinSolveLocalSystem(A, x, b, m_options);
+}
+
+void
+LinearSolver::_builtinSolveLocalSystem(const SmallMatrix<double>& A,
+ SmallVector<double>& x,
+ const SmallVector<const double>& b)
+{
+ Internals::_builtinSolveLocalSystem(A, x, b, m_options);
+}
+
+void
+LinearSolver::_eigen3SolveLocalSystem(const Eigen3SparseMatrixEmbedder& A,
+ VectorWrapper<double> x,
+ const VectorWrapper<const double>& b)
+{
+ Internals::_eigen3SolveLocalSystem(A, x, b, m_options);
+}
+
+void
+LinearSolver::_eigen3SolveLocalSystem(const DenseMatrixWrapper<double>& A,
+ VectorWrapper<double> x,
+ const VectorWrapper<const double>& b)
+{
+ Eigen3DenseMatrixEmbedder A_wrapper{A};
+ Internals::_eigen3SolveLocalSystem(A_wrapper, x, b, m_options);
+}
+
+void
+LinearSolver::_petscSolveLocalSystem(const CRSMatrix<double, int>& A,
+ VectorWrapper<double> x,
+ const VectorWrapper<const double>& b)
{
- Internals::solveLocalSystem(A, x, b, m_options);
+ Internals::_petscSolveLocalSystem(A, x, b, m_options);
}
void
-LinearSolver::solveLocalSystem(const SmallMatrix<double>& A, SmallVector<double>& x, const SmallVector<const double>& b)
+LinearSolver::_petscSolveLocalSystem(const DenseMatrixWrapper<double>& A,
+ VectorWrapper<double> x,
+ const VectorWrapper<const double>& b)
{
- Internals::solveLocalSystem(A, x, b, m_options);
+ Internals::_petscSolveLocalSystem(A, x, b, m_options);
}
LinearSolver::LinearSolver(const LinearSolverOptions& options) : m_options{options}
diff --git a/src/algebra/LinearSolver.hpp b/src/algebra/LinearSolver.hpp
index e90402c8ddc78c9426f97300ab6b10730b4c503b..2d03a9cf5e7cced7f99b6707c042395f48054199 100644
--- a/src/algebra/LinearSolver.hpp
+++ b/src/algebra/LinearSolver.hpp
@@ -2,11 +2,14 @@
#define LINEAR_SOLVER_HPP
#include <algebra/CRSMatrix.hpp>
+#include <algebra/DenseMatrixWrapper.hpp>
+#include <algebra/Eigen3Utils.hpp>
#include <algebra/LinearSolverOptions.hpp>
#include <algebra/SmallMatrix.hpp>
#include <algebra/TinyMatrix.hpp>
#include <algebra/TinyVector.hpp>
#include <algebra/Vector.hpp>
+#include <algebra/VectorWrapper.hpp>
#include <utils/Exceptions.hpp>
class LinearSolver
@@ -16,29 +19,139 @@ class LinearSolver
const LinearSolverOptions m_options;
+ void _builtinSolveLocalSystem(const CRSMatrix<double, int>& A, Vector<double>& x, const Vector<const double>& b);
+
+ void _builtinSolveLocalSystem(const SmallMatrix<double>& A,
+ SmallVector<double>& x,
+ const SmallVector<const double>& b);
+
+ void _eigen3SolveLocalSystem(const Eigen3SparseMatrixEmbedder& A,
+ VectorWrapper<double> x,
+ const VectorWrapper<const double>& b);
+
+ void _eigen3SolveLocalSystem(const DenseMatrixWrapper<double>& A,
+ VectorWrapper<double> x,
+ const VectorWrapper<const double>& b);
+
+ void _petscSolveLocalSystem(const CRSMatrix<double, int>& A,
+ VectorWrapper<double> x,
+ const VectorWrapper<const double>& b);
+
+ void _petscSolveLocalSystem(const DenseMatrixWrapper<double>& A,
+ VectorWrapper<double> x,
+ const VectorWrapper<const double>& b);
+
public:
bool hasLibrary(LSLibrary library) const;
void checkOptions(const LinearSolverOptions& options) const;
+ void
+ solveLocalSystem(const CRSMatrix<double, int>& A, Vector<double>& x, const Vector<const double>& b)
+ {
+ Assert((x.size() == b.size()) and (x.size() - A.numberOfColumns() == 0) and A.isSquare());
+ switch (m_options.library()) {
+ case LSLibrary::builtin: {
+ _builtinSolveLocalSystem(A, x, b);
+ break;
+ }
+ // LCOV_EXCL_START
+ case LSLibrary::eigen3: {
+ // not covered since if Eigen3 is not linked this point is
+ // unreachable: LinearSolver throws an exception at construction
+ // in this case.
+ _eigen3SolveLocalSystem(A, x, b);
+ break;
+ }
+ case LSLibrary::petsc: {
+ // not covered since if PETSc is not linked this point is
+ // unreachable: LinearSolver throws an exception at construction
+ // in this case.
+ _petscSolveLocalSystem(A, x, b);
+ break;
+ }
+ default: {
+ throw UnexpectedError(::name(m_options.library()) + " cannot solve local systems");
+ }
+ // LCOV_EXCL_STOP
+ }
+ }
+
+ void
+ solveLocalSystem(const SmallMatrix<double>& A, SmallVector<double>& x, const SmallVector<const double>& b)
+ {
+ Assert((x.size() == b.size()) and (x.size() - A.numberOfColumns() == 0) and A.isSquare());
+ switch (m_options.library()) {
+ case LSLibrary::builtin: {
+ _builtinSolveLocalSystem(A, x, b);
+ break;
+ }
+ // LCOV_EXCL_START
+ case LSLibrary::eigen3: {
+ // not covered since if Eigen3 is not linked this point is
+ // unreachable: LinearSolver throws an exception at construction
+ // in this case.
+ _eigen3SolveLocalSystem(A, x, b);
+ break;
+ }
+ case LSLibrary::petsc: {
+ // not covered since if PETSc is not linked this point is
+ // unreachable: LinearSolver throws an exception at construction
+ // in this case.
+ _petscSolveLocalSystem(A, x, b);
+ break;
+ }
+ default: {
+ throw UnexpectedError(::name(m_options.library()) + " cannot solve local systems");
+ }
+ // LCOV_EXCL_STOP
+ }
+ }
+
template <size_t N>
void
solveLocalSystem(const TinyMatrix<N>& A, TinyVector<N>& x, const TinyVector<N>& b)
{
- SmallMatrix A_dense{A};
+ Assert((x.size() == b.size()) and (x.size() - A.numberOfColumns() == 0) and A.isSquare());
+ switch (m_options.library()) {
+ case LSLibrary::builtin: {
+ // Probably expensive but builtin is not the preferred way to
+ // solve linear systems
+ SmallMatrix A_dense{A};
+
+ SmallVector x_vector{x};
+ SmallVector<const double> b_vector{b};
- SmallVector x_vector{x};
- SmallVector b_vector{b};
+ this->solveLocalSystem(A_dense, x_vector, b_vector);
- this->solveLocalSystem(A_dense, x_vector, b_vector);
+ for (size_t i = 0; i < N; ++i) {
+ x[i] = x_vector[i];
+ }
- for (size_t i = 0; i < N; ++i) {
- x[i] = x_vector[i];
+ _builtinSolveLocalSystem(A, x, b);
+ break;
+ }
+ // LCOV_EXCL_START
+ case LSLibrary::eigen3: {
+ // not covered since if Eigen3 is not linked this point is
+ // unreachable: LinearSolver throws an exception at construction
+ // in this case.
+ _eigen3SolveLocalSystem(A, x, b);
+ break;
+ }
+ case LSLibrary::petsc: {
+ // not covered since if PETSc is not linked this point is
+ // unreachable: LinearSolver throws an exception at construction
+ // in this case.
+ _petscSolveLocalSystem(A, x, b);
+ break;
+ }
+ default: {
+ throw UnexpectedError(::name(m_options.library()) + " cannot solve local systems");
+ }
+ // LCOV_EXCL_STOP
}
}
- void solveLocalSystem(const CRSMatrix<double, int>& A, Vector<double>& x, const Vector<const double>& b);
- void solveLocalSystem(const SmallMatrix<double>& A, SmallVector<double>& x, const SmallVector<const double>& b);
-
LinearSolver();
LinearSolver(const LinearSolverOptions& options);
};
diff --git a/src/algebra/PETScUtils.hpp b/src/algebra/PETScUtils.hpp
index c4db9ee5910ef8e58333dbfb6dc20c8a1987bb46..c6d2e8be88c6098e2f187cda9177c120facc6f23 100644
--- a/src/algebra/PETScUtils.hpp
+++ b/src/algebra/PETScUtils.hpp
@@ -6,6 +6,7 @@
#ifdef PUGS_HAS_PETSC
#include <algebra/CRSMatrix.hpp>
+#include <algebra/DenseMatrixWrapper.hpp>
#include <algebra/SmallMatrix.hpp>
#include <algebra/TinyMatrix.hpp>
@@ -52,12 +53,8 @@ class PETScAijMatrixEmbedder
return m_petscMat;
}
- template <size_t N>
- PETScAijMatrixEmbedder(const TinyMatrix<N>& A) : PETScAijMatrixEmbedder{N, N, &A(0, 0)}
- {}
-
- PETScAijMatrixEmbedder(const SmallMatrix<double>& A)
- : PETScAijMatrixEmbedder{A.numberOfRows(), A.numberOfColumns(), &A(0, 0)}
+ PETScAijMatrixEmbedder(const DenseMatrixWrapper<double>& A)
+ : PETScAijMatrixEmbedder{A.numberOfRows(), A.numberOfColumns(), A.ptr()}
{}
PETScAijMatrixEmbedder(const CRSMatrix<double, int>& A);
diff --git a/src/algebra/VectorWrapper.hpp b/src/algebra/VectorWrapper.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..0bdeefbfd3f10407358613e7f1d2487bba22b26d
--- /dev/null
+++ b/src/algebra/VectorWrapper.hpp
@@ -0,0 +1,47 @@
+#ifndef VECTOR_WRAPPER_HPP
+#define VECTOR_WRAPPER_HPP
+
+#include <algebra/SmallVector.hpp>
+#include <algebra/TinyVector.hpp>
+#include <algebra/Vector.hpp>
+
+template <typename DataType>
+class VectorWrapper
+{
+ private:
+ const size_t m_size;
+ DataType* const m_vector_ptr;
+
+ public:
+ PUGS_INLINE
+ size_t
+ size() const
+ {
+ return m_size;
+ }
+
+ PUGS_INLINE
+ DataType* const&
+ ptr() const
+ {
+ return m_vector_ptr;
+ }
+
+ VectorWrapper(const Vector<DataType>& x) : m_size{x.size()}, m_vector_ptr{(m_size > 0) ? (&x[0]) : nullptr} {}
+ VectorWrapper(Vector<DataType>& x) : m_size{x.size()}, m_vector_ptr{(m_size > 0) ? (&x[0]) : nullptr} {}
+
+ VectorWrapper(const SmallVector<DataType>& x) : m_size{x.size()}, m_vector_ptr{(m_size > 0) ? (&x[0]) : nullptr} {}
+
+ template <size_t N>
+ VectorWrapper(const TinyVector<N, DataType>& x) : m_size{N}, m_vector_ptr{(N > 0) ? (&x[0]) : nullptr}
+ {}
+
+ template <size_t N>
+ VectorWrapper(TinyVector<N, DataType>& x) : m_size{N}, m_vector_ptr{(N > 0) ? (&x[0]) : nullptr}
+ {}
+
+ VectorWrapper(const VectorWrapper&) = delete;
+ VectorWrapper(VectorWrapper&&) = delete;
+};
+
+#endif // VECTOR_WRAPPER_HPP
diff --git a/src/utils/PugsTraits.hpp b/src/utils/PugsTraits.hpp
index 2e1d409578c9295a1ad032466842c2f336ba272f..3c7a4e40a217d957dea388223d99574e6b1ecf74 100644
--- a/src/utils/PugsTraits.hpp
+++ b/src/utils/PugsTraits.hpp
@@ -15,6 +15,12 @@ class TinyVector;
template <size_t M, size_t N, typename T>
class TinyMatrix;
+template <typename DataType>
+class SmallMatrix;
+
+template <typename DataType, typename IndexType>
+class CRSMatrix;
+
template <typename DataType, ItemType item_type, typename ConnectivityPtr>
class ItemValue;
template <typename DataType, ItemType item_type, typename ConnectivityPtr>
@@ -114,6 +120,22 @@ inline constexpr bool is_tiny_matrix_v = false;
template <size_t M, size_t N, typename T>
inline constexpr bool is_tiny_matrix_v<TinyMatrix<M, N, T>> = true;
+// Traits is_small_matrix
+
+template <typename T>
+inline constexpr bool is_small_matrix_v = false;
+
+template <typename DataType>
+inline constexpr bool is_small_matrix_v<SmallMatrix<DataType>> = true;
+
+// Traits is_crs_matrix
+
+template <typename T>
+inline constexpr bool is_crs_matrix_v = false;
+
+template <typename DataType, typename IndexType>
+inline constexpr bool is_crs_matrix_v<CRSMatrix<DataType, IndexType>> = true;
+
// Trais is ItemValue
template <typename T>
diff --git a/tests/test_LinearSolver.cpp b/tests/test_LinearSolver.cpp
index 4d01df2f0605e2bd141c1de3d95f28776f3178a0..76ef6aa2bb75c3ea94fcf392f59ddb74350dc50a 100644
--- a/tests/test_LinearSolver.cpp
+++ b/tests/test_LinearSolver.cpp
@@ -752,46 +752,10 @@ TEST_CASE("LinearSolver", "[algebra]")
}
}
- SECTION("BICGStab2")
- {
- LinearSolverOptions options;
- options.library() = LSLibrary::petsc;
- options.method() = LSMethod::bicgstab2;
- options.precond() = LSPrecond::none;
-
- SECTION("BICGStab2 no preconditioner")
- {
- options.precond() = LSPrecond::none;
-
- Vector<double> x{5};
- x = zero;
-
- LinearSolver solver{options};
-
- solver.solveLocalSystem(A, x, b);
- Vector error = x - x_exact;
- REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
- }
-
- SECTION("BICGStab2 Diagonal")
- {
- options.precond() = LSPrecond::diagonal;
-
- Vector<double> x{5};
- x = zero;
-
- LinearSolver solver{options};
-
- solver.solveLocalSystem(A, x, b);
- Vector error = x - x_exact;
- REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
- }
- }
-
SECTION("GMRES")
{
LinearSolverOptions options;
- options.library() = LSLibrary::petsc;
+ options.library() = LSLibrary::eigen3;
options.method() = LSMethod::gmres;
options.precond() = LSPrecond::none;
@@ -841,7 +805,7 @@ TEST_CASE("LinearSolver", "[algebra]")
SECTION("LU")
{
LinearSolverOptions options;
- options.library() = LSLibrary::petsc;
+ options.library() = LSLibrary::eigen3;
options.method() = LSMethod::lu;
options.precond() = LSPrecond::none;
@@ -859,7 +823,7 @@ TEST_CASE("LinearSolver", "[algebra]")
SECTION("Eigen3 not linked")
{
LinearSolverOptions options;
- options.library() = LSLibrary::petsc;
+ options.library() = LSLibrary::eigen3;
options.method() = LSMethod::cg;
options.precond() = LSPrecond::none;