DOLFINx 0.10.0.0
DOLFINx C++ interface
Loading...
Searching...
No Matches
MatrixCSR.h
1// Copyright (C) 2021-2022 Garth N. Wells and Chris N. Richardson
2//
3// This file is part of DOLFINx (https://www.fenicsproject.org)
4//
5// SPDX-License-Identifier: LGPL-3.0-or-later
6
7#pragma once
8
9#include "SparsityPattern.h"
10#include "Vector.h"
11#include "matrix_csr_impl.h"
12#include <algorithm>
13#include <dolfinx/common/IndexMap.h>
14#include <dolfinx/common/MPI.h>
15#include <dolfinx/graph/AdjacencyList.h>
16#include <mpi.h>
17#include <numeric>
18#include <span>
19#include <utility>
20#include <vector>
21
22namespace dolfinx::la
23{
25enum class BlockMode : int
26{
27 compact = 0,
32};
33
48template <class Scalar, class Container = std::vector<Scalar>,
49 class ColContainer = std::vector<std::int32_t>,
50 class RowPtrContainer = std::vector<std::int64_t>>
52{
53 static_assert(std::is_same_v<typename Container::value_type, Scalar>);
54
55public:
57 using value_type = Scalar;
58
60 using container_type = Container;
61
63 using column_container_type = ColContainer;
64
66 using rowptr_container_type = RowPtrContainer;
67
68 static_assert(std::is_same_v<value_type, typename container_type::value_type>,
69 "Scalar type and container value type must be the same.");
70
94 template <int BS0 = 1, int BS1 = 1>
96 {
97 if ((BS0 != _bs[0] and BS0 > 1 and _bs[0] > 1)
98 or (BS1 != _bs[1] and BS1 > 1 and _bs[1] > 1))
99 {
100 throw std::runtime_error(
101 "Cannot insert blocks of different size than matrix block size");
102 }
103
104 return [&](std::span<const std::int32_t> rows,
105 std::span<const std::int32_t> cols,
106 std::span<const value_type> data) -> int
107 {
108 this->set<BS0, BS1>(data, rows, cols);
109 return 0;
110 };
111 }
112
136 template <int BS0 = 1, int BS1 = 1>
138 {
139 if ((BS0 != _bs[0] and BS0 > 1 and _bs[0] > 1)
140 or (BS1 != _bs[1] and BS1 > 1 and _bs[1] > 1))
141 {
142 throw std::runtime_error(
143 "Cannot insert blocks of different size than matrix block size");
144 }
145
146 return [&](std::span<const std::int32_t> rows,
147 std::span<const std::int32_t> cols,
148 std::span<const value_type> data) -> int
149 {
150 this->add<BS0, BS1>(data, rows, cols);
151 return 0;
152 };
153 }
154
178 MatrixCSR(const SparsityPattern& p, BlockMode mode = BlockMode::compact);
179
182 MatrixCSR(MatrixCSR&& A) = default;
183
187 void set(value_type x) { std::ranges::fill(_data, x); }
188
205 template <int BS0, int BS1>
206 void set(std::span<const value_type> x, std::span<const std::int32_t> rows,
207 std::span<const std::int32_t> cols)
208 {
209 auto set_fn = [](value_type& y, const value_type& x) { y = x; };
210
211 std::int32_t num_rows
212 = _index_maps[0]->size_local() + _index_maps[0]->num_ghosts();
213 assert(x.size() == rows.size() * cols.size() * BS0 * BS1);
214 if (_bs[0] == BS0 and _bs[1] == BS1)
215 {
216 impl::insert_csr<BS0, BS1>(_data, _cols, _row_ptr, x, rows, cols, set_fn,
217 num_rows);
218 }
219 else if (_bs[0] == 1 and _bs[1] == 1)
220 {
221 // Set blocked data in a regular CSR matrix (_bs[0]=1, _bs[1]=1)
222 // with correct sparsity
223 impl::insert_blocked_csr<BS0, BS1>(_data, _cols, _row_ptr, x, rows, cols,
224 set_fn, num_rows);
225 }
226 else
227 {
228 assert(BS0 == 1 and BS1 == 1);
229 // Set non-blocked data in a blocked CSR matrix (BS0=1, BS1=1)
230 impl::insert_nonblocked_csr(_data, _cols, _row_ptr, x, rows, cols, set_fn,
231 num_rows, _bs[0], _bs[1]);
232 }
233 }
234
250 template <int BS0 = 1, int BS1 = 1>
251 void add(std::span<const value_type> x, std::span<const std::int32_t> rows,
252 std::span<const std::int32_t> cols)
253 {
254 auto add_fn = [](value_type& y, const value_type& x) { y += x; };
255
256 assert(x.size() == rows.size() * cols.size() * BS0 * BS1);
257 if (_bs[0] == BS0 and _bs[1] == BS1)
258 {
259 impl::insert_csr<BS0, BS1>(_data, _cols, _row_ptr, x, rows, cols, add_fn,
260 _row_ptr.size());
261 }
262 else if (_bs[0] == 1 and _bs[1] == 1)
263 {
264 // Add blocked data to a regular CSR matrix (_bs[0]=1, _bs[1]=1)
265 impl::insert_blocked_csr<BS0, BS1>(_data, _cols, _row_ptr, x, rows, cols,
266 add_fn, _row_ptr.size());
267 }
268 else
269 {
270 assert(BS0 == 1 and BS1 == 1);
271 // Add non-blocked data to a blocked CSR matrix (BS0=1, BS1=1)
272 impl::insert_nonblocked_csr(_data, _cols, _row_ptr, x, rows, cols, add_fn,
273 _row_ptr.size(), _bs[0], _bs[1]);
274 }
275 }
276
278 std::int32_t num_owned_rows() const { return _index_maps[0]->size_local(); }
279
281 std::int32_t num_all_rows() const { return _row_ptr.size() - 1; }
282
292 std::vector<value_type> to_dense() const
293 {
294 const std::size_t nrows = num_all_rows();
295 const std::size_t ncols = _index_maps[1]->size_global();
296 std::vector<value_type> A(nrows * ncols * _bs[0] * _bs[1], 0.0);
297 for (std::size_t r = 0; r < nrows; ++r)
298 {
299 for (std::int32_t j = _row_ptr[r]; j < _row_ptr[r + 1]; ++j)
300 {
301 for (int i0 = 0; i0 < _bs[0]; ++i0)
302 {
303 for (int i1 = 0; i1 < _bs[1]; ++i1)
304 {
305 std::array<std::int32_t, 1> local_col{_cols[j]};
306 std::array<std::int64_t, 1> global_col{0};
307 _index_maps[1]->local_to_global(local_col, global_col);
308 A[(r * _bs[1] + i0) * ncols * _bs[0] + global_col[0] * _bs[1] + i1]
309 = _data[j * _bs[0] * _bs[1] + i0 * _bs[1] + i1];
310 }
311 }
312 }
313 }
314
315 return A;
316 }
317
325 {
328 }
329
338 {
339 const std::int32_t local_size0 = _index_maps[0]->size_local();
340 const std::int32_t num_ghosts0 = _index_maps[0]->num_ghosts();
341 const int bs2 = _bs[0] * _bs[1];
342
343 // For each ghost row, pack and send values to send to neighborhood
344 std::vector<int> insert_pos = _val_send_disp;
345 _ghost_value_data.resize(_val_send_disp.back());
346 for (int i = 0; i < num_ghosts0; ++i)
347 {
348 const int rank = _ghost_row_to_rank[i];
349
350 // Get position in send buffer to place data to send to this
351 // neighbour
352 const std::int32_t val_pos = insert_pos[rank];
353 std::copy(std::next(_data.data(), _row_ptr[local_size0 + i] * bs2),
354 std::next(_data.data(), _row_ptr[local_size0 + i + 1] * bs2),
355 std::next(_ghost_value_data.begin(), val_pos));
356 insert_pos[rank]
357 += bs2 * (_row_ptr[local_size0 + i + 1] - _row_ptr[local_size0 + i]);
358 }
359
360 _ghost_value_data_in.resize(_val_recv_disp.back());
361
362 // Compute data sizes for send and receive from displacements
363 std::vector<int> val_send_count(_val_send_disp.size() - 1);
364 std::adjacent_difference(std::next(_val_send_disp.begin()),
365 _val_send_disp.end(), val_send_count.begin());
366
367 std::vector<int> val_recv_count(_val_recv_disp.size() - 1);
368 std::adjacent_difference(std::next(_val_recv_disp.begin()),
369 _val_recv_disp.end(), val_recv_count.begin());
370
371 int status = MPI_Ineighbor_alltoallv(
372 _ghost_value_data.data(), val_send_count.data(), _val_send_disp.data(),
373 dolfinx::MPI::mpi_t<value_type>, _ghost_value_data_in.data(),
374 val_recv_count.data(), _val_recv_disp.data(),
375 dolfinx::MPI::mpi_t<value_type>, _comm.comm(), &_request);
376 dolfinx::MPI::check_error(_comm.comm(), status);
377 }
378
385 {
386 int status = MPI_Wait(&_request, MPI_STATUS_IGNORE);
387 dolfinx::MPI::check_error(_comm.comm(), status);
388
389 _ghost_value_data.clear();
390 _ghost_value_data.shrink_to_fit();
391
392 // Add to local rows
393 const int bs2 = _bs[0] * _bs[1];
394 assert(_ghost_value_data_in.size() == _unpack_pos.size() * bs2);
395 for (std::size_t i = 0; i < _unpack_pos.size(); ++i)
396 for (int j = 0; j < bs2; ++j)
397 _data[_unpack_pos[i] * bs2 + j] += _ghost_value_data_in[i * bs2 + j];
398
399 _ghost_value_data_in.clear();
400 _ghost_value_data_in.shrink_to_fit();
401
402 // Set ghost row data to zero
403 const std::int32_t local_size0 = _index_maps[0]->size_local();
404 std::fill(std::next(_data.begin(), _row_ptr[local_size0] * bs2),
405 _data.end(), 0);
406 }
407
410 double squared_norm() const
411 {
412 const std::size_t num_owned_rows = _index_maps[0]->size_local();
413 const int bs2 = _bs[0] * _bs[1];
414 assert(num_owned_rows < _row_ptr.size());
415 double norm_sq_local = std::accumulate(
416 _data.cbegin(),
417 std::next(_data.cbegin(), _row_ptr[num_owned_rows] * bs2), double(0),
418 [](auto norm, value_type y) { return norm + std::norm(y); });
419 double norm_sq;
420 MPI_Allreduce(&norm_sq_local, &norm_sq, 1, MPI_DOUBLE, MPI_SUM,
421 _comm.comm());
422 return norm_sq;
423 }
424
433
441 std::shared_ptr<const common::IndexMap> index_map(int dim) const
442 {
443 return _index_maps.at(dim);
444 }
445
448 container_type& values() { return _data; }
449
452 const container_type& values() const { return _data; }
453
456 const rowptr_container_type& row_ptr() const { return _row_ptr; }
457
460 const column_container_type& cols() const { return _cols; }
461
470 {
471 return _off_diagonal_offset;
472 }
473
476 std::array<int, 2> block_size() const { return _bs; }
477
478private:
479 // Maps for the distribution of the ows and columns
480 std::array<std::shared_ptr<const common::IndexMap>, 2> _index_maps;
481
482 // Block mode (compact or expanded)
483 BlockMode _block_mode;
484
485 // Block sizes
486 std::array<int, 2> _bs;
487
488 // Matrix data
489 container_type _data;
491 rowptr_container_type _row_ptr;
492
493 // Start of off-diagonal (unowned columns) on each row
494 rowptr_container_type _off_diagonal_offset;
495
496 // Neighborhood communicator (ghost->owner communicator for rows)
497 dolfinx::MPI::Comm _comm;
498
499 // -- Precomputed data for scatter_rev/update
500
501 // Request in non-blocking communication
502 MPI_Request _request;
503
504 // Position in _data to add received data
505 std::vector<int> _unpack_pos;
506
507 // Displacements for alltoall for each neighbor when sending and
508 // receiving
509 std::vector<int> _val_send_disp, _val_recv_disp;
510
511 // Ownership of each row, by neighbor (for the neighbourhood defined
512 // on _comm)
513 std::vector<int> _ghost_row_to_rank;
514
515 // Temporary stores for data during non-blocking communication
516 container_type _ghost_value_data;
517 container_type _ghost_value_data_in;
518};
519//-----------------------------------------------------------------------------
520template <class U, class V, class W, class X>
522 : _index_maps({p.index_map(0),
523 std::make_shared<common::IndexMap>(p.column_index_map())}),
524 _block_mode(mode), _bs({p.block_size(0), p.block_size(1)}),
525 _data(p.num_nonzeros() * _bs[0] * _bs[1], 0),
526 _cols(p.graph().first.begin(), p.graph().first.end()),
527 _row_ptr(p.graph().second.begin(), p.graph().second.end()),
528 _comm(MPI_COMM_NULL)
529{
530 if (_block_mode == BlockMode::expanded)
531 {
532 // Rebuild IndexMaps
533 for (int i = 0; i < 2; ++i)
534 {
535 const auto im = _index_maps[i];
536 const std::int32_t size_local = im->size_local() * _bs[i];
537 std::span ghost_i = im->ghosts();
538 std::vector<std::int64_t> ghosts;
539 const std::vector<int> ghost_owner_i(im->owners().begin(),
540 im->owners().end());
541 std::vector<int> src_rank;
542 for (std::size_t j = 0; j < ghost_i.size(); ++j)
543 {
544 for (int k = 0; k < _bs[i]; ++k)
545 {
546 ghosts.push_back(ghost_i[j] * _bs[i] + k);
547 src_rank.push_back(ghost_owner_i[j]);
548 }
549 }
550
551 std::array<std::vector<int>, 2> src_dest0
552 = {std::vector(_index_maps[i]->src().begin(),
553 _index_maps[i]->src().end()),
554 std::vector(_index_maps[i]->dest().begin(),
555 _index_maps[i]->dest().end())};
556 _index_maps[i] = std::make_shared<common::IndexMap>(
557 _index_maps[i]->comm(), size_local, src_dest0, ghosts, src_rank);
558 }
559
560 // Convert sparsity pattern and set _bs to 1
561
562 column_container_type new_cols;
563 new_cols.reserve(_data.size());
564 rowptr_container_type new_row_ptr = {0};
565 new_row_ptr.reserve(_row_ptr.size() * _bs[0]);
566 std::span<const std::int32_t> num_diag_nnz = p.off_diagonal_offsets();
567 for (std::size_t i = 0; i < _row_ptr.size() - 1; ++i)
568 {
569 // Repeat row _bs[0] times
570 for (int q0 = 0; q0 < _bs[0]; ++q0)
571 {
572 _off_diagonal_offset.push_back(new_row_ptr.back()
573 + num_diag_nnz[i] * _bs[1]);
574 for (auto j = _row_ptr[i]; j < _row_ptr[i + 1]; ++j)
575 {
576 for (int q1 = 0; q1 < _bs[1]; ++q1)
577 new_cols.push_back(_cols[j] * _bs[1] + q1);
578 }
579 new_row_ptr.push_back(new_cols.size());
580 }
581 }
582 _cols = new_cols;
583 _row_ptr = new_row_ptr;
584 _bs[0] = 1;
585 _bs[1] = 1;
586 }
587 else
588 {
589 // Compute off-diagonal offset for each row (compact)
590 std::span<const std::int32_t> num_diag_nnz = p.off_diagonal_offsets();
591 _off_diagonal_offset.reserve(num_diag_nnz.size());
592 std::ranges::transform(num_diag_nnz, _row_ptr,
593 std::back_inserter(_off_diagonal_offset),
594 std::plus{});
595 }
596
597 // Some short-hand
598 const std::array local_size
599 = {_index_maps[0]->size_local(), _index_maps[1]->size_local()};
600 const std::array local_range
601 = {_index_maps[0]->local_range(), _index_maps[1]->local_range()};
602 std::span ghosts1 = _index_maps[1]->ghosts();
603
604 std::span ghosts0 = _index_maps[0]->ghosts();
605 std::span src_ranks = _index_maps[0]->src();
606 std::span dest_ranks = _index_maps[0]->dest();
607
608 // Create neighbourhood communicator (owner <- ghost)
609 MPI_Comm comm;
610 MPI_Dist_graph_create_adjacent(_index_maps[0]->comm(), dest_ranks.size(),
611 dest_ranks.data(), MPI_UNWEIGHTED,
612 src_ranks.size(), src_ranks.data(),
613 MPI_UNWEIGHTED, MPI_INFO_NULL, false, &comm);
614 _comm = dolfinx::MPI::Comm(comm, false);
615
616 // Build map from ghost row index position to owning (neighborhood)
617 // rank
618 _ghost_row_to_rank.reserve(_index_maps[0]->owners().size());
619 for (int r : _index_maps[0]->owners())
620 {
621 auto it = std::ranges::lower_bound(src_ranks, r);
622 assert(it != src_ranks.end() and *it == r);
623 std::size_t pos = std::distance(src_ranks.begin(), it);
624 _ghost_row_to_rank.push_back(pos);
625 }
626
627 // Compute size of data to send to each neighbor
628 std::vector<std::int32_t> data_per_proc(src_ranks.size(), 0);
629 for (std::size_t i = 0; i < _ghost_row_to_rank.size(); ++i)
630 {
631 assert(_ghost_row_to_rank[i] < (int)data_per_proc.size());
632 std::size_t pos = local_size[0] + i;
633 data_per_proc[_ghost_row_to_rank[i]] += _row_ptr[pos + 1] - _row_ptr[pos];
634 }
635
636 // Compute send displacements
637 _val_send_disp.resize(src_ranks.size() + 1, 0);
638 std::partial_sum(data_per_proc.begin(), data_per_proc.end(),
639 std::next(_val_send_disp.begin()));
640
641 // For each ghost row, pack and send indices to neighborhood
642 std::vector<std::int64_t> ghost_index_data(2 * _val_send_disp.back());
643 {
644 std::vector<int> insert_pos = _val_send_disp;
645 for (std::size_t i = 0; i < _ghost_row_to_rank.size(); ++i)
646 {
647 const int rank = _ghost_row_to_rank[i];
648 std::int32_t row_id = local_size[0] + i;
649 for (int j = _row_ptr[row_id]; j < _row_ptr[row_id + 1]; ++j)
650 {
651 // Get position in send buffer
652 const std::int32_t idx_pos = 2 * insert_pos[rank];
653
654 // Pack send data (row, col) as global indices
655 ghost_index_data[idx_pos] = ghosts0[i];
656 if (std::int32_t col_local = _cols[j]; col_local < local_size[1])
657 ghost_index_data[idx_pos + 1] = col_local + local_range[1][0];
658 else
659 ghost_index_data[idx_pos + 1] = ghosts1[col_local - local_size[1]];
660
661 insert_pos[rank] += 1;
662 }
663 }
664 }
665
666 // Communicate data with neighborhood
667 std::vector<std::int64_t> ghost_index_array;
668 std::vector<int> recv_disp;
669 {
670 std::vector<int> send_sizes;
671 std::ranges::transform(data_per_proc, std::back_inserter(send_sizes),
672 [](auto x) { return 2 * x; });
673
674 std::vector<int> recv_sizes(dest_ranks.size());
675 send_sizes.reserve(1);
676 recv_sizes.reserve(1);
677 MPI_Neighbor_alltoall(send_sizes.data(), 1, MPI_INT, recv_sizes.data(), 1,
678 MPI_INT, _comm.comm());
679
680 // Build send/recv displacement
681 std::vector<int> send_disp = {0};
682 std::partial_sum(send_sizes.begin(), send_sizes.end(),
683 std::back_inserter(send_disp));
684 recv_disp = {0};
685 std::partial_sum(recv_sizes.begin(), recv_sizes.end(),
686 std::back_inserter(recv_disp));
687
688 ghost_index_array.resize(recv_disp.back());
689 MPI_Neighbor_alltoallv(ghost_index_data.data(), send_sizes.data(),
690 send_disp.data(), MPI_INT64_T,
691 ghost_index_array.data(), recv_sizes.data(),
692 recv_disp.data(), MPI_INT64_T, _comm.comm());
693 }
694
695 // Store receive displacements for future use, when transferring
696 // data values
697 _val_recv_disp.resize(recv_disp.size());
698 const int bs2 = _bs[0] * _bs[1];
699 std::ranges::transform(recv_disp, _val_recv_disp.begin(),
700 [&bs2](auto d) { return bs2 * d / 2; });
701 std::ranges::transform(_val_send_disp, _val_send_disp.begin(),
702 [&bs2](auto d) { return d * bs2; });
703
704 // Global-to-local map for ghost columns
705 std::vector<std::pair<std::int64_t, std::int32_t>> global_to_local;
706 global_to_local.reserve(ghosts1.size());
707 for (std::int64_t idx : ghosts1)
708 global_to_local.push_back({idx, global_to_local.size() + local_size[1]});
709 std::ranges::sort(global_to_local);
710
711 // Compute location in which data for each index should be stored
712 // when received
713 for (std::size_t i = 0; i < ghost_index_array.size(); i += 2)
714 {
715 // Row must be on this process
716 const std::int32_t local_row = ghost_index_array[i] - local_range[0][0];
717 assert(local_row >= 0 and local_row < local_size[0]);
718
719 // Column may be owned or unowned
720 std::int32_t local_col = ghost_index_array[i + 1] - local_range[1][0];
721 if (local_col < 0 or local_col >= local_size[1])
722 {
723 auto it = std::ranges::lower_bound(
724 global_to_local, std::pair(ghost_index_array[i + 1], -1),
725 [](auto& a, auto& b) { return a.first < b.first; });
726 assert(it != global_to_local.end()
727 and it->first == ghost_index_array[i + 1]);
728 local_col = it->second;
729 }
730 auto cit0 = std::next(_cols.begin(), _row_ptr[local_row]);
731 auto cit1 = std::next(_cols.begin(), _row_ptr[local_row + 1]);
732
733 // Find position of column index and insert data
734 auto cit = std::lower_bound(cit0, cit1, local_col);
735 assert(cit != cit1);
736 assert(*cit == local_col);
737 std::size_t d = std::distance(_cols.begin(), cit);
738 _unpack_pos.push_back(d);
739 }
740}
741//-----------------------------------------------------------------------------
742
743// The matrix A is distributed across P processes by blocks of rows:
744// A = | A_0 |
745// | A_1 |
746// | ... |
747// | A_P-1 |
748//
749// Each submatrix A_i is owned by a single process "i" and can be further
750// decomposed into diagonal (Ai[0]) and off diagonal (Ai[1]) blocks:
751// Ai = |Ai[0] Ai[1]|
752//
753// If A is square, the diagonal block Ai[0] is also square and contains
754// only owned columns and rows. The block Ai[1] contains ghost columns
755// (unowned dofs).
756
757// Likewise, a local vector x can be decomposed into owned and ghost blocks:
758// xi = | x[0] |
759// | x[1] |
760//
761// So the product y = Ax can be computed into two separate steps:
762// y[0] = |Ai[0] Ai[1]| | x[0] | = Ai[0] x[0] + Ai[1] x[1]
763// | x[1] |
764//
767template <typename Scalar, typename V, typename W, typename X>
770{
771 // start communication (update ghosts)
773
774 const std::int32_t nrowslocal = num_owned_rows();
775 std::span<const std::int64_t> Arow_ptr(row_ptr().data(), nrowslocal + 1);
776 std::span<const std::int32_t> Acols(cols().data(), Arow_ptr[nrowslocal]);
777 std::span<const std::int64_t> Aoff_diag_offset(off_diag_offset().data(),
778 nrowslocal);
779 std::span<const Scalar> Avalues(values().data(), Arow_ptr[nrowslocal]);
780
781 std::span<const Scalar> _x = x.array();
782 std::span<Scalar> _y = y.mutable_array();
783
784 std::span<const std::int64_t> Arow_begin(Arow_ptr.data(), nrowslocal);
785 std::span<const std::int64_t> Arow_end(Arow_ptr.data() + 1, nrowslocal);
786
787 // First stage: spmv - diagonal
788 // yi[0] += Ai[0] * xi[0]
789 if (_bs[1] == 1)
790 {
791 impl::spmv<Scalar, 1>(Avalues, Arow_begin, Aoff_diag_offset, Acols, _x, _y,
792 _bs[0], 1);
793 }
794 else
795 {
796 impl::spmv<Scalar, -1>(Avalues, Arow_begin, Aoff_diag_offset, Acols, _x, _y,
797 _bs[0], _bs[1]);
798 }
799
800 // finalize ghost update
801 x.scatter_fwd_end();
802
803 // Second stage: spmv - off-diagonal
804 // yi[0] += Ai[1] * xi[1]
805 if (_bs[1] == 1)
806 {
807 impl::spmv<Scalar, 1>(Avalues, Aoff_diag_offset, Arow_end, Acols, _x, _y,
808 _bs[0], 1);
809 }
810 else
811 {
812 impl::spmv<Scalar, -1>(Avalues, Aoff_diag_offset, Arow_end, Acols, _x, _y,
813 _bs[0], _bs[1]);
814 }
815}
816
817} // namespace dolfinx::la
A duplicate MPI communicator and manage lifetime of the communicator.
Definition MPI.h:43
const container_type & values() const
Definition MatrixCSR.h:452
std::shared_ptr< const common::IndexMap > index_map(int dim) const
Index maps for the row and column space.
Definition MatrixCSR.h:441
const rowptr_container_type & off_diag_offset() const
Definition MatrixCSR.h:469
void set(std::span< const value_type > x, std::span< const std::int32_t > rows, std::span< const std::int32_t > cols)
Set values in the matrix.
Definition MatrixCSR.h:206
RowPtrContainer rowptr_container_type
Row pointer container type.
Definition MatrixCSR.h:66
void scatter_rev_end()
End transfer of ghost row data to owning ranks.
Definition MatrixCSR.h:384
container_type & values()
Definition MatrixCSR.h:448
auto mat_add_values()
Insertion functor for adding values to a matrix. It is typically used in finite element assembly func...
Definition MatrixCSR.h:137
void add(std::span< const value_type > x, std::span< const std::int32_t > rows, std::span< const std::int32_t > cols)
Accumulate values in the matrix.
Definition MatrixCSR.h:251
std::int32_t num_owned_rows() const
Number of local rows excluding ghost rows.
Definition MatrixCSR.h:278
MatrixCSR(const SparsityPattern &p, BlockMode mode=BlockMode::compact)
Create a distributed matrix.
Definition MatrixCSR.h:521
ColContainer column_container_type
Column index container type.
Definition MatrixCSR.h:63
MatrixCSR(MatrixCSR &&A)=default
void mult(Vector< value_type > &x, Vector< value_type > &y)
Compute the product y += Ax.
Definition MatrixCSR.h:768
double squared_norm() const
Compute the Frobenius norm squared across all processes.
Definition MatrixCSR.h:410
void scatter_rev()
Transfer ghost row data to the owning ranks accumulating received values on the owned rows,...
Definition MatrixCSR.h:324
Container container_type
Matrix entries container type.
Definition MatrixCSR.h:60
Scalar value_type
Scalar type.
Definition MatrixCSR.h:57
void scatter_rev_begin()
Begin transfer of ghost row data to owning ranks, where it will be accumulated into existing owned ro...
Definition MatrixCSR.h:337
const column_container_type & cols() const
Definition MatrixCSR.h:460
void set(value_type x)
Set all non-zero local entries to a value including entries in ghost rows.
Definition MatrixCSR.h:187
std::array< int, 2 > block_size() const
Definition MatrixCSR.h:476
std::int32_t num_all_rows() const
Number of local rows including ghost rows.
Definition MatrixCSR.h:281
const rowptr_container_type & row_ptr() const
Definition MatrixCSR.h:456
std::vector< value_type > to_dense() const
Copy to a dense matrix.
Definition MatrixCSR.h:292
auto mat_set_values()
Insertion functor for setting values in a matrix. It is typically used in finite element assembly fun...
Definition MatrixCSR.h:95
Definition SparsityPattern.h:26
Definition Vector.h:32
void scatter_fwd_begin()
Definition Vector.h:86
std::span< const value_type > array() const
Get local part of the vector (const version)
Definition Vector.h:190
void scatter_fwd_end()
Definition Vector.h:105
std::span< value_type > mutable_array()
Get local part of the vector.
Definition Vector.h:196
MPI_Datatype mpi_t
Retrieves the MPI data type associated to the provided type.
Definition MPI.h:281
void check_error(MPI_Comm comm, int code)
Check MPI error code. If the error code is not equal to MPI_SUCCESS, then std::abort is called.
Definition MPI.cpp:80
int size(MPI_Comm comm)
Definition MPI.cpp:72
int rank(MPI_Comm comm)
Return process rank for the communicator.
Definition MPI.cpp:64
constexpr std::array< std::int64_t, 2 > local_range(int rank, std::int64_t N, int size)
Return local range for the calling process, partitioning the global [0, N - 1] range across all ranks...
Definition MPI.h:90
Linear algebra interface.
Definition sparsitybuild.h:15
BlockMode
Modes for representing block structured matrices.
Definition MatrixCSR.h:26
@ expanded
Definition MatrixCSR.h:29
auto norm(const V &x, Norm type=Norm::l2)
Definition Vector.h:268