Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Sparsity pattern insert #3338

Merged
merged 15 commits into from
Aug 27, 2024
Merged

Sparsity pattern insert #3338

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
15 commits
Select commit Hold shift + click to select a range
a01b318
Add single index insert to sparsity pattern.
schnellerhase Aug 7, 2024
d8ea430
Remove some superfluous {}
schnellerhase Aug 7, 2024
9982c7e
export new insert function
schnellerhase Aug 7, 2024
c946c0f
Make use of new insert functionality to simplify code
schnellerhase Aug 7, 2024
c1b5746
Extend insert doc to make clear how the block insertion works
schnellerhase Aug 7, 2024
01c9236
Fix doc
schnellerhase Aug 7, 2024
f5ca9d5
Merge branch 'main' into sparsity_pattern_insert
schnellerhase Aug 7, 2024
b50640e
to -> into
schnellerhase Aug 12, 2024
eae7cb9
Remove duplicate comment
schnellerhase Aug 12, 2024
51675dc
Merge branch 'main' into sparsity_pattern_insert
schnellerhase Aug 12, 2024
0766807
Merge branch 'main' into sparsity_pattern_insert
chrisrichardson Aug 27, 2024
7dee1fa
Remove dense matrix example
schnellerhase Aug 27, 2024
2791e0a
{} style adapted
schnellerhase Aug 27, 2024
b7c310e
Missed some
schnellerhase Aug 27, 2024
28ec172
Update doctring
schnellerhase Aug 27, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
22 changes: 21 additions & 1 deletion cpp/dolfinx/la/SparsityPattern.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -140,6 +140,27 @@ SparsityPattern::SparsityPattern(
}
}
//-----------------------------------------------------------------------------
void SparsityPattern::insert(std::int32_t row, std::int32_t col)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Isn't this covered by passing rows and cols of length 1?

Copy link
Contributor Author

@schnellerhase schnellerhase Aug 7, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It it, but the interface of the caller becomes cleaner. Instead of needing to call

sp.insert(std::vector<std::int32_t>{row}, std::vector<std::int32_t>{col});

we can now simply call

sp.insert(row, col);

same holds for the python interface with numpy arrays.

Also, in combination with the span based version it makes clearer how the sparsity pattern is to be used, i.e. loop over entries and dense matrices and do not pre-compute the entries and then pass it in bulk. Sparsity pattern already performs the caching and associated memory management.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Well, it's clear from the improved docstring. I don't see the point of this new method, the clarity of overloading functions is debatable.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

To add, the interface is vectorised from the Python end to encourage users to pass in multiple entries, avoiding the overhead of the function call.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

How/where is the interface vectorized? I can not find that.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@chrisrichardson mentioned the other day, we might instead of the integer version consider to create a non owning const std::vector<std::int32_t>& version, this should allow for the interface to be insert({row}, {col}), i.e. it allows for conversion from an initializer list.

{
if (!_offsets.empty())
{
throw std::runtime_error(
"Cannot insert into sparsity pattern. It has already been finalized");
}

assert(_index_maps[0]);
const std::int32_t max_row
= _index_maps[0]->size_local() + _index_maps[0]->num_ghosts() - 1;

if (row > max_row or row < 0)
{
throw std::runtime_error(
"Cannot insert rows that do not exist in the IndexMap.");
}

_row_cache[row].push_back(col);
}
//-----------------------------------------------------------------------------
void SparsityPattern::insert(std::span<const std::int32_t> rows,
std::span<const std::int32_t> cols)
{
Expand All @@ -160,7 +181,6 @@ void SparsityPattern::insert(std::span<const std::int32_t> rows,
throw std::runtime_error(
"Cannot insert rows that do not exist in the IndexMap.");
}

_row_cache[row].insert(_row_cache[row].end(), cols.begin(), cols.end());
}
}
Expand Down
13 changes: 13 additions & 0 deletions cpp/dolfinx/la/SparsityPattern.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -67,6 +67,19 @@ class SparsityPattern

/// @brief Insert non-zero locations using local (process-wise)
/// indices.
/// @param[in] row local row index
/// @param[in] col local column index
void insert(std::int32_t row, std::int32_t col);

/// @brief Insert non-zero locations using local (process-wise)
/// indices.
///
/// This routine inserts non-zero locations at the outer product of rows and
/// cols into the sparsity pattern, i.e. adds the matrix entries at
/// A[row[i], col[j]] for all i, j.
///
/// @param[in] rows list of the local row indices
/// @param[in] cols list of the local column indices
void insert(std::span<const std::int32_t> rows,
std::span<const std::int32_t> cols);

Expand Down
6 changes: 3 additions & 3 deletions cpp/test/matrix.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -201,9 +201,9 @@ void test_matrix()
{
auto map0 = std::make_shared<common::IndexMap>(MPI_COMM_SELF, 8);
la::SparsityPattern p(MPI_COMM_SELF, {map0, map0}, {1, 1});
p.insert(std::vector{0}, std::vector{0});
p.insert(std::vector{4}, std::vector{5});
p.insert(std::vector{5}, std::vector{4});
p.insert(0, 0);
p.insert(4, 5);
p.insert(5, 4);
p.finalize();

using T = float;
Expand Down
5 changes: 5 additions & 0 deletions python/dolfinx/wrappers/la.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@
#include "caster_mpi.h"
#include "numpy_dtype.h"
#include <complex>
#include <cstdint>
#include <dolfinx/common/IndexMap.h>
#include <dolfinx/la/MatrixCSR.h>
#include <dolfinx/la/SparsityPattern.h>
Expand Down Expand Up @@ -265,6 +266,10 @@ void la(nb::module_& m)
std::span(cols.data(), cols.size()));
},
nb::arg("rows"), nb::arg("cols"))
.def("insert",
nb::overload_cast<int32_t, int32_t>(
&dolfinx::la::SparsityPattern::insert),
nb::arg("row"), nb::arg("col"))
.def(
"insert_diagonal",
[](dolfinx::la::SparsityPattern& self,
Expand Down
8 changes: 4 additions & 4 deletions python/test/unit/la/test_matrix_csr.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,9 +25,9 @@ def create_test_sparsity(n, bs):
if bs == 1:
for i in range(2):
for j in range(2):
sp.insert(np.array([2 + i]), np.array([4 + j]))
sp.insert(2 + i, 4 + j)
elif bs == 2:
sp.insert(np.array([1]), np.array([2]))
sp.insert(1, 2)
sp.finalize()
return sp

Expand Down Expand Up @@ -126,10 +126,10 @@ def test_distributed_csr(dtype):
sp = SparsityPattern(MPI.COMM_WORLD, [im, im], [1, 1])
for i in range(n):
for j in range(n + nghost):
sp.insert(np.array([i]), np.array([j]))
sp.insert(i, j)
for i in range(n, n + nghost):
for j in range(n, n + nghost):
sp.insert(np.array([i]), np.array([j]))
sp.insert(i, j)
sp.finalize()

mat = matrix_csr(sp, dtype=dtype)
Expand Down
Loading