From 344591351adf4b47decc426f7c32c2727298be7c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Thu, 1 Jul 2021 16:36:32 +0200
Subject: [PATCH] Replace SparseMatrixDescriptor by CRSMatrixDescriptor
The new class is meant to be much faster
- it requires an estimation of the number values per row
- if the number of stored values per line exceeds the reserved one, an
a special storage zone (per row) is used.
It is important to give a precise estimation of the number of values
per rows:
- if the provided expected number of values is too small with regard
to the real number of values stored for a given row, then building
performances can be noticeably affected
- if the provided expected number of values per row corresponds
exactly to the effective number of stored values, then the
construction of the CRSMatrix costs nothing.
In the later case, one benefits from a significant CPU and memory
costs reduction.
---
src/algebra/CRSMatrix.hpp | 84 +++++---
src/algebra/CRSMatrixDescriptor.hpp | 265 +++++++++++++++++++++++++
src/algebra/LinearSolver.cpp | 4 +-
src/algebra/LinearSolver.hpp | 2 +-
src/algebra/PETScUtils.cpp | 38 ++--
src/algebra/PETScUtils.hpp | 7 +-
src/algebra/SparseMatrixDescriptor.hpp | 170 ----------------
tests/CMakeLists.txt | 2 +-
tests/test_BiCGStab.cpp | 17 +-
tests/test_CG.cpp | 23 ++-
tests/test_CRSMatrix.cpp | 118 ++++++++---
tests/test_CRSMatrixDescriptor.cpp | 166 ++++++++++++++++
tests/test_EigenvalueSolver.cpp | 9 +-
tests/test_LinearSolver.cpp | 18 +-
tests/test_SparseMatrixDescriptor.cpp | 201 -------------------
15 files changed, 640 insertions(+), 484 deletions(-)
create mode 100644 src/algebra/CRSMatrixDescriptor.hpp
delete mode 100644 src/algebra/SparseMatrixDescriptor.hpp
create mode 100644 tests/test_CRSMatrixDescriptor.cpp
delete mode 100644 tests/test_SparseMatrixDescriptor.cpp
diff --git a/src/algebra/CRSMatrix.hpp b/src/algebra/CRSMatrix.hpp
index edb85d5e8..bf650afe2 100644
--- a/src/algebra/CRSMatrix.hpp
+++ b/src/algebra/CRSMatrix.hpp
@@ -1,28 +1,31 @@
#ifndef CRS_MATRIX_HPP
#define CRS_MATRIX_HPP
-#include <Kokkos_StaticCrsGraph.hpp>
-
#include <utils/Array.hpp>
#include <utils/PugsAssert.hpp>
-#include <algebra/SparseMatrixDescriptor.hpp>
#include <algebra/Vector.hpp>
#include <iostream>
-template <typename DataType = double, typename IndexType = size_t>
+template <typename DataType, typename IndexType>
+class CRSMatrixDescriptor;
+
+template <typename DataType = double, typename IndexType = int>
class CRSMatrix
{
+ public:
using MutableDataType = std::remove_const_t<DataType>;
+ using index_type = IndexType;
+ using data_type = DataType;
private:
- using HostMatrix = Kokkos::StaticCrsGraph<const IndexType, Kokkos::HostSpace>;
+ const IndexType m_nb_rows;
+ const IndexType m_nb_columns;
- HostMatrix m_host_matrix;
+ Array<const IndexType> m_row_map;
Array<const DataType> m_values;
- size_t m_nb_rows;
- size_t m_nb_columns;
+ Array<const IndexType> m_column_indices;
public:
PUGS_INLINE
@@ -33,14 +36,14 @@ class CRSMatrix
}
PUGS_INLINE
- size_t
+ IndexType
numberOfRows() const
{
return m_nb_rows;
}
PUGS_INLINE
- size_t
+ IndexType
numberOfColumns() const
{
return m_nb_columns;
@@ -53,15 +56,15 @@ class CRSMatrix
}
auto
- rowIndices() const
+ rowMap() const
{
- return encapsulate(m_host_matrix.row_map);
+ return m_row_map;
}
auto
- row(size_t i) const
+ columnIndices() const
{
- return m_host_matrix.rowConst(i);
+ return m_column_indices;
}
template <typename DataType2>
@@ -71,18 +74,17 @@ class CRSMatrix
static_assert(std::is_same<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>(),
"Cannot multiply matrix and vector of different type");
- Assert(x.size() == m_host_matrix.numRows());
+ Assert(x.size() - m_nb_columns == 0);
- Vector<MutableDataType> Ax{m_host_matrix.numRows()};
- auto host_row_map = m_host_matrix.row_map;
+ Vector<MutableDataType> Ax(m_nb_rows);
parallel_for(
- m_host_matrix.numRows(), PUGS_LAMBDA(const IndexType& i_row) {
- const auto row_begin = host_row_map(i_row);
- const auto row_end = host_row_map(i_row + 1);
+ m_nb_rows, PUGS_LAMBDA(const IndexType& i_row) {
+ const auto row_begin = m_row_map[i_row];
+ const auto row_end = m_row_map[i_row + 1];
MutableDataType sum{0};
for (IndexType j = row_begin; j < row_end; ++j) {
- sum += m_values[j] * x[m_host_matrix.entries(j)];
+ sum += m_values[j] * x[m_column_indices[j]];
}
Ax[i_row] = sum;
});
@@ -90,21 +92,37 @@ class CRSMatrix
return Ax;
}
- CRSMatrix(const SparseMatrixDescriptor<DataType, IndexType>& M)
- : m_nb_rows{M.numberOfRows()}, m_nb_columns{M.numberOfColumns()}
+ friend PUGS_INLINE std::ostream&
+ operator<<(std::ostream& os, const CRSMatrix& A)
{
- {
- auto host_matrix =
- Kokkos::create_staticcrsgraph<Kokkos::StaticCrsGraph<IndexType, Kokkos::HostSpace>>("connectivity_matrix",
- M.graphVector());
-
- // This is a bit crappy but it is the price to pay to avoid
- m_host_matrix.entries = host_matrix.entries;
- m_host_matrix.row_map = host_matrix.row_map;
- m_host_matrix.row_block_offsets = host_matrix.row_block_offsets;
+ for (IndexType i = 0; i < A.m_nb_rows; ++i) {
+ const auto row_begin = A.m_row_map[i];
+ const auto row_end = A.m_row_map[i + 1];
+ os << i << "|";
+ for (IndexType j = row_begin; j < row_end; ++j) {
+ os << ' ' << A.m_column_indices[j] << ':' << A.m_values[j];
+ }
+ os << '\n';
}
- m_values = M.valueArray();
+ return os;
}
+
+ private:
+ friend class CRSMatrixDescriptor<DataType, IndexType>;
+
+ CRSMatrix(const IndexType nb_rows,
+ const IndexType nb_columns,
+ const Array<const IndexType>& row_map,
+ const Array<const DataType>& values,
+ const Array<const IndexType>& column_indices)
+ : m_nb_rows{nb_rows},
+ m_nb_columns{nb_columns},
+ m_row_map{row_map},
+ m_values{values},
+ m_column_indices{column_indices}
+ {}
+
+ public:
~CRSMatrix() = default;
};
diff --git a/src/algebra/CRSMatrixDescriptor.hpp b/src/algebra/CRSMatrixDescriptor.hpp
new file mode 100644
index 000000000..cbf82a904
--- /dev/null
+++ b/src/algebra/CRSMatrixDescriptor.hpp
@@ -0,0 +1,265 @@
+#ifndef CRS_MATRIX_DESCRIPTOR_HPP
+#define CRS_MATRIX_DESCRIPTOR_HPP
+
+#include <algebra/CRSMatrix.hpp>
+#include <utils/Array.hpp>
+#include <utils/Exceptions.hpp>
+
+#include <map>
+
+template <typename DataType = double, typename IndexType = int>
+class CRSMatrixDescriptor
+{
+ public:
+ using index_type = IndexType;
+ using data_type = DataType;
+
+ static_assert(std::is_integral_v<IndexType>);
+
+ private:
+ const IndexType m_nb_rows;
+ const IndexType m_nb_columns;
+
+ Array<const IndexType> m_row_map;
+ Array<data_type> m_values;
+ Array<IndexType> m_column_indices;
+
+ Array<IndexType> m_nb_values_per_row;
+ Array<std::map<IndexType, DataType>> m_overflow_per_row;
+
+ Array<const IndexType>
+ _computeRowMap(const Array<IndexType>& non_zeros) const
+ {
+ Assert(non_zeros.size() - m_nb_rows == 0);
+
+ Array<IndexType> row_map(m_nb_rows + 1);
+ row_map[0] = 0;
+ for (IndexType i = 0; i < m_nb_rows; ++i) {
+ row_map[i + 1] = row_map[i] + non_zeros[i];
+ }
+
+ return row_map;
+ }
+
+ public:
+ IndexType
+ numberOfRows() const
+ {
+ return m_nb_rows;
+ }
+
+ IndexType
+ numberOfColumns() const
+ {
+ return m_nb_columns;
+ }
+
+ bool PUGS_INLINE
+ hasOverflow() const
+ {
+ for (IndexType i = 0; i < m_nb_rows; ++i) {
+ if (m_overflow_per_row[i].size() > 0) {
+ return true;
+ }
+ }
+ return false;
+ }
+
+ bool PUGS_INLINE
+ isFilled() const
+ {
+ for (IndexType i = 0; i < m_nb_rows; ++i) {
+ if (m_nb_values_per_row[i] < m_row_map[i + 1] - m_row_map[i]) {
+ return false;
+ }
+ }
+ return true;
+ }
+
+ friend PUGS_INLINE std::ostream&
+ operator<<(std::ostream& os, const CRSMatrixDescriptor& A)
+ {
+ for (IndexType i = 0; i < A.m_nb_rows; ++i) {
+ const auto& overflow_row = A.m_overflow_per_row[i];
+ os << i << "|";
+
+ if (A.m_nb_values_per_row[i] > 0) {
+ IndexType j = 0;
+ auto i_overflow = overflow_row.begin();
+ while (j < A.m_nb_values_per_row[i] or i_overflow != overflow_row.end()) {
+ const IndexType j_index = [&] {
+ if (j < A.m_nb_values_per_row[i]) {
+ return A.m_column_indices[A.m_row_map[i] + j];
+ } else {
+ return std::numeric_limits<IndexType>::max();
+ }
+ }();
+
+ const IndexType overflow_index = [&] {
+ if (i_overflow != overflow_row.end()) {
+ return i_overflow->first;
+ } else {
+ return std::numeric_limits<IndexType>::max();
+ }
+ }();
+
+ if (j_index < overflow_index) {
+ os << ' ' << A.m_column_indices[A.m_row_map[i] + j] << ':' << A.m_values[A.m_row_map[i] + j];
+ ++j;
+ } else {
+ os << ' ' << i_overflow->first << ':' << i_overflow->second;
+ ++i_overflow;
+ }
+ }
+ }
+
+ os << '\n';
+ }
+ return os;
+ }
+
+ [[nodiscard]] PUGS_INLINE DataType&
+ operator()(const IndexType i, const IndexType j)
+ {
+ Assert(i < m_nb_rows, "invalid row index");
+ Assert(j < m_nb_columns, "invalid column index");
+
+ const IndexType row_start = m_row_map[i];
+ const IndexType row_end = m_row_map[i + 1];
+
+ IndexType position_candidate = 0;
+ bool found = false;
+ for (IndexType i_row = 0; i_row < m_nb_values_per_row[i]; ++i_row) {
+ if (m_column_indices[row_start + i_row] > j) {
+ position_candidate = i_row;
+ break;
+ } else {
+ if (m_column_indices[row_start + i_row] == j) {
+ position_candidate = i_row;
+ found = true;
+ break;
+ } else {
+ ++position_candidate;
+ }
+ }
+ }
+
+ if (not found) {
+ if (m_nb_values_per_row[i] < row_end - row_start) {
+ for (IndexType i_shift = 0; i_shift < m_nb_values_per_row[i] - position_candidate; ++i_shift) {
+ Assert(std::make_signed_t<IndexType>(row_start + m_nb_values_per_row[i]) -
+ std::make_signed_t<IndexType>(i_shift + 1) >=
+ 0);
+ const IndexType i_destination = row_start + m_nb_values_per_row[i] - i_shift;
+ const IndexType i_source = i_destination - 1;
+
+ m_column_indices[i_destination] = m_column_indices[i_source];
+ m_values[i_destination] = m_values[i_source];
+ }
+
+ m_column_indices[row_start + position_candidate] = j;
+ m_values[row_start + position_candidate] = 0;
+ ++m_nb_values_per_row[i];
+
+ return m_values[row_start + position_candidate];
+ } else {
+ auto& overflow_row = m_overflow_per_row[i];
+
+ auto iterator = overflow_row.insert(std::make_pair(j, DataType{0})).first;
+
+ return iterator->second;
+ }
+ } else {
+ return m_values[row_start + position_candidate];
+ }
+ }
+
+ CRSMatrix<DataType, IndexType>
+ getCRSMatrix() const
+ {
+ const bool is_filled = this->isFilled();
+ const bool has_overflow = this->hasOverflow();
+ if (is_filled and not has_overflow) {
+ return CRSMatrix<DataType, IndexType>{m_nb_rows, m_nb_columns, m_row_map, m_values, m_column_indices};
+ } else {
+ std::cout << rang::fgB::yellow << "warning:" << rang::style::reset
+ << " CRSMatrix profile was not properly set:\n";
+ if (not is_filled) {
+ std::cout << "- some rows are " << rang::style::bold << "too long" << rang::style::reset << '\n';
+ }
+ if (has_overflow) {
+ std::cout << "- some rows are " << rang::style::bold << "too short" << rang::style::reset << '\n';
+ }
+
+ Array<IndexType> row_map(m_nb_rows + 1);
+
+ row_map[0] = 0;
+ for (IndexType i = 0; i < m_nb_rows; ++i) {
+ row_map[i + 1] = row_map[i] + m_nb_values_per_row[i] + m_overflow_per_row[i].size();
+ }
+
+ IndexType nb_values = row_map[row_map.size() - 1];
+
+ Array<data_type> values(nb_values);
+ Array<IndexType> column_indices(nb_values);
+
+ IndexType l = 0;
+ for (IndexType i = 0; i < m_nb_rows; ++i) {
+ const auto& overflow_row = m_overflow_per_row[i];
+
+ if (m_nb_values_per_row[i] > 0) {
+ IndexType j = 0;
+ auto i_overflow = overflow_row.begin();
+ while (j < m_nb_values_per_row[i] or i_overflow != overflow_row.end()) {
+ const IndexType j_index = [&] {
+ if (j < m_nb_values_per_row[i]) {
+ return m_column_indices[m_row_map[i] + j];
+ } else {
+ return std::numeric_limits<IndexType>::max();
+ }
+ }();
+
+ const IndexType overflow_index = [&] {
+ if (i_overflow != overflow_row.end()) {
+ return i_overflow->first;
+ } else {
+ return std::numeric_limits<IndexType>::max();
+ }
+ }();
+
+ if (j_index < overflow_index) {
+ column_indices[l] = m_column_indices[m_row_map[i] + j];
+ values[l] = m_values[m_row_map[i] + j];
+ ++j;
+
+ } else {
+ column_indices[l] = i_overflow->first;
+ values[l] = i_overflow->second;
+ ++i_overflow;
+ }
+ ++l;
+ }
+ }
+ }
+
+ return CRSMatrix<DataType, IndexType>{m_nb_rows, m_nb_columns, row_map, values, column_indices};
+ }
+ }
+
+ CRSMatrixDescriptor(const IndexType nb_rows, const IndexType nb_columns, const Array<IndexType>& non_zeros)
+ : m_nb_rows{nb_rows},
+ m_nb_columns{nb_columns},
+ m_row_map{this->_computeRowMap(non_zeros)},
+ m_values(m_row_map[m_row_map.size() - 1]),
+ m_column_indices(m_row_map[m_row_map.size() - 1]),
+ m_nb_values_per_row(m_nb_rows),
+ m_overflow_per_row(m_nb_rows)
+ {
+ m_nb_values_per_row.fill(0);
+ m_values.fill(0);
+ }
+
+ ~CRSMatrixDescriptor() = default;
+};
+
+#endif // CRS_MATRIX_DESCRIPTOR_HPP
diff --git a/src/algebra/LinearSolver.cpp b/src/algebra/LinearSolver.cpp
index 98e2a3a1c..3cffca729 100644
--- a/src/algebra/LinearSolver.cpp
+++ b/src/algebra/LinearSolver.cpp
@@ -175,7 +175,7 @@ struct LinearSolver::Internals
const RHSVectorType& b,
const LinearSolverOptions& options)
{
- Assert(x.size() == b.size() and x.size() == A.numberOfColumns() and A.isSquare());
+ Assert((x.size() == b.size()) and (x.size() - A.numberOfColumns() == 0) and A.isSquare());
Vec petscB;
VecCreateMPIWithArray(PETSC_COMM_SELF, 1, b.size(), b.size(), &b[0], &petscB);
@@ -329,7 +329,7 @@ LinearSolver::checkOptions(const LinearSolverOptions& options) const
}
void
-LinearSolver::solveLocalSystem(const CRSMatrix<double, size_t>& A, Vector<double>& x, const Vector<const double>& b)
+LinearSolver::solveLocalSystem(const CRSMatrix<double, int>& A, Vector<double>& x, const Vector<const double>& b)
{
Internals::solveLocalSystem(A, x, b, m_options);
}
diff --git a/src/algebra/LinearSolver.hpp b/src/algebra/LinearSolver.hpp
index 387337400..1ab43a140 100644
--- a/src/algebra/LinearSolver.hpp
+++ b/src/algebra/LinearSolver.hpp
@@ -39,7 +39,7 @@ class LinearSolver
}
}
- void solveLocalSystem(const CRSMatrix<double, size_t>& A, Vector<double>& x, const Vector<const double>& b);
+ void solveLocalSystem(const CRSMatrix<double, int>& A, Vector<double>& x, const Vector<const double>& b);
void solveLocalSystem(const DenseMatrix<double>& A, Vector<double>& x, const Vector<const double>& b);
LinearSolver();
diff --git a/src/algebra/PETScUtils.cpp b/src/algebra/PETScUtils.cpp
index 5ed261353..63c5a68ac 100644
--- a/src/algebra/PETScUtils.cpp
+++ b/src/algebra/PETScUtils.cpp
@@ -35,32 +35,22 @@ PETScAijMatrixEmbedder::PETScAijMatrixEmbedder(const size_t nb_rows, const size_
MatAssemblyEnd(m_petscMat, MAT_FINAL_ASSEMBLY);
}
-PETScAijMatrixEmbedder::PETScAijMatrixEmbedder(const CRSMatrix<double, size_t>& A)
- : m_nb_rows{A.numberOfRows()}, m_nb_columns{A.numberOfColumns()}
+PETScAijMatrixEmbedder::PETScAijMatrixEmbedder(const CRSMatrix<double, int>& A)
+ : m_nb_rows(A.numberOfRows()), m_nb_columns(A.numberOfColumns())
{
- const Array<PetscReal>& values = copy(A.values());
+ static_assert(std::is_same_v<PetscInt, typename std::decay_t<decltype(A)>::index_type>,
+ "index type conversion is not implemented yet");
+ static_assert(std::is_same_v<PetscScalar, typename std::decay_t<decltype(A)>::data_type>,
+ "value type conversion is not implemented yet");
+ m_values = copy(A.values());
+ m_row_indices = copy(A.rowMap());
+ m_column_indices = copy(A.columnIndices());
- {
- 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);
+ MatCreateSeqAIJWithArrays(PETSC_COMM_SELF, A.numberOfRows(), A.numberOfColumns(),
+ const_cast<PetscInt*>(&m_row_indices[0]), // These ugly const_cast
+ const_cast<PetscInt*>(&m_column_indices[0]), // permit CRSMatrix to
+ const_cast<PetscScalar*>(&m_values[0]), // be non-mutable
+ &m_petscMat);
MatAssemblyBegin(m_petscMat, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(m_petscMat, MAT_FINAL_ASSEMBLY);
diff --git a/src/algebra/PETScUtils.hpp b/src/algebra/PETScUtils.hpp
index 68748c02f..3e7cff5b8 100644
--- a/src/algebra/PETScUtils.hpp
+++ b/src/algebra/PETScUtils.hpp
@@ -15,8 +15,9 @@ class PETScAijMatrixEmbedder
{
private:
Mat m_petscMat;
- Array<PetscInt> m_row_indices;
- Array<PetscInt> m_column_indices;
+ Array<const PetscScalar> m_values;
+ Array<const PetscInt> m_row_indices;
+ Array<const PetscInt> m_column_indices;
const size_t m_nb_rows;
const size_t m_nb_columns;
@@ -57,7 +58,7 @@ class PETScAijMatrixEmbedder
: PETScAijMatrixEmbedder{A.numberOfRows(), A.numberOfColumns(), &A(0, 0)}
{}
- PETScAijMatrixEmbedder(const CRSMatrix<double, size_t>& A);
+ PETScAijMatrixEmbedder(const CRSMatrix<double, int>& A);
~PETScAijMatrixEmbedder();
};
diff --git a/src/algebra/SparseMatrixDescriptor.hpp b/src/algebra/SparseMatrixDescriptor.hpp
deleted file mode 100644
index d94ed4686..000000000
--- a/src/algebra/SparseMatrixDescriptor.hpp
+++ /dev/null
@@ -1,170 +0,0 @@
-#ifndef SPARSE_MATRIX_DESCRIPTOR_HPP
-#define SPARSE_MATRIX_DESCRIPTOR_HPP
-
-#include <utils/Array.hpp>
-
-#include <map>
-#include <type_traits>
-
-template <typename DataType = double, typename IndexType = size_t>
-class SparseMatrixDescriptor
-{
- static_assert(std::is_integral_v<IndexType>, "Index type must be an integral type");
- static_assert(std::is_unsigned_v<IndexType>, "Index type must be unsigned");
-
- public:
- using data_type = DataType;
- using index_type = IndexType;
-
- class SparseRowDescriptor
- {
- friend class SparseMatrixDescriptor;
- std::map<IndexType, DataType> m_id_value_map;
-
- public:
- IndexType
- numberOfValues() const noexcept
- {
- return m_id_value_map.size();
- }
-
- const DataType&
- operator[](const IndexType& j) const
- {
- auto i_value = m_id_value_map.find(j);
- Assert(i_value != m_id_value_map.end());
- return i_value->second;
- }
-
- DataType&
- operator[](const IndexType& j)
- {
- auto i_value = m_id_value_map.find(j);
- if (i_value != m_id_value_map.end()) {
- return i_value->second;
- } else {
- auto [i_inserted, success] = m_id_value_map.insert(std::make_pair(j, 0));
- Assert(success);
- return i_inserted->second;
- }
- }
-
- friend std::ostream&
- operator<<(std::ostream& os, const SparseRowDescriptor& row)
- {
- for (auto [j, value] : row.m_id_value_map) {
- os << ' ' << static_cast<size_t>(j) << ':' << value;
- }
- return os;
- }
-
- SparseRowDescriptor() = default;
- };
-
- private:
- Array<SparseRowDescriptor> m_row_array;
- const size_t m_nb_columns;
-
- public:
- PUGS_INLINE
- size_t
- numberOfRows() const noexcept
- {
- return m_row_array.size();
- }
-
- PUGS_INLINE
- size_t
- numberOfColumns() const noexcept
- {
- return m_nb_columns;
- }
-
- SparseRowDescriptor&
- row(const IndexType i)
- {
- Assert(i < m_row_array.size());
- return m_row_array[i];
- }
-
- const SparseRowDescriptor&
- row(const IndexType i) const
- {
- Assert(i < m_row_array.size());
- return m_row_array[i];
- }
-
- DataType&
- operator()(const IndexType& i, const IndexType& j)
- {
- Assert(i < m_row_array.size());
- Assert(j < m_nb_columns);
- return m_row_array[i][j];
- }
-
- const DataType&
- operator()(const IndexType& i, const IndexType& j) const
- {
- Assert(i < 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];
- }
-
- friend std::ostream&
- operator<<(std::ostream& os, const SparseMatrixDescriptor& M)
- {
- for (IndexType i = 0; i < M.m_row_array.size(); ++i) {
- os << static_cast<size_t>(i) << " |" << M.m_row_array[i] << '\n';
- }
- return os;
- }
-
- auto
- graphVector() const
- {
- std::vector<std::vector<IndexType>> graph_vector(m_row_array.size());
- for (IndexType i_row = 0; i_row < m_row_array.size(); ++i_row) {
- const SparseRowDescriptor& row = m_row_array[i_row];
- for (auto [id, value] : row.m_id_value_map) {
- graph_vector[i_row].push_back(id);
- }
- }
- return graph_vector;
- }
-
- Array<DataType>
- valueArray() const
- {
- IndexType size = 0;
- for (IndexType i_row = 0; i_row < m_row_array.size(); ++i_row) {
- size += m_row_array[i_row].m_id_value_map.size();
- }
-
- Array<DataType> values(size);
-
- IndexType cpt = 0;
- for (IndexType i_row = 0; i_row < m_row_array.size(); ++i_row) {
- const SparseRowDescriptor& row = m_row_array[i_row];
- for (auto [id, value] : row.m_id_value_map) {
- values[cpt++] = value;
- }
- }
- return values;
- }
-
- SparseMatrixDescriptor(size_t nb_rows, size_t nb_columns) : m_row_array{nb_rows}, m_nb_columns{nb_columns}
- {
- 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;
-};
-
-#endif // SPARSE_MATRIX_DESCRIPTOR_HPP
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index a56fb1f9d..7904eedc0 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -56,6 +56,7 @@ add_executable (unit_tests
test_ConcatExpressionProcessor.cpp
test_CRSGraph.cpp
test_CRSMatrix.cpp
+ test_CRSMatrixDescriptor.cpp
test_DataVariant.cpp
test_DenseMatrix.cpp
test_Demangle.cpp
@@ -89,7 +90,6 @@ add_executable (unit_tests
test_PugsFunctionAdapter.cpp
test_PugsUtils.cpp
test_RevisionInfo.cpp
- test_SparseMatrixDescriptor.cpp
test_SymbolTable.cpp
test_Table.cpp
test_Timer.cpp
diff --git a/tests/test_BiCGStab.cpp b/tests/test_BiCGStab.cpp
index e1be80698..231783090 100644
--- a/tests/test_BiCGStab.cpp
+++ b/tests/test_BiCGStab.cpp
@@ -3,6 +3,7 @@
#include <algebra/BiCGStab.hpp>
#include <algebra/CRSMatrix.hpp>
+#include <algebra/CRSMatrixDescriptor.hpp>
// clazy:excludeall=non-pod-global-static
@@ -10,7 +11,13 @@ TEST_CASE("BiCGStab", "[algebra]")
{
SECTION("no preconditionner")
{
- SparseMatrixDescriptor S{5};
+ Array<int> non_zeros(5);
+ non_zeros[0] = 2;
+ non_zeros[1] = 3;
+ non_zeros[2] = 3;
+ non_zeros[3] = 3;
+ non_zeros[4] = 2;
+ CRSMatrixDescriptor<double> S{5, 5, non_zeros};
S(0, 0) = 2;
S(0, 1) = -1;
@@ -29,7 +36,7 @@ TEST_CASE("BiCGStab", "[algebra]")
S(4, 3) = 1;
S(4, 4) = 3;
- CRSMatrix A{S};
+ CRSMatrix A{S.getCRSMatrix()};
Vector<const double> x_exact = [] {
Vector<double> y{5};
@@ -53,7 +60,9 @@ TEST_CASE("BiCGStab", "[algebra]")
SECTION("no preconditionner non-converged")
{
- SparseMatrixDescriptor S{5};
+ Array<int> non_zeros(5);
+ non_zeros.fill(2);
+ CRSMatrixDescriptor<double> S{5, 5, non_zeros};
S(0, 0) = 2;
S(0, 1) = -1;
@@ -72,7 +81,7 @@ TEST_CASE("BiCGStab", "[algebra]")
S(4, 3) = 1;
S(4, 4) = 3;
- CRSMatrix A{S};
+ CRSMatrix A{S.getCRSMatrix()};
Vector<const double> x_exact = [] {
Vector<double> y{5};
diff --git a/tests/test_CG.cpp b/tests/test_CG.cpp
index 05919858f..b76d90a04 100644
--- a/tests/test_CG.cpp
+++ b/tests/test_CG.cpp
@@ -3,6 +3,7 @@
#include <algebra/CG.hpp>
#include <algebra/CRSMatrix.hpp>
+#include <algebra/CRSMatrixDescriptor.hpp>
// clazy:excludeall=non-pod-global-static
@@ -10,7 +11,11 @@ TEST_CASE("CG", "[algebra]")
{
SECTION("no preconditionner")
{
- SparseMatrixDescriptor S{5};
+ Array<int> non_zeros{5};
+ non_zeros.fill(3);
+ non_zeros[0] = 2;
+ non_zeros[4] = 2;
+ CRSMatrixDescriptor S{5, 5, non_zeros};
S(0, 0) = 2;
S(0, 1) = -1;
@@ -29,7 +34,7 @@ TEST_CASE("CG", "[algebra]")
S(4, 3) = -1;
S(4, 4) = 2;
- CRSMatrix A{S};
+ CRSMatrix A{S.getCRSMatrix()};
Vector<const double> x_exact = [] {
Vector<double> y{5};
@@ -53,9 +58,11 @@ TEST_CASE("CG", "[algebra]")
SECTION("singular matrix with RHS in its image")
{
- SparseMatrixDescriptor<double, uint8_t> S{5};
+ Array<int> non_zeros(5);
+ non_zeros.fill(0);
+ CRSMatrixDescriptor<double> S{5, 5, non_zeros};
- CRSMatrix<double, uint8_t> A{S};
+ CRSMatrix<double> A{S.getCRSMatrix()};
Vector<double> b{5};
b = 0;
@@ -69,7 +76,11 @@ TEST_CASE("CG", "[algebra]")
SECTION("no preconditionner non-converged")
{
- SparseMatrixDescriptor S{5};
+ Array<int> non_zeros(5);
+ non_zeros.fill(3);
+ non_zeros[0] = 2;
+ non_zeros[4] = 2;
+ CRSMatrixDescriptor<double> S{5, 5, non_zeros};
S(0, 0) = 2;
S(0, 1) = -1;
@@ -88,7 +99,7 @@ TEST_CASE("CG", "[algebra]")
S(4, 3) = -1;
S(4, 4) = 2;
- CRSMatrix A{S};
+ CRSMatrix A{S.getCRSMatrix()};
Vector<const double> x_exact = [] {
Vector<double> y{5};
diff --git a/tests/test_CRSMatrix.cpp b/tests/test_CRSMatrix.cpp
index d6bfdaba9..e0aa27b57 100644
--- a/tests/test_CRSMatrix.cpp
+++ b/tests/test_CRSMatrix.cpp
@@ -2,9 +2,7 @@
#include <catch2/matchers/catch_matchers_all.hpp>
#include <algebra/CRSMatrix.hpp>
-
-// Instantiate to ensure full coverage is performed
-template class SparseMatrixDescriptor<int, uint8_t>;
+#include <algebra/CRSMatrixDescriptor.hpp>
// clazy:excludeall=non-pod-global-static
@@ -12,7 +10,10 @@ TEST_CASE("CRSMatrix", "[algebra]")
{
SECTION("matrix size")
{
- SparseMatrixDescriptor<int, uint8_t> S{5};
+ Array<int> non_zeros{6};
+ non_zeros.fill(3);
+
+ CRSMatrixDescriptor<int> S{6, 5, non_zeros};
S(0, 2) = 3;
S(1, 2) = 11;
S(1, 1) = 1;
@@ -22,16 +23,19 @@ TEST_CASE("CRSMatrix", "[algebra]")
S(2, 2) = 4;
S(4, 3) = 2;
S(4, 4) = -2;
+ S(5, 2) = 1;
- CRSMatrix<int, uint8_t> A{S};
+ CRSMatrix<int> A{S.getCRSMatrix()};
REQUIRE(A.numberOfRows() == S.numberOfRows());
REQUIRE(A.numberOfColumns() == S.numberOfColumns());
}
- SECTION("matrix vector product (simple)")
+ SECTION("square matrix vector product (simple)")
{
- SparseMatrixDescriptor<int, uint8_t> S{5};
+ Array<int> non_zeros{5};
+ non_zeros.fill(3);
+ CRSMatrixDescriptor<int> S{5, 5, non_zeros};
S(0, 2) = 3;
S(1, 2) = 11;
S(1, 1) = 1;
@@ -42,9 +46,9 @@ TEST_CASE("CRSMatrix", "[algebra]")
S(4, 3) = 2;
S(4, 4) = -2;
- CRSMatrix<int, uint8_t> A{S};
+ CRSMatrix<int> A{S.getCRSMatrix()};
- Vector<int> x{A.numberOfColumns()};
+ Vector<int> x(A.numberOfColumns());
x = 0;
x[0] = 1;
@@ -96,9 +100,39 @@ TEST_CASE("CRSMatrix", "[algebra]")
REQUIRE(y[4] == -2);
}
+ SECTION("rectangular matrix vector product (simple)")
+ {
+ Array<int> non_zeros{3};
+ non_zeros.fill(3);
+ CRSMatrixDescriptor<int> S{3, 5, non_zeros};
+ S(0, 2) = 3;
+ S(0, 4) = 2;
+ S(1, 2) = 11;
+ S(1, 1) = 2;
+ S(1, 0) = 1;
+ S(2, 3) = 6;
+ S(2, 2) = 4;
+
+ CRSMatrix<int> A{S.getCRSMatrix()};
+
+ Vector<int> x(A.numberOfColumns());
+ x[0] = 1;
+ x[1] = 3;
+ x[2] = 2;
+ x[3] = -3;
+ x[4] = 5;
+
+ Vector<int> y = A * x;
+ REQUIRE(y[0] == 16);
+ REQUIRE(y[1] == 29);
+ REQUIRE(y[2] == -10);
+ }
+
SECTION("matrix vector product (complete)")
{
- SparseMatrixDescriptor<int, uint8_t> S{4};
+ Array<int> non_zeros{4};
+ non_zeros.fill(4);
+ CRSMatrixDescriptor<int> S{4, 4, non_zeros};
S(0, 0) = 1;
S(0, 1) = 2;
S(0, 2) = 3;
@@ -116,9 +150,9 @@ TEST_CASE("CRSMatrix", "[algebra]")
S(3, 2) = 15;
S(3, 3) = 16;
- CRSMatrix<int, uint8_t> A{S};
+ CRSMatrix<int> A{S.getCRSMatrix()};
- Vector<int> x{A.numberOfColumns()};
+ Vector<int> x(A.numberOfColumns());
x[0] = 1;
x[1] = 2;
x[2] = 3;
@@ -133,7 +167,10 @@ TEST_CASE("CRSMatrix", "[algebra]")
SECTION("check values")
{
- SparseMatrixDescriptor<int, uint8_t> S{4};
+ Array<int> non_zeros{4};
+ non_zeros.fill(3);
+ non_zeros[2] = 4;
+ CRSMatrixDescriptor<int> S{4, 4, non_zeros};
S(3, 0) = 13;
S(0, 0) = 1;
S(0, 1) = 2;
@@ -148,7 +185,7 @@ TEST_CASE("CRSMatrix", "[algebra]")
S(1, 0) = 5;
S(2, 1) = 10;
- CRSMatrix<int, uint8_t> A{S};
+ CRSMatrix<int> A{S.getCRSMatrix()};
auto values = A.values();
REQUIRE(values.size() == 13);
@@ -166,29 +203,46 @@ TEST_CASE("CRSMatrix", "[algebra]")
REQUIRE(values[11] == 15);
REQUIRE(values[12] == 16);
- auto row_indices = A.rowIndices();
-
- REQUIRE(row_indices.size() == 5);
-
- REQUIRE(A.row(0).colidx(0) == 0);
- REQUIRE(A.row(0).colidx(1) == 1);
- REQUIRE(A.row(0).colidx(2) == 3);
- REQUIRE(A.row(1).colidx(0) == 0);
- REQUIRE(A.row(1).colidx(1) == 1);
- REQUIRE(A.row(1).colidx(2) == 2);
- REQUIRE(A.row(2).colidx(0) == 0);
- REQUIRE(A.row(2).colidx(1) == 1);
- REQUIRE(A.row(2).colidx(2) == 2);
- REQUIRE(A.row(2).colidx(3) == 3);
- REQUIRE(A.row(3).colidx(0) == 0);
- REQUIRE(A.row(3).colidx(1) == 2);
- REQUIRE(A.row(3).colidx(2) == 3);
+ auto row_map = A.rowMap();
+ REQUIRE(row_map.size() == 5);
+ REQUIRE(row_map[0] == 0);
+ REQUIRE(row_map[1] == 3);
+ REQUIRE(row_map[2] == 6);
+ REQUIRE(row_map[3] == 10);
+ REQUIRE(row_map[4] == 13);
+
+ auto column_indices = A.columnIndices();
+ REQUIRE(column_indices.size() == 13);
+ REQUIRE(column_indices[0] == 0);
+ REQUIRE(column_indices[1] == 1);
+ REQUIRE(column_indices[2] == 3);
+ REQUIRE(column_indices[3] == 0);
+ REQUIRE(column_indices[4] == 1);
+ REQUIRE(column_indices[5] == 2);
+ REQUIRE(column_indices[6] == 0);
+ REQUIRE(column_indices[7] == 1);
+ REQUIRE(column_indices[8] == 2);
+ REQUIRE(column_indices[9] == 3);
+ REQUIRE(column_indices[10] == 0);
+ REQUIRE(column_indices[11] == 2);
+ REQUIRE(column_indices[12] == 3);
+
+ std::ostringstream crs_descriptor_os;
+ std::ostringstream crs_os;
+
+ crs_descriptor_os << S;
+ crs_os << A;
+
+ REQUIRE(crs_descriptor_os.str() == crs_os.str());
}
#ifndef NDEBUG
SECTION("incompatible runtime matrix/vector product")
{
- CRSMatrix<int, uint8_t> A{SparseMatrixDescriptor<int, uint8_t>{4}};
+ Array<int> non_zeros(2);
+ non_zeros.fill(0);
+ CRSMatrixDescriptor<int> S{2, 4, non_zeros};
+ CRSMatrix<int> A{S.getCRSMatrix()};
Vector<int> x{2};
REQUIRE_THROWS_AS(A * x, AssertError);
}
diff --git a/tests/test_CRSMatrixDescriptor.cpp b/tests/test_CRSMatrixDescriptor.cpp
new file mode 100644
index 000000000..2bf59af22
--- /dev/null
+++ b/tests/test_CRSMatrixDescriptor.cpp
@@ -0,0 +1,166 @@
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_all.hpp>
+
+#include <utils/PugsAssert.hpp>
+
+#include <algebra/CRSMatrixDescriptor.hpp>
+
+// Instantiate to ensure full coverage is performed
+template class CRSMatrixDescriptor<int>;
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("CRSMatrixDescriptor", "[algebra]")
+{
+ SECTION("has overflow / not filled")
+ {
+ const size_t nb_lines = 2;
+ const size_t nb_columns = 5;
+
+ Array<int> non_zeros{nb_lines};
+ non_zeros.fill(2);
+ CRSMatrixDescriptor<int> S(nb_lines, nb_columns, non_zeros);
+
+ S(1, 3) = 5;
+ S(1, 1) += 3;
+ S(0, 1) = 4;
+ S(1, 0) = 2;
+ S(1, 2) = 7;
+
+ REQUIRE(S.numberOfRows() == 2);
+ REQUIRE(S.numberOfColumns() == 5);
+
+ REQUIRE(S.isFilled() == false);
+ REQUIRE(S.hasOverflow() == true);
+
+ std::ostringstream os;
+ os << S;
+
+ std::string_view output =
+ R"(0| 1:4
+1| 0:2 1:3 2:7 3:5
+)";
+ REQUIRE(os.str() == output);
+ }
+
+ SECTION("no overflow / not filled")
+ {
+ const size_t nb_lines = 2;
+ const size_t nb_columns = 5;
+
+ Array<int> non_zeros{nb_lines};
+ non_zeros.fill(2);
+ CRSMatrixDescriptor<int> S(nb_lines, nb_columns, non_zeros);
+
+ S(1, 3) = 5;
+ S(1, 1) += 3;
+ S(0, 1) = 4;
+
+ REQUIRE(S.numberOfRows() == 2);
+ REQUIRE(S.numberOfColumns() == 5);
+
+ REQUIRE(S.isFilled() == false);
+ REQUIRE(S.hasOverflow() == false);
+
+ std::ostringstream os;
+ os << S;
+
+ std::string_view output =
+ R"(0| 1:4
+1| 1:3 3:5
+)";
+ REQUIRE(os.str() == output);
+ }
+
+ SECTION("no overflow / filled")
+ {
+ const size_t nb_lines = 2;
+ const size_t nb_columns = 5;
+
+ Array<int> non_zeros{nb_lines};
+ non_zeros[0] = 1;
+ non_zeros[1] = 5;
+ CRSMatrixDescriptor<int> S(nb_lines, nb_columns, non_zeros);
+
+ S(1, 3) = 2;
+ S(1, 3) += 3;
+ S(1, 1) += 3;
+ S(0, 1) = 4;
+ S(1, 2) = 2;
+ S(1, 0) = 1;
+ S(1, 4) = 13;
+
+ REQUIRE(S.numberOfRows() == 2);
+ REQUIRE(S.numberOfColumns() == 5);
+
+ REQUIRE(S.isFilled() == true);
+ REQUIRE(S.hasOverflow() == false);
+
+ std::ostringstream os;
+ os << S;
+
+ std::string_view output =
+ R"(0| 1:4
+1| 0:1 1:3 2:2 3:5 4:13
+)";
+ REQUIRE(os.str() == output);
+ }
+
+ SECTION("has overflow / filled")
+ {
+ const size_t nb_lines = 2;
+ const size_t nb_columns = 5;
+
+ Array<int> non_zeros{nb_lines};
+ non_zeros.fill(2);
+ CRSMatrixDescriptor<int> S(nb_lines, nb_columns, non_zeros);
+
+ S(1, 3) = 3;
+ S(1, 0) = -2;
+ S(1, 1) += 1;
+ S(1, 1) += 1;
+ S(1, 1) = 3;
+ S(0, 1) = 4;
+ S(0, 0) = 1;
+ S(0, 3) = 5;
+ S(0, 2) = 1;
+
+ REQUIRE(S.numberOfRows() == 2);
+ REQUIRE(S.numberOfColumns() == 5);
+
+ REQUIRE(S.isFilled() == true);
+ REQUIRE(S.hasOverflow() == true);
+
+ std::ostringstream os;
+ os << S;
+
+ std::string_view output =
+ R"(0| 0:1 1:4 2:1 3:5
+1| 0:-2 1:3 3:3
+)";
+ REQUIRE(os.str() == output);
+ }
+
+#ifndef NDEBUG
+ SECTION("incompatible non zero vector size and number of columns")
+ {
+ Array<int> non_zeros(3);
+ REQUIRE_THROWS_AS((CRSMatrixDescriptor<int>{2, 4, non_zeros}), AssertError);
+ }
+
+ SECTION("bad row number")
+ {
+ Array<int> non_zeros(2);
+ CRSMatrixDescriptor<int> S{2, 4, non_zeros};
+ REQUIRE_THROWS_AS(S(2, 3) = 1, AssertError);
+ }
+
+ SECTION("bad column number")
+ {
+ Array<int> non_zeros(2);
+ CRSMatrixDescriptor<int> S{2, 4, non_zeros};
+ REQUIRE_THROWS_AS(S(0, 5) = 1, AssertError);
+ }
+
+#endif // NDEBUG
+}
diff --git a/tests/test_EigenvalueSolver.cpp b/tests/test_EigenvalueSolver.cpp
index e779bce56..3ca8f7f51 100644
--- a/tests/test_EigenvalueSolver.cpp
+++ b/tests/test_EigenvalueSolver.cpp
@@ -2,9 +2,9 @@
#include <catch2/matchers/catch_matchers_all.hpp>
#include <algebra/CRSMatrix.hpp>
+#include <algebra/CRSMatrixDescriptor.hpp>
#include <algebra/DenseMatrix.hpp>
#include <algebra/EigenvalueSolver.hpp>
-#include <algebra/SparseMatrixDescriptor.hpp>
#include <utils/pugs_config.hpp>
@@ -16,7 +16,10 @@ TEST_CASE("EigenvalueSolver", "[algebra]")
{
SECTION("symmetric system")
{
- SparseMatrixDescriptor S{3};
+ Array<int> non_zeros(3);
+ non_zeros.fill(3);
+ CRSMatrixDescriptor<double> S{3, 3, non_zeros};
+
S(0, 0) = 3;
S(0, 1) = 2;
S(0, 2) = 4;
@@ -29,7 +32,7 @@ TEST_CASE("EigenvalueSolver", "[algebra]")
S(2, 1) = 2;
S(2, 2) = 3;
- CRSMatrix A{S};
+ CRSMatrix A{S.getCRSMatrix()};
SECTION("eigenvalues")
{
diff --git a/tests/test_LinearSolver.cpp b/tests/test_LinearSolver.cpp
index 1219c48e6..a9b4cab1f 100644
--- a/tests/test_LinearSolver.cpp
+++ b/tests/test_LinearSolver.cpp
@@ -3,6 +3,7 @@
#include <utils/pugs_config.hpp>
+#include <algebra/CRSMatrixDescriptor.hpp>
#include <algebra/LinearSolver.hpp>
// clazy:excludeall=non-pod-global-static
@@ -147,7 +148,11 @@ TEST_CASE("LinearSolver", "[algebra]")
{
SECTION("symmetric system")
{
- SparseMatrixDescriptor S{5};
+ Array<int> non_zeros{5};
+ non_zeros.fill(3);
+ non_zeros[0] = 2;
+ non_zeros[4] = 2;
+ CRSMatrixDescriptor S{5, 5, non_zeros};
S(0, 0) = 2;
S(0, 1) = -1;
@@ -166,7 +171,7 @@ TEST_CASE("LinearSolver", "[algebra]")
S(4, 3) = -1;
S(4, 4) = 2;
- CRSMatrix A{S};
+ CRSMatrix A{S.getCRSMatrix()};
Vector<const double> x_exact = [] {
Vector<double> y{5};
@@ -302,7 +307,12 @@ TEST_CASE("LinearSolver", "[algebra]")
SECTION("none symmetric system")
{
- SparseMatrixDescriptor S{5};
+ Array<int> non_zeros{5};
+ non_zeros.fill(3);
+ non_zeros[0] = 2;
+ non_zeros[4] = 2;
+ CRSMatrixDescriptor S{5, 5, non_zeros};
+
S(0, 0) = 2;
S(0, 1) = -1;
@@ -321,7 +331,7 @@ TEST_CASE("LinearSolver", "[algebra]")
S(4, 3) = 1;
S(4, 4) = 3;
- CRSMatrix A{S};
+ CRSMatrix A{S.getCRSMatrix()};
Vector<const double> x_exact = [] {
Vector<double> y{5};
diff --git a/tests/test_SparseMatrixDescriptor.cpp b/tests/test_SparseMatrixDescriptor.cpp
deleted file mode 100644
index 77bc890f7..000000000
--- a/tests/test_SparseMatrixDescriptor.cpp
+++ /dev/null
@@ -1,201 +0,0 @@
-#include <catch2/catch_test_macros.hpp>
-#include <catch2/matchers/catch_matchers_all.hpp>
-
-#include <utils/PugsAssert.hpp>
-
-#include <algebra/SparseMatrixDescriptor.hpp>
-
-// Instantiate to ensure full coverage is performed
-template class SparseMatrixDescriptor<int, uint8_t>;
-
-// clazy:excludeall=non-pod-global-static
-
-TEST_CASE("SparseMatrixDescriptor", "[algebra]")
-{
- SECTION("SparseRowDescriptor subclass")
- {
- SECTION("number of values")
- {
- SparseMatrixDescriptor<int, uint8_t>::SparseRowDescriptor r;
- REQUIRE(r.numberOfValues() == 0);
-
- r[0] = 3;
- r[0] += 2;
-
- REQUIRE(r.numberOfValues() == 1);
-
- r[1] = -2;
-
- REQUIRE(r.numberOfValues() == 2);
- }
-
- SECTION("values access")
- {
- SparseMatrixDescriptor<int, uint8_t>::SparseRowDescriptor r;
-
- r[0] = 3;
- r[0] += 2;
-
- r[3] = 8;
-
- REQUIRE(r[0] == 5);
- REQUIRE(r[3] == 8);
-
- const auto& const_r = r;
- REQUIRE(const_r[0] == 5);
- REQUIRE(const_r[3] == 8);
-
-#ifndef NDEBUG
- REQUIRE_THROWS_AS(const_r[2], AssertError);
-#endif // NDEBUG
- }
- }
-
- SECTION("matrix size")
- {
- SparseMatrixDescriptor<int, uint8_t> S{5};
- REQUIRE(S.numberOfRows() == 5);
- REQUIRE(S.numberOfColumns() == 5);
- }
-
- SECTION("location operators")
- {
- SparseMatrixDescriptor<int, uint8_t> S{5};
- S(0, 2) = 3;
- S(1, 1) = 1;
- S(1, 2) = 11;
- S(0, 2) += 2;
- S(2, 2) = 4;
- S(3, 1) = 5;
- S(4, 1) = 1;
- S(4, 4) = -2;
-
- REQUIRE(S.row(0).numberOfValues() == 2);
- REQUIRE(S.row(1).numberOfValues() == 2);
- REQUIRE(S.row(2).numberOfValues() == 1);
- REQUIRE(S.row(3).numberOfValues() == 2);
- REQUIRE(S.row(4).numberOfValues() == 2);
-
- const auto& const_S = S;
-
- REQUIRE(const_S.row(0).numberOfValues() == 2);
- REQUIRE(const_S.row(1).numberOfValues() == 2);
- REQUIRE(const_S.row(2).numberOfValues() == 1);
- REQUIRE(const_S.row(3).numberOfValues() == 2);
- REQUIRE(const_S.row(4).numberOfValues() == 2);
-
-#ifndef NDEBUG
- REQUIRE_THROWS_AS(S.row(5), AssertError);
- REQUIRE_THROWS_AS(const_S.row(5), AssertError);
-#endif // NDEBUG
-
- REQUIRE(S(0, 2) == 5);
- REQUIRE(S(1, 1) == 1);
- REQUIRE(S(1, 2) == 11);
- REQUIRE(S(2, 2) == 4);
- REQUIRE(S(3, 1) == 5);
- REQUIRE(S(4, 1) == 1);
- REQUIRE(S(4, 4) == -2);
-
- REQUIRE(const_S(0, 2) == 5);
- REQUIRE(const_S(1, 1) == 1);
- REQUIRE(const_S(1, 2) == 11);
- REQUIRE(const_S(2, 2) == 4);
- REQUIRE(const_S(3, 1) == 5);
- REQUIRE(const_S(4, 1) == 1);
- REQUIRE(const_S(4, 4) == -2);
-
-#ifndef NDEBUG
- REQUIRE_THROWS_AS(S(5, 0), AssertError);
- REQUIRE_THROWS_AS(const_S(6, 1), AssertError);
- REQUIRE_THROWS_AS(const_S(0, 1), AssertError);
-#endif // NDEBUG
- }
-
- SECTION("vector-graph")
- {
- SparseMatrixDescriptor<int, uint8_t> S{5};
- S(0, 2) = 3;
- S(1, 2) = 11;
- S(1, 1) = 1;
- S(0, 2) += 2;
- S(3, 3) = 5;
- S(3, 1) = 5;
- S(4, 1) = 1;
- S(2, 2) = 4;
- S(4, 4) = -2;
-
- const auto graph = S.graphVector();
-
- REQUIRE(graph.size() == S.numberOfRows());
- REQUIRE(graph[0].size() == 2);
- REQUIRE(graph[1].size() == 2);
- REQUIRE(graph[2].size() == 1);
- REQUIRE(graph[3].size() == 2);
- REQUIRE(graph[4].size() == 2);
-
- REQUIRE(graph[0][0] == 0);
- REQUIRE(graph[0][1] == 2);
- REQUIRE(graph[1][0] == 1);
- REQUIRE(graph[1][1] == 2);
- REQUIRE(graph[2][0] == 2);
- REQUIRE(graph[3][0] == 1);
- REQUIRE(graph[3][1] == 3);
- REQUIRE(graph[4][0] == 1);
- REQUIRE(graph[4][1] == 4);
- }
-
- SECTION("value array")
- {
- SparseMatrixDescriptor<int, uint8_t> S{5};
- S(0, 2) = 3;
- S(1, 2) = 11;
- S(1, 1) = 1;
- S(0, 2) += 2;
- S(3, 3) = 5;
- S(3, 1) = -3;
- S(4, 1) = 1;
- S(2, 2) = 4;
- S(4, 4) = -2;
-
- const auto value_array = S.valueArray();
-
- REQUIRE(value_array.size() == 9);
-
- REQUIRE(value_array[0] == 0);
- REQUIRE(value_array[1] == 5);
- REQUIRE(value_array[2] == 1);
- REQUIRE(value_array[3] == 11);
- REQUIRE(value_array[4] == 4);
- REQUIRE(value_array[5] == -3);
- REQUIRE(value_array[6] == 5);
- REQUIRE(value_array[7] == 1);
- REQUIRE(value_array[8] == -2);
- }
-
- SECTION("output")
- {
- SparseMatrixDescriptor<int, uint8_t> S{5};
- S(0, 2) = 3;
- S(1, 2) = 11;
- S(1, 1) = 1;
- S(0, 2) += 2;
- S(3, 3) = 5;
- S(3, 1) = -3;
- S(4, 1) = 1;
- S(2, 2) = 4;
- S(4, 4) = -2;
-
- std::ostringstream output;
- output << '\n' << S;
-
- std::string expected_output = R"(
-0 | 0:0 2:5
-1 | 1:1 2:11
-2 | 2:4
-3 | 1:-3 3:5
-4 | 1:1 4:-2
-)";
- REQUIRE(output.str() == expected_output);
- }
-}
--
GitLab