diff --git a/src/algebra/CRSMatrix.hpp b/src/algebra/CRSMatrix.hpp
index edb85d5e87b6a347ad5bc4f82987fd9255154484..bf650afe2784035023fbb5311c1c3a12c68b2a3b 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 0000000000000000000000000000000000000000..cbf82a904e85dcb8aa735e978612a4baec0affa7
--- /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 98e2a3a1c2caf8a041b6e69e71502bb6441647e7..3cffca7298b0a7377996c130195ac9e43a66ca44 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 3873374002a146c79a447a815d4ef6a57669883e..1ab43a1405ee312fdb9255c6b1f4508ade4e168e 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 5ed261353d93673267c3f4107c35f855bf643211..63c5a68ac5f3a248ba7e5bc51bf381f5418dbe4a 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 68748c02f5b36092cba6dce8ba3281a6d1cbb93b..3e7cff5b87b843f84575ed9b2af4a425e6ed4f5a 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 d94ed468606a10ee15857a971abab20fb4be5b1a..0000000000000000000000000000000000000000
--- 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 a56fb1f9d0b49108f1770f186cec453758cfd793..7904eedc03c7a203667e5c604096a5cded6e136b 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 e1be8069885b80a3e26ede233dac3941b7525521..231783090756e52187b46da02a97d44cf2279ea9 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 05919858f3ada681eb371883c6fcf6fd7d15a724..b76d90a047f2aab41e0576ab6f69ed5b133b474e 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 d6bfdaba966c0c36a5245af457954a6f53a0919d..e0aa27b57c1a9fd773c99d7d96c6d5e47d79f74e 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 0000000000000000000000000000000000000000..2bf59af22201ee9a8f164d2d6344b16357c56303
--- /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 e779bce5600e32cf26b0fde7339693f8da69bcb1..3ca8f7f5136d72f7003357113b7b597d1b84ef54 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 1219c48e62d28a749383e115aef8b6e51a17e5ec..a9b4cab1fe49efafa6a8ec4d38e34b7b947c6eff 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 77bc890f7c4a53ea170de08bebb83de9b7e0fbed..0000000000000000000000000000000000000000
--- 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);
- }
-}