Skip to content
Snippets Groups Projects
Commit fcedbaf6 authored by Stéphane Del Pino's avatar Stéphane Del Pino
Browse files

Displace SparseMatrixDescriptor in its own class

SparseMatrixDescriptor is also renamed from FlexibleMatrix. This is more clear
since it is actually not a matrix but a tool to allow simple construction of
matrices.
parent f7b3ec60
No related branches found
No related tags found
1 merge request!26Feature/linear systems
......@@ -5,129 +5,12 @@
#include <Kokkos_StaticCrsGraph.hpp>
#include <PugsAssert.hpp>
#include <SparseMatrixDescriptor.hpp>
#include <Vector.hpp>
#include <iostream>
template <typename DataType = double, typename IndexType = size_t>
class FlexibleMatrix
{
public:
using data_type = DataType;
using index_type = IndexType;
class FlexibleRow
{
friend class FlexibleMatrix;
std::map<IndexType, DataType> m_id_value_map;
public:
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 FlexibleRow& row)
{
for (auto [j, value] : row.m_id_value_map) {
os << ' ' << j << ':' << value;
}
return os;
}
FlexibleRow() = default;
};
private:
Array<FlexibleRow> m_row_array;
public:
FlexibleRow&
row(const IndexType i)
{
return m_row_array[i];
}
const FlexibleRow&
row(const IndexType i) const
{
return m_row_array[i];
}
DataType&
operator()(const IndexType& i, const IndexType& j)
{
return m_row_array[i][j];
}
const DataType&
operator()(const IndexType& i, const IndexType& j) const
{
return m_row_array[i][j];
}
friend std::ostream&
operator<<(std::ostream& os, const FlexibleMatrix& M)
{
for (IndexType i = 0; i < M.m_row_array.size(); ++i) {
os << 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 FlexibleRow& 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 FlexibleRow& row = m_row_array[i_row];
for (auto [id, value] : row.m_id_value_map) {
values[cpt++] = value;
}
}
return values;
}
FlexibleMatrix(IndexType nb_row) : m_row_array{nb_row} {}
~FlexibleMatrix() = default;
};
template <typename DataType = double, typename IndexType = size_t>
class CRSMatrix
{
......@@ -173,7 +56,7 @@ class CRSMatrix
return os;
}
CRSMatrix(const FlexibleMatrix<DataType, IndexType>& M)
CRSMatrix(const SparseMatrixDescriptor<DataType, IndexType>& M)
{
m_host_matrix = Kokkos::create_staticcrsgraph<HostMatrix>("connectivity_matrix", M.graphVector());
m_values = M.valueArray();
......
#ifndef SPARSE_MATRIX_DESCRIPTOR_HPP
#define SPARSE_MATRIX_DESCRIPTOR_HPP
#include <Array.hpp>
#include <map>
template <typename DataType = double, typename IndexType = size_t>
class SparseMatrixDescriptor
{
public:
using data_type = DataType;
using index_type = IndexType;
class SparseRowDescriptor
{
friend class SparseMatrixDescriptor;
std::map<IndexType, DataType> m_id_value_map;
public:
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 << ' ' << j << ':' << value;
}
return os;
}
SparseRowDescriptor() = default;
};
private:
Array<SparseRowDescriptor> m_row_array;
public:
SparseRowDescriptor&
row(const IndexType i)
{
return m_row_array[i];
}
const SparseRowDescriptor&
row(const IndexType i) const
{
return m_row_array[i];
}
DataType&
operator()(const IndexType& i, const IndexType& j)
{
return m_row_array[i][j];
}
const DataType&
operator()(const IndexType& i, const IndexType& j) const
{
return m_row_array[i][j];
}
friend std::ostream&
operator<<(std::ostream& os, const SparseMatrixDescriptor& M)
{
for (IndexType i = 0; i < M.m_row_array.size(); ++i) {
os << 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(IndexType nb_row) : m_row_array{nb_row} {}
~SparseMatrixDescriptor() = default;
};
#endif // SPARSE_MATRIX_DESCRIPTOR_HPP
......@@ -44,17 +44,17 @@ main(int argc, char* argv[])
{
{ // PCG simple test
size_t matrix_size = 10;
FlexibleMatrix F(matrix_size);
SparseMatrixDescriptor D(matrix_size);
for (size_t i = 0; i < matrix_size; ++i) {
F(i, i) = 4;
D(i, i) = 4;
if (i + 1 < matrix_size) {
F(i, i + 1) = -1;
F(i + 1, i) = -1;
D(i, i + 1) = -1;
D(i + 1, i) = -1;
}
}
CRSMatrix A{F};
CRSMatrix A{D};
Vector<double> x{10};
for (size_t i = 0; i < x.size(); ++i) {
......@@ -77,17 +77,17 @@ main(int argc, char* argv[])
{ // BiCGStab simple test
size_t matrix_size = 10;
FlexibleMatrix F(matrix_size);
SparseMatrixDescriptor D(matrix_size);
for (size_t i = 0; i < matrix_size; ++i) {
F(i, i) = 4;
D(i, i) = 4;
if (i + 1 < matrix_size) {
F(i, i + 1) = -1;
F(i + 1, i) = -0.3;
D(i, i + 1) = -1;
D(i + 1, i) = -0.3;
}
}
CRSMatrix A{F};
CRSMatrix A{D};
Vector<double> x{10};
for (size_t i = 0; i < x.size(); ++i) {
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment