diff --git a/CMakeLists.txt b/CMakeLists.txt
index a4aa863b266890293cb48ceaf21990120f4d303b..e9255bff6ac523e8bfdf61fe42044e4366d3993f 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -187,6 +187,34 @@ else()
endif()
endif()
+#------------------------------------------------------
+# Check for SLEPc
+# defaults use SLEPc
+set(PUGS_ENABLE_SLEPC AUTO CACHE STRING
+ "Choose one of: AUTO ON OFF")
+
+if (PUGS_ENABLE_SLEPC MATCHES "^(AUTO|ON)$")
+ if (PETSC_FOUND)
+ # SLEPc support is deactivated if PETSc is not found
+ pkg_check_modules(SLEPC SLEPc)
+ else()
+ message(STATUS "SLEPc support is deactivated since pugs will not be build with PETSc support")
+ set(SLEPc_FOUND FALSE)
+ unset(PUGS_HAS_SLEPC)
+ endif()
+ set(PUGS_HAS_SLEPC ${SLEPC_FOUND})
+else()
+ unset(PUGS_HAS_SLEPC)
+endif()
+
+if (${SLEPC_FOUND})
+ include_directories(SYSTEM ${SLEPC_INCLUDE_DIRS})
+else()
+ if (PUGS_ENABLE_SLEPC MATCHES "^ON$")
+ message(FATAL_ERROR "Could not find SLEPc!")
+ endif()
+endif()
+
# -----------------------------------------------------
if (${MPI_FOUND})
@@ -558,6 +586,7 @@ target_link_libraries(
PugsLanguageUtils
kokkos
${PETSC_LIBRARIES}
+ ${SLEPC_LIBRARIES}
${PARMETIS_LIBRARIES}
${MPI_CXX_LINK_FLAGS} ${MPI_CXX_LIBRARIES}
${KOKKOS_CXX_FLAGS}
@@ -615,6 +644,16 @@ else()
endif()
endif()
+if (SLEPC_FOUND)
+ message(" SLEPc: ${SLEPC_VERSION}")
+else()
+ if (PUGS_ENABLE_SLEPC MATCHES "^(AUTO|ON)$")
+ message(" SLEPc: not found!")
+ else()
+ message(" SLEPc: explicitly deactivated!")
+ endif()
+endif()
+
if(CLANG_FORMAT)
message(" clang-format: ${CLANG_FORMAT}")
else()
diff --git a/src/algebra/BiCGStab.hpp b/src/algebra/BiCGStab.hpp
index dc62e370e6ec304cb2e4fd44d220218d9bd94283..82d3446791cd0a65065169d4f8d15f87120931a1 100644
--- a/src/algebra/BiCGStab.hpp
+++ b/src/algebra/BiCGStab.hpp
@@ -9,10 +9,10 @@
struct BiCGStab
{
- template <typename VectorType, typename MatrixType>
+ template <typename MatrixType, typename VectorType, typename RHSVectorType>
BiCGStab(const MatrixType& A,
VectorType& x,
- const VectorType& b,
+ const RHSVectorType& b,
const double epsilon,
const size_t maximum_iteration,
const bool verbose)
diff --git a/src/algebra/CG.hpp b/src/algebra/CG.hpp
index 7e653d8f95cddab4d5d1b7eb36481fd7690ef170..4ab52012a08732f21d57cf114bd873406bfd3898 100644
--- a/src/algebra/CG.hpp
+++ b/src/algebra/CG.hpp
@@ -8,10 +8,10 @@
struct CG
{
- template <typename VectorType, typename MatrixType>
+ template <typename MatrixType, typename VectorType, typename RHSVectorType>
CG(const MatrixType& A,
VectorType& x,
- const VectorType& f,
+ const RHSVectorType& f,
const double epsilon,
const size_t maximum_iteration,
const bool verbose = false)
diff --git a/src/algebra/CMakeLists.txt b/src/algebra/CMakeLists.txt
index d1f01ec310858be0e775514f934ce0a7c0e9da28..6310f48dcf69eb93e2c1178e05d32dbc7fe5afad 100644
--- a/src/algebra/CMakeLists.txt
+++ b/src/algebra/CMakeLists.txt
@@ -2,6 +2,13 @@
add_library(
PugsAlgebra
+ EigenvalueSolver.cpp
LinearSolver.cpp
LinearSolverOptions.cpp
- PETScWrapper.cpp)
+ PETScUtils.cpp)
+
+target_link_libraries(
+ PugsAlgebra
+ ${PETSC_LIBRARIES}
+ ${SLEPC_LIBRARIES}
+)
diff --git a/src/algebra/CRSMatrix.hpp b/src/algebra/CRSMatrix.hpp
index 98297d759675fdb7fdbf8be1e638f19cba4cc562..edb85d5e87b6a347ad5bc4f82987fd9255154484 100644
--- a/src/algebra/CRSMatrix.hpp
+++ b/src/algebra/CRSMatrix.hpp
@@ -21,13 +21,29 @@ class CRSMatrix
HostMatrix m_host_matrix;
Array<const DataType> m_values;
+ size_t m_nb_rows;
+ size_t m_nb_columns;
public:
+ PUGS_INLINE
+ bool
+ isSquare() const noexcept
+ {
+ return m_nb_rows == m_nb_columns;
+ }
+
PUGS_INLINE
size_t
numberOfRows() const
{
- return m_host_matrix.numRows();
+ return m_nb_rows;
+ }
+
+ PUGS_INLINE
+ size_t
+ numberOfColumns() const
+ {
+ return m_nb_columns;
}
auto
@@ -75,6 +91,7 @@ class CRSMatrix
}
CRSMatrix(const SparseMatrixDescriptor<DataType, IndexType>& M)
+ : m_nb_rows{M.numberOfRows()}, m_nb_columns{M.numberOfColumns()}
{
{
auto host_matrix =
diff --git a/src/algebra/DenseMatrix.hpp b/src/algebra/DenseMatrix.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..8f15f125861d0cc9cce0aef716e63959f736d36b
--- /dev/null
+++ b/src/algebra/DenseMatrix.hpp
@@ -0,0 +1,296 @@
+#ifndef DENSE_MATRIX_HPP
+#define DENSE_MATRIX_HPP
+
+#include <algebra/TinyMatrix.hpp>
+#include <algebra/Vector.hpp>
+#include <utils/Array.hpp>
+#include <utils/PugsAssert.hpp>
+#include <utils/PugsMacros.hpp>
+#include <utils/PugsUtils.hpp>
+#include <utils/Types.hpp>
+
+template <typename DataType>
+class DenseMatrix // LCOV_EXCL_LINE
+{
+ public:
+ using data_type = DataType;
+ using index_type = size_t;
+
+ private:
+ size_t m_nb_rows;
+ size_t m_nb_columns;
+ Array<DataType> m_values;
+
+ static_assert(std::is_same_v<typename decltype(m_values)::index_type, index_type>);
+ static_assert(std::is_arithmetic_v<DataType>, "Dense matrices expect arithmetic data");
+
+ // Allows const version to access our data
+ friend DenseMatrix<std::add_const_t<DataType>>;
+
+ public:
+ PUGS_INLINE
+ bool
+ isSquare() const noexcept
+ {
+ return m_nb_rows == m_nb_columns;
+ }
+
+ friend PUGS_INLINE DenseMatrix<std::remove_const_t<DataType>>
+ copy(const DenseMatrix& A) noexcept
+ {
+ return {A.m_nb_rows, A.m_nb_columns, copy(A.m_values)};
+ }
+
+ friend PUGS_INLINE DenseMatrix<std::remove_const_t<DataType>>
+ transpose(const DenseMatrix& A)
+ {
+ DenseMatrix<std::remove_const_t<DataType>> A_transpose{A.m_nb_columns, A.m_nb_rows};
+ for (size_t i = 0; i < A.m_nb_rows; ++i) {
+ for (size_t j = 0; j < A.m_nb_columns; ++j) {
+ A_transpose(j, i) = A(i, j);
+ }
+ }
+ return A_transpose;
+ }
+
+ friend PUGS_INLINE DenseMatrix
+ operator*(const DataType& a, const DenseMatrix& A)
+ {
+ DenseMatrix<std::remove_const_t<DataType>> aA = copy(A);
+ return aA *= a;
+ }
+
+ template <typename DataType2>
+ PUGS_INLINE Vector<std::remove_const_t<DataType>>
+ operator*(const Vector<DataType2>& x) const
+ {
+ static_assert(std::is_same_v<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>,
+ "incompatible data types");
+ Assert(m_nb_columns == x.size());
+ const DenseMatrix& A = *this;
+ Vector<std::remove_const_t<DataType>> Ax{m_nb_rows};
+ for (size_t i = 0; i < m_nb_rows; ++i) {
+ std::remove_const_t<DataType> Axi = A(i, 0) * x[0];
+ for (size_t j = 1; j < m_nb_columns; ++j) {
+ Axi += A(i, j) * x[j];
+ }
+ Ax[i] = Axi;
+ }
+ return Ax;
+ }
+
+ template <typename DataType2>
+ PUGS_INLINE DenseMatrix<std::remove_const_t<DataType>>
+ operator*(const DenseMatrix<DataType2>& B) const
+ {
+ static_assert(std::is_same_v<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>,
+ "incompatible data types");
+ Assert(m_nb_columns == B.numberOfRows());
+ const DenseMatrix& A = *this;
+ DenseMatrix<std::remove_const_t<DataType>> AB{m_nb_rows, B.numberOfColumns()};
+
+ for (size_t i = 0; i < m_nb_rows; ++i) {
+ for (size_t j = 0; j < B.numberOfColumns(); ++j) {
+ std::remove_const_t<DataType> ABij = 0;
+ for (size_t k = 0; k < m_nb_columns; ++k) {
+ ABij += A(i, k) * B(k, j);
+ }
+ AB(i, j) = ABij;
+ }
+ }
+ return AB;
+ }
+
+ template <typename DataType2>
+ PUGS_INLINE DenseMatrix&
+ operator/=(const DataType2& a)
+ {
+ const auto inv_a = 1. / a;
+ return (*this) *= inv_a;
+ }
+
+ template <typename DataType2>
+ PUGS_INLINE DenseMatrix&
+ operator*=(const DataType2& a)
+ {
+ parallel_for(
+ m_values.size(), PUGS_LAMBDA(index_type i) { m_values[i] *= a; });
+ return *this;
+ }
+
+ template <typename DataType2>
+ PUGS_INLINE DenseMatrix&
+ operator-=(const DenseMatrix<DataType2>& B)
+ {
+ static_assert(std::is_same_v<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>,
+ "incompatible data types");
+ Assert(m_nb_rows == B.numberOfRows());
+ Assert(m_nb_columns == B.numberOfColumns());
+
+ parallel_for(
+ m_values.size(), PUGS_LAMBDA(index_type i) { m_values[i] -= B.m_values[i]; });
+ return *this;
+ }
+
+ template <typename DataType2>
+ PUGS_INLINE DenseMatrix&
+ operator+=(const DenseMatrix<DataType2>& B)
+ {
+ static_assert(std::is_same_v<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>,
+ "incompatible data types");
+ Assert(m_nb_rows == B.numberOfRows());
+ Assert(m_nb_columns == B.numberOfColumns());
+
+ parallel_for(
+ m_values.size(), PUGS_LAMBDA(index_type i) { m_values[i] += B.m_values[i]; });
+ return *this;
+ }
+
+ template <typename DataType2>
+ PUGS_INLINE DenseMatrix<std::remove_const_t<DataType>>
+ operator+(const DenseMatrix<DataType2>& B) const
+ {
+ static_assert(std::is_same_v<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>,
+ "incompatible data types");
+ Assert(m_nb_rows == B.numberOfRows());
+ Assert(m_nb_columns == B.numberOfColumns());
+ DenseMatrix<std::remove_const_t<DataType>> sum{B.numberOfRows(), B.numberOfColumns()};
+
+ parallel_for(
+ m_values.size(), PUGS_LAMBDA(index_type i) { sum.m_values[i] = m_values[i] + B.m_values[i]; });
+
+ return sum;
+ }
+
+ template <typename DataType2>
+ PUGS_INLINE DenseMatrix<std::remove_const_t<DataType>>
+ operator-(const DenseMatrix<DataType2>& B) const
+ {
+ static_assert(std::is_same_v<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>,
+ "incompatible data types");
+ Assert(m_nb_rows == B.numberOfRows());
+ Assert(m_nb_columns == B.numberOfColumns());
+ DenseMatrix<std::remove_const_t<DataType>> difference{B.numberOfRows(), B.numberOfColumns()};
+
+ parallel_for(
+ m_values.size(), PUGS_LAMBDA(index_type i) { difference.m_values[i] = m_values[i] - B.m_values[i]; });
+
+ return difference;
+ }
+
+ PUGS_INLINE
+ DataType&
+ operator()(index_type i, index_type j) const noexcept(NO_ASSERT)
+ {
+ Assert((i < m_nb_rows and j < m_nb_columns), "invalid indices");
+ return m_values[i * m_nb_columns + j];
+ }
+
+ PUGS_INLINE
+ size_t
+ numberOfRows() const noexcept
+ {
+ return m_nb_rows;
+ }
+
+ PUGS_INLINE
+ size_t
+ numberOfColumns() const noexcept
+ {
+ return m_nb_columns;
+ }
+
+ PUGS_INLINE void
+ fill(const DataType& value) noexcept
+ {
+ m_values.fill(value);
+ }
+
+ PUGS_INLINE DenseMatrix& operator=(ZeroType) noexcept
+ {
+ m_values.fill(0);
+ return *this;
+ }
+
+ PUGS_INLINE DenseMatrix& operator=(IdentityType) noexcept(NO_ASSERT)
+ {
+ Assert(m_nb_rows == m_nb_columns, "Identity must be a square matrix");
+
+ m_values.fill(0);
+ parallel_for(
+ m_nb_rows, PUGS_LAMBDA(const index_type i) { m_values[i * m_nb_rows + i] = 1; });
+ return *this;
+ }
+
+ template <typename DataType2>
+ PUGS_INLINE DenseMatrix&
+ operator=(const DenseMatrix<DataType2>& A) noexcept
+ {
+ // ensures that DataType is the same as source DataType2
+ static_assert(std::is_same<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>(),
+ "Cannot assign DenseMatrix of different type");
+ // ensures that const is not lost through copy
+ static_assert(((std::is_const<DataType2>() and std::is_const<DataType>()) or not std::is_const<DataType2>()),
+ "Cannot assign DenseMatrix of const to DenseMatrix of non-const");
+
+ m_nb_rows = A.m_nb_rows;
+ m_nb_columns = A.m_nb_columns;
+ m_values = A.m_values;
+ return *this;
+ }
+
+ PUGS_INLINE
+ DenseMatrix& operator=(const DenseMatrix&) = default;
+
+ PUGS_INLINE
+ DenseMatrix& operator=(DenseMatrix&&) = default;
+
+ template <typename DataType2>
+ DenseMatrix(const DenseMatrix<DataType2>& A)
+ {
+ // ensures that DataType is the same as source DataType2
+ static_assert(std::is_same<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>(),
+ "Cannot assign DenseMatrix of different type");
+ // ensures that const is not lost through copy
+ static_assert(((std::is_const<DataType2>() and std::is_const<DataType>()) or not std::is_const<DataType2>()),
+ "Cannot assign DenseMatrix of const to DenseMatrix of non-const");
+
+ this->operator=(A);
+ }
+
+ DenseMatrix(const DenseMatrix&) = default;
+
+ DenseMatrix(DenseMatrix&&) = default;
+
+ explicit DenseMatrix(size_t nb_rows, size_t nb_columns) noexcept
+ : m_nb_rows{nb_rows}, m_nb_columns{nb_columns}, m_values{nb_rows * nb_columns}
+ {}
+
+ explicit DenseMatrix(size_t nb_rows) noexcept : m_nb_rows{nb_rows}, m_nb_columns{nb_rows}, m_values{nb_rows * nb_rows}
+ {}
+
+ template <size_t N>
+ explicit DenseMatrix(const TinyMatrix<N, DataType>& M) noexcept : m_nb_rows{N}, m_nb_columns{N}, m_values{N * N}
+ {
+ parallel_for(
+ N, PUGS_LAMBDA(const index_type i) {
+ for (size_t j = 0; j < N; ++j) {
+ m_values[i * N + j] = M(i, j);
+ }
+ });
+ }
+
+ DenseMatrix() noexcept : m_nb_rows{0}, m_nb_columns{0} {}
+
+ private:
+ DenseMatrix(size_t nb_rows, size_t nb_columns, const Array<DataType> values) noexcept(NO_ASSERT)
+ : m_nb_rows{nb_rows}, m_nb_columns{nb_columns}, m_values{values}
+ {
+ Assert(m_values.size() == m_nb_rows * m_nb_columns, "incompatible sizes");
+ }
+
+ public:
+ ~DenseMatrix() = default;
+};
+
+#endif // DENSE_MATRIX_HPP
diff --git a/src/algebra/EigenvalueSolver.cpp b/src/algebra/EigenvalueSolver.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..1d96d01f992d22a8346417dce85750d9b186749b
--- /dev/null
+++ b/src/algebra/EigenvalueSolver.cpp
@@ -0,0 +1,152 @@
+#include <algebra/EigenvalueSolver.hpp>
+#include <utils/pugs_config.hpp>
+
+#ifdef PUGS_HAS_SLEPC
+#include <slepc.h>
+
+struct EigenvalueSolver::Internals
+{
+ static PetscReal
+ computeSmallestRealEigenvalueOfSymmetricMatrix(EPS& eps)
+ {
+ EPSSetWhichEigenpairs(eps, EPS_SMALLEST_REAL);
+ EPSSolve(eps);
+
+ PetscReal smallest_eigenvalue;
+ EPSGetEigenpair(eps, 0, &smallest_eigenvalue, nullptr, nullptr, nullptr);
+ return smallest_eigenvalue;
+ }
+
+ static PetscReal
+ computeLargestRealEigenvalueOfSymmetricMatrix(EPS& eps)
+ {
+ EPSSetWhichEigenpairs(eps, EPS_LARGEST_REAL);
+ EPSSolve(eps);
+
+ PetscReal largest_eigenvalue;
+ EPSGetEigenpair(eps, 0, &largest_eigenvalue, nullptr, nullptr, nullptr);
+ return largest_eigenvalue;
+ }
+
+ static void
+ computeAllEigenvaluesOfSymmetricMatrixInInterval(EPS& eps, const PetscReal left_bound, const PetscReal right_bound)
+ {
+ Assert(left_bound < right_bound);
+ EPSSetWhichEigenpairs(eps, EPS_ALL);
+ EPSSetInterval(eps, left_bound - 0.01 * std::abs(left_bound), right_bound + 0.01 * std::abs(right_bound));
+
+ ST st;
+ EPSGetST(eps, &st);
+ STSetType(st, STSINVERT);
+
+ KSP ksp;
+ STGetKSP(st, &ksp);
+ KSPSetType(ksp, KSPPREONLY);
+
+ PC pc;
+ KSPGetPC(ksp, &pc);
+ PCSetType(pc, PCCHOLESKY);
+ EPSSetFromOptions(eps);
+
+ EPSSolve(eps);
+ }
+
+ static void
+ computeAllEigenvaluesOfSymmetricMatrix(EPS& eps)
+ {
+ const PetscReal left_bound = computeSmallestRealEigenvalueOfSymmetricMatrix(eps);
+ const PetscReal right_bound = computeLargestRealEigenvalueOfSymmetricMatrix(eps);
+
+ computeAllEigenvaluesOfSymmetricMatrixInInterval(eps, left_bound - 0.01 * std::abs(left_bound),
+ right_bound + 0.01 * std::abs(right_bound));
+ }
+};
+
+void
+EigenvalueSolver::computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A, Array<double>& eigenvalues)
+{
+ EPS eps;
+
+ EPSCreate(PETSC_COMM_SELF, &eps);
+ EPSSetOperators(eps, A, nullptr);
+ EPSSetProblemType(eps, EPS_HEP);
+
+ Internals::computeAllEigenvaluesOfSymmetricMatrix(eps);
+
+ PetscInt nb_eigenvalues;
+ EPSGetDimensions(eps, &nb_eigenvalues, nullptr, nullptr);
+
+ eigenvalues = Array<double>(nb_eigenvalues);
+ for (PetscInt i = 0; i < nb_eigenvalues; ++i) {
+ EPSGetEigenpair(eps, i, &(eigenvalues[i]), nullptr, nullptr, nullptr);
+ }
+
+ EPSDestroy(&eps);
+}
+
+void
+EigenvalueSolver::computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A,
+ Array<double>& eigenvalues,
+ std::vector<Vector<double>>& eigenvector_list)
+{
+ EPS eps;
+
+ EPSCreate(PETSC_COMM_SELF, &eps);
+ EPSSetOperators(eps, A, nullptr);
+ EPSSetProblemType(eps, EPS_HEP);
+
+ Internals::computeAllEigenvaluesOfSymmetricMatrix(eps);
+
+ PetscInt nb_eigenvalues;
+ EPSGetDimensions(eps, &nb_eigenvalues, nullptr, nullptr);
+
+ eigenvalues = Array<double>(nb_eigenvalues);
+ eigenvector_list.reserve(nb_eigenvalues);
+ for (PetscInt i = 0; i < nb_eigenvalues; ++i) {
+ Vec Vr;
+ Vector<double> eigenvector{A.numberOfRows()};
+ VecCreateSeqWithArray(PETSC_COMM_SELF, 1, A.numberOfRows(), &(eigenvector[0]), &Vr);
+ EPSGetEigenpair(eps, i, &(eigenvalues[i]), nullptr, Vr, nullptr);
+ VecDestroy(&Vr);
+ eigenvector_list.push_back(eigenvector);
+ }
+
+ EPSDestroy(&eps);
+}
+
+void
+EigenvalueSolver::computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A,
+ Array<double>& eigenvalues,
+ DenseMatrix<double>& P)
+{
+ EPS eps;
+
+ EPSCreate(PETSC_COMM_SELF, &eps);
+ EPSSetOperators(eps, A, nullptr);
+ EPSSetProblemType(eps, EPS_HEP);
+
+ Internals::computeAllEigenvaluesOfSymmetricMatrix(eps);
+
+ PetscInt nb_eigenvalues;
+ EPSGetDimensions(eps, &nb_eigenvalues, nullptr, nullptr);
+
+ eigenvalues = Array<double>(nb_eigenvalues);
+ P = DenseMatrix<double>(nb_eigenvalues, nb_eigenvalues);
+
+ Array<double> eigenvector(nb_eigenvalues);
+ for (PetscInt i = 0; i < nb_eigenvalues; ++i) {
+ Vec Vr;
+ VecCreateSeqWithArray(PETSC_COMM_SELF, 1, A.numberOfRows(), &(eigenvector[0]), &Vr);
+ EPSGetEigenpair(eps, i, &(eigenvalues[i]), nullptr, Vr, nullptr);
+ VecDestroy(&Vr);
+ for (size_t j = 0; j < eigenvector.size(); ++j) {
+ P(j, i) = eigenvector[j];
+ }
+ }
+
+ EPSDestroy(&eps);
+}
+
+#endif // PUGS_HAS_SLEPC
+
+EigenvalueSolver::EigenvalueSolver() {}
diff --git a/src/algebra/EigenvalueSolver.hpp b/src/algebra/EigenvalueSolver.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..a9dd941b9bd87d3f0a79be5f36e5683ae8a1afb6
--- /dev/null
+++ b/src/algebra/EigenvalueSolver.hpp
@@ -0,0 +1,68 @@
+#ifndef EIGENVALUE_SOLVER_HPP
+#define EIGENVALUE_SOLVER_HPP
+
+#include <algebra/PETScUtils.hpp>
+
+#include <algebra/DenseMatrix.hpp>
+#include <algebra/Vector.hpp>
+#include <utils/Array.hpp>
+#include <utils/Exceptions.hpp>
+
+class EigenvalueSolver
+{
+ private:
+ struct Internals;
+
+#ifdef PUGS_HAS_SLEPC
+ void computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A, Array<double>& eigenvalues);
+
+ void computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A,
+ Array<double>& eigenvalues,
+ std::vector<Vector<double>>& eigenvectors);
+
+ void computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A, Array<double>& eigenvalues, DenseMatrix<double>& P);
+#endif // PUGS_HAS_SLEPC
+
+ public:
+ template <typename MatrixType>
+ void
+ computeForSymmetricMatrix([[maybe_unused]] const MatrixType& A, [[maybe_unused]] Array<double>& eigenvalues)
+ {
+#ifdef PUGS_HAS_SLEPC
+ this->computeForSymmetricMatrix(PETScAijMatrixEmbedder{A}, eigenvalues);
+#else // PUGS_HAS_SLEPC
+ throw NotImplementedError("SLEPc is required to solve eigenvalue problems");
+#endif // PUGS_HAS_SLEPC
+ }
+
+ template <typename MatrixType>
+ void
+ computeForSymmetricMatrix([[maybe_unused]] const MatrixType& A,
+ [[maybe_unused]] Array<double>& eigenvalues,
+ [[maybe_unused]] std::vector<Vector<double>>& eigenvectors)
+ {
+#ifdef PUGS_HAS_SLEPC
+ this->computeForSymmetricMatrix(PETScAijMatrixEmbedder{A}, eigenvalues, eigenvectors);
+#else // PUGS_HAS_SLEPC
+ throw NotImplementedError("SLEPc is required to solve eigenvalue problems");
+#endif // PUGS_HAS_SLEPC
+ }
+
+ template <typename MatrixType>
+ void
+ computeForSymmetricMatrix([[maybe_unused]] const MatrixType& A,
+ [[maybe_unused]] Array<double>& eigenvalues,
+ [[maybe_unused]] DenseMatrix<double>& P)
+ {
+#ifdef PUGS_HAS_SLEPC
+ this->computeForSymmetricMatrix(PETScAijMatrixEmbedder{A}, eigenvalues, P);
+#else // PUGS_HAS_SLEPC
+ throw NotImplementedError("SLEPc is required to solve eigenvalue problems");
+#endif // PUGS_HAS_SLEPC
+ }
+
+ EigenvalueSolver();
+ ~EigenvalueSolver() = default;
+};
+
+#endif // EIGENVALUE_SOLVER_HPP
diff --git a/src/algebra/LinearSolver.cpp b/src/algebra/LinearSolver.cpp
index b3a4d9816f6bc6634f73a80b70ba304f82ed60b1..98e2a3a1c2caf8a041b6e69e71502bb6441647e7 100644
--- a/src/algebra/LinearSolver.cpp
+++ b/src/algebra/LinearSolver.cpp
@@ -3,6 +3,7 @@
#include <algebra/BiCGStab.hpp>
#include <algebra/CG.hpp>
+#include <algebra/PETScUtils.hpp>
#ifdef PUGS_HAS_PETSC
#include <petsc.h>
@@ -130,10 +131,11 @@ struct LinearSolver::Internals
}
}
+ template <typename MatrixType, typename SolutionVectorType, typename RHSVectorType>
static void
- builtinSolveLocalSystem(const CRSMatrix<double, size_t>& A,
- Vector<double>& x,
- const Vector<double>& b,
+ builtinSolveLocalSystem(const MatrixType& A,
+ SolutionVectorType& x,
+ const RHSVectorType& b,
const LinearSolverOptions& options)
{
if (options.precond() != LSPrecond::none) {
@@ -166,48 +168,27 @@ struct LinearSolver::Internals
return 0;
}
+ template <typename MatrixType, typename SolutionVectorType, typename RHSVectorType>
static void
- petscSolveLocalSystem(const CRSMatrix<double, size_t>& A,
- Vector<double>& x,
- const Vector<double>& b,
+ petscSolveLocalSystem(const MatrixType& A,
+ SolutionVectorType& x,
+ const RHSVectorType& b,
const LinearSolverOptions& options)
{
- Assert(x.size() == b.size() and x.size() == A.numberOfRows());
+ Assert(x.size() == b.size() and x.size() == A.numberOfColumns() and A.isSquare());
Vec petscB;
- VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, b.size(), b.size(), &b[0], &petscB);
+ VecCreateMPIWithArray(PETSC_COMM_SELF, 1, b.size(), b.size(), &b[0], &petscB);
Vec petscX;
- VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, x.size(), x.size(), &x[0], &petscX);
+ VecCreateMPIWithArray(PETSC_COMM_SELF, 1, x.size(), x.size(), &x[0], &petscX);
- Array<PetscScalar> values = copy(A.values());
-
- const auto A_row_indices = A.rowIndices();
- Array<PetscInt> row_indices{A_row_indices.size()};
- for (size_t i = 0; i < row_indices.size(); ++i) {
- row_indices[i] = A_row_indices[i];
- }
-
- Array<PetscInt> column_indices{values.size()};
- size_t l = 0;
- for (size_t i = 0; i < A.numberOfRows(); ++i) {
- const auto row_i = A.row(i);
- for (size_t j = 0; j < row_i.length; ++j) {
- column_indices[l++] = row_i.colidx(j);
- }
- }
-
- Mat petscMat;
- 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);
+ PETScAijMatrixEmbedder petscA(A);
KSP ksp;
- KSPCreate(PETSC_COMM_WORLD, &ksp);
+ KSPCreate(PETSC_COMM_SELF, &ksp);
KSPSetTolerances(ksp, options.epsilon(), 1E-100, 1E5, options.maximumIteration());
- KSPSetOperators(ksp, petscMat, petscMat);
+ KSPSetOperators(ksp, petscA, petscA);
PC pc;
KSPGetPC(ksp, &pc);
@@ -288,8 +269,6 @@ struct LinearSolver::Internals
KSPSolve(ksp, petscB, petscX);
- // free used memory
- MatDestroy(&petscMat);
VecDestroy(&petscB);
VecDestroy(&petscX);
KSPDestroy(&ksp);
@@ -298,11 +277,9 @@ struct LinearSolver::Internals
#else // PUGS_HAS_PETSC
// LCOV_EXCL_START
+ template <typename MatrixType, typename SolutionVectorType, typename RHSVectorType>
static void
- petscSolveLocalSystem(const CRSMatrix<double, size_t>&,
- Vector<double>&,
- const Vector<double>&,
- const LinearSolverOptions&)
+ petscSolveLocalSystem(const MatrixType&, SolutionVectorType&, const RHSVectorType&, const LinearSolverOptions&)
{
checkHasLibrary(LSLibrary::petsc);
throw UnexpectedError("unexpected situation should not reach this point!");
@@ -310,6 +287,33 @@ 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::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
@@ -325,26 +329,15 @@ LinearSolver::checkOptions(const LinearSolverOptions& options) const
}
void
-LinearSolver::solveLocalSystem(const CRSMatrix<double, size_t>& A, Vector<double>& x, const Vector<double>& b)
+LinearSolver::solveLocalSystem(const CRSMatrix<double, size_t>& A, Vector<double>& x, const Vector<const double>& b)
{
- switch (m_options.library()) {
- case LSLibrary::builtin: {
- Internals::builtinSolveLocalSystem(A, x, b, m_options);
- break;
- }
- // LCOV_EXCL_START
- case LSLibrary::petsc: {
- // not covered since if PETSc is not linked this point is
- // unreachable: LinearSolver throws an exception at construction
- // in this case.
- Internals::petscSolveLocalSystem(A, x, b, m_options);
- break;
- }
- default: {
- throw UnexpectedError(::name(m_options.library()) + " cannot solve local systems for sparse matrices");
- }
- // LCOV_EXCL_STOP
- }
+ Internals::solveLocalSystem(A, x, b, m_options);
+}
+
+void
+LinearSolver::solveLocalSystem(const DenseMatrix<double>& A, Vector<double>& x, const Vector<const double>& b)
+{
+ Internals::solveLocalSystem(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 3ed5660bd9ce51c1e5ba5de3fe00c595aa96f1e7..3873374002a146c79a447a815d4ef6a57669883e 100644
--- a/src/algebra/LinearSolver.hpp
+++ b/src/algebra/LinearSolver.hpp
@@ -2,6 +2,7 @@
#define LINEAR_SOLVER_HPP
#include <algebra/CRSMatrix.hpp>
+#include <algebra/DenseMatrix.hpp>
#include <algebra/LinearSolverOptions.hpp>
#include <algebra/TinyMatrix.hpp>
#include <algebra/TinyVector.hpp>
@@ -15,21 +16,32 @@ class LinearSolver
const LinearSolverOptions m_options;
- void _solveLocalDense(size_t N, const double* A, double* x, const double* b);
-
public:
bool hasLibrary(LSLibrary library) const;
void checkOptions(const LinearSolverOptions& options) 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]);
+ DenseMatrix A_dense{A};
+ Vector<double> x_vector{N};
+ Vector<double> b_vector{N};
+ for (size_t i = 0; i < N; ++i) {
+ x_vector[i] = x[i];
+ b_vector[i] = b[i];
+ }
+
+ this->solveLocalSystem(A_dense, x_vector, b_vector);
+
+ for (size_t i = 0; i < N; ++i) {
+ x[i] = x_vector[i];
+ }
}
+ void solveLocalSystem(const CRSMatrix<double, size_t>& A, Vector<double>& x, const Vector<const double>& b);
+ void solveLocalSystem(const DenseMatrix<double>& A, Vector<double>& x, const Vector<const double>& b);
+
LinearSolver();
LinearSolver(const LinearSolverOptions& options);
};
diff --git a/src/algebra/PETScUtils.cpp b/src/algebra/PETScUtils.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..5ed261353d93673267c3f4107c35f855bf643211
--- /dev/null
+++ b/src/algebra/PETScUtils.cpp
@@ -0,0 +1,74 @@
+#include <algebra/PETScUtils.hpp>
+
+#ifdef PUGS_HAS_PETSC
+
+PETScAijMatrixEmbedder::PETScAijMatrixEmbedder(const size_t nb_rows, const size_t nb_columns, const double* A)
+ : m_nb_rows{nb_rows}, m_nb_columns{nb_columns}
+{
+ MatCreate(PETSC_COMM_SELF, &m_petscMat);
+ MatSetSizes(m_petscMat, PETSC_DECIDE, PETSC_DECIDE, nb_rows, nb_columns);
+ MatSetFromOptions(m_petscMat);
+ MatSetType(m_petscMat, MATAIJ);
+ MatSetUp(m_petscMat);
+
+ {
+ Array<PetscInt> row_indices(nb_rows);
+ for (size_t i = 0; i < nb_rows; ++i) {
+ row_indices[i] = i;
+ }
+ m_row_indices = row_indices;
+ }
+
+ if (nb_rows == nb_columns) {
+ m_column_indices = m_row_indices;
+ } else {
+ Array<PetscInt> column_indices(nb_columns);
+ for (size_t i = 0; i < nb_columns; ++i) {
+ column_indices[i] = i;
+ }
+ m_column_indices = column_indices;
+ }
+
+ MatSetValuesBlocked(m_petscMat, nb_rows, &(m_row_indices[0]), nb_columns, &(m_column_indices[0]), A, INSERT_VALUES);
+
+ MatAssemblyBegin(m_petscMat, MAT_FINAL_ASSEMBLY);
+ MatAssemblyEnd(m_petscMat, MAT_FINAL_ASSEMBLY);
+}
+
+PETScAijMatrixEmbedder::PETScAijMatrixEmbedder(const CRSMatrix<double, size_t>& A)
+ : m_nb_rows{A.numberOfRows()}, m_nb_columns{A.numberOfColumns()}
+{
+ const Array<PetscReal>& values = copy(A.values());
+
+ {
+ const auto& A_row_indices = A.rowIndices();
+ Array<PetscInt> row_indices{A_row_indices.size()};
+ for (size_t i = 0; i < row_indices.size(); ++i) {
+ row_indices[i] = A_row_indices[i];
+ }
+ m_row_indices = row_indices;
+
+ Array<PetscInt> column_indices{values.size()};
+ size_t l = 0;
+ for (size_t i = 0; i < A.numberOfRows(); ++i) {
+ const auto row_i = A.row(i);
+ for (size_t j = 0; j < row_i.length; ++j) {
+ column_indices[l++] = row_i.colidx(j);
+ }
+ }
+ m_column_indices = column_indices;
+ }
+
+ MatCreateSeqAIJWithArrays(PETSC_COMM_SELF, A.numberOfRows(), A.numberOfColumns(), &m_row_indices[0],
+ &m_column_indices[0], &values[0], &m_petscMat);
+
+ MatAssemblyBegin(m_petscMat, MAT_FINAL_ASSEMBLY);
+ MatAssemblyEnd(m_petscMat, MAT_FINAL_ASSEMBLY);
+}
+
+PETScAijMatrixEmbedder::~PETScAijMatrixEmbedder()
+{
+ MatDestroy(&m_petscMat);
+}
+
+#endif // PUGS_HAS_PETSC
diff --git a/src/algebra/PETScUtils.hpp b/src/algebra/PETScUtils.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..68748c02f5b36092cba6dce8ba3281a6d1cbb93b
--- /dev/null
+++ b/src/algebra/PETScUtils.hpp
@@ -0,0 +1,67 @@
+#ifndef PETSC_UTILS_HPP
+#define PETSC_UTILS_HPP
+
+#include <utils/pugs_config.hpp>
+
+#ifdef PUGS_HAS_PETSC
+
+#include <algebra/CRSMatrix.hpp>
+#include <algebra/DenseMatrix.hpp>
+#include <algebra/TinyMatrix.hpp>
+
+#include <petsc.h>
+
+class PETScAijMatrixEmbedder
+{
+ private:
+ Mat m_petscMat;
+ Array<PetscInt> m_row_indices;
+ Array<PetscInt> m_column_indices;
+ const size_t m_nb_rows;
+ const size_t m_nb_columns;
+
+ PETScAijMatrixEmbedder(const size_t nb_rows, const size_t nb_columns, const double* A);
+
+ public:
+ PUGS_INLINE
+ size_t
+ numberOfRows() const
+ {
+ return m_nb_rows;
+ }
+
+ PUGS_INLINE
+ size_t
+ numberOfColumns() const
+ {
+ return m_nb_columns;
+ }
+
+ PUGS_INLINE
+ operator Mat&()
+ {
+ return m_petscMat;
+ }
+
+ PUGS_INLINE
+ operator const Mat&() const
+ {
+ return m_petscMat;
+ }
+
+ template <size_t N>
+ PETScAijMatrixEmbedder(const TinyMatrix<N>& A) : PETScAijMatrixEmbedder{N, N, &A(0, 0)}
+ {}
+
+ PETScAijMatrixEmbedder(const DenseMatrix<double>& A)
+ : PETScAijMatrixEmbedder{A.numberOfRows(), A.numberOfColumns(), &A(0, 0)}
+ {}
+
+ PETScAijMatrixEmbedder(const CRSMatrix<double, size_t>& A);
+
+ ~PETScAijMatrixEmbedder();
+};
+
+#endif // PUGS_HAS_PETSC
+
+#endif // PETSC_UTILS_HPP
diff --git a/src/algebra/SparseMatrixDescriptor.hpp b/src/algebra/SparseMatrixDescriptor.hpp
index f63a17576d619bf692f36a00e5c0c4e373effdbf..d94ed468606a10ee15857a971abab20fb4be5b1a 100644
--- a/src/algebra/SparseMatrixDescriptor.hpp
+++ b/src/algebra/SparseMatrixDescriptor.hpp
@@ -63,6 +63,7 @@ class SparseMatrixDescriptor
private:
Array<SparseRowDescriptor> m_row_array;
+ const size_t m_nb_columns;
public:
PUGS_INLINE
@@ -72,6 +73,13 @@ class SparseMatrixDescriptor
return m_row_array.size();
}
+ PUGS_INLINE
+ size_t
+ numberOfColumns() const noexcept
+ {
+ return m_nb_columns;
+ }
+
SparseRowDescriptor&
row(const IndexType i)
{
@@ -90,7 +98,7 @@ class SparseMatrixDescriptor
operator()(const IndexType& i, const IndexType& j)
{
Assert(i < m_row_array.size());
- Assert(j < m_row_array.size());
+ Assert(j < m_nb_columns);
return m_row_array[i][j];
}
@@ -98,7 +106,7 @@ class SparseMatrixDescriptor
operator()(const IndexType& i, const IndexType& j) const
{
Assert(i < m_row_array.size());
- Assert(j < m_row_array.size());
+ Assert(j < m_nb_columns);
const auto& r = m_row_array[i]; // split to ensure const-ness of call
return r[j];
}
@@ -145,13 +153,17 @@ class SparseMatrixDescriptor
return values;
}
- SparseMatrixDescriptor(size_t nb_row) : m_row_array{nb_row}
+ SparseMatrixDescriptor(size_t nb_rows, size_t nb_columns) : m_row_array{nb_rows}, m_nb_columns{nb_columns}
{
- for (size_t i = 0; i < nb_row; ++i) {
- m_row_array[i][i] = 0;
+ if (nb_rows == nb_columns) {
+ for (size_t i = 0; i < nb_rows; ++i) {
+ m_row_array[i][i] = 0;
+ }
}
}
+ SparseMatrixDescriptor(size_t nb_rows) : SparseMatrixDescriptor{nb_rows, nb_rows} {}
+
~SparseMatrixDescriptor() = default;
};
diff --git a/src/algebra/TinyMatrix.hpp b/src/algebra/TinyMatrix.hpp
index 636ff19c18112c582c8968e8dc6f4df5ebdee726..6c23b068bf1a8c8744e24ffbf60ec88b42b1317b 100644
--- a/src/algebra/TinyMatrix.hpp
+++ b/src/algebra/TinyMatrix.hpp
@@ -36,6 +36,12 @@ class [[nodiscard]] TinyMatrix
}
public:
+ PUGS_INLINE
+ constexpr bool isSquare() const noexcept
+ {
+ return true;
+ }
+
PUGS_INLINE
constexpr size_t dimension() const
{
@@ -43,13 +49,19 @@ class [[nodiscard]] TinyMatrix
}
PUGS_INLINE
- constexpr size_t nbRows() const
+ constexpr size_t numberOfValues() const
+ {
+ return this->dimension();
+ }
+
+ PUGS_INLINE
+ constexpr size_t numberOfRows() const
{
return N;
}
PUGS_INLINE
- constexpr size_t nbColumns() const
+ constexpr size_t numberOfColumns() const
{
return N;
}
diff --git a/src/algebra/Vector.hpp b/src/algebra/Vector.hpp
index 089efa633161f4d3e83d16a5fc666ad5f9fc1696..e9eace62afa018c31c9f727e779a01c9ce39ccd6 100644
--- a/src/algebra/Vector.hpp
+++ b/src/algebra/Vector.hpp
@@ -17,12 +17,23 @@ class Vector // LCOV_EXCL_LINE
private:
Array<DataType> m_values;
+ const bool m_deep_copies;
+
static_assert(std::is_same_v<typename decltype(m_values)::index_type, index_type>);
// Allows const version to access our data
friend Vector<std::add_const_t<DataType>>;
public:
+ friend std::ostream&
+ operator<<(std::ostream& os, const Vector& x)
+ {
+ for (size_t i = 0; i < x.size(); ++i) {
+ os << ' ' << x[i];
+ }
+ return os;
+ }
+
friend Vector<std::remove_const_t<DataType>>
copy(const Vector& x)
{
@@ -56,7 +67,7 @@ class Vector // LCOV_EXCL_LINE
// Does not use parallel_for to preserve sum order
for (index_type i = 0; i < x.size(); ++i) {
- sum += x.m_values[i] * y.m_values[i];
+ sum += x[i] * y[i];
}
return sum;
@@ -152,20 +163,49 @@ class Vector // LCOV_EXCL_LINE
template <typename DataType2>
PUGS_INLINE Vector&
- operator=(const Vector<DataType2>& vector) noexcept
+ operator=(const Vector<DataType2>& x) noexcept
{
- m_values = vector.m_values;
+ // ensures that DataType is the same as source DataType2
+ static_assert(std::is_same<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>(),
+ "Cannot assign Vector of different type");
+ // ensures that const is not lost through copy
+ static_assert(((std::is_const<DataType2>() and std::is_const<DataType>()) or not std::is_const<DataType2>()),
+ "Cannot assign Vector of const to Vector of non-const");
+
+ if (m_deep_copies) {
+ copy_to(x.m_values, m_values);
+ } else {
+ m_values = x.m_values;
+ }
return *this;
}
PUGS_INLINE
- Vector& operator=(const Vector&) = default;
+ Vector&
+ operator=(const Vector& x)
+ {
+ if (m_deep_copies) {
+ copy_to(x.m_values, m_values);
+ } else {
+ m_values = x.m_values;
+ }
+ return *this;
+ }
PUGS_INLINE
- Vector& operator=(Vector&&) = default;
+ Vector&
+ operator=(Vector&& x)
+ {
+ if (m_deep_copies) {
+ copy_to(x.m_values, m_values);
+ } else {
+ m_values = x.m_values;
+ }
+ return *this;
+ }
template <typename DataType2>
- Vector(const Vector<DataType2>& x)
+ Vector(const Vector<DataType2>& x) : m_deep_copies{x.m_deep_copies}
{
// ensures that DataType is the same as source DataType2
static_assert(std::is_same<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>(),
@@ -174,14 +214,19 @@ class Vector // LCOV_EXCL_LINE
static_assert(((std::is_const<DataType2>() and std::is_const<DataType>()) or not std::is_const<DataType2>()),
"Cannot assign Vector of const to Vector of non-const");
- this->operator=(x);
+ m_values = x.m_values;
+ }
+
+ explicit Vector(const Array<DataType>& values) : m_deep_copies{true}
+ {
+ m_values = values;
}
Vector(const Vector&) = default;
Vector(Vector&&) = default;
- Vector(size_t size) : m_values{size} {}
+ Vector(size_t size) : m_values{size}, m_deep_copies{false} {}
~Vector() = default;
};
diff --git a/src/language/ast/ASTNodeArraySubscriptExpressionBuilder.cpp b/src/language/ast/ASTNodeArraySubscriptExpressionBuilder.cpp
index e717044d1dea5fe6b85aec6c41941699c1e5ef65..dcb2441c8bf14452c1f7af56a6a55bcdf340ede0 100644
--- a/src/language/ast/ASTNodeArraySubscriptExpressionBuilder.cpp
+++ b/src/language/ast/ASTNodeArraySubscriptExpressionBuilder.cpp
@@ -29,9 +29,9 @@ ASTNodeArraySubscriptExpressionBuilder::ASTNodeArraySubscriptExpressionBuilder(A
}
}
} else if (array_expression.m_data_type == ASTNodeDataType::matrix_t) {
- Assert(array_expression.m_data_type.nbRows() == array_expression.m_data_type.nbColumns());
+ Assert(array_expression.m_data_type.numberOfRows() == array_expression.m_data_type.numberOfColumns());
- switch (array_expression.m_data_type.nbRows()) {
+ switch (array_expression.m_data_type.numberOfRows()) {
case 1: {
node.m_node_processor = std::make_unique<ArraySubscriptProcessor<TinyMatrix<1>>>(node);
break;
diff --git a/src/language/ast/ASTNodeBuiltinFunctionExpressionBuilder.cpp b/src/language/ast/ASTNodeBuiltinFunctionExpressionBuilder.cpp
index 412906b9d67291ed7901d3ac00943cc8a9107449..70d8e3b3ac10456eb28c0036fdfa31eb1630bcab 100644
--- a/src/language/ast/ASTNodeBuiltinFunctionExpressionBuilder.cpp
+++ b/src/language/ast/ASTNodeBuiltinFunctionExpressionBuilder.cpp
@@ -120,8 +120,8 @@ ASTNodeBuiltinFunctionExpressionBuilder::_getArgumentConverter(const ASTNodeData
if constexpr (std::is_same_v<ParameterT, TinyMatrix<1>>) {
switch (argument_node_sub_data_type.m_data_type) {
case ASTNodeDataType::matrix_t: {
- if ((argument_node_sub_data_type.m_data_type.nbRows() == 1) and
- (argument_node_sub_data_type.m_data_type.nbColumns() == 1)) {
+ if ((argument_node_sub_data_type.m_data_type.numberOfRows() == 1) and
+ (argument_node_sub_data_type.m_data_type.numberOfColumns() == 1)) {
return std::make_unique<FunctionTinyMatrixArgumentConverter<ParameterT, ParameterT>>(argument_number);
} else {
// LCOV_EXCL_START
@@ -152,8 +152,8 @@ ASTNodeBuiltinFunctionExpressionBuilder::_getArgumentConverter(const ASTNodeData
} else {
switch (argument_node_sub_data_type.m_data_type) {
case ASTNodeDataType::matrix_t: {
- if ((argument_node_sub_data_type.m_data_type.nbRows() == parameter_v.nbRows()) and
- (argument_node_sub_data_type.m_data_type.nbColumns() == parameter_v.nbColumns())) {
+ if ((argument_node_sub_data_type.m_data_type.numberOfRows() == parameter_v.numberOfRows()) and
+ (argument_node_sub_data_type.m_data_type.numberOfColumns() == parameter_v.numberOfColumns())) {
return std::make_unique<FunctionTinyMatrixArgumentConverter<ParameterT, ParameterT>>(argument_number);
} else {
// LCOV_EXCL_START
@@ -164,7 +164,7 @@ ASTNodeBuiltinFunctionExpressionBuilder::_getArgumentConverter(const ASTNodeData
}
case ASTNodeDataType::list_t: {
if (argument_node_sub_data_type.m_parent_node.children.size() ==
- (parameter_v.nbRows() * parameter_v.nbColumns())) {
+ (parameter_v.numberOfRows() * parameter_v.numberOfColumns())) {
return std::make_unique<FunctionTinyMatrixArgumentConverter<ParameterT, ParameterT>>(argument_number);
} else {
// LCOV_EXCL_START
@@ -293,8 +293,8 @@ ASTNodeBuiltinFunctionExpressionBuilder::_getArgumentConverter(const ASTNodeData
}
}
case ASTNodeDataType::matrix_t: {
- Assert(arg_data_type.nbRows() == arg_data_type.nbColumns());
- switch (arg_data_type.nbRows()) {
+ Assert(arg_data_type.numberOfRows() == arg_data_type.numberOfColumns());
+ switch (arg_data_type.numberOfRows()) {
case 1: {
return std::make_unique<FunctionTupleArgumentConverter<ParameterContentT, TinyMatrix<1>>>(argument_number);
}
@@ -370,8 +370,8 @@ ASTNodeBuiltinFunctionExpressionBuilder::_getArgumentConverter(const ASTNodeData
}
}
case ASTNodeDataType::matrix_t: {
- Assert(parameter_type.nbRows() == parameter_type.nbColumns());
- switch (parameter_type.nbRows()) {
+ Assert(parameter_type.numberOfRows() == parameter_type.numberOfColumns());
+ switch (parameter_type.numberOfRows()) {
case 1: {
return get_function_argument_converter_for_matrix(TinyMatrix<1>{});
}
@@ -439,8 +439,8 @@ ASTNodeBuiltinFunctionExpressionBuilder::_getArgumentConverter(const ASTNodeData
}
}
case ASTNodeDataType::matrix_t: {
- Assert(parameter_type.contentType().nbRows() == parameter_type.contentType().nbColumns());
- switch (parameter_type.contentType().nbRows()) {
+ Assert(parameter_type.contentType().numberOfRows() == parameter_type.contentType().numberOfColumns());
+ switch (parameter_type.contentType().numberOfRows()) {
case 1: {
return get_function_argument_to_tuple_converter(TinyMatrix<1>{});
}
diff --git a/src/language/ast/ASTNodeDataTypeBuilder.cpp b/src/language/ast/ASTNodeDataTypeBuilder.cpp
index d8cf06824ea5cf070a07a84d4bf3bcc711f71a79..eb307e2c16893be52128166c973f4992e86f3cb4 100644
--- a/src/language/ast/ASTNodeDataTypeBuilder.cpp
+++ b/src/language/ast/ASTNodeDataTypeBuilder.cpp
@@ -266,10 +266,11 @@ ASTNodeDataTypeBuilder::_buildNodeDataTypes(ASTNode& n) const
}
} else if (image_domain_node.is_type<language::matrix_type>()) {
ASTNodeDataType image_type = getMatrixDataType(image_domain_node);
- if (image_type.nbRows() * image_type.nbColumns() != nb_image_expressions) {
+ if (image_type.numberOfRows() * image_type.numberOfColumns() != nb_image_expressions) {
std::ostringstream message;
- message << "expecting " << image_type.nbRows() * image_type.nbColumns() << " scalar expressions or an "
- << dataTypeName(image_type) << ", found " << nb_image_expressions << " scalar expressions";
+ message << "expecting " << image_type.numberOfRows() * image_type.numberOfColumns()
+ << " scalar expressions or an " << dataTypeName(image_type) << ", found " << nb_image_expressions
+ << " scalar expressions";
throw ParseError(message.str(), image_expression_node.begin());
}
} else {
diff --git a/src/language/ast/ASTNodeFunctionExpressionBuilder.cpp b/src/language/ast/ASTNodeFunctionExpressionBuilder.cpp
index 15720a9f645616ad08588c5d00e5185d06482276..6c1db0335b2f46a8410baed7f02e59efa6ae49c8 100644
--- a/src/language/ast/ASTNodeFunctionExpressionBuilder.cpp
+++ b/src/language/ast/ASTNodeFunctionExpressionBuilder.cpp
@@ -112,8 +112,8 @@ ASTNodeFunctionExpressionBuilder::_getArgumentConverter(SymbolType& parameter_sy
using ParameterT = std::decay_t<decltype(parameter_v)>;
switch (node_sub_data_type.m_data_type) {
case ASTNodeDataType::matrix_t: {
- if ((node_sub_data_type.m_data_type.nbRows() == parameter_v.nbRows()) and
- (node_sub_data_type.m_data_type.nbColumns() == parameter_v.nbColumns())) {
+ if ((node_sub_data_type.m_data_type.numberOfRows() == parameter_v.numberOfRows()) and
+ (node_sub_data_type.m_data_type.numberOfColumns() == parameter_v.numberOfColumns())) {
return std::make_unique<FunctionTinyMatrixArgumentConverter<ParameterT, ParameterT>>(parameter_id);
} else {
// LCOV_EXCL_START
@@ -187,9 +187,10 @@ ASTNodeFunctionExpressionBuilder::_getArgumentConverter(SymbolType& parameter_sy
// LCOV_EXCL_STOP
}
case ASTNodeDataType::matrix_t: {
- Assert(parameter_symbol.attributes().dataType().nbRows() == parameter_symbol.attributes().dataType().nbColumns());
+ Assert(parameter_symbol.attributes().dataType().numberOfRows() ==
+ parameter_symbol.attributes().dataType().numberOfColumns());
- switch (parameter_symbol.attributes().dataType().nbRows()) {
+ switch (parameter_symbol.attributes().dataType().numberOfRows()) {
case 1: {
return get_function_argument_converter_for_matrix(TinyMatrix<1>{});
}
@@ -361,8 +362,8 @@ ASTNodeFunctionExpressionBuilder::_getFunctionProcessor(const ASTNodeDataType& r
using ReturnT = std::decay_t<decltype(return_v)>;
switch (function_component_expression.m_data_type) {
case ASTNodeDataType::matrix_t: {
- if ((function_component_expression.m_data_type.nbRows() == return_v.nbRows()) and
- (function_component_expression.m_data_type.nbColumns() == return_v.nbColumns())) {
+ if ((function_component_expression.m_data_type.numberOfRows() == return_v.numberOfRows()) and
+ (function_component_expression.m_data_type.numberOfColumns() == return_v.numberOfColumns())) {
return std::make_unique<FunctionExpressionProcessor<ReturnT, ReturnT>>(function_component_expression);
} else {
// LCOV_EXCL_START
@@ -461,9 +462,9 @@ ASTNodeFunctionExpressionBuilder::_getFunctionProcessor(const ASTNodeDataType& r
}
}
case ASTNodeDataType::matrix_t: {
- Assert(return_value_type.nbRows() == return_value_type.nbColumns());
+ Assert(return_value_type.numberOfRows() == return_value_type.numberOfColumns());
- switch (return_value_type.nbRows()) {
+ switch (return_value_type.numberOfRows()) {
case 1: {
if (function_component_expression.m_data_type == ASTNodeDataType::matrix_t) {
return get_function_processor_for_expression_matrix(TinyMatrix<1>{});
@@ -624,15 +625,15 @@ ASTNodeFunctionExpressionBuilder::ASTNodeFunctionExpressionBuilder(ASTNode& node
ASTNodeNaturalConversionChecker<AllowRToR1Conversion>{function_expression, matrix_type};
if (function_expression.is_type<language::expression_list>()) {
- Assert(matrix_type.nbRows() * matrix_type.nbColumns() == function_expression.children.size());
+ Assert(matrix_type.numberOfRows() * matrix_type.numberOfColumns() == function_expression.children.size());
- for (size_t i = 0; i < matrix_type.nbRows() * matrix_type.nbColumns(); ++i) {
+ for (size_t i = 0; i < matrix_type.numberOfRows() * matrix_type.numberOfColumns(); ++i) {
function_processor->addFunctionExpressionProcessor(
this->_getFunctionProcessor(ASTNodeDataType::build<ASTNodeDataType::double_t>(), node,
*function_expression.children[i]));
}
- switch (matrix_type.nbRows()) {
+ switch (matrix_type.numberOfRows()) {
case 2: {
node.m_node_processor =
std::make_unique<TupleToTinyMatrixProcessor<FunctionProcessor, 2>>(node, std::move(function_processor));
@@ -651,7 +652,7 @@ ASTNodeFunctionExpressionBuilder::ASTNodeFunctionExpressionBuilder(ASTNode& node
}
} else if (function_expression.is_type<language::integer>()) {
if (std::stoi(function_expression.string()) == 0) {
- switch (matrix_type.nbRows()) {
+ switch (matrix_type.numberOfRows()) {
case 1: {
node.m_node_processor =
std::make_unique<FunctionExpressionProcessor<TinyMatrix<1>, ZeroType>>(function_expression);
diff --git a/src/language/ast/ASTNodeListAffectationExpressionBuilder.cpp b/src/language/ast/ASTNodeListAffectationExpressionBuilder.cpp
index 6b53db72a1a257f74196a56f6e2c83694531fd5f..dc438e61f5aa8ea2ad0f4d3d2b0879c0e2ec0f0e 100644
--- a/src/language/ast/ASTNodeListAffectationExpressionBuilder.cpp
+++ b/src/language/ast/ASTNodeListAffectationExpressionBuilder.cpp
@@ -85,19 +85,19 @@ ASTNodeListAffectationExpressionBuilder::_buildAffectationProcessor(
using ValueT = std::decay_t<decltype(value)>;
if constexpr (std::is_same_v<ValueT, TinyMatrix<1>>) {
if ((node_sub_data_type.m_data_type == ASTNodeDataType::matrix_t) and
- (node_sub_data_type.m_data_type.nbRows() == value.nbRows()) and
- (node_sub_data_type.m_data_type.nbColumns() == value.nbColumns())) {
+ (node_sub_data_type.m_data_type.numberOfRows() == value.numberOfRows()) and
+ (node_sub_data_type.m_data_type.numberOfColumns() == value.numberOfColumns())) {
list_affectation_processor->template add<ValueT, ValueT>(value_node);
} else {
add_affectation_processor_for_data(value, node_sub_data_type);
}
} else if constexpr (std::is_same_v<ValueT, TinyMatrix<2>> or std::is_same_v<ValueT, TinyMatrix<3>>) {
if ((node_sub_data_type.m_data_type == ASTNodeDataType::matrix_t) and
- (node_sub_data_type.m_data_type.nbRows() == value.nbRows()) and
- (node_sub_data_type.m_data_type.nbColumns() == value.nbColumns())) {
+ (node_sub_data_type.m_data_type.numberOfRows() == value.numberOfRows()) and
+ (node_sub_data_type.m_data_type.numberOfColumns() == value.numberOfColumns())) {
list_affectation_processor->template add<ValueT, ValueT>(value_node);
} else if ((node_sub_data_type.m_data_type == ASTNodeDataType::list_t) and
- (node_sub_data_type.m_parent_node.children.size() == value.nbRows() * value.nbColumns())) {
+ (node_sub_data_type.m_parent_node.children.size() == value.numberOfRows() * value.numberOfColumns())) {
list_affectation_processor->template add<ValueT, AggregateDataVariant>(value_node);
} else if (node_sub_data_type.m_parent_node.is_type<language::integer>()) {
if (std::stoi(node_sub_data_type.m_parent_node.string()) == 0) {
@@ -241,8 +241,8 @@ ASTNodeListAffectationExpressionBuilder::_buildAffectationProcessor(
break;
}
case ASTNodeDataType::matrix_t: {
- Assert(value_type.nbRows() == value_type.nbColumns());
- switch (value_type.nbRows()) {
+ Assert(value_type.numberOfRows() == value_type.numberOfColumns());
+ switch (value_type.numberOfRows()) {
case 1: {
add_affectation_processor_for_matrix_data(TinyMatrix<1>{}, node_sub_data_type);
break;
diff --git a/src/language/node_processor/AffectationProcessor.hpp b/src/language/node_processor/AffectationProcessor.hpp
index 07a52de44e09cba57cf5b60aba65e844ff7faa94..f24e8473be7717f4ddd2e8f20d3d2ab666f19001 100644
--- a/src/language/node_processor/AffectationProcessor.hpp
+++ b/src/language/node_processor/AffectationProcessor.hpp
@@ -141,8 +141,8 @@ class AffectationExecutor final : public IAffectationExecutor
v[i]);
}
} else if constexpr (is_tiny_matrix_v<ValueT>) {
- for (size_t i = 0, l = 0; i < m_lhs.nbRows(); ++i) {
- for (size_t j = 0; j < m_lhs.nbColumns(); ++j, ++l) {
+ for (size_t i = 0, l = 0; i < m_lhs.numberOfRows(); ++i) {
+ for (size_t j = 0; j < m_lhs.numberOfColumns(); ++j, ++l) {
std::visit(
[&](auto&& Aij) {
using Aij_T = std::decay_t<decltype(Aij)>;
@@ -435,12 +435,12 @@ class AffectationProcessor final : public INodeProcessor
}
} else if (array_expression.m_data_type == ASTNodeDataType::matrix_t) {
Assert(lhs_node.children.size() == 3);
- Assert(array_expression.m_data_type.nbRows() == array_expression.m_data_type.nbColumns());
+ Assert(array_expression.m_data_type.numberOfRows() == array_expression.m_data_type.numberOfColumns());
auto& index0_expression = *lhs_node.children[1];
auto& index1_expression = *lhs_node.children[2];
- switch (array_expression.m_data_type.nbRows()) {
+ switch (array_expression.m_data_type.numberOfRows()) {
case 1: {
using ArrayTypeT = TinyMatrix<1>;
if (not std::holds_alternative<ArrayTypeT>(value)) {
@@ -576,8 +576,8 @@ class AffectationToTinyMatrixFromListProcessor final : public AffectationToDataV
static_assert(std::is_same_v<OperatorT, language::eq_op>, "forbidden affection operator for list to vectors");
ValueT v;
- for (size_t i = 0, l = 0; i < v.nbRows(); ++i) {
- for (size_t j = 0; j < v.nbColumns(); ++j, ++l) {
+ for (size_t i = 0, l = 0; i < v.numberOfRows(); ++i) {
+ for (size_t j = 0; j < v.numberOfColumns(); ++j, ++l) {
std::visit(
[&](auto&& child_value) {
using T = std::decay_t<decltype(child_value)>;
@@ -721,9 +721,9 @@ class AffectationToTupleFromListProcessor final : public AffectationToDataVarian
} else if constexpr (is_tiny_matrix_v<ValueT>) {
if constexpr (std::is_same_v<T, AggregateDataVariant>) {
ValueT& A = tuple_value[i];
- Assert(A.nbRows() * A.nbColumns() == child_value.size());
- for (size_t j = 0, l = 0; j < A.nbRows(); ++j) {
- for (size_t k = 0; k < A.nbColumns(); ++k, ++l) {
+ Assert(A.numberOfRows() * A.numberOfColumns() == child_value.size());
+ for (size_t j = 0, l = 0; j < A.numberOfRows(); ++j) {
+ for (size_t k = 0; k < A.numberOfColumns(); ++k, ++l) {
std::visit(
[&](auto&& Ajk) {
using Ti = std::decay_t<decltype(Ajk)>;
@@ -927,9 +927,9 @@ class ListAffectationProcessor final : public INodeProcessor
auto& index0_expression = *array_subscript_expression.children[1];
auto& index1_expression = *array_subscript_expression.children[2];
- Assert(array_expression.m_data_type.nbRows() == array_expression.m_data_type.nbColumns());
+ Assert(array_expression.m_data_type.numberOfRows() == array_expression.m_data_type.numberOfColumns());
- switch (array_expression.m_data_type.nbRows()) {
+ switch (array_expression.m_data_type.numberOfRows()) {
case 1: {
using ArrayTypeT = TinyMatrix<1>;
if (not std::holds_alternative<ArrayTypeT>(value)) {
diff --git a/src/language/node_processor/ExecutionPolicy.hpp b/src/language/node_processor/ExecutionPolicy.hpp
index 47fc3174ca5236c067bde9045a97f4d1c0fe160a..672daf5915808232285d76b2f6ea8492a3bdb7c2 100644
--- a/src/language/node_processor/ExecutionPolicy.hpp
+++ b/src/language/node_processor/ExecutionPolicy.hpp
@@ -34,12 +34,14 @@ class ExecutionPolicy
return m_shared_values->size();
}
- DataVariant& operator[](size_t i)
+ DataVariant&
+ operator[](size_t i)
{
return (*m_shared_values)[i];
}
- const DataVariant& operator[](size_t i) const
+ const DataVariant&
+ operator[](size_t i) const
{
return (*m_shared_values)[i];
}
diff --git a/src/language/node_processor/FunctionArgumentConverter.hpp b/src/language/node_processor/FunctionArgumentConverter.hpp
index b13eb222b08cd2c0f5dde5e72d4de5626ef7758b..f747672e2c036ab59402610cb7170218f96c15e3 100644
--- a/src/language/node_processor/FunctionArgumentConverter.hpp
+++ b/src/language/node_processor/FunctionArgumentConverter.hpp
@@ -139,8 +139,8 @@ class FunctionTinyMatrixArgumentConverter final : public IFunctionArgumentConver
exec_policy.currentContext()[m_argument_id] = std::move(value);
} else if constexpr (std::is_same_v<ValueT, AggregateDataVariant>) {
ExpectedValueType matrix_value{};
- for (size_t i = 0, l = 0; i < matrix_value.nbRows(); ++i) {
- for (size_t j = 0; j < matrix_value.nbColumns(); ++j, ++l) {
+ for (size_t i = 0, l = 0; i < matrix_value.numberOfRows(); ++i) {
+ for (size_t j = 0; j < matrix_value.numberOfColumns(); ++j, ++l) {
std::visit(
[&](auto&& A_ij) {
using Vi_T = std::decay_t<decltype(A_ij)>;
diff --git a/src/language/node_processor/FunctionProcessor.hpp b/src/language/node_processor/FunctionProcessor.hpp
index f8c86ea9e5cb51e169e48ecadc0b4770ef873cc8..cead55c20fd2e6b34edb4df510840d683fc8b5b8 100644
--- a/src/language/node_processor/FunctionProcessor.hpp
+++ b/src/language/node_processor/FunctionProcessor.hpp
@@ -38,8 +38,8 @@ class FunctionExpressionProcessor final : public INodeProcessor
} else {
static_assert(is_tiny_matrix_v<ReturnType>);
- for (size_t i = 0, l = 0; i < return_value.nbRows(); ++i) {
- for (size_t j = 0; j < return_value.nbColumns(); ++j, ++l) {
+ for (size_t i = 0, l = 0; i < return_value.numberOfRows(); ++i) {
+ for (size_t j = 0; j < return_value.numberOfColumns(); ++j, ++l) {
std::visit(
[&](auto&& Aij) {
using Vi_T = std::decay_t<decltype(Aij)>;
diff --git a/src/language/utils/ASTNodeDataType.cpp b/src/language/utils/ASTNodeDataType.cpp
index 3550025dd36e05c46e8ca20ce49a45d6b0fff30b..ddea3a6dd9a6a79fda474bb93e05645d55dd7503 100644
--- a/src/language/utils/ASTNodeDataType.cpp
+++ b/src/language/utils/ASTNodeDataType.cpp
@@ -76,7 +76,7 @@ dataTypeName(const ASTNodeDataType& data_type)
name = "R^" + std::to_string(data_type.dimension());
break;
case ASTNodeDataType::matrix_t:
- name = "R^" + std::to_string(data_type.nbRows()) + "x" + std::to_string(data_type.nbColumns());
+ name = "R^" + std::to_string(data_type.numberOfRows()) + "x" + std::to_string(data_type.numberOfColumns());
break;
case ASTNodeDataType::tuple_t:
name = "(" + dataTypeName(data_type.contentType()) + "...)";
@@ -173,8 +173,8 @@ isNaturalConversion(const ASTNodeDataType& data_type, const ASTNodeDataType& tar
} else if (data_type == ASTNodeDataType::vector_t) {
return (data_type.dimension() == target_data_type.dimension());
} else if (data_type == ASTNodeDataType::matrix_t) {
- return ((data_type.nbRows() == target_data_type.nbRows()) and
- (data_type.nbColumns() == target_data_type.nbColumns()));
+ return ((data_type.numberOfRows() == target_data_type.numberOfRows()) and
+ (data_type.numberOfColumns() == target_data_type.numberOfColumns()));
} else {
return true;
}
diff --git a/src/language/utils/ASTNodeDataType.hpp b/src/language/utils/ASTNodeDataType.hpp
index c31c9e15b6e8df0055aa4854bb7ede0a73318abd..ea744757cf68b5055605111d35acc14564b06b30 100644
--- a/src/language/utils/ASTNodeDataType.hpp
+++ b/src/language/utils/ASTNodeDataType.hpp
@@ -70,7 +70,7 @@ class ASTNodeDataType
PUGS_INLINE
size_t
- nbRows() const
+ numberOfRows() const
{
Assert(std::holds_alternative<std::array<size_t, 2>>(m_details));
return std::get<std::array<size_t, 2>>(m_details)[0];
@@ -78,7 +78,7 @@ class ASTNodeDataType
PUGS_INLINE
size_t
- nbColumns() const
+ numberOfColumns() const
{
Assert(std::holds_alternative<std::array<size_t, 2>>(m_details));
return std::get<std::array<size_t, 2>>(m_details)[1];
diff --git a/src/language/utils/ASTNodeNaturalConversionChecker.cpp b/src/language/utils/ASTNodeNaturalConversionChecker.cpp
index 81a277574753dc3e8cc75871a07d8d8a3b1a0cc6..a05c5aa1b4fcc79d531ea51544fb47dc1b40698e 100644
--- a/src/language/utils/ASTNodeNaturalConversionChecker.cpp
+++ b/src/language/utils/ASTNodeNaturalConversionChecker.cpp
@@ -14,8 +14,8 @@ ASTNodeNaturalConversionChecker<RToR1Conversion>::_checkIsNaturalTypeConversion(
if (not isNaturalConversion(data_type, target_data_type)) {
if constexpr (std::is_same_v<RToR1ConversionStrategy, AllowRToR1Conversion>) {
if (((target_data_type == ASTNodeDataType::vector_t) and (target_data_type.dimension() == 1)) or
- ((target_data_type == ASTNodeDataType::matrix_t) and (target_data_type.nbRows() == 1) and
- (target_data_type.nbColumns() == 1))) {
+ ((target_data_type == ASTNodeDataType::matrix_t) and (target_data_type.numberOfRows() == 1) and
+ (target_data_type.numberOfColumns() == 1))) {
if (isNaturalConversion(data_type, ASTNodeDataType::build<ASTNodeDataType::double_t>())) {
return;
}
@@ -93,10 +93,11 @@ ASTNodeNaturalConversionChecker<RToR1Conversion>::_checkIsNaturalExpressionConve
switch (data_type) {
case ASTNodeDataType::list_t: {
const auto& content_type_list = data_type.contentTypeList();
- if (content_type_list.size() != (target_data_type.nbRows() * target_data_type.nbColumns())) {
+ if (content_type_list.size() != (target_data_type.numberOfRows() * target_data_type.numberOfColumns())) {
std::ostringstream os;
os << "incompatible dimensions in affectation: expecting "
- << target_data_type.nbRows() * target_data_type.nbColumns() << ", but provided " << content_type_list.size();
+ << target_data_type.numberOfRows() * target_data_type.numberOfColumns() << ", but provided "
+ << content_type_list.size();
throw ParseError(os.str(), std::vector{node.begin()});
}
@@ -112,8 +113,8 @@ ASTNodeNaturalConversionChecker<RToR1Conversion>::_checkIsNaturalExpressionConve
break;
}
case ASTNodeDataType::matrix_t: {
- if ((data_type.nbRows() != target_data_type.nbRows()) or
- (data_type.nbColumns() != target_data_type.nbColumns())) {
+ if ((data_type.numberOfRows() != target_data_type.numberOfRows()) or
+ (data_type.numberOfColumns() != target_data_type.numberOfColumns())) {
std::ostringstream error_message;
error_message << "invalid implicit conversion: ";
error_message << rang::fgB::red << dataTypeName(data_type) << " -> " << dataTypeName(target_data_type)
diff --git a/src/language/utils/BuiltinFunctionEmbedderUtils.cpp b/src/language/utils/BuiltinFunctionEmbedderUtils.cpp
index ffe1112d0543fa644aa952a9599d9fb12eb5e383..1c7723ff1c3f2b510db742dd9b647f643b83d303 100644
--- a/src/language/utils/BuiltinFunctionEmbedderUtils.cpp
+++ b/src/language/utils/BuiltinFunctionEmbedderUtils.cpp
@@ -109,13 +109,13 @@ getBuiltinFunctionEmbedder(ASTNode& n)
}
bool is_castable = true;
- if (target_type.nbRows() > 1) {
+ if (target_type.numberOfRows() > 1) {
switch (arg_type) {
case ASTNodeDataType::int_t: {
break;
}
case ASTNodeDataType::list_t: {
- if (arg_type.contentTypeList().size() != target_type.nbRows() * target_type.nbColumns()) {
+ if (arg_type.contentTypeList().size() != target_type.numberOfRows() * target_type.numberOfColumns()) {
is_castable = false;
break;
}
diff --git a/src/language/utils/EmbeddedIDiscreteFunctionOperators.cpp b/src/language/utils/EmbeddedIDiscreteFunctionOperators.cpp
index 0e5a789ce31e645b8ccdde6c48009eca4275d76b..13def86d07990370bb7dce658dde09f8d06b565a 100644
--- a/src/language/utils/EmbeddedIDiscreteFunctionOperators.cpp
+++ b/src/language/utils/EmbeddedIDiscreteFunctionOperators.cpp
@@ -58,8 +58,8 @@ applyUnaryOperation(const std::shared_ptr<const IDiscreteFunction>& f)
}
}
case ASTNodeDataType::matrix_t: {
- Assert(f->dataType().nbRows() == f->dataType().nbColumns());
- switch (f->dataType().nbRows()) {
+ Assert(f->dataType().numberOfRows() == f->dataType().numberOfColumns());
+ switch (f->dataType().numberOfRows()) {
case 1: {
auto fh = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>&>(*f);
return applyUnaryOperation<UnaryOperatorT>(fh);
@@ -205,8 +205,8 @@ innerCompositionLaw(const std::shared_ptr<const IDiscreteFunction>& f,
}
case ASTNodeDataType::matrix_t: {
Assert(f->descriptor().type() == DiscreteFunctionType::P0);
- Assert(f->dataType().nbRows() == f->dataType().nbColumns());
- switch (f->dataType().nbRows()) {
+ Assert(f->dataType().numberOfRows() == f->dataType().numberOfColumns());
+ switch (f->dataType().numberOfRows()) {
case 1: {
auto fh = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>&>(*f);
auto gh = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>&>(*g);
@@ -350,9 +350,9 @@ applyBinaryOperation(const DiscreteFunctionT& fh, const std::shared_ptr<const ID
}
}
case ASTNodeDataType::matrix_t: {
- Assert(g->dataType().nbRows() == g->dataType().nbColumns());
+ Assert(g->dataType().numberOfRows() == g->dataType().numberOfColumns());
if constexpr (std::is_same_v<lhs_data_type, double> and std::is_same_v<language::multiply_op, BinOperatorT>) {
- switch (g->dataType().nbRows()) {
+ switch (g->dataType().numberOfRows()) {
case 1: {
auto gh = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>&>(*g);
@@ -397,8 +397,8 @@ applyBinaryOperation(const std::shared_ptr<const IDiscreteFunction>& f,
return applyBinaryOperation<BinOperatorT, Dimension>(fh, g);
}
case ASTNodeDataType::matrix_t: {
- Assert(f->dataType().nbRows() == f->dataType().nbColumns());
- switch (f->dataType().nbRows()) {
+ Assert(f->dataType().numberOfRows() == f->dataType().numberOfColumns());
+ switch (f->dataType().numberOfRows()) {
case 1: {
auto fh = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>&>(*f);
@@ -621,9 +621,9 @@ applyBinaryOperationWithLeftConstant(const DataType& a, const std::shared_ptr<co
}
}
case ASTNodeDataType::matrix_t: {
- Assert(f->dataType().nbRows() == f->dataType().nbColumns());
+ Assert(f->dataType().numberOfRows() == f->dataType().numberOfColumns());
if constexpr (is_tiny_matrix_v<DataType>) {
- switch (f->dataType().nbRows()) {
+ switch (f->dataType().numberOfRows()) {
case 1: {
if constexpr (std::is_same_v<DataType, TinyMatrix<1>>) {
auto fh = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>&>(*f);
@@ -653,7 +653,7 @@ applyBinaryOperationWithLeftConstant(const DataType& a, const std::shared_ptr<co
}
}
} else {
- switch (f->dataType().nbRows()) {
+ switch (f->dataType().numberOfRows()) {
case 1: {
auto fh = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>&>(*f);
return applyBinaryOperationWithLeftConstant<BinOperatorT>(a, fh);
@@ -747,9 +747,9 @@ applyBinaryOperationWithRightConstant(const std::shared_ptr<const IDiscreteFunct
return applyBinaryOperationWithRightConstant<BinOperatorT>(fh, a);
}
case ASTNodeDataType::matrix_t: {
- Assert(f->dataType().nbRows() == f->dataType().nbColumns());
+ Assert(f->dataType().numberOfRows() == f->dataType().numberOfColumns());
if constexpr (is_tiny_matrix_v<DataType>) {
- switch (f->dataType().nbRows()) {
+ switch (f->dataType().numberOfRows()) {
case 1: {
if constexpr (std::is_same_v<DataType, TinyMatrix<1>>) {
auto fh = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>&>(*f);
@@ -779,7 +779,7 @@ applyBinaryOperationWithRightConstant(const std::shared_ptr<const IDiscreteFunct
}
}
} else {
- switch (f->dataType().nbRows()) {
+ switch (f->dataType().numberOfRows()) {
case 1: {
auto fh = dynamic_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>&>(*f);
return applyBinaryOperationWithRightConstant<BinOperatorT>(fh, a);
diff --git a/src/language/utils/EmbeddedIDiscreteFunctionUtils.hpp b/src/language/utils/EmbeddedIDiscreteFunctionUtils.hpp
index 9321c3273d7955e77934f9f4233ae3603cf14990..981cd8be9ca166038b5aa4b913dae63488868b64 100644
--- a/src/language/utils/EmbeddedIDiscreteFunctionUtils.hpp
+++ b/src/language/utils/EmbeddedIDiscreteFunctionUtils.hpp
@@ -35,8 +35,8 @@ isSameDiscretization(const IDiscreteFunction& f, const IDiscreteFunction& g)
return f.dataType().dimension() == g.dataType().dimension();
}
case ASTNodeDataType::matrix_t: {
- return (f.dataType().nbRows() == g.dataType().nbRows()) and
- (f.dataType().nbColumns() == g.dataType().nbColumns());
+ return (f.dataType().numberOfRows() == g.dataType().numberOfRows()) and
+ (f.dataType().numberOfColumns() == g.dataType().numberOfColumns());
}
default: {
throw UnexpectedError("invalid data type " + operand_type_name(f));
diff --git a/src/language/utils/EmbedderTable.hpp b/src/language/utils/EmbedderTable.hpp
index 60441cec257aee69afc0caa72faf9a487fc736d0..68e1ddb5243e7c954745842e9e7d1e13677fcbca 100644
--- a/src/language/utils/EmbedderTable.hpp
+++ b/src/language/utils/EmbedderTable.hpp
@@ -21,7 +21,8 @@ class EmbedderTable
}
PUGS_INLINE
- const std::shared_ptr<EmbedderT>& operator[](size_t embedder_id) const
+ const std::shared_ptr<EmbedderT>&
+ operator[](size_t embedder_id) const
{
Assert(embedder_id < m_embedder_list.size());
return m_embedder_list[embedder_id];
diff --git a/src/language/utils/PugsFunctionAdapter.hpp b/src/language/utils/PugsFunctionAdapter.hpp
index b6d3a259dbb676d5f1259c058e0ea1f36fb925fd..3047652e3017fb2cd193b799a8edf454e3e4c22b 100644
--- a/src/language/utils/PugsFunctionAdapter.hpp
+++ b/src/language/utils/PugsFunctionAdapter.hpp
@@ -209,7 +209,7 @@ class PugsFunctionAdapter<OutputType(InputType...)>
};
} else {
// If this point is reached must be a 0 vector
- return [](DataVariant&&) -> OutputType { return OutputType{ZeroType{}}; };
+ return [](DataVariant &&) -> OutputType { return OutputType{ZeroType{}}; };
}
}
case ASTNodeDataType::double_t: {
@@ -236,8 +236,8 @@ class PugsFunctionAdapter<OutputType(InputType...)>
AggregateDataVariant& v = std::get<AggregateDataVariant>(result);
OutputType x;
- for (size_t i = 0, l = 0; i < x.nbRows(); ++i) {
- for (size_t j = 0; j < x.nbColumns(); ++j, ++l) {
+ for (size_t i = 0, l = 0; i < x.numberOfRows(); ++i) {
+ for (size_t j = 0; j < x.numberOfColumns(); ++j, ++l) {
std::visit(
[&](auto&& Aij) {
using Aij_T = std::decay_t<decltype(Aij)>;
@@ -288,7 +288,7 @@ class PugsFunctionAdapter<OutputType(InputType...)>
};
} else {
// If this point is reached must be a 0 matrix
- return [](DataVariant&&) -> OutputType { return OutputType{ZeroType{}}; };
+ return [](DataVariant &&) -> OutputType { return OutputType{ZeroType{}}; };
}
}
case ASTNodeDataType::double_t: {
diff --git a/src/mesh/ItemArray.hpp b/src/mesh/ItemArray.hpp
index 1a2b004c9533921ab6531b6fa73c9a955badff57..5a16861a37b5642ba10977d8c9e2d5ecd25fd194 100644
--- a/src/mesh/ItemArray.hpp
+++ b/src/mesh/ItemArray.hpp
@@ -109,7 +109,7 @@ class ItemArray
numberOfItems() const noexcept(NO_ASSERT)
{
Assert(this->isBuilt());
- return m_values.nbRows();
+ return m_values.numberOfRows();
}
PUGS_INLINE
@@ -117,7 +117,7 @@ class ItemArray
sizeOfArrays() const
{
Assert(this->isBuilt());
- return m_values.nbColumns();
+ return m_values.numberOfColumns();
}
template <typename DataType2>
@@ -131,11 +131,11 @@ class ItemArray
static_assert(((std::is_const_v<DataType2> and std::is_const_v<DataType>) or not std::is_const_v<DataType2>),
"Cannot assign ItemArray of const to ItemArray of non-const");
- Assert((table.nbRows() * table.nbColumns() == 0) or this->isBuilt(),
+ Assert((table.numberOfRows() * table.numberOfColumns() == 0) or this->isBuilt(),
"Cannot assign array of arrays to a non-built ItemArray\n");
- Assert(this->numberOfItems() == table.nbRows(), "Cannot assign a table of a different dimensions\n");
- Assert(this->sizeOfArrays() == table.nbColumns(), "Cannot assign a table of a different dimensions\n");
+ Assert(this->numberOfItems() == table.numberOfRows(), "Cannot assign a table of a different dimensions\n");
+ Assert(this->sizeOfArrays() == table.numberOfColumns(), "Cannot assign a table of a different dimensions\n");
m_values = table;
diff --git a/src/mesh/SubItemArrayPerItem.hpp b/src/mesh/SubItemArrayPerItem.hpp
index 501762ea6270b7ebf59b26c4afdea6edd552cc3e..9fd488bc5e6ebe3ea189368de8eeb67073e590ec 100644
--- a/src/mesh/SubItemArrayPerItem.hpp
+++ b/src/mesh/SubItemArrayPerItem.hpp
@@ -100,7 +100,7 @@ class SubItemArrayPerItem
{
static_assert(std::is_integral_v<ArrayIndexType>, "index must be an integral type");
Assert(this->isBuilt());
- Assert(static_cast<size_t>(i) < m_values.nbRows());
+ Assert(static_cast<size_t>(i) < m_values.numberOfRows());
return m_values[i];
}
@@ -109,7 +109,7 @@ class SubItemArrayPerItem
numberOfArrays() const noexcept(NO_ASSERT)
{
Assert(this->isBuilt());
- return m_values.nbRows();
+ return m_values.numberOfRows();
}
PUGS_INLINE
@@ -125,7 +125,7 @@ class SubItemArrayPerItem
sizeOfArrays() const noexcept(NO_ASSERT)
{
Assert(this->isBuilt());
- return m_values.nbColumns();
+ return m_values.numberOfColumns();
}
template <typename IndexType>
diff --git a/src/output/GnuplotWriter.cpp b/src/output/GnuplotWriter.cpp
index 5cb59553a54df6f1b7c173066079110058c84330..324620fd67f3330bd92958fd8df8a8ccadfc5088 100644
--- a/src/output/GnuplotWriter.cpp
+++ b/src/output/GnuplotWriter.cpp
@@ -64,8 +64,8 @@ GnuplotWriter::_writePreamble(const OutputNamedItemDataSet& output_named_item_da
fout << ' ' << i++ << ':' << name << '[' << j << ']';
}
} else if constexpr (is_tiny_matrix_v<DataType>) {
- for (size_t j = 0; j < DataType{}.nbRows(); ++j) {
- for (size_t k = 0; k < DataType{}.nbColumns(); ++k) {
+ for (size_t j = 0; j < DataType{}.numberOfRows(); ++j) {
+ for (size_t k = 0; k < DataType{}.numberOfColumns(); ++k) {
fout << ' ' << i++ << ':' << name << '(' << j << ',' << k << ')';
}
}
@@ -101,8 +101,8 @@ GnuplotWriter::_writeCellData(const CellValue<DataType>& cell_value, CellId cell
fout << ' ' << value[i];
}
} else if constexpr (is_tiny_matrix_v<std::decay_t<DataType>>) {
- for (size_t i = 0; i < value.nbRows(); ++i) {
- for (size_t j = 0; j < value.nbColumns(); ++j) {
+ for (size_t i = 0; i < value.numberOfRows(); ++i) {
+ for (size_t j = 0; j < value.numberOfColumns(); ++j) {
fout << ' ' << value(i, j);
}
}
@@ -209,8 +209,8 @@ GnuplotWriter::_writeNodeData(const NodeValue<DataType>& node_value, NodeId node
fout << ' ' << value[i];
}
} else if constexpr (is_tiny_matrix_v<std::decay_t<DataType>>) {
- for (size_t i = 0; i < value.nbRows(); ++i) {
- for (size_t j = 0; j < value.nbColumns(); ++j) {
+ for (size_t i = 0; i < value.numberOfRows(); ++i) {
+ for (size_t j = 0; j < value.numberOfColumns(); ++j) {
fout << ' ' << value(i, j);
}
}
diff --git a/src/output/GnuplotWriter1D.cpp b/src/output/GnuplotWriter1D.cpp
index 31320054e6f244e9e0771190c539eecd733d66ee..a57405d3aabc4ef23aefd6bd2b5f74f2f1681fa3 100644
--- a/src/output/GnuplotWriter1D.cpp
+++ b/src/output/GnuplotWriter1D.cpp
@@ -74,8 +74,8 @@ GnuplotWriter1D::_writePreamble(const OutputNamedItemDataSet& output_named_item_
fout << ' ' << i++ << ':' << name << '[' << j << ']';
}
} else if constexpr (is_tiny_matrix_v<DataType>) {
- for (size_t j = 0; j < DataType{}.nbRows(); ++j) {
- for (size_t k = 0; k < DataType{}.nbColumns(); ++k) {
+ for (size_t j = 0; j < DataType{}.numberOfRows(); ++j) {
+ for (size_t k = 0; k < DataType{}.numberOfColumns(); ++k) {
fout << ' ' << i++ << ':' << name << '(' << j << ',' << k << ')';
}
}
@@ -199,8 +199,8 @@ GnuplotWriter1D::_writeItemDatas(const std::shared_ptr<const MeshType>& mesh,
}
} else if constexpr (is_tiny_matrix_v<DataT>) {
size_t k = number_of_columns * index + column_number;
- for (size_t i = 0; i < DataT{}.nbRows(); ++i) {
- for (size_t j = 0; j < DataT{}.nbColumns(); ++j) {
+ for (size_t i = 0; i < DataT{}.numberOfRows(); ++i) {
+ for (size_t j = 0; j < DataT{}.numberOfColumns(); ++j) {
values[k++] = item_data[item_id](i, j);
}
}
diff --git a/src/output/WriterBase.cpp b/src/output/WriterBase.cpp
index 40d08b305094c6c33576181a9d730e579656497e..21589e732ecde0a85491438ae7dba259b8ed3060 100644
--- a/src/output/WriterBase.cpp
+++ b/src/output/WriterBase.cpp
@@ -85,8 +85,8 @@ WriterBase::_register(const std::string& name,
IDiscreteFunction::HandledItemDataType::vector) {
throw UnexpectedError("invalid data type for vector data");
} else {
- Assert(data_type.nbRows() == data_type.nbColumns(), "invalid matrix dimensions");
- switch (data_type.nbRows()) {
+ Assert(data_type.numberOfRows() == data_type.numberOfColumns(), "invalid matrix dimensions");
+ switch (data_type.numberOfRows()) {
case 1: {
_register(name,
dynamic_cast<const DiscreteFunctionType<Dimension, TinyMatrix<1, double>>&>(i_discrete_function),
diff --git a/src/scheme/DiscreteFunctionInterpoler.cpp b/src/scheme/DiscreteFunctionInterpoler.cpp
index c8eeaf94cb8d4f1037b28ebf65d377660464c093..11b19835eab3af91c0ea125359364831916cbc18 100644
--- a/src/scheme/DiscreteFunctionInterpoler.cpp
+++ b/src/scheme/DiscreteFunctionInterpoler.cpp
@@ -60,8 +60,8 @@ DiscreteFunctionInterpoler::_interpolate() const
}
}
case ASTNodeDataType::matrix_t: {
- Assert(data_type.nbColumns() == data_type.nbRows(), "undefined matrix type");
- switch (data_type.nbColumns()) {
+ Assert(data_type.numberOfColumns() == data_type.numberOfRows(), "undefined matrix type");
+ switch (data_type.numberOfColumns()) {
case 1: {
return this->_interpolate<Dimension, TinyMatrix<1>>();
}
diff --git a/src/scheme/DiscreteFunctionUtils.cpp b/src/scheme/DiscreteFunctionUtils.cpp
index 56453a746a19b9ee999fe1548d7067250e2884eb..6d2e62da4fcdd9d1c785b01894ca3d4776dfae8c 100644
--- a/src/scheme/DiscreteFunctionUtils.cpp
+++ b/src/scheme/DiscreteFunctionUtils.cpp
@@ -53,12 +53,12 @@ shallowCopy(const std::shared_ptr<const Mesh<Connectivity<Dimension>>>& mesh,
}
}
case ASTNodeDataType::matrix_t: {
- if (discrete_function->dataType().nbRows() != discrete_function->dataType().nbColumns()) {
+ if (discrete_function->dataType().numberOfRows() != discrete_function->dataType().numberOfColumns()) {
throw UnexpectedError(
- "invalid data matrix dimensions: " + std::to_string(discrete_function->dataType().nbRows()) + "x" +
- std::to_string(discrete_function->dataType().nbColumns()));
+ "invalid data matrix dimensions: " + std::to_string(discrete_function->dataType().numberOfRows()) + "x" +
+ std::to_string(discrete_function->dataType().numberOfColumns()));
}
- switch (discrete_function->dataType().nbRows()) {
+ switch (discrete_function->dataType().numberOfRows()) {
case 1: {
return shallowCopy(mesh, std::dynamic_pointer_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>>(
discrete_function));
@@ -73,8 +73,8 @@ shallowCopy(const std::shared_ptr<const Mesh<Connectivity<Dimension>>>& mesh,
}
default: {
throw UnexpectedError(
- "invalid data matrix dimensions: " + std::to_string(discrete_function->dataType().nbRows()) + "x" +
- std::to_string(discrete_function->dataType().nbColumns()));
+ "invalid data matrix dimensions: " + std::to_string(discrete_function->dataType().numberOfRows()) + "x" +
+ std::to_string(discrete_function->dataType().numberOfColumns()));
}
}
}
diff --git a/src/utils/Array.hpp b/src/utils/Array.hpp
index 5eb834b85fb8f0679636531e740863286e912c3d..9fa02e95a41a5200260fab84cda211412c4747e2 100644
--- a/src/utils/Array.hpp
+++ b/src/utils/Array.hpp
@@ -85,10 +85,15 @@ class [[nodiscard]] Array
Array& operator=(Array&&) = default;
PUGS_INLINE
- explicit Array(size_t size) : m_values("anonymous", size)
+ explicit Array(size_t size)
{
static_assert(not std::is_const<DataType>(), "Cannot allocate Array of const data: only view is "
"supported");
+ if constexpr (std::is_arithmetic_v<DataType>) {
+ m_values = Kokkos::View<DataType*>{Kokkos::view_alloc(Kokkos::WithoutInitializing, "anonymous"), size};
+ } else {
+ m_values = Kokkos::View<DataType*>{"anonymous", size};
+ }
}
PUGS_INLINE
diff --git a/src/utils/BuildInfo.cpp b/src/utils/BuildInfo.cpp
index 38d46007ef6cfcf470051e8546714db7041a7852..71d702c99ab86f7737b23a8ef42fa53845875cef 100644
--- a/src/utils/BuildInfo.cpp
+++ b/src/utils/BuildInfo.cpp
@@ -12,6 +12,10 @@
#include <petsc.h>
#endif // PUGS_HAS_PETSC
+#ifdef PUGS_HAS_SLEPC
+#include <slepc.h>
+#endif // PUGS_HAS_PETSC
+
std::string
BuildInfo::type()
{
@@ -57,3 +61,14 @@ BuildInfo::petscLibrary()
return "none";
#endif // PUGS_HAS_PETSC
}
+
+std::string
+BuildInfo::slepcLibrary()
+{
+#ifdef PUGS_HAS_SLEPC
+ return std::to_string(SLEPC_VERSION_MAJOR) + "." + std::to_string(SLEPC_VERSION_MINOR) + "." +
+ std::to_string(SLEPC_VERSION_SUBMINOR);
+#else
+ return "none";
+#endif // PUGS_HAS_SLEPC
+}
diff --git a/src/utils/BuildInfo.hpp b/src/utils/BuildInfo.hpp
index 67134a782f595ccf68dc89cd29a6f273c7ce6140..bc83cf3f1426bcab2972ebb8466f790cadf8a9ef 100644
--- a/src/utils/BuildInfo.hpp
+++ b/src/utils/BuildInfo.hpp
@@ -10,6 +10,7 @@ struct BuildInfo
static std::string kokkosDevices();
static std::string mpiLibrary();
static std::string petscLibrary();
+ static std::string slepcLibrary();
};
#endif // BUILD_INFO_HPP
diff --git a/src/utils/CMakeLists.txt b/src/utils/CMakeLists.txt
index fe9b5e6fed3ff7ddd796ab0081ce144e8b535b90..7639963d4fc73094f2e79f64e9f9c9383426ac43 100644
--- a/src/utils/CMakeLists.txt
+++ b/src/utils/CMakeLists.txt
@@ -10,9 +10,17 @@ add_library(
FPEManager.cpp
Messenger.cpp
Partitioner.cpp
+ PETScWrapper.cpp
PugsUtils.cpp
RevisionInfo.cpp
- SignalManager.cpp)
+ SignalManager.cpp
+ SLEPcWrapper.cpp)
+
+target_link_libraries(
+ PugsUtils
+ ${PETSC_LIBRARIES}
+ ${SLEPC_LIBRARIES}
+)
# --------------- get git revision info ---------------
diff --git a/src/algebra/PETScWrapper.cpp b/src/utils/PETScWrapper.cpp
similarity index 92%
rename from src/algebra/PETScWrapper.cpp
rename to src/utils/PETScWrapper.cpp
index dd11dc9772f11ad29ca4970c5b88656c623a3508..f88bbfd8a6bfd9de92d153e54c68e111f0b692b6 100644
--- a/src/algebra/PETScWrapper.cpp
+++ b/src/utils/PETScWrapper.cpp
@@ -1,4 +1,4 @@
-#include <algebra/PETScWrapper.hpp>
+#include <utils/PETScWrapper.hpp>
#include <utils/pugs_config.hpp>
diff --git a/src/algebra/PETScWrapper.hpp b/src/utils/PETScWrapper.hpp
similarity index 100%
rename from src/algebra/PETScWrapper.hpp
rename to src/utils/PETScWrapper.hpp
diff --git a/src/utils/PugsUtils.cpp b/src/utils/PugsUtils.cpp
index ad9e12e7ec4cdd01305cf8701c327894c7942304..d35fa95f43cb03810814ec6b62ed039842de1289 100644
--- a/src/utils/PugsUtils.cpp
+++ b/src/utils/PugsUtils.cpp
@@ -1,11 +1,12 @@
#include <utils/PugsUtils.hpp>
-#include <algebra/PETScWrapper.hpp>
#include <utils/BuildInfo.hpp>
#include <utils/ConsoleManager.hpp>
#include <utils/FPEManager.hpp>
#include <utils/Messenger.hpp>
+#include <utils/PETScWrapper.hpp>
#include <utils/RevisionInfo.hpp>
+#include <utils/SLEPcWrapper.hpp>
#include <utils/SignalManager.hpp>
#include <utils/pugs_build_info.hpp>
@@ -55,6 +56,7 @@ pugsBuildInfo()
os << "kokkos: " << rang::style::bold << BuildInfo::kokkosDevices() << rang::style::reset << '\n';
os << "MPI: " << rang::style::bold << BuildInfo::mpiLibrary() << rang::style::reset << '\n';
os << "PETSc: " << rang::style::bold << BuildInfo::petscLibrary() << rang::style::reset << '\n';
+ os << "SLEPc: " << rang::style::bold << BuildInfo::slepcLibrary() << rang::style::reset << '\n';
os << "-------------------------------------------------------";
return os.str();
@@ -123,6 +125,7 @@ initialize(int& argc, char* argv[])
}
PETScWrapper::initialize(argc, argv);
+ SLEPcWrapper::initialize(argc, argv);
setDefaultOMPEnvironment();
Kokkos::initialize(argc, argv);
@@ -152,6 +155,7 @@ void
finalize()
{
Kokkos::finalize();
+ SLEPcWrapper::finalize();
PETScWrapper::finalize();
parallel::Messenger::destroy();
}
diff --git a/src/utils/SLEPcWrapper.cpp b/src/utils/SLEPcWrapper.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..2c7076d54f9cf438e34e54302041cd6e76633689
--- /dev/null
+++ b/src/utils/SLEPcWrapper.cpp
@@ -0,0 +1,26 @@
+#include <utils/SLEPcWrapper.hpp>
+
+#include <utils/pugs_config.hpp>
+
+#ifdef PUGS_HAS_SLEPC
+#include <slepc.h>
+#endif // PUGS_HAS_SLEPC
+
+namespace SLEPcWrapper
+{
+void
+initialize([[maybe_unused]] int& argc, [[maybe_unused]] char* argv[])
+{
+#ifdef PUGS_HAS_SLEPC
+ SlepcInitialize(&argc, &argv, 0, 0);
+#endif // PUGS_HAS_SLEPC
+}
+
+void
+finalize()
+{
+#ifdef PUGS_HAS_SLEPC
+ SlepcFinalize();
+#endif // PUGS_HAS_SLEPC
+}
+} // namespace SLEPcWrapper
diff --git a/src/utils/SLEPcWrapper.hpp b/src/utils/SLEPcWrapper.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..8b64802339ba2c6d8c78aa3c7a86005e3fd013bf
--- /dev/null
+++ b/src/utils/SLEPcWrapper.hpp
@@ -0,0 +1,10 @@
+#ifndef SLEPC_WRAPPER_HPP
+#define SLEPC_WRAPPER_HPP
+
+namespace SLEPcWrapper
+{
+void initialize(int& argc, char* argv[]);
+void finalize();
+} // namespace SLEPcWrapper
+
+#endif // SLEPC_WRAPPER_HPP
diff --git a/src/utils/Table.hpp b/src/utils/Table.hpp
index eedc5e7d72c4c2f935f0ecafe0d24747c3b2d00e..8dae32b97b0d15a99256fd38c4f4f7c5ecf29588 100644
--- a/src/utils/Table.hpp
+++ b/src/utils/Table.hpp
@@ -21,12 +21,12 @@ class [[nodiscard]] Table
friend Table<std::add_const_t<DataType>>;
public:
- PUGS_INLINE size_t nbRows() const noexcept
+ PUGS_INLINE size_t numberOfRows() const noexcept
{
return m_values.extent(0);
}
- PUGS_INLINE size_t nbColumns() const noexcept
+ PUGS_INLINE size_t numberOfColumns() const noexcept
{
return m_values.extent(1);
}
@@ -34,13 +34,13 @@ class [[nodiscard]] Table
PUGS_INLINE
Array<DataType> operator[](index_type i) const
{
- Assert(i < this->nbRows());
+ Assert(i < this->numberOfRows());
return encapsulate(Kokkos::View<DataType*>(m_values, i, Kokkos::ALL));
}
friend PUGS_INLINE Table<std::remove_const_t<DataType>> copy(const Table<DataType>& source)
{
- Table<std::remove_const_t<DataType>> image(source.nbRows(), source.nbColumns());
+ Table<std::remove_const_t<DataType>> image(source.numberOfRows(), source.numberOfColumns());
Kokkos::deep_copy(image.m_values, source.m_values);
return image;
@@ -49,8 +49,8 @@ class [[nodiscard]] Table
friend PUGS_INLINE void copy_to(const Table<DataType>& source,
const Table<std::remove_const_t<DataType>>& destination)
{
- Assert(source.nbRows() == destination.nbRows());
- Assert(source.nbColumns() == destination.nbColumns());
+ Assert(source.numberOfRows() == destination.numberOfRows());
+ Assert(source.numberOfColumns() == destination.numberOfColumns());
Kokkos::deep_copy(destination.m_values, source.m_values);
}
@@ -66,8 +66,8 @@ class [[nodiscard]] Table
PUGS_INLINE DataType& operator()(index_type i, index_type j) const noexcept(NO_ASSERT)
{
- Assert(i < this->nbRows());
- Assert(j < this->nbColumns());
+ Assert(i < this->numberOfRows());
+ Assert(j < this->numberOfColumns());
return m_values(i, j);
}
@@ -99,10 +99,16 @@ class [[nodiscard]] Table
Table& operator=(Table&&) = default;
PUGS_INLINE
- explicit Table(size_t nb_lines, size_t nb_columns) : m_values("anonymous", nb_lines, nb_columns)
+ explicit Table(size_t nb_lines, size_t nb_columns)
{
static_assert(not std::is_const<DataType>(), "Cannot allocate Table of const data: only view is "
"supported");
+ if constexpr (std::is_arithmetic_v<DataType>) {
+ m_values =
+ Kokkos::View<DataType**>{Kokkos::view_alloc(Kokkos::WithoutInitializing, "anonymous"), nb_lines, nb_columns};
+ } else {
+ m_values = Kokkos::View<DataType**>{"anonymous", nb_lines, nb_columns};
+ }
}
PUGS_INLINE
@@ -141,10 +147,10 @@ subTable(const Table<DataType>& table,
typename Table<DataType>::index_type column_begin,
typename Table<DataType>::index_type column_size)
{
- Assert(row_begin < table.nbRows());
- Assert(row_begin + row_size <= table.nbRows());
- Assert(column_begin < table.nbColumns());
- Assert(column_begin + column_size <= table.nbColumns());
+ Assert(row_begin < table.numberOfRows());
+ Assert(row_begin + row_size <= table.numberOfRows());
+ Assert(column_begin < table.numberOfColumns());
+ Assert(column_begin + column_size <= table.numberOfColumns());
return encapsulate(Kokkos::View<DataType**>(table.m_values, std::make_pair(row_begin, row_begin + row_size),
std::make_pair(column_begin, column_begin + column_size)));
}
diff --git a/src/utils/pugs_config.hpp.in b/src/utils/pugs_config.hpp.in
index 84a6d8fa0c2d50eb623d5d1c2b52765f5610de24..0736003b9f40c4f71feea089072c19d9238e4033 100644
--- a/src/utils/pugs_config.hpp.in
+++ b/src/utils/pugs_config.hpp.in
@@ -4,6 +4,7 @@
#cmakedefine PUGS_HAS_FENV_H
#cmakedefine PUGS_HAS_MPI
#cmakedefine PUGS_HAS_PETSC
+#cmakedefine PUGS_HAS_SLEPC
#cmakedefine SYSTEM_IS_LINUX
#cmakedefine SYSTEM_IS_DARWIN
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 4bb227aac16a433471caa45bdfb8b8bc48f4f2bb..a56fb1f9d0b49108f1770f186cec453758cfd793 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -57,11 +57,13 @@ add_executable (unit_tests
test_CRSGraph.cpp
test_CRSMatrix.cpp
test_DataVariant.cpp
+ test_DenseMatrix.cpp
test_Demangle.cpp
test_DiscreteFunctionDescriptorP0.cpp
test_DiscreteFunctionDescriptorP0Vector.cpp
test_DiscreteFunctionType.cpp
test_DoWhileProcessor.cpp
+ test_EigenvalueSolver.cpp
test_EmbeddedData.cpp
test_EscapedString.cpp
test_Exceptions.cpp
diff --git a/tests/mpi_test_main.cpp b/tests/mpi_test_main.cpp
index 68736ebe3c04ff1c171be3a96dcaca4fea58c17b..51a01b0cfa99b36995b8849fc1fdb94360d971ad 100644
--- a/tests/mpi_test_main.cpp
+++ b/tests/mpi_test_main.cpp
@@ -2,13 +2,13 @@
#include <Kokkos_Core.hpp>
-#include <algebra/PETScWrapper.hpp>
#include <language/utils/OperatorRepository.hpp>
#include <mesh/DiamondDualConnectivityManager.hpp>
#include <mesh/DiamondDualMeshManager.hpp>
#include <mesh/MeshDataManager.hpp>
#include <mesh/SynchronizerManager.hpp>
#include <utils/Messenger.hpp>
+#include <utils/PETScWrapper.hpp>
#include <utils/pugs_config.hpp>
#include <MeshDataBaseForTests.hpp>
diff --git a/tests/test_ASTNodeDataType.cpp b/tests/test_ASTNodeDataType.cpp
index 9316c27dbcf09f0f114500a6c04ca3dad1c8806b..023ec9963bcd66286c9ffcfcadbc4fb3b2abdbee 100644
--- a/tests/test_ASTNodeDataType.cpp
+++ b/tests/test_ASTNodeDataType.cpp
@@ -213,8 +213,8 @@ TEST_CASE("ASTNodeDataType", "[language]")
SECTION("good node")
{
REQUIRE(getMatrixDataType(*type_node) == ASTNodeDataType::build<ASTNodeDataType::matrix_t>(3, 3));
- REQUIRE(getMatrixDataType(*type_node).nbRows() == 3);
- REQUIRE(getMatrixDataType(*type_node).nbColumns() == 3);
+ REQUIRE(getMatrixDataType(*type_node).numberOfRows() == 3);
+ REQUIRE(getMatrixDataType(*type_node).numberOfColumns() == 3);
}
SECTION("bad node type")
diff --git a/tests/test_CRSMatrix.cpp b/tests/test_CRSMatrix.cpp
index 0f748fb37971fefb53cdc881d7452607c6dd3e6e..d6bfdaba966c0c36a5245af457954a6f53a0919d 100644
--- a/tests/test_CRSMatrix.cpp
+++ b/tests/test_CRSMatrix.cpp
@@ -26,6 +26,7 @@ TEST_CASE("CRSMatrix", "[algebra]")
CRSMatrix<int, uint8_t> A{S};
REQUIRE(A.numberOfRows() == S.numberOfRows());
+ REQUIRE(A.numberOfColumns() == S.numberOfColumns());
}
SECTION("matrix vector product (simple)")
@@ -43,7 +44,7 @@ TEST_CASE("CRSMatrix", "[algebra]")
CRSMatrix<int, uint8_t> A{S};
- Vector<int> x{A.numberOfRows()};
+ Vector<int> x{A.numberOfColumns()};
x = 0;
x[0] = 1;
@@ -117,7 +118,7 @@ TEST_CASE("CRSMatrix", "[algebra]")
CRSMatrix<int, uint8_t> A{S};
- Vector<int> x{A.numberOfRows()};
+ Vector<int> x{A.numberOfColumns()};
x[0] = 1;
x[1] = 2;
x[2] = 3;
diff --git a/tests/test_DenseMatrix.cpp b/tests/test_DenseMatrix.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..6e2470fa0f325fab1e17eb0f51c2ad34a61c318b
--- /dev/null
+++ b/tests/test_DenseMatrix.cpp
@@ -0,0 +1,344 @@
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_all.hpp>
+
+#include <utils/PugsAssert.hpp>
+
+#include <algebra/DenseMatrix.hpp>
+#include <algebra/Vector.hpp>
+
+// Instantiate to ensure full coverage is performed
+template class DenseMatrix<int>;
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("DenseMatrix", "[algebra]")
+{
+ SECTION("size")
+ {
+ DenseMatrix<int> A{2, 3};
+ REQUIRE(A.numberOfRows() == 2);
+ REQUIRE(A.numberOfColumns() == 3);
+ }
+
+ SECTION("write access")
+ {
+ DenseMatrix<int> A{2, 3};
+ A(0, 0) = 0;
+ A(0, 1) = 1;
+ A(0, 2) = 2;
+ A(1, 0) = 3;
+ A(1, 1) = 4;
+ A(1, 2) = 5;
+ REQUIRE(A(0, 0) == 0);
+ REQUIRE(A(0, 1) == 1);
+ REQUIRE(A(0, 2) == 2);
+ REQUIRE(A(1, 0) == 3);
+ REQUIRE(A(1, 1) == 4);
+ REQUIRE(A(1, 2) == 5);
+ DenseMatrix<const int> const_A = A;
+ REQUIRE(const_A(0, 0) == 0);
+ REQUIRE(const_A(0, 1) == 1);
+ REQUIRE(const_A(0, 2) == 2);
+ REQUIRE(const_A(1, 0) == 3);
+ REQUIRE(const_A(1, 1) == 4);
+ REQUIRE(const_A(1, 2) == 5);
+ }
+
+ SECTION("fill")
+ {
+ DenseMatrix<int> A{2, 3};
+ A.fill(2);
+ REQUIRE(A(0, 0) == 2);
+ REQUIRE(A(0, 1) == 2);
+ REQUIRE(A(0, 2) == 2);
+ REQUIRE(A(1, 0) == 2);
+ REQUIRE(A(1, 1) == 2);
+ REQUIRE(A(1, 2) == 2);
+ }
+
+ SECTION("copy constructor (shallow)")
+ {
+ DenseMatrix<int> A{2, 3};
+ A(0, 0) = 0;
+ A(0, 1) = 1;
+ A(0, 2) = 2;
+ A(1, 0) = 3;
+ A(1, 1) = 4;
+ A(1, 2) = 5;
+
+ const DenseMatrix<int> B = A;
+ REQUIRE(B(0, 0) == 0);
+ REQUIRE(B(0, 1) == 1);
+ REQUIRE(B(0, 2) == 2);
+ REQUIRE(B(1, 0) == 3);
+ REQUIRE(B(1, 1) == 4);
+ REQUIRE(B(1, 2) == 5);
+ }
+
+ SECTION("copy (deep)")
+ {
+ DenseMatrix<int> A{2, 3};
+ A(0, 0) = 0;
+ A(0, 1) = 1;
+ A(0, 2) = 2;
+ A(1, 0) = 3;
+ A(1, 1) = 4;
+ A(1, 2) = 5;
+
+ const DenseMatrix<int> B = copy(A);
+
+ A(0, 0) = 10;
+ A(0, 1) = 11;
+ A(0, 2) = 12;
+ A(1, 0) = 13;
+ A(1, 1) = 14;
+ A(1, 2) = 15;
+
+ REQUIRE(B(0, 0) == 0);
+ REQUIRE(B(0, 1) == 1);
+ REQUIRE(B(0, 2) == 2);
+ REQUIRE(B(1, 0) == 3);
+ REQUIRE(B(1, 1) == 4);
+ REQUIRE(B(1, 2) == 5);
+
+ REQUIRE(A(0, 0) == 10);
+ REQUIRE(A(0, 1) == 11);
+ REQUIRE(A(0, 2) == 12);
+ REQUIRE(A(1, 0) == 13);
+ REQUIRE(A(1, 1) == 14);
+ REQUIRE(A(1, 2) == 15);
+ }
+
+ SECTION("self scalar multiplication")
+ {
+ DenseMatrix<int> A{2, 3};
+ A(0, 0) = 0;
+ A(0, 1) = 1;
+ A(0, 2) = 2;
+ A(1, 0) = 3;
+ A(1, 1) = 4;
+ A(1, 2) = 5;
+
+ A *= 2;
+ REQUIRE(A(0, 0) == 0);
+ REQUIRE(A(0, 1) == 2);
+ REQUIRE(A(0, 2) == 4);
+ REQUIRE(A(1, 0) == 6);
+ REQUIRE(A(1, 1) == 8);
+ REQUIRE(A(1, 2) == 10);
+ }
+
+ SECTION("left scalar multiplication")
+ {
+ DenseMatrix<int> A{2, 3};
+ A(0, 0) = 0;
+ A(0, 1) = 1;
+ A(0, 2) = 2;
+ A(1, 0) = 3;
+ A(1, 1) = 4;
+ A(1, 2) = 5;
+
+ const DenseMatrix<int> B = 2 * A;
+
+ REQUIRE(B(0, 0) == 0);
+ REQUIRE(B(0, 1) == 2);
+ REQUIRE(B(0, 2) == 4);
+ REQUIRE(B(1, 0) == 6);
+ REQUIRE(B(1, 1) == 8);
+ REQUIRE(B(1, 2) == 10);
+ }
+
+ SECTION("product matrix vector")
+ {
+ DenseMatrix<int> A{2, 3};
+ A(0, 0) = 6;
+ A(0, 1) = 1;
+ A(0, 2) = 2;
+ A(1, 0) = 3;
+ A(1, 1) = 4;
+ A(1, 2) = 5;
+
+ Vector<int> x{3};
+ x[0] = 7;
+ x[1] = 3;
+ x[2] = 4;
+
+ Vector y = A * x;
+ REQUIRE(y[0] == 53);
+ REQUIRE(y[1] == 53);
+ }
+
+ SECTION("self scalar division")
+ {
+ DenseMatrix<int> A{2, 3};
+ A(0, 0) = 0;
+ A(0, 1) = 2;
+ A(0, 2) = 4;
+ A(1, 0) = 6;
+ A(1, 1) = 8;
+ A(1, 2) = 10;
+
+ A /= 2;
+
+ REQUIRE(A(0, 0) == 0);
+ REQUIRE(A(0, 1) == 1);
+ REQUIRE(A(0, 2) == 2);
+ REQUIRE(A(1, 0) == 3);
+ REQUIRE(A(1, 1) == 4);
+ REQUIRE(A(1, 2) == 5);
+ }
+
+ SECTION("self minus")
+ {
+ DenseMatrix<int> A{2, 3};
+ A(0, 0) = 0;
+ A(0, 1) = 1;
+ A(0, 2) = 2;
+ A(1, 0) = 3;
+ A(1, 1) = 4;
+ A(1, 2) = 5;
+
+ DenseMatrix<int> B{2, 3};
+ B(0, 0) = 5;
+ B(0, 1) = 6;
+ B(0, 2) = 4;
+ B(1, 0) = 2;
+ B(1, 1) = 1;
+ B(1, 2) = 3;
+
+ A -= B;
+
+ REQUIRE(A(0, 0) == -5);
+ REQUIRE(A(0, 1) == -5);
+ REQUIRE(A(0, 2) == -2);
+ REQUIRE(A(1, 0) == 1);
+ REQUIRE(A(1, 1) == 3);
+ REQUIRE(A(1, 2) == 2);
+ }
+
+ SECTION("self sum")
+ {
+ DenseMatrix<int> A{2, 3};
+ A(0, 0) = 0;
+ A(0, 1) = 1;
+ A(0, 2) = 2;
+ A(1, 0) = 3;
+ A(1, 1) = 4;
+ A(1, 2) = 5;
+
+ DenseMatrix<int> B{2, 3};
+ B(0, 0) = 5;
+ B(0, 1) = 6;
+ B(0, 2) = 4;
+ B(1, 0) = 2;
+ B(1, 1) = 1;
+ B(1, 2) = 3;
+
+ A += B;
+
+ REQUIRE(A(0, 0) == 5);
+ REQUIRE(A(0, 1) == 7);
+ REQUIRE(A(0, 2) == 6);
+ REQUIRE(A(1, 0) == 5);
+ REQUIRE(A(1, 1) == 5);
+ REQUIRE(A(1, 2) == 8);
+ }
+
+ SECTION("sum")
+ {
+ DenseMatrix<int> A{2, 3};
+ A(0, 0) = 6;
+ A(0, 1) = 5;
+ A(0, 2) = 4;
+ A(1, 0) = 3;
+ A(1, 1) = 2;
+ A(1, 2) = 1;
+
+ DenseMatrix<int> B{2, 3};
+ B(0, 0) = 0;
+ B(0, 1) = 1;
+ B(0, 2) = 2;
+ B(1, 0) = 3;
+ B(1, 1) = 4;
+ B(1, 2) = 5;
+
+ DenseMatrix C = A + B;
+ REQUIRE(C(0, 0) == 6);
+ REQUIRE(C(0, 1) == 6);
+ REQUIRE(C(0, 2) == 6);
+ REQUIRE(C(1, 0) == 6);
+ REQUIRE(C(1, 1) == 6);
+ REQUIRE(C(1, 2) == 6);
+ }
+
+ SECTION("difference")
+ {
+ DenseMatrix<int> A{2, 3};
+ A(0, 0) = 6;
+ A(0, 1) = 5;
+ A(0, 2) = 4;
+ A(1, 0) = 3;
+ A(1, 1) = 2;
+ A(1, 2) = 1;
+
+ DenseMatrix<int> B{2, 3};
+ B(0, 0) = 0;
+ B(0, 1) = 1;
+ B(0, 2) = 2;
+ B(1, 0) = 3;
+ B(1, 1) = 4;
+ B(1, 2) = 5;
+
+ DenseMatrix C = A - B;
+ REQUIRE(C(0, 0) == 6);
+ REQUIRE(C(0, 1) == 4);
+ REQUIRE(C(0, 2) == 2);
+ REQUIRE(C(1, 0) == 0);
+ REQUIRE(C(1, 1) == -2);
+ REQUIRE(C(1, 2) == -4);
+ }
+
+ SECTION("transpose")
+ {
+ DenseMatrix<int> A{2, 3};
+ A(0, 0) = 0;
+ A(0, 1) = 1;
+ A(0, 2) = 2;
+ A(1, 0) = 3;
+ A(1, 1) = 4;
+ A(1, 2) = 5;
+
+ DenseMatrix B = transpose(A);
+ REQUIRE(B(0, 0) == 0);
+ REQUIRE(B(0, 1) == 3);
+ REQUIRE(B(1, 0) == 1);
+ REQUIRE(B(1, 1) == 4);
+ REQUIRE(B(2, 0) == 2);
+ REQUIRE(B(2, 1) == 5);
+ }
+
+ SECTION("product matrix vector")
+ {
+ DenseMatrix<int> A{2, 3};
+ A(0, 0) = 1;
+ A(0, 1) = 2;
+ A(0, 2) = 3;
+ A(1, 0) = 4;
+ A(1, 1) = 5;
+ A(1, 2) = 6;
+
+ DenseMatrix<int> B{3, 2};
+ B(0, 0) = 2;
+ B(0, 1) = 8;
+ B(1, 0) = 4;
+ B(1, 1) = 9;
+ B(2, 0) = 6;
+ B(2, 1) = 10;
+
+ DenseMatrix C = A * B;
+ REQUIRE(C(0, 0) == 28);
+ REQUIRE(C(0, 1) == 56);
+ REQUIRE(C(1, 0) == 64);
+ REQUIRE(C(1, 1) == 137);
+ }
+}
diff --git a/tests/test_EigenvalueSolver.cpp b/tests/test_EigenvalueSolver.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e779bce5600e32cf26b0fde7339693f8da69bcb1
--- /dev/null
+++ b/tests/test_EigenvalueSolver.cpp
@@ -0,0 +1,118 @@
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_all.hpp>
+
+#include <algebra/CRSMatrix.hpp>
+#include <algebra/DenseMatrix.hpp>
+#include <algebra/EigenvalueSolver.hpp>
+#include <algebra/SparseMatrixDescriptor.hpp>
+
+#include <utils/pugs_config.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("EigenvalueSolver", "[algebra]")
+{
+ SECTION("Sparse Matrices")
+ {
+ SECTION("symmetric system")
+ {
+ SparseMatrixDescriptor S{3};
+ S(0, 0) = 3;
+ S(0, 1) = 2;
+ S(0, 2) = 4;
+
+ S(1, 0) = 2;
+ S(1, 1) = 0;
+ S(1, 2) = 2;
+
+ S(2, 0) = 4;
+ S(2, 1) = 2;
+ S(2, 2) = 3;
+
+ CRSMatrix A{S};
+
+ SECTION("eigenvalues")
+ {
+ Array<double> eigenvalues;
+
+#ifdef PUGS_HAS_SLEPC
+ EigenvalueSolver{}.computeForSymmetricMatrix(A, eigenvalues);
+ REQUIRE(eigenvalues[0] == Catch::Approx(-1));
+ REQUIRE(eigenvalues[1] == Catch::Approx(-1));
+ REQUIRE(eigenvalues[2] == Catch::Approx(8));
+#else // PUGS_HAS_SLEPC
+ REQUIRE_THROWS_WITH(EigenvalueSolver{}.computeForSymmetricMatrix(A, eigenvalues),
+ "not implemented yet: SLEPc is required to solve eigenvalue problems");
+#endif // PUGS_HAS_SLEPC
+ }
+
+ SECTION("eigenvalues and eigenvectors")
+ {
+ Array<double> eigenvalue_list;
+ std::vector<Vector<double>> eigenvector_list;
+
+#ifdef PUGS_HAS_SLEPC
+ EigenvalueSolver{}.computeForSymmetricMatrix(A, eigenvalue_list, eigenvector_list);
+
+ REQUIRE(eigenvalue_list[0] == Catch::Approx(-1));
+ REQUIRE(eigenvalue_list[1] == Catch::Approx(-1));
+ REQUIRE(eigenvalue_list[2] == Catch::Approx(8));
+
+ DenseMatrix<double> P{3};
+ DenseMatrix<double> PT{3};
+ DenseMatrix<double> D{3};
+ D = zero;
+
+ for (size_t i = 0; i < 3; ++i) {
+ for (size_t j = 0; j < 3; ++j) {
+ P(i, j) = eigenvector_list[j][i];
+ PT(i, j) = eigenvector_list[i][j];
+ }
+ D(i, i) = eigenvalue_list[i];
+ }
+
+ DenseMatrix PDPT = P * D * PT;
+ for (size_t i = 0; i < 3; ++i) {
+ for (size_t j = 0; j < 3; ++j) {
+ REQUIRE(PDPT(i, j) - S(i, j) == Catch::Approx(0).margin(1E-13));
+ }
+ }
+#else // PUGS_HAS_SLEPC
+ REQUIRE_THROWS_WITH(EigenvalueSolver{}.computeForSymmetricMatrix(A, eigenvalue_list, eigenvector_list),
+ "not implemented yet: SLEPc is required to solve eigenvalue problems");
+#endif // PUGS_HAS_SLEPC
+ }
+
+ SECTION("eigenvalues and passage matrix")
+ {
+ Array<double> eigenvalue_list;
+ DenseMatrix<double> P{};
+
+#ifdef PUGS_HAS_SLEPC
+ EigenvalueSolver{}.computeForSymmetricMatrix(A, eigenvalue_list, P);
+
+ REQUIRE(eigenvalue_list[0] == Catch::Approx(-1));
+ REQUIRE(eigenvalue_list[1] == Catch::Approx(-1));
+ REQUIRE(eigenvalue_list[2] == Catch::Approx(8));
+
+ DenseMatrix<double> D{3};
+ D = zero;
+ for (size_t i = 0; i < 3; ++i) {
+ D(i, i) = eigenvalue_list[i];
+ }
+ DenseMatrix PT = transpose(P);
+
+ DenseMatrix PDPT = P * D * PT;
+ for (size_t i = 0; i < 3; ++i) {
+ for (size_t j = 0; j < 3; ++j) {
+ REQUIRE(PDPT(i, j) - S(i, j) == Catch::Approx(0).margin(1E-13));
+ }
+ }
+#else // PUGS_HAS_SLEPC
+ REQUIRE_THROWS_WITH(EigenvalueSolver{}.computeForSymmetricMatrix(A, eigenvalue_list, P),
+ "not implemented yet: SLEPc is required to solve eigenvalue problems");
+#endif // PUGS_HAS_SLEPC
+ }
+ }
+ }
+}
diff --git a/tests/test_InterpolateItemArray.cpp b/tests/test_InterpolateItemArray.cpp
index cf65afdf42c803b15900be7fc96fdad23e9c230f..5e0bbd962bd9c24165956b9e17a223f80a187055 100644
--- a/tests/test_InterpolateItemArray.cpp
+++ b/tests/test_InterpolateItemArray.cpp
@@ -241,8 +241,8 @@ let scalar_non_linear_3d: R^3 -> R, x -> 2 * exp(x[0]) * sin(x[1]) * x[2] + 3;
SECTION("interpolate on items list")
{
auto same_cell_value = [](auto interpolated, auto reference) -> bool {
- for (size_t i = 0; i < interpolated.nbRows(); ++i) {
- for (size_t j = 0; j < interpolated.nbColumns(); ++j) {
+ for (size_t i = 0; i < interpolated.numberOfRows(); ++i) {
+ for (size_t j = 0; j < interpolated.numberOfColumns(); ++j) {
if (interpolated[i][j] != reference[i][j]) {
return false;
}
diff --git a/tests/test_ItemArray.cpp b/tests/test_ItemArray.cpp
index 2a76275324c4922ed0bad182ebcfa106c3d9b462..17337200786737dcc7b98c75eb18021cfd62a82d 100644
--- a/tests/test_ItemArray.cpp
+++ b/tests/test_ItemArray.cpp
@@ -120,8 +120,8 @@ TEST_CASE("ItemArray", "[mesh]")
Table<size_t> table{cell_array.numberOfItems(), cell_array.sizeOfArrays()};
{
size_t k = 0;
- for (size_t i = 0; i < table.nbRows(); ++i) {
- for (size_t j = 0; j < table.nbColumns(); ++j) {
+ for (size_t i = 0; i < table.numberOfRows(); ++i) {
+ for (size_t j = 0; j < table.numberOfColumns(); ++j) {
table(i, j) = k++;
}
}
diff --git a/tests/test_LinearSolver.cpp b/tests/test_LinearSolver.cpp
index 5504ae390fd259eaa700a911b158bf6a6e46a625..1219c48e62d28a749383e115aef8b6e51a17e5ec 100644
--- a/tests/test_LinearSolver.cpp
+++ b/tests/test_LinearSolver.cpp
@@ -141,77 +141,52 @@ TEST_CASE("LinearSolver", "[algebra]")
#endif // PUGS_HAS_PETSC
}
- SECTION("check linear solvers")
+ SECTION("Sparse Matrices")
{
- SECTION("symmetric system")
+ SECTION("check linear solvers")
{
- SparseMatrixDescriptor S{5};
- S(0, 0) = 2;
- S(0, 1) = -1;
+ SECTION("symmetric system")
+ {
+ SparseMatrixDescriptor S{5};
+ S(0, 0) = 2;
+ S(0, 1) = -1;
- S(1, 0) = -1;
- S(1, 1) = 2;
- S(1, 2) = -1;
+ S(1, 0) = -1;
+ S(1, 1) = 2;
+ S(1, 2) = -1;
- S(2, 1) = -1;
- S(2, 2) = 2;
- S(2, 3) = -1;
+ S(2, 1) = -1;
+ S(2, 2) = 2;
+ S(2, 3) = -1;
- S(3, 2) = -1;
- S(3, 3) = 2;
- S(3, 4) = -1;
+ S(3, 2) = -1;
+ S(3, 3) = 2;
+ S(3, 4) = -1;
- S(4, 3) = -1;
- S(4, 4) = 2;
+ S(4, 3) = -1;
+ S(4, 4) = 2;
- CRSMatrix A{S};
+ CRSMatrix A{S};
- Vector<const double> x_exact = [] {
- Vector<double> y{5};
- y[0] = 1;
- y[1] = 3;
- y[2] = 2;
- y[3] = 4;
- y[4] = 5;
- return y;
- }();
+ Vector<const double> x_exact = [] {
+ Vector<double> y{5};
+ y[0] = 1;
+ y[1] = 3;
+ y[2] = 2;
+ y[3] = 4;
+ y[4] = 5;
+ return y;
+ }();
- Vector<double> b = A * x_exact;
+ Vector<double> b = A * x_exact;
- SECTION("builtin")
- {
- SECTION("CG no preconditioner")
+ SECTION("builtin")
{
- LinearSolverOptions options;
- options.library() = LSLibrary::builtin;
- options.method() = LSMethod::cg;
- options.precond() = LSPrecond::none;
-
- Vector<double> x{5};
- x = 0;
-
- 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("PETSc")
- {
-#ifdef PUGS_HAS_PETSC
-
- SECTION("CG")
- {
- LinearSolverOptions options;
- options.library() = LSLibrary::petsc;
- options.method() = LSMethod::cg;
- options.precond() = LSPrecond::none;
- options.verbose() = true;
-
SECTION("CG no preconditioner")
{
+ LinearSolverOptions options;
+ options.library() = LSLibrary::builtin;
+ options.method() = LSMethod::cg;
options.precond() = LSPrecond::none;
Vector<double> x{5};
@@ -223,24 +198,83 @@ TEST_CASE("LinearSolver", "[algebra]")
Vector error = x - x_exact;
REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
}
+ }
- SECTION("CG Diagonal")
+ SECTION("PETSc")
+ {
+#ifdef PUGS_HAS_PETSC
+
+ SECTION("CG")
{
- options.precond() = LSPrecond::diagonal;
+ LinearSolverOptions options;
+ options.library() = LSLibrary::petsc;
+ options.method() = LSMethod::cg;
+ options.precond() = LSPrecond::none;
+ options.verbose() = true;
- Vector<double> x{5};
- x = 0;
+ SECTION("CG no preconditioner")
+ {
+ options.precond() = LSPrecond::none;
- LinearSolver solver{options};
+ Vector<double> x{5};
+ x = 0;
- 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)));
+ 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("CG Diagonal")
+ {
+ options.precond() = LSPrecond::diagonal;
+
+ Vector<double> x{5};
+ x = 0;
+
+ 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("CG ICholeski")
+ {
+ options.precond() = LSPrecond::incomplete_choleski;
+
+ Vector<double> x{5};
+ x = 0;
+
+ 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("CG AMG")
+ {
+ options.precond() = LSPrecond::amg;
+
+ Vector<double> x{5};
+ x = 0;
+
+ 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("CG ICholeski")
+ SECTION("Choleski")
{
- options.precond() = LSPrecond::incomplete_choleski;
+ LinearSolverOptions options;
+ options.library() = LSLibrary::petsc;
+ options.method() = LSMethod::choleski;
+ options.precond() = LSPrecond::none;
Vector<double> x{5};
x = 0;
@@ -252,9 +286,63 @@ TEST_CASE("LinearSolver", "[algebra]")
REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
}
- SECTION("CG AMG")
+#else // PUGS_HAS_PETSC
+ SECTION("PETSc not linked")
{
- options.precond() = LSPrecond::amg;
+ LinearSolverOptions options;
+ options.library() = LSLibrary::petsc;
+ options.method() = LSMethod::cg;
+ options.precond() = LSPrecond::none;
+
+ REQUIRE_THROWS_WITH(LinearSolver{options}, "error: PETSc is not linked to pugs. Cannot use it!");
+ }
+#endif // PUGS_HAS_PETSC
+ }
+ }
+
+ SECTION("none symmetric system")
+ {
+ SparseMatrixDescriptor S{5};
+ S(0, 0) = 2;
+ S(0, 1) = -1;
+
+ S(1, 0) = -0.2;
+ S(1, 1) = 2;
+ S(1, 2) = -1;
+
+ S(2, 1) = -1;
+ S(2, 2) = 4;
+ S(2, 3) = -2;
+
+ S(3, 2) = -1;
+ S(3, 3) = 2;
+ S(3, 4) = -0.1;
+
+ S(4, 3) = 1;
+ S(4, 4) = 3;
+
+ CRSMatrix A{S};
+
+ Vector<const double> x_exact = [] {
+ Vector<double> y{5};
+ y[0] = 1;
+ y[1] = 3;
+ y[2] = 2;
+ y[3] = 4;
+ y[4] = 5;
+ return y;
+ }();
+
+ Vector<double> b = A * x_exact;
+
+ SECTION("builtin")
+ {
+ SECTION("BICGStab no preconditioner")
+ {
+ LinearSolverOptions options;
+ options.library() = LSLibrary::builtin;
+ options.method() = LSMethod::bicgstab;
+ options.precond() = LSPrecond::none;
Vector<double> x{5};
x = 0;
@@ -267,121 +355,153 @@ TEST_CASE("LinearSolver", "[algebra]")
}
}
- SECTION("Choleski")
+ SECTION("PETSc")
{
- LinearSolverOptions options;
- options.library() = LSLibrary::petsc;
- options.method() = LSMethod::choleski;
- options.precond() = LSPrecond::none;
+#ifdef PUGS_HAS_PETSC
- Vector<double> x{5};
- x = 0;
+ SECTION("BICGStab")
+ {
+ LinearSolverOptions options;
+ options.library() = LSLibrary::petsc;
+ options.method() = LSMethod::bicgstab;
+ options.precond() = LSPrecond::none;
+ options.verbose() = true;
- LinearSolver solver{options};
+ SECTION("BICGStab no preconditioner")
+ {
+ options.precond() = LSPrecond::none;
- 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)));
- }
+ Vector<double> x{5};
+ x = 0;
-#else // PUGS_HAS_PETSC
- SECTION("PETSc not linked")
- {
- LinearSolverOptions options;
- options.library() = LSLibrary::petsc;
- options.method() = LSMethod::cg;
- options.precond() = LSPrecond::none;
+ LinearSolver solver{options};
- REQUIRE_THROWS_WITH(LinearSolver{options}, "error: PETSc is not linked to pugs. Cannot use it!");
- }
-#endif // PUGS_HAS_PETSC
- }
- }
+ 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("none symmetric system")
- {
- SparseMatrixDescriptor S{5};
- S(0, 0) = 2;
- S(0, 1) = -1;
+ SECTION("BICGStab Diagonal")
+ {
+ options.precond() = LSPrecond::diagonal;
+
+ Vector<double> x{5};
+ x = 0;
- S(1, 0) = -0.2;
- S(1, 1) = 2;
- S(1, 2) = -1;
+ LinearSolver solver{options};
- S(2, 1) = -1;
- S(2, 2) = 4;
- S(2, 3) = -2;
+ 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)));
+ }
- S(3, 2) = -1;
- S(3, 3) = 2;
- S(3, 4) = -0.1;
+ SECTION("BICGStab ILU")
+ {
+ options.precond() = LSPrecond::incomplete_LU;
- S(4, 3) = 1;
- S(4, 4) = 3;
+ Vector<double> x{5};
+ x = 0;
- CRSMatrix A{S};
+ LinearSolver solver{options};
- Vector<const double> x_exact = [] {
- Vector<double> y{5};
- y[0] = 1;
- y[1] = 3;
- y[2] = 2;
- y[3] = 4;
- y[4] = 5;
- return y;
- }();
+ 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)));
+ }
+ }
- Vector<double> b = A * x_exact;
+ SECTION("BICGStab2")
+ {
+ LinearSolverOptions options;
+ options.library() = LSLibrary::petsc;
+ options.method() = LSMethod::bicgstab2;
+ options.precond() = LSPrecond::none;
- SECTION("builtin")
- {
- SECTION("BICGStab no preconditioner")
- {
- LinearSolverOptions options;
- options.library() = LSLibrary::builtin;
- options.method() = LSMethod::bicgstab;
- options.precond() = LSPrecond::none;
+ SECTION("BICGStab2 no preconditioner")
+ {
+ options.precond() = LSPrecond::none;
- Vector<double> x{5};
- x = 0;
+ Vector<double> x{5};
+ x = 0;
- LinearSolver solver{options};
+ 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)));
- }
- }
+ 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("PETSc")
- {
-#ifdef PUGS_HAS_PETSC
+ SECTION("BICGStab2 Diagonal")
+ {
+ options.precond() = LSPrecond::diagonal;
- SECTION("BICGStab")
- {
- LinearSolverOptions options;
- options.library() = LSLibrary::petsc;
- options.method() = LSMethod::bicgstab;
- options.precond() = LSPrecond::none;
- options.verbose() = true;
+ Vector<double> x{5};
+ x = 0;
- SECTION("BICGStab no preconditioner")
+ 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.method() = LSMethod::gmres;
options.precond() = LSPrecond::none;
- Vector<double> x{5};
- x = 0;
+ SECTION("GMRES no preconditioner")
+ {
+ options.precond() = LSPrecond::none;
- LinearSolver solver{options};
+ Vector<double> x{5};
+ x = 0;
- 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)));
+ 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 Diagonal")
+ {
+ options.precond() = LSPrecond::diagonal;
+
+ Vector<double> x{5};
+ x = 0;
+
+ 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 ILU")
+ {
+ options.precond() = LSPrecond::incomplete_LU;
+
+ Vector<double> x{5};
+ x = 0;
+
+ 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("BICGStab Diagonal")
+ SECTION("LU")
{
- options.precond() = LSPrecond::diagonal;
+ LinearSolverOptions options;
+ options.library() = LSLibrary::petsc;
+ options.method() = LSMethod::lu;
+ options.precond() = LSPrecond::none;
Vector<double> x{5};
x = 0;
@@ -393,9 +513,69 @@ TEST_CASE("LinearSolver", "[algebra]")
REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
}
- SECTION("BICGStab ILU")
+#else // PUGS_HAS_PETSC
+ SECTION("PETSc not linked")
+ {
+ LinearSolverOptions options;
+ options.library() = LSLibrary::petsc;
+ options.method() = LSMethod::cg;
+ options.precond() = LSPrecond::none;
+
+ REQUIRE_THROWS_WITH(LinearSolver{options}, "error: PETSc is not linked to pugs. Cannot use it!");
+ }
+#endif // PUGS_HAS_PETSC
+ }
+ }
+ }
+ }
+
+ SECTION("Dense Matrices")
+ {
+ SECTION("check linear solvers")
+ {
+ SECTION("symmetric system")
+ {
+ DenseMatrix<double> A{5};
+ A = zero;
+
+ A(0, 0) = 2;
+ A(0, 1) = -1;
+
+ A(1, 0) = -1;
+ A(1, 1) = 2;
+ A(1, 2) = -1;
+
+ A(2, 1) = -1;
+ A(2, 2) = 2;
+ A(2, 3) = -1;
+
+ A(3, 2) = -1;
+ A(3, 3) = 2;
+ A(3, 4) = -1;
+
+ A(4, 3) = -1;
+ A(4, 4) = 2;
+
+ Vector<const double> x_exact = [] {
+ Vector<double> y{5};
+ y[0] = 1;
+ y[1] = 3;
+ y[2] = 2;
+ y[3] = 4;
+ y[4] = 5;
+ return y;
+ }();
+
+ Vector<double> b = A * x_exact;
+
+ SECTION("builtin")
+ {
+ SECTION("CG no preconditioner")
{
- options.precond() = LSPrecond::incomplete_LU;
+ LinearSolverOptions options;
+ options.library() = LSLibrary::builtin;
+ options.method() = LSMethod::cg;
+ options.precond() = LSPrecond::none;
Vector<double> x{5};
x = 0;
@@ -408,15 +588,80 @@ TEST_CASE("LinearSolver", "[algebra]")
}
}
- SECTION("BICGStab2")
+ SECTION("PETSc")
{
- LinearSolverOptions options;
- options.library() = LSLibrary::petsc;
- options.method() = LSMethod::bicgstab2;
- options.precond() = LSPrecond::none;
+#ifdef PUGS_HAS_PETSC
- SECTION("BICGStab2 no preconditioner")
+ SECTION("CG")
{
+ LinearSolverOptions options;
+ options.library() = LSLibrary::petsc;
+ options.method() = LSMethod::cg;
+ options.precond() = LSPrecond::none;
+ options.verbose() = true;
+
+ SECTION("CG no preconditioner")
+ {
+ options.precond() = LSPrecond::none;
+
+ Vector<double> x{5};
+ x = 0;
+
+ 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("CG Diagonal")
+ {
+ options.precond() = LSPrecond::diagonal;
+
+ Vector<double> x{5};
+ x = 0;
+
+ 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("CG ICholeski")
+ {
+ options.precond() = LSPrecond::incomplete_choleski;
+
+ Vector<double> x{5};
+ x = 0;
+
+ 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("CG AMG")
+ {
+ options.precond() = LSPrecond::amg;
+
+ Vector<double> x{5};
+ x = 0;
+
+ 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("Choleski")
+ {
+ LinearSolverOptions options;
+ options.library() = LSLibrary::petsc;
+ options.method() = LSMethod::choleski;
options.precond() = LSPrecond::none;
Vector<double> x{5};
@@ -429,9 +674,63 @@ TEST_CASE("LinearSolver", "[algebra]")
REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
}
- SECTION("BICGStab2 Diagonal")
+#else // PUGS_HAS_PETSC
+ SECTION("PETSc not linked")
+ {
+ LinearSolverOptions options;
+ options.library() = LSLibrary::petsc;
+ options.method() = LSMethod::cg;
+ options.precond() = LSPrecond::none;
+
+ REQUIRE_THROWS_WITH(LinearSolver{options}, "error: PETSc is not linked to pugs. Cannot use it!");
+ }
+#endif // PUGS_HAS_PETSC
+ }
+ }
+
+ SECTION("none symmetric system")
+ {
+ DenseMatrix<double> A{5};
+ A = zero;
+
+ A(0, 0) = 2;
+ A(0, 1) = -1;
+
+ A(1, 0) = -0.2;
+ A(1, 1) = 2;
+ A(1, 2) = -1;
+
+ A(2, 1) = -1;
+ A(2, 2) = 4;
+ A(2, 3) = -2;
+
+ A(3, 2) = -1;
+ A(3, 3) = 2;
+ A(3, 4) = -0.1;
+
+ A(4, 3) = 1;
+ A(4, 4) = 3;
+
+ Vector<const double> x_exact = [] {
+ Vector<double> y{5};
+ y[0] = 1;
+ y[1] = 3;
+ y[2] = 2;
+ y[3] = 4;
+ y[4] = 5;
+ return y;
+ }();
+
+ Vector<double> b = A * x_exact;
+
+ SECTION("builtin")
+ {
+ SECTION("BICGStab no preconditioner")
{
- options.precond() = LSPrecond::diagonal;
+ LinearSolverOptions options;
+ options.library() = LSLibrary::builtin;
+ options.method() = LSMethod::bicgstab;
+ options.precond() = LSPrecond::none;
Vector<double> x{5};
x = 0;
@@ -444,44 +743,153 @@ TEST_CASE("LinearSolver", "[algebra]")
}
}
- SECTION("GMRES")
+ SECTION("PETSc")
{
- LinearSolverOptions options;
- options.library() = LSLibrary::petsc;
- options.method() = LSMethod::gmres;
- options.precond() = LSPrecond::none;
+#ifdef PUGS_HAS_PETSC
- SECTION("GMRES no preconditioner")
+ SECTION("BICGStab")
{
+ LinearSolverOptions options;
+ options.library() = LSLibrary::petsc;
+ options.method() = LSMethod::bicgstab;
options.precond() = LSPrecond::none;
+ options.verbose() = true;
- Vector<double> x{5};
- x = 0;
+ SECTION("BICGStab no preconditioner")
+ {
+ options.precond() = LSPrecond::none;
- LinearSolver solver{options};
+ Vector<double> x{5};
+ x = 0;
- 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)));
+ 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("BICGStab Diagonal")
+ {
+ options.precond() = LSPrecond::diagonal;
+
+ Vector<double> x{5};
+ x = 0;
+
+ 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("BICGStab ILU")
+ {
+ options.precond() = LSPrecond::incomplete_LU;
+
+ Vector<double> x{5};
+ x = 0;
+
+ 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 Diagonal")
+ SECTION("BICGStab2")
{
- options.precond() = LSPrecond::diagonal;
+ LinearSolverOptions options;
+ options.library() = LSLibrary::petsc;
+ options.method() = LSMethod::bicgstab2;
+ options.precond() = LSPrecond::none;
- Vector<double> x{5};
- x = 0;
+ SECTION("BICGStab2 no preconditioner")
+ {
+ options.precond() = LSPrecond::none;
- LinearSolver solver{options};
+ Vector<double> x{5};
+ x = 0;
- 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)));
+ 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 = 0;
+
+ 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.method() = LSMethod::gmres;
+ options.precond() = LSPrecond::none;
+
+ SECTION("GMRES no preconditioner")
+ {
+ options.precond() = LSPrecond::none;
+
+ Vector<double> x{5};
+ x = 0;
+
+ 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 Diagonal")
+ {
+ options.precond() = LSPrecond::diagonal;
+
+ Vector<double> x{5};
+ x = 0;
+
+ 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 ILU")
+ {
+ options.precond() = LSPrecond::incomplete_LU;
+
+ Vector<double> x{5};
+ x = 0;
+
+ 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 ILU")
+ SECTION("LU")
{
- options.precond() = LSPrecond::incomplete_LU;
+ LinearSolverOptions options;
+ options.library() = LSLibrary::petsc;
+ options.method() = LSMethod::lu;
+ options.precond() = LSPrecond::none;
Vector<double> x{5};
x = 0;
@@ -492,36 +900,19 @@ TEST_CASE("LinearSolver", "[algebra]")
Vector error = x - x_exact;
REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
}
- }
-
- SECTION("LU")
- {
- LinearSolverOptions options;
- options.library() = LSLibrary::petsc;
- options.method() = LSMethod::lu;
- options.precond() = LSPrecond::none;
-
- Vector<double> x{5};
- x = 0;
-
- 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)));
- }
#else // PUGS_HAS_PETSC
- SECTION("PETSc not linked")
- {
- LinearSolverOptions options;
- options.library() = LSLibrary::petsc;
- options.method() = LSMethod::cg;
- options.precond() = LSPrecond::none;
+ SECTION("PETSc not linked")
+ {
+ LinearSolverOptions options;
+ options.library() = LSLibrary::petsc;
+ options.method() = LSMethod::cg;
+ options.precond() = LSPrecond::none;
- REQUIRE_THROWS_WITH(LinearSolver{options}, "error: PETSc is not linked to pugs. Cannot use it!");
- }
+ REQUIRE_THROWS_WITH(LinearSolver{options}, "error: PETSc is not linked to pugs. Cannot use it!");
+ }
#endif // PUGS_HAS_PETSC
+ }
}
}
}
diff --git a/tests/test_PugsUtils.cpp b/tests/test_PugsUtils.cpp
index f32e8014a76561a07aa3d70756641e0a9231238d..4632ae825bf9089834f90ba082d7481746885d28 100644
--- a/tests/test_PugsUtils.cpp
+++ b/tests/test_PugsUtils.cpp
@@ -50,6 +50,7 @@ TEST_CASE("PugsUtils", "[utils]")
os << "kokkos: " << rang::style::bold << BuildInfo::kokkosDevices() << rang::style::reset << '\n';
os << "MPI: " << rang::style::bold << BuildInfo::mpiLibrary() << rang::style::reset << '\n';
os << "PETSc: " << rang::style::bold << BuildInfo::petscLibrary() << rang::style::reset << '\n';
+ os << "SLEPc: " << rang::style::bold << BuildInfo::slepcLibrary() << rang::style::reset << '\n';
os << "-------------------------------------------------------";
return os.str();
diff --git a/tests/test_SparseMatrixDescriptor.cpp b/tests/test_SparseMatrixDescriptor.cpp
index 8bd0c37c5cc9e294462fc92a759980dbc7a4a336..77bc890f7c4a53ea170de08bebb83de9b7e0fbed 100644
--- a/tests/test_SparseMatrixDescriptor.cpp
+++ b/tests/test_SparseMatrixDescriptor.cpp
@@ -51,10 +51,11 @@ TEST_CASE("SparseMatrixDescriptor", "[algebra]")
}
}
- SECTION("number of columns")
+ SECTION("matrix size")
{
SparseMatrixDescriptor<int, uint8_t> S{5};
REQUIRE(S.numberOfRows() == 5);
+ REQUIRE(S.numberOfColumns() == 5);
}
SECTION("location operators")
diff --git a/tests/test_SubItemArrayPerItem.cpp b/tests/test_SubItemArrayPerItem.cpp
index b03353f0ebbff8fe3d6c5634020d1ea739668f57..d1cb315046b5450e975b7ced5e795103bb6453ad 100644
--- a/tests/test_SubItemArrayPerItem.cpp
+++ b/tests/test_SubItemArrayPerItem.cpp
@@ -80,7 +80,7 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
for (CellId cell_id = 0; cell_id < connectivity.numberOfCells(); ++cell_id) {
is_correct &=
(cell_to_node_matrix[cell_id].size() == node_array_per_cell.numberOfSubArrays(cell_id)) and
- (node_array_per_cell.itemTable(cell_id).nbRows() == node_array_per_cell.numberOfSubArrays(cell_id));
+ (node_array_per_cell.itemTable(cell_id).numberOfRows() == node_array_per_cell.numberOfSubArrays(cell_id));
}
REQUIRE(is_correct);
}
@@ -89,8 +89,8 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
{
bool is_correct = true;
for (CellId cell_id = 0; cell_id < connectivity.numberOfCells(); ++cell_id) {
- is_correct &=
- (const_node_array_per_cell.itemTable(cell_id).nbRows() == node_array_per_cell.numberOfSubArrays(cell_id));
+ is_correct &= (const_node_array_per_cell.itemTable(cell_id).numberOfRows() ==
+ node_array_per_cell.numberOfSubArrays(cell_id));
}
REQUIRE(is_correct);
}
@@ -576,7 +576,7 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
size_t i = 0;
for (CellId cell_id = 0; cell_id < connectivity.numberOfCells(); ++cell_id) {
const auto& cell_table = const_edge_arrays_per_cell.itemTable(cell_id);
- for (size_t i_edge = 0; i_edge < cell_table.nbRows(); ++i_edge) {
+ for (size_t i_edge = 0; i_edge < cell_table.numberOfRows(); ++i_edge) {
const auto& array = cell_table[i_edge];
for (size_t l = 0; l < array.size(); ++l, ++i) {
is_same &= (array[l] == i * i + 1);
@@ -640,7 +640,7 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
size_t i = 0;
for (FaceId face_id = 0; face_id < connectivity.numberOfFaces(); ++face_id) {
const auto& face_table = const_cell_arrays_per_face.itemTable(face_id);
- for (size_t i_cell = 0; i_cell < face_table.nbRows(); ++i_cell) {
+ for (size_t i_cell = 0; i_cell < face_table.numberOfRows(); ++i_cell) {
const auto& array = face_table[i_cell];
for (size_t l = 0; l < array.size(); ++l, ++i) {
is_same &= (array[l] == 3 * i + 1);
@@ -703,7 +703,7 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
size_t i = 0;
for (NodeId node_id = 0; node_id < connectivity.numberOfNodes(); ++node_id) {
const auto& node_table = const_face_arrays_per_node.itemTable(node_id);
- for (size_t i_face = 0; i_face < node_table.nbRows(); ++i_face) {
+ for (size_t i_face = 0; i_face < node_table.numberOfRows(); ++i_face) {
const auto& array = node_table[i_face];
for (size_t l = 0; l < array.size(); ++l, ++i) {
is_same &= (array[l] == 3 + i * i);
@@ -900,9 +900,9 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
if (connectivity.numberOfFaces() > 0) {
FaceId face_id = 0;
const auto& face_table = cell_array_per_face.itemTable(face_id);
- REQUIRE_THROWS_AS(cell_array_per_face(face_id, face_table.nbRows()), AssertError);
- REQUIRE_THROWS_AS(face_table[face_table.nbRows()], AssertError);
- REQUIRE_THROWS_AS(cell_array_per_face.itemTable(face_id)[face_table.nbRows()], AssertError);
+ REQUIRE_THROWS_AS(cell_array_per_face(face_id, face_table.numberOfRows()), AssertError);
+ REQUIRE_THROWS_AS(face_table[face_table.numberOfRows()], AssertError);
+ REQUIRE_THROWS_AS(cell_array_per_face.itemTable(face_id)[face_table.numberOfRows()], AssertError);
REQUIRE_THROWS_AS(cell_array_per_face.itemTable(face_id)[0][cell_array_per_face.sizeOfArrays()], AssertError);
REQUIRE_THROWS_AS(cell_array_per_face.itemTable(face_id)[0][cell_array_per_face.sizeOfArrays()] = 2,
AssertError);
@@ -916,9 +916,9 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
if (connectivity.numberOfNodes() > 0) {
NodeId node_id = 0;
const auto& node_table = face_array_per_node.itemTable(node_id);
- REQUIRE_THROWS_AS(face_array_per_node(node_id, node_table.nbRows()), AssertError);
- REQUIRE_THROWS_AS(node_table[node_table.nbRows()], AssertError);
- REQUIRE_THROWS_AS(face_array_per_node.itemTable(node_id)[node_table.nbRows()], AssertError);
+ REQUIRE_THROWS_AS(face_array_per_node(node_id, node_table.numberOfRows()), AssertError);
+ REQUIRE_THROWS_AS(node_table[node_table.numberOfRows()], AssertError);
+ REQUIRE_THROWS_AS(face_array_per_node.itemTable(node_id)[node_table.numberOfRows()], AssertError);
REQUIRE_THROWS_AS(face_array_per_node.itemTable(node_id)[0][face_array_per_node.sizeOfArrays()], AssertError);
REQUIRE_THROWS_AS(face_array_per_node.itemTable(node_id)[0][face_array_per_node.sizeOfArrays()] = 2,
AssertError);
@@ -932,9 +932,9 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
if (connectivity.numberOfCells() > 0) {
CellId cell_id = 0;
const auto& cell_table = edge_array_per_cell.itemTable(cell_id);
- REQUIRE_THROWS_AS(edge_array_per_cell(cell_id, cell_table.nbRows()), AssertError);
- REQUIRE_THROWS_AS(cell_table[cell_table.nbRows()], AssertError);
- REQUIRE_THROWS_AS(edge_array_per_cell.itemTable(cell_id)[cell_table.nbRows()], AssertError);
+ REQUIRE_THROWS_AS(edge_array_per_cell(cell_id, cell_table.numberOfRows()), AssertError);
+ REQUIRE_THROWS_AS(cell_table[cell_table.numberOfRows()], AssertError);
+ REQUIRE_THROWS_AS(edge_array_per_cell.itemTable(cell_id)[cell_table.numberOfRows()], AssertError);
REQUIRE_THROWS_AS(edge_array_per_cell.itemTable(cell_id)[0][edge_array_per_cell.sizeOfArrays()], AssertError);
REQUIRE_THROWS_AS(edge_array_per_cell.itemTable(cell_id)[0][edge_array_per_cell.sizeOfArrays()] == 2,
AssertError);
@@ -948,9 +948,9 @@ TEST_CASE("SubItemArrayPerItem", "[mesh]")
if (connectivity.numberOfEdges() > 0) {
EdgeId edge_id = 0;
const auto& edge_table = node_array_per_edge.itemTable(edge_id);
- REQUIRE_THROWS_AS(node_array_per_edge(edge_id, edge_table.nbRows()), AssertError);
- REQUIRE_THROWS_AS(edge_table[edge_table.nbRows()], AssertError);
- REQUIRE_THROWS_AS(node_array_per_edge.itemTable(edge_id)[edge_table.nbRows()], AssertError);
+ REQUIRE_THROWS_AS(node_array_per_edge(edge_id, edge_table.numberOfRows()), AssertError);
+ REQUIRE_THROWS_AS(edge_table[edge_table.numberOfRows()], AssertError);
+ REQUIRE_THROWS_AS(node_array_per_edge.itemTable(edge_id)[edge_table.numberOfRows()], AssertError);
REQUIRE_THROWS_AS(node_array_per_edge.itemTable(edge_id)[0][node_array_per_edge.sizeOfArrays()], AssertError);
REQUIRE_THROWS_AS(node_array_per_edge.itemTable(edge_id)[0][node_array_per_edge.sizeOfArrays()] = 2,
AssertError);
diff --git a/tests/test_Table.cpp b/tests/test_Table.cpp
index 3521476c375dba6a31e57dd2f55ad09430026d8a..5e738867a592c3f06a9fe8ac7bd958dbcf390e4e 100644
--- a/tests/test_Table.cpp
+++ b/tests/test_Table.cpp
@@ -13,15 +13,15 @@ template class Table<int>;
TEST_CASE("Table", "[utils]")
{
Table<int> a(4, 3);
- REQUIRE(a.nbRows() == 4);
- REQUIRE(a.nbColumns() == 3);
+ REQUIRE(a.numberOfRows() == 4);
+ REQUIRE(a.numberOfColumns() == 3);
REQUIRE(a[0].size() == 3);
REQUIRE(a[1].size() == 3);
REQUIRE(a[2].size() == 3);
REQUIRE(a[3].size() == 3);
- for (size_t i = 0; i < a.nbRows(); ++i) {
- for (size_t j = 0; j < a.nbColumns(); ++j) {
+ for (size_t i = 0; i < a.numberOfRows(); ++i) {
+ for (size_t j = 0; j < a.numberOfColumns(); ++j) {
a(i, j) = 2 * i + j;
}
}
@@ -138,7 +138,7 @@ TEST_CASE("Table", "[utils]")
(c(0, 1) == 2) and (c(1, 1) == 2) and (c(2, 1) == 2) and (c(3, 1) == 2) and //
(c(0, 2) == 2) and (c(1, 2) == 2) and (c(2, 2) == 2) and (c(3, 2) == 2)));
- Table<int> d{a.nbRows(), a.nbColumns()};
+ Table<int> d{a.numberOfRows(), a.numberOfColumns()};
copy_to(a, d);
REQUIRE(((a(0, 0) == 0) and (a(1, 0) == 2) and (a(2, 0) == 4) and (a(3, 0) == 6) and //
@@ -168,10 +168,10 @@ TEST_CASE("Table", "[utils]")
Table table = encapsulate(kokkos_view);
- REQUIRE(table.nbRows() == kokkos_view.extent(0));
- REQUIRE(table.nbColumns() == kokkos_view.extent(1));
- for (size_t i = 0; i < table.nbRows(); ++i) {
- for (size_t j = 0; j < table.nbColumns(); ++j) {
+ REQUIRE(table.numberOfRows() == kokkos_view.extent(0));
+ REQUIRE(table.numberOfColumns() == kokkos_view.extent(1));
+ for (size_t i = 0; i < table.numberOfRows(); ++i) {
+ for (size_t j = 0; j < table.numberOfColumns(); ++j) {
REQUIRE(&table(i, j) == &kokkos_view(i, j));
}
}
@@ -189,13 +189,13 @@ TEST_CASE("Table", "[utils]")
{
SECTION("wrong row number")
{
- Table<int> b{2 * a.nbRows(), a.nbColumns()};
+ Table<int> b{2 * a.numberOfRows(), a.numberOfColumns()};
REQUIRE_THROWS_AS(copy_to(a, b), AssertError);
}
SECTION("wrong column number")
{
- Table<int> c{a.nbRows(), 2 * a.nbColumns()};
+ Table<int> c{a.numberOfRows(), 2 * a.numberOfColumns()};
REQUIRE_THROWS_AS(copy_to(a, c), AssertError);
}
}
diff --git a/tests/test_TinyMatrix.cpp b/tests/test_TinyMatrix.cpp
index 158512dbf67e80751c362c67385b454a1bd1e5c9..464e78448192f3dacc248d6c7fde8c63338449bb 100644
--- a/tests/test_TinyMatrix.cpp
+++ b/tests/test_TinyMatrix.cpp
@@ -206,16 +206,16 @@ TEST_CASE("TinyMatrix", "[algebra]")
SECTION("checking for sizes")
{
- REQUIRE(TinyMatrix<1>{}.nbRows() == 1);
- REQUIRE(TinyMatrix<1>{}.nbColumns() == 1);
+ REQUIRE(TinyMatrix<1>{}.numberOfRows() == 1);
+ REQUIRE(TinyMatrix<1>{}.numberOfColumns() == 1);
REQUIRE(TinyMatrix<1>{}.dimension() == 1);
- REQUIRE(TinyMatrix<2>{}.nbRows() == 2);
- REQUIRE(TinyMatrix<2>{}.nbColumns() == 2);
+ REQUIRE(TinyMatrix<2>{}.numberOfRows() == 2);
+ REQUIRE(TinyMatrix<2>{}.numberOfColumns() == 2);
REQUIRE(TinyMatrix<2>{}.dimension() == 4);
- REQUIRE(TinyMatrix<3>{}.nbRows() == 3);
- REQUIRE(TinyMatrix<3>{}.nbColumns() == 3);
+ REQUIRE(TinyMatrix<3>{}.numberOfRows() == 3);
+ REQUIRE(TinyMatrix<3>{}.numberOfColumns() == 3);
REQUIRE(TinyMatrix<3>{}.dimension() == 9);
}
diff --git a/tests/test_main.cpp b/tests/test_main.cpp
index ee1d0769b152c5b483214e0edc7d9313775cf416..b5316614d10c21cbe13bb5540cc0a13d7e4b4b24 100644
--- a/tests/test_main.cpp
+++ b/tests/test_main.cpp
@@ -2,13 +2,14 @@
#include <Kokkos_Core.hpp>
-#include <algebra/PETScWrapper.hpp>
#include <language/utils/OperatorRepository.hpp>
#include <mesh/DiamondDualConnectivityManager.hpp>
#include <mesh/DiamondDualMeshManager.hpp>
#include <mesh/MeshDataManager.hpp>
#include <mesh/SynchronizerManager.hpp>
#include <utils/Messenger.hpp>
+#include <utils/PETScWrapper.hpp>
+#include <utils/SLEPcWrapper.hpp>
#include <MeshDataBaseForTests.hpp>
@@ -19,6 +20,7 @@ main(int argc, char* argv[])
Kokkos::initialize({4, -1, -1, true});
PETScWrapper::initialize(argc, argv);
+ SLEPcWrapper::initialize(argc, argv);
Catch::Session session;
int result = session.applyCommandLine(argc, argv);
@@ -53,6 +55,7 @@ main(int argc, char* argv[])
}
}
+ SLEPcWrapper::finalize();
PETScWrapper::finalize();
Kokkos::finalize();