ryujin 2.1.1 revision ee5cbcbf2346c1299c942d0e1f13b46449973c18
Loading...
Searching...
No Matches
sparse_matrix.h
Go to the documentation of this file.
1//
2// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
3// Copyright (C) 2020 - 2026 by the ryujin authors
4//
5
6#pragma once
7
8#include <compile_time_options.h>
9
10#include "gpu.h"
11#include "sparsity_pattern.h"
12
13#include <deal.II/base/exceptions.h>
14#include <deal.II/base/partitioner.h>
15#include <deal.II/lac/affine_constraints.h>
16#include <deal.II/lac/dynamic_sparsity_pattern.h>
17
18#include <type_traits>
19
20// #define DEBUG_MPI_EXCHANGE
21
22#ifdef DEBUG_MPI_EXCHANGE
23#include <chrono>
24#endif
25
26namespace ryujin
27{
28 template <typename Number,
29 int n_comp,
30 int warp_size,
31 int simd_length,
32 typename MemorySpace = dealii::MemorySpace::Host,
33 bool writable = true>
34 class SparseMatrixView;
35
54 template <typename Number,
55 int n_comp = 1,
56 int warp_size = ryujin::warp_size,
57 int simd_length = dealii::VectorizedArray<Number>::size()>
58 class SparseMatrix : public MirroredStorage<
59 SparseMatrix<Number, n_comp, warp_size, simd_length>>
60 {
61 public:
62 static_assert(warp_size % simd_length == 0,
63 "The warp size must be an integer multiple of the SIMD "
64 "length");
65
70
74 SparseMatrix() = default;
75
79 SparseMatrix(const SparsityPattern<warp_size> &sparsity,
80 const TransferPolicy transfer_policy =
81 TransferPolicy::explicit_transfers);
82
97 void reinit(const SparsityPattern<warp_size> &sparsity,
98 const TransferPolicy transfer_policy =
99 TransferPolicy::explicit_transfers);
100
104 ACCESSOR_READ_ONLY(sparsity_pattern);
105
107
111
120 template <typename MemorySpace = dealii::MemorySpace::Host>
121 SparseMatrixView<Number, n_comp, warp_size, simd_length, MemorySpace, true>
122 view();
123
131 template <typename MemorySpace = dealii::MemorySpace::Host>
132 SparseMatrixView<Number, n_comp, warp_size, simd_length, MemorySpace, false>
133 view() const;
134
135 /*
136 * The is_resident(), copy_to_memory_space(), move_to_memory_space(),
137 * transfer_policy(), and set_transfer_policy() methods are inherited
138 * from the MirroredStorage base class.
139 */
140
142
146
150 template <typename MemorySpace>
151 void zero_out_ghost_rows_on_memory_space();
152
157 template <typename MemorySpace>
158 void update_ghost_rows_on_memory_space();
159
165 template <typename MemorySpace>
166 void compress_on_memory_space(dealii::VectorOperation::values operation);
167
177 template <typename MemorySpace>
178 void populate_exchange_buffer_on_memory_space();
179
180 private:
182
186
187 const SparsityPattern<warp_size> *sparsity_pattern_ =
188 nullptr; // FIXME shared_ptr
189
190 /*
191 * Note: The storage is marked mutable so that the (logically const)
192 * copy_to_memory_space() operation can populate a mirror from within
193 * a const view() under the implicit_transfers policy.
194 */
195
196 using KokkosHost = dealii::MemorySpace::Host::kokkos_space;
197 mutable Kokkos::View<Number *, KokkosHost> data_host_;
198
199 using KokkosDefault = dealii::MemorySpace::Default::kokkos_space;
200 mutable Kokkos::View<Number *, KokkosDefault> data_default_;
201
202 /*
203 * Note: For the time being, let's assume MPI only operates on host
204 * memory. Meaning, when we populate exchange_buffer_default_ on the
205 * device, we will copy over to the host prior to sending.
206 */
207
208 mutable Kokkos::View<Number *, KokkosHost> exchange_buffer_host_;
209 mutable Kokkos::View<Number *, KokkosHost> ghost_buffer_host_;
210 mutable Kokkos::View<Number *, KokkosDefault> exchange_buffer_default_;
211
212 std::vector<MPI_Request> requests_;
213
214 /*
215 * Storage primitives used by the MirroredStorage base class:
216 */
217
218 template <typename MemorySpace>
219 void allocate_storage() const;
220
221 template <typename To, typename From>
222 void deep_copy_storage() const;
223
224 template <typename MemorySpace>
225 void deallocate_storage();
226
227
228 friend class MirroredStorage<
229 SparseMatrix<Number, n_comp, warp_size, simd_length>>;
230
231 template <typename, int, int, int, typename, bool>
232 friend class SparseMatrixView;
233
235 };
236
237
252 template <typename Number,
253 int n_comp,
254 int warp_size,
255 int simd_length,
256 typename MemorySpace,
257 bool writable>
258 class SparseMatrixView
259 {
260 public:
261 static_assert(warp_size % simd_length == 0,
262 "The warp size must be an integer multiple of the SIMD "
263 "length");
264
269
270 SparseMatrixView() = default;
271
272 SparseMatrixView(
273 SparseMatrix<Number, n_comp, warp_size, simd_length> &sparse_matrix)
274 requires(writable);
275
276 SparseMatrixView(const SparseMatrix<Number, n_comp, warp_size, simd_length>
277 &sparse_matrix)
278 requires(!writable);
279
288 template <bool other_writable>
289 DEAL_II_HOST_DEVICE
290 SparseMatrixView(const SparseMatrixView<Number,
291 n_comp,
292 warp_size,
293 simd_length,
294 MemorySpace,
295 other_writable> &other)
296 requires(!writable && other_writable);
297
298 template <typename SparseMatrix>
299 void reinit(SparseMatrix &sparse_matrix)
300 requires(writable != std::is_const_v<SparseMatrix>);
301
305 ACCESSOR_READ_ONLY(sparsity_pattern);
306
308
312
313 /* Get scalar or tensor-valued entry: */
314
325 template <typename Number2 = Number>
326 DEAL_II_HOST_DEVICE Number2
327 read_entry(const unsigned int row, const unsigned int column_index) const;
328
340 template <typename Number2 = Number,
341 typename Tensor = dealii::Tensor<1, n_comp, Number2>>
342 DEAL_II_HOST_DEVICE Tensor
343 read_tensor(const unsigned int row, const unsigned int column_index) const;
344
345 /* Get transposed scalar or tensor-valued entry: */
346
358 template <typename Number2 = Number>
359 DEAL_II_HOST_DEVICE Number2 read_transposed_entry(
360 const unsigned int row, const unsigned int column_index) const;
361
373 template <typename Number2 = Number,
374 typename Tensor = dealii::Tensor<1, n_comp, Number2>>
375 DEAL_II_HOST_DEVICE Tensor read_transposed_tensor(
376 const unsigned int row, const unsigned int column_index) const;
377
378 /* Write scalar or tensor entry: */
379
391 template <typename Number2 = Number>
392 DEAL_II_HOST_DEVICE void
393 write_entry(const Number2 entry,
394 const unsigned int row,
395 const unsigned int column_index,
396 const bool do_streaming_store = false) const
397 requires(writable);
398
408 template <typename Number2 = Number,
409 typename Tensor = dealii::Tensor<1, n_comp, Number2>>
410 DEAL_II_HOST_DEVICE void
411 write_tensor(const Tensor &tensor,
412 const unsigned int row,
413 const unsigned int column_index,
414 const bool do_streaming_store = false) const
415 requires(writable);
416
428 template <typename Number2 = Number>
429 DEAL_II_HOST_DEVICE void add_entry(const Number2 entry,
430 const unsigned int row,
431 const unsigned int column_index) const
432 requires(writable);
433
443 template <typename Number2 = Number,
444 typename Tensor = dealii::Tensor<1, n_comp, Number2>>
445 DEAL_II_HOST_DEVICE void add_tensor(const Tensor &tensor,
446 const unsigned int row,
447 const unsigned int column_index) const
448 requires(writable);
449
451
455
456 void zero_out_ghost_rows() const
457 requires(writable);
458
459 void update_ghost_rows() const
460 requires(writable);
461
462 void compress(dealii::VectorOperation::values operation) const
463 requires(writable);
464
465 private:
467
471
472 using SM = SparseMatrix<Number, n_comp, warp_size, simd_length>;
473 std::conditional_t<writable, SM *, const SM *> sparse_matrix_;
474
475 SparsityPatternView<warp_size, MemorySpace> sparsity_pattern_;
476
477 using KokkosSpace = typename MemorySpace::kokkos_space;
478 Kokkos::View<Number *, KokkosSpace> data_;
479
480 template <typename, int, int, int, typename, bool>
481 friend class SparseMatrixView;
482
484 };
485
486
487 /*
488 * Given a matrix contribution and row and column indices in (deal.II
489 * typical) global numbering, add the contribution to the sparse matrix.
490 * The function takes an @p affine_constraints object in (deal.II
491 * typical) global numbering and resolves constrained degrees of freedom
492 * prior to distributing to the matrix.
493 *
494 * @note The method will not modify the diagonal entry of constrained
495 * degrees of freedom, in contrast to the deal.II version
496 * AffineConstraints<Number>::distribute_local_to_global().
497 *
498 * @note For a vector-valued matrix with n_comp > 1 the @p cell_matrix
499 * must be a container with a subscript operator[] returning a matrix for
500 * each component.
501 */
502 template <typename Number,
503 int n_comp,
504 int warp_size,
505 int simd_length,
506 typename FullMatrix>
507 void distribute_local_to_global(
508 const FullMatrix &cell_matrix,
509 const std::vector<dealii::types::global_dof_index> &dof_indices_row,
510 const std::vector<dealii::types::global_dof_index> &dof_indices_column,
511 const dealii::AffineConstraints<Number> &affine_constraints,
512 SparseMatrix<Number, n_comp, warp_size, simd_length> &sparse_matrix);
513
514
515 /*
516 * Variant of the function above that takes a symmetric matrix
517 * constribution where the column and row indices are the same.
518 */
519 template <typename Number,
520 int n_comp,
521 int warp_size,
522 int simd_length,
523 typename FullMatrix,
524 typename Vector>
525 void distribute_local_to_global(
526 const FullMatrix &cell_matrix,
527 const std::vector<dealii::types::global_dof_index> &dof_indices,
528 const dealii::AffineConstraints<Number> &affine_constraints,
529 SparseMatrix<Number, n_comp, warp_size, simd_length> &sparse_matrix);
530
531
532#ifndef DOXYGEN
533 /*
534 * -------------------------------------------------------------------------
535 * Inline function definitions
536 * -------------------------------------------------------------------------
537 */
538
539
540 template <typename Number, int n_components, int warp_size, int simd_length>
541 SparseMatrix<Number, n_components, warp_size, simd_length>::SparseMatrix(
542 const SparsityPattern<warp_size> &sparsity,
543 const TransferPolicy transfer_policy)
544 {
545 reinit(sparsity, transfer_policy);
546 }
547
548
549 template <typename Number, int n_components, int warp_size, int simd_length>
550 void SparseMatrix<Number, n_components, warp_size, simd_length>::reinit(
551 const SparsityPattern<warp_size> &sparsity,
552 const TransferPolicy transfer_policy)
553 {
554 /*
555 * Drop the transfer policy for the duration of the reinit so that a
556 * pinned memory space of a previous policy does not interfere:
557 */
558 this->set_transfer_policy(TransferPolicy::explicit_transfers);
559
560 this->sparsity_pattern_ = &sparsity;
561
562 const auto sparsity_view = sparsity.view();
563
564 using KokkosHost = dealii::MemorySpace::Host::kokkos_space;
565 using Aligned = Kokkos::MemoryTraits<Kokkos::Aligned>;
566
567 data_host_ = Kokkos::View<Number *, KokkosHost, Aligned>(
568 "sparse_matrix_data",
569 sparsity_view.n_nonzero_elements() * n_components);
570
571 const std::size_t n_indices = sparsity.entries_to_be_sent().size();
572
573 exchange_buffer_host_ = Kokkos::View<Number *, KokkosHost, Aligned>(
574 "sparse_matrix_exchange_buffer", n_components * n_indices);
575
576 if constexpr (have_separate_memory_spaces) {
577 const auto ghost_offset =
578 sparsity_view.template ghost_offset<n_components>();
579 const auto end_offset = sparsity_view.n_nonzero_elements() * n_components;
580
581 ghost_buffer_host_ = Kokkos::View<Number *, KokkosHost, Aligned>(
582 "sparse_matrix_ghost_buffer", end_offset - ghost_offset);
583 }
584
585 /*
586 * The matrix is resident on the host memory space only. Device
587 * storage is allocated lazily on the first copy_to_memory_space() /
588 * move_to_memory_space(); drop possibly stale device storage from a
589 * previous reinit():
590 */
591 data_default_ = {};
592 exchange_buffer_default_ = {};
593
594 this->reset_residency(/*host*/ true,
595 /*default*/ !have_separate_memory_spaces);
596
597 this->set_transfer_policy(transfer_policy);
598 }
599
600
601 template <typename Number, int n_comp, int warp_size, int simd_length>
602 template <typename MemorySpace>
603 SparseMatrixView<Number, n_comp, warp_size, simd_length, MemorySpace, true>
604 SparseMatrix<Number, n_comp, warp_size, simd_length>::view()
605 {
606 this->template prepare_write_access<MemorySpace>();
607
608 return SparseMatrixView<Number,
609 n_comp,
610 warp_size,
611 simd_length,
612 MemorySpace,
613 true>(*this);
614 }
615
616
617 template <typename Number, int n_comp, int warp_size, int simd_length>
618 template <typename MemorySpace>
619 SparseMatrixView<Number, n_comp, warp_size, simd_length, MemorySpace, false>
620 SparseMatrix<Number, n_comp, warp_size, simd_length>::view() const
621 {
622 this->template prepare_read_access<MemorySpace>();
623
624 return SparseMatrixView<Number,
625 n_comp,
626 warp_size,
627 simd_length,
628 MemorySpace,
629 false>(*this);
630 }
631
632
633 template <typename Number, int n_components, int warp_size, int simd_length>
634 template <typename MemorySpace>
635 void
636 SparseMatrix<Number, n_components, warp_size, simd_length>::allocate_storage()
637 const
638 {
639 using HostSpace = dealii::MemorySpace::Host;
640 using Aligned = Kokkos::MemoryTraits<Kokkos::Aligned>;
641
642 Assert(sparsity_pattern_ != nullptr, dealii::ExcNotInitialized());
643
644 const auto sparsity_view = sparsity_pattern_->view();
645
646 /*
647 * Note: We allocate without initializing because a deep_copy_storage()
648 * always follows.
649 */
650
651 const std::size_t n_data =
652 sparsity_view.n_nonzero_elements() * n_components;
653 const std::size_t n_exchange =
654 n_components * sparsity_pattern_->entries_to_be_sent().size();
655
656 if constexpr (std::is_same_v<MemorySpace, HostSpace>) {
657 data_host_ = Kokkos::View<Number *, KokkosHost, Aligned>(
658 Kokkos::view_alloc(Kokkos::WithoutInitializing, "sparse_matrix_data"),
659 n_data);
660
661 } else {
662 data_default_ = Kokkos::View<Number *, KokkosDefault>(
663 Kokkos::view_alloc(Kokkos::WithoutInitializing, "sparse_matrix_data"),
664 n_data);
665
666 exchange_buffer_default_ = Kokkos::View<Number *, KokkosDefault>(
667 Kokkos::view_alloc(Kokkos::WithoutInitializing,
668 "sparse_matrix_exchange_buffer"),
669 n_exchange);
670 }
671 }
672
673
674 template <typename Number, int n_components, int warp_size, int simd_length>
675 template <typename To, typename From>
676 void SparseMatrix<Number, n_components, warp_size, simd_length>::
677 deep_copy_storage() const
678 {
679 using HostSpace = dealii::MemorySpace::Host;
680
681 if constexpr (std::is_same_v<To, HostSpace>) {
682 Kokkos::deep_copy(/*dst*/ data_host_, /*src*/ data_default_);
683 } else {
684 Kokkos::deep_copy(/*dst*/ data_default_, /*src*/ data_host_);
685 }
686 }
687
688
689 template <typename Number, int n_components, int warp_size, int simd_length>
690 template <typename MemorySpace>
691 void SparseMatrix<Number, n_components, warp_size, simd_length>::
692 deallocate_storage()
693 {
694 using HostSpace = dealii::MemorySpace::Host;
695
696 if constexpr (std::is_same_v<MemorySpace, HostSpace>) {
697 data_host_ = {};
698
699 } else {
700 data_default_ = {};
701 exchange_buffer_default_ = {};
702 }
703 }
704
705
706 template <typename Number, int n_components, int warp_size, int simd_length>
707 template <typename MemorySpace>
708 void SparseMatrix<Number, n_components, warp_size, simd_length>::
709 zero_out_ghost_rows_on_memory_space()
710 {
711 using HostSpace = dealii::MemorySpace::Host;
712 using DefaultSpace = dealii::MemorySpace::Default;
713
714 Assert(this->template is_resident<MemorySpace>(),
715 dealii::ExcMessage("The chosen memory space is not resident."));
716
717 AssertThrow((std::is_same_v<MemorySpace, HostSpace>),
718 dealii::ExcNotImplemented());
719
720 const auto sparsity_view = sparsity_pattern_->view();
721
722 const auto ghost_offset =
723 sparsity_view.template ghost_offset<n_components>();
724 const auto end_offset = sparsity_view.n_nonzero_elements() * n_components;
725 std::fill(data_host_.data() + ghost_offset,
726 data_host_.data() + end_offset,
727 Number{});
728 }
729
730
731 template <typename Number, int n_components, int warp_size, int simd_length>
732 template <typename MemorySpace>
733 void SparseMatrix<Number, n_components, warp_size, simd_length>::
734 populate_exchange_buffer_on_memory_space()
735 {
736 using HostSpace = dealii::MemorySpace::Host;
737
738 const auto sparsity_view = sparsity_pattern_->template view<MemorySpace>();
739
740 const auto *entries_to_be_sent =
741 sparsity_pattern_->entries_to_be_sent().template view<MemorySpace>();
742 const std::size_t n_entries_to_be_sent =
743 sparsity_pattern_->entries_to_be_sent().size();
744
745 /*
746 * Note: If the host and default memory spaces coincide all views
747 * reference the host storage.
748 */
749 constexpr bool on_host =
750 !have_separate_memory_spaces || std::is_same_v<MemorySpace, HostSpace>;
751
752 const auto data = [&]() {
753 if constexpr (on_host)
754 return data_host_;
755 else
756 return data_default_;
757 }();
758
759 const auto exchange_buffer = [&]() {
760 if constexpr (on_host)
761 return exchange_buffer_host_;
762 else
763 return exchange_buffer_default_;
764 }();
765
766 using ExecutionSpace = typename MemorySpace::kokkos_space::execution_space;
767 using Policy =
768 Kokkos::RangePolicy<ExecutionSpace, Kokkos::IndexType<std::size_t>>;
769
770 const auto exec = ExecutionSpace{};
771
772 Kokkos::parallel_for(
773 "sparse_matrix_populate_exchange_buffer",
774 Policy(exec, 0, n_entries_to_be_sent),
775 KOKKOS_LAMBDA(const std::size_t c) {
776 const auto &[row, column_index] = entries_to_be_sent[c];
777 for (unsigned int d = 0; d < n_components; ++d) {
778 const auto offset = sparsity_view.template offset<n_components>(
779 row, column_index, d);
780 exchange_buffer(n_components * c + d) = data(offset);
781 }
782 });
783
784 exec.fence();
785 }
786
787
788 template <typename Number, int n_components, int warp_size, int simd_length>
789 template <typename MemorySpace>
790 void SparseMatrix<Number, n_components, warp_size, simd_length>::
791 update_ghost_rows_on_memory_space()
792 {
793 using HostSpace = dealii::MemorySpace::Host;
794
795 Assert(this->template is_resident<MemorySpace>(),
796 dealii::ExcMessage("The chosen memory space is not resident."));
797
798 constexpr bool on_host =
799 !have_separate_memory_spaces || std::is_same_v<MemorySpace, HostSpace>;
800
801 const auto sparsity_view = sparsity_pattern_->view();
802
803 const auto &receive_targets = sparsity_pattern_->receive_targets();
804 const auto &send_targets = sparsity_pattern_->send_targets();
805
806 const unsigned int mpi_tag =
807 dealii::Utilities::MPI::internal::Tags::partitioner_export_start + 0;
808 Assert(mpi_tag <=
809 dealii::Utilities::MPI::internal::Tags::partitioner_export_end,
810 dealii::ExcInternalError());
811
812 const unsigned int n_requests =
813 receive_targets.size() + send_targets.size();
814 std::vector<MPI_Request> requests(n_requests);
815
816 const auto ghost_offset =
817 sparsity_view.template ghost_offset<n_components>();
818 const auto end_offset = sparsity_view.n_nonzero_elements() * n_components;
819
820 Number *const receive_pointer =
821 on_host ? data_host_.data() + ghost_offset : ghost_buffer_host_.data();
822 Number *const send_pointer = exchange_buffer_host_.data();
823
824 for (unsigned int p = 0; p < receive_targets.size(); ++p) {
825 const auto receive_offset =
826 n_components * (p == 0 ? 0 : receive_targets[p - 1].second);
827 const auto receive_size =
828 (receive_targets[p].second * n_components - receive_offset);
829
830#ifdef DEBUG_MPI_EXCHANGE
831 const auto mpi_rank =
832 dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
833 std::cout << "Rank " << mpi_rank << " receive from "
834 << receive_targets[p].first << " offset = " << receive_offset
835 << " size = " << receive_size << std::endl;
836#endif
837
838 const int ierr =
839 MPI_Irecv(receive_pointer + receive_offset,
840 receive_size,
841 dealii::Utilities::MPI::mpi_type_id_for_type<Number>,
842 receive_targets[p].first,
843 mpi_tag,
844 sparsity_pattern_->partitioner()->get_mpi_communicator(),
845 &requests[p]);
846 AssertThrowMPI(ierr);
847 }
848
849 populate_exchange_buffer_on_memory_space<MemorySpace>();
850
851 if constexpr (!on_host)
852 Kokkos::deep_copy(/*dst*/ exchange_buffer_host_,
853 /*src*/ exchange_buffer_default_);
854
855 for (unsigned int p = 0; p < send_targets.size(); ++p) {
856 const auto send_offset =
857 n_components * (p == 0 ? 0 : send_targets[p - 1].second);
858 const auto send_size =
859 (send_targets[p].second * n_components - send_offset);
860
861#ifdef DEBUG_MPI_EXCHANGE
862 const auto mpi_rank =
863 dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
864 std::cout << "Rank " << mpi_rank << " send to " << send_targets[p].first
865 << " offset = " << send_offset << " size = " << send_size
866 << std::endl;
867#endif
868
869 const int ierr =
870 MPI_Isend(send_pointer + send_offset,
871 send_size,
872 dealii::Utilities::MPI::mpi_type_id_for_type<Number>,
873 send_targets[p].first,
874 mpi_tag,
875 sparsity_pattern_->partitioner()->get_mpi_communicator(),
876 &requests[receive_targets.size() + p]);
877 AssertThrowMPI(ierr);
878 }
879
880#ifdef DEBUG_MPI_EXCHANGE
881 using namespace std::chrono_literals;
882 std::this_thread::sleep_for(200ms);
883#endif
884
885 const int ierr =
886 MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
887 AssertThrowMPI(ierr);
888
889 if constexpr (!on_host) {
890 /* Copy the received ghost range back into device memory: */
891 const auto ghost_range =
892 Kokkos::subview(data_default_,
893 Kokkos::make_pair(std::size_t(ghost_offset),
894 std::size_t(end_offset)));
895 Kokkos::deep_copy(/*dst*/ ghost_range, /*src*/ ghost_buffer_host_);
896 }
897 }
898
899
900 template <typename Number, int n_components, int warp_size, int simd_length>
901 template <typename MemorySpace>
902 void SparseMatrix<Number, n_components, warp_size, simd_length>::
903 compress_on_memory_space(dealii::VectorOperation::values operation
904 [[maybe_unused]])
905 {
906 AssertThrow(operation == dealii::VectorOperation::add,
907 dealii::ExcNotImplemented());
908
909 using HostSpace = dealii::MemorySpace::Host;
910 using DefaultSpace = dealii::MemorySpace::Default;
911
912 AssertThrow((std::is_same_v<MemorySpace, HostSpace>),
913 dealii::ExcNotImplemented());
914
915 Assert(this->template is_resident<MemorySpace>(),
916 dealii::ExcMessage("The chosen memory space is not resident."));
917
918 const auto sparsity_view = sparsity_pattern_->view();
919
920 const auto &receive_targets = sparsity_pattern_->receive_targets();
921 const auto &send_targets = sparsity_pattern_->send_targets();
922 const auto *entries_to_be_sent =
923 sparsity_pattern_->entries_to_be_sent().view();
924 const auto n_entries_to_be_sent =
925 sparsity_pattern_->entries_to_be_sent().size();
926
927 const unsigned int mpi_tag =
928 dealii::Utilities::MPI::internal::Tags::partitioner_export_start + 0;
929 Assert(mpi_tag <=
930 dealii::Utilities::MPI::internal::Tags::partitioner_export_end,
931 dealii::ExcInternalError());
932
933 const unsigned int n_requests =
934 receive_targets.size() + send_targets.size();
935 std::vector<MPI_Request> requests(n_requests);
936
937 /*
938 * Note: For the compress() operation we receive from the "send
939 * targets" and store in the exchange buffer.
940 */
941
942 for (unsigned int p = 0; p < send_targets.size(); ++p) {
943 const auto receive_offset =
944 n_components * (p == 0 ? 0 : send_targets[p - 1].second);
945 const auto receive_size =
946 (send_targets[p].second * n_components - receive_offset);
947
948#ifdef DEBUG_MPI_EXCHANGE
949 const auto mpi_rank =
950 dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
951 std::cout << "Rank " << mpi_rank << " receive from "
952 << send_targets[p].first << " offset = " << receive_offset
953 << " size = " << receive_size << std::endl;
954#endif
955
956 const int ierr =
957 MPI_Irecv(exchange_buffer_host_.data() + receive_offset,
958 receive_size,
959 dealii::Utilities::MPI::mpi_type_id_for_type<Number>,
960 send_targets[p].first,
961 mpi_tag,
962 sparsity_pattern_->partitioner()->get_mpi_communicator(),
963 &requests[p]);
964 AssertThrowMPI(ierr);
965 }
966
967 const auto ghost_offset =
968 sparsity_view.template ghost_offset<n_components>();
969
970 /*
971 * Note: For the compress() operation we send our ghost range to the
972 * "receive targets".
973 */
974
975 for (unsigned int p = 0; p < receive_targets.size(); ++p) {
976 const auto send_offset =
977 n_components * (p == 0 ? 0 : receive_targets[p - 1].second);
978 const auto send_size =
979 (receive_targets[p].second * n_components - send_offset);
980
981#ifdef DEBUG_MPI_EXCHANGE
982 const auto mpi_rank =
983 dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
984 std::cout << "Rank " << mpi_rank << " send to "
985 << receive_targets[p].first << " offset = " << send_offset
986 << " size = " << send_size << std::endl;
987#endif
988
989 const int ierr =
990 MPI_Isend(data_host_.data() + ghost_offset + send_offset,
991 send_size,
992 dealii::Utilities::MPI::mpi_type_id_for_type<Number>,
993 receive_targets[p].first,
994 mpi_tag,
995 sparsity_pattern_->partitioner()->get_mpi_communicator(),
996 &requests[send_targets.size() + p]);
997 AssertThrowMPI(ierr);
998 }
999
1000#ifdef DEBUG_MPI_EXCHANGE
1001 using namespace std::chrono_literals;
1002 std::this_thread::sleep_for(200ms);
1003#endif
1004
1005 const int ierr =
1006 MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
1007 AssertThrowMPI(ierr);
1008
1009 /* Add back contributions and clear ghost range: */
1010
1011 for (std::size_t c = 0; c < n_entries_to_be_sent; ++c) {
1012 const auto &[row, column_index] = entries_to_be_sent[c];
1013 for (unsigned int d = 0; d < n_components; ++d) {
1014 const auto offset =
1015 sparsity_view.template offset<n_components>(row, column_index, d);
1016 data_host_(offset) += exchange_buffer_host_(n_components * c + d);
1017 }
1018 }
1019
1020 zero_out_ghost_rows_on_memory_space<MemorySpace>();
1021 }
1022
1023
1024 template <typename Number,
1025 int n_comp,
1026 int warp_size,
1027 int simd_length,
1028 typename MemorySpace,
1029 bool writable>
1030 SparseMatrixView<Number,
1031 n_comp,
1032 warp_size,
1033 simd_length,
1034 MemorySpace,
1035 writable>::
1036 SparseMatrixView(
1037 SparseMatrix<Number, n_comp, warp_size, simd_length> &sparse_matrix)
1038 requires(writable)
1039 {
1040 reinit(sparse_matrix);
1041 }
1042
1043
1044 template <typename Number,
1045 int n_comp,
1046 int warp_size,
1047 int simd_length,
1048 typename MemorySpace,
1049 bool writable>
1050 SparseMatrixView<Number,
1051 n_comp,
1052 warp_size,
1053 simd_length,
1054 MemorySpace,
1055 writable>::
1056 SparseMatrixView(
1057 const SparseMatrix<Number, n_comp, warp_size, simd_length>
1058 &sparse_matrix)
1059 requires(!writable)
1060 {
1061 reinit(sparse_matrix);
1062 }
1063
1064
1065 template <typename Number,
1066 int n_comp,
1067 int warp_size,
1068 int simd_length,
1069 typename MemorySpace,
1070 bool writable>
1071 template <bool other_writable>
1072 DEAL_II_HOST_DEVICE SparseMatrixView<
1073 Number,
1074 n_comp,
1075 warp_size,
1076 simd_length,
1077 MemorySpace,
1078 writable>::SparseMatrixView(const SparseMatrixView<Number,
1079 n_comp,
1080 warp_size,
1081 simd_length,
1082 MemorySpace,
1083 other_writable> &other)
1084 requires(!writable && other_writable)
1085 : sparse_matrix_(other.sparse_matrix_)
1086 , sparsity_pattern_(other.sparsity_pattern_)
1087 , data_(other.data_)
1088 {
1089 }
1090
1091
1092 template <typename Number,
1093 int n_comp,
1094 int warp_size,
1095 int simd_length,
1096 typename MemorySpace,
1097 bool writable>
1098 template <typename SparseMatrix>
1099 void SparseMatrixView<Number,
1100 n_comp,
1101 warp_size,
1102 simd_length,
1103 MemorySpace,
1104 writable>::reinit(SparseMatrix &sparse_matrix)
1105 requires(writable != std::is_const_v<SparseMatrix>)
1106 {
1107 using HostSpace = dealii::MemorySpace::Host;
1108 using DefaultSpace = dealii::MemorySpace::Default;
1109
1110 static_assert(std::is_same_v<MemorySpace, HostSpace> ||
1111 std::is_same_v<MemorySpace, DefaultSpace>,
1112 "Unexpected Kokkos memory space");
1113
1114 sparse_matrix_ = &sparse_matrix;
1115
1116 /*
1117 * Note: If the host and default memory spaces coincide all views
1118 * reference the host storage.
1119 */
1120 if constexpr (have_separate_memory_spaces &&
1121 !std::is_same_v<MemorySpace, HostSpace>) {
1122 data_ = sparse_matrix.data_default_;
1123 } else {
1124 data_ = sparse_matrix.data_host_;
1125 }
1126
1127 sparsity_pattern_ =
1128 sparse_matrix.sparsity_pattern_->template view<MemorySpace>();
1129 }
1130
1131
1132 template <typename Number,
1133 int n_comp,
1134 int warp_size,
1135 int simd_length,
1136 typename MemorySpace,
1137 bool writable>
1138 template <typename Number2>
1139 DEAL_II_HOST_DEVICE_ALWAYS_INLINE Number2
1140 SparseMatrixView<Number,
1141 n_comp,
1142 warp_size,
1143 simd_length,
1144 MemorySpace,
1145 writable>::read_entry(const unsigned int row,
1146 const unsigned int column_index) const
1147 {
1148 static_assert(
1149 n_comp == 1,
1150 "Attempted to write a scalar value into a tensor-valued matrix entry");
1151
1152 const auto result = read_tensor<Number2>(row, column_index);
1153 return result[0];
1154 }
1155
1156
1157 template <typename Number,
1158 int n_comp,
1159 int warp_size,
1160 int simd_length,
1161 typename MemorySpace,
1162 bool writable>
1163 template <typename Number2, typename Tensor>
1164 DEAL_II_HOST_DEVICE_ALWAYS_INLINE Tensor
1165 SparseMatrixView<Number,
1166 n_comp,
1167 warp_size,
1168 simd_length,
1169 MemorySpace,
1170 writable>::read_tensor(const unsigned int row,
1171 const unsigned int column_index) const
1172 {
1173 static_assert(std::is_same_v<Number2, typename Tensor::value_type>,
1174 "type mismatch");
1175
1176 AssertIndexRange(row, sparsity_pattern_.n_rows());
1177 AssertIndexRange(column_index, sparsity_pattern_.row_length(row));
1178
1179 Tensor result;
1180
1181 using VA = dealii::VectorizedArray<Number, simd_length>;
1182 if constexpr (std::is_same_v<VA, Number2>) {
1183 /*
1184 * Vectorized fast access. Indices must be in the range
1185 * [0,n_internal), index must be divisible by simd_length
1186 */
1187
1188 Assert(row < sparsity_pattern_.n_internal_dofs(),
1189 dealii::ExcMessage(
1190 "Vectorized access only possible in vectorized part"));
1191 Assert(row % simd_length == 0,
1192 dealii::ExcMessage(
1193 "Access only supported for rows at the SIMD granularity"));
1194
1195 const Number *load_pos = data_.data();
1196 load_pos +=
1197 sparsity_pattern_.template offset_internal<n_comp>(row, column_index);
1198
1199 for (unsigned int d = 0; d < n_comp; ++d)
1200 result[d].load(load_pos + d * warp_size);
1201
1202 } else {
1203 /*
1204 * Non-vectorized slow access. Supports all row indices in [0,n_owned):
1205 */
1206
1207 for (unsigned int d = 0; d < n_comp; ++d) {
1208 const auto offset =
1209 sparsity_pattern_.template offset<n_comp>(row, column_index, d);
1210 result[d] = data_(offset);
1211 }
1212 }
1213
1214 return result;
1215 }
1216
1217
1218 template <typename Number,
1219 int n_comp,
1220 int warp_size,
1221 int simd_length,
1222 typename MemorySpace,
1223 bool writable>
1224 template <typename Number2>
1225 DEAL_II_HOST_DEVICE_ALWAYS_INLINE Number2
1226 SparseMatrixView<Number,
1227 n_comp,
1228 warp_size,
1229 simd_length,
1230 MemorySpace,
1231 writable>::read_transposed_entry(const unsigned int row,
1232 const unsigned int
1233 column_index) const
1234 {
1235 static_assert(
1236 n_comp == 1,
1237 "Attempted to write a scalar value into a tensor-valued matrix entry");
1238
1239 const auto result = read_transposed_tensor<Number2>(row, column_index);
1240 return result[0];
1241 }
1242
1243
1244 template <typename Number,
1245 int n_comp,
1246 int warp_size,
1247 int simd_length,
1248 typename MemorySpace,
1249 bool writable>
1250 template <typename Number2, typename Tensor>
1251 DEAL_II_HOST_DEVICE_ALWAYS_INLINE Tensor
1252 SparseMatrixView<Number,
1253 n_comp,
1254 warp_size,
1255 simd_length,
1256 MemorySpace,
1257 writable>::read_transposed_tensor(const unsigned int row,
1258 const unsigned int
1259 column_index) const
1260 {
1261 static_assert(std::is_same_v<Number2, typename Tensor::value_type>,
1262 "type mismatch");
1263
1264 AssertIndexRange(row, sparsity_pattern_.n_rows());
1265 AssertIndexRange(column_index, sparsity_pattern_.row_length(row));
1266
1267 dealii::Tensor<1, n_comp, Number2> result;
1268
1269 using VA = dealii::VectorizedArray<Number, simd_length>;
1270 if constexpr (std::is_same_v<VA, Number2> && (n_comp == 1)) {
1271 /*
1272 * Vectorized fast access. Indices must be in the range
1273 * [0,n_internal), index must be divisible by simd_length
1274 */
1275
1276 Assert(row < sparsity_pattern_.n_internal_dofs(),
1277 dealii::ExcMessage(
1278 "Vectorized access only possible in vectorized part"));
1279 Assert(row % simd_length == 0,
1280 dealii::ExcMessage(
1281 "Access only supported for rows at the SIMD granularity"));
1282
1283 const auto offsets =
1284 sparsity_pattern_.template transposed_offset_internal<1>(
1285 row, column_index);
1286 result[0].gather(data_.data(), offsets);
1287
1288 } else if constexpr (std::is_same_v<VA, Number2> && (n_comp != 1)) {
1289
1290 /* not implemented */
1291 Assert(false,
1292 dealii::ExcMessage("Vectorized transposed access to multiple "
1293 "components is not implemented."));
1294 __builtin_trap();
1295
1296 } else {
1297 /*
1298 * Non-vectorized slow access. Supports all row indices in [0,n_owned):
1299 */
1300
1301 for (unsigned int d = 0; d < n_comp; ++d) {
1302 const auto offset =
1303 sparsity_pattern_.template transposed_offset<n_comp>(
1304 row, column_index, d);
1305 result[d] = data_(offset);
1306 }
1307 }
1308
1309 return result;
1310 }
1311
1312
1313 template <typename Number,
1314 int n_comp,
1315 int warp_size,
1316 int simd_length,
1317 typename MemorySpace,
1318 bool writable>
1319 template <typename Number2>
1320 DEAL_II_HOST_DEVICE_ALWAYS_INLINE void
1321 SparseMatrixView<Number,
1322 n_comp,
1323 warp_size,
1324 simd_length,
1325 MemorySpace,
1326 writable>::write_entry(const Number2 entry,
1327 const unsigned int row,
1328 const unsigned int column_index,
1329 const bool do_streaming_store) const
1330 requires(writable)
1331 {
1332 static_assert(
1333 n_comp == 1,
1334 "Attempted to write a scalar value into a tensor-valued matrix entry");
1335
1336 AssertIndexRange(row, sparsity_pattern_.n_rows());
1337 AssertIndexRange(column_index, sparsity_pattern_.row_length(row));
1338
1339 dealii::Tensor<1, n_comp, Number2> tensor;
1340 tensor[0] = entry;
1341
1342 write_tensor<Number2>(tensor, row, column_index, do_streaming_store);
1343 }
1344
1345
1346 template <typename Number,
1347 int n_comp,
1348 int warp_size,
1349 int simd_length,
1350 typename MemorySpace,
1351 bool writable>
1352 template <typename Number2, typename Tensor>
1353 DEAL_II_HOST_DEVICE_ALWAYS_INLINE void
1354 SparseMatrixView<Number,
1355 n_comp,
1356 warp_size,
1357 simd_length,
1358 MemorySpace,
1359 writable>::write_tensor(const Tensor &tensor,
1360 const unsigned int row,
1361 const unsigned int column_index,
1362 const bool do_streaming_store) const
1363 requires(writable)
1364 {
1365 AssertIndexRange(row, sparsity_pattern_.n_rows());
1366 AssertIndexRange(column_index, sparsity_pattern_.row_length(row));
1367
1368 using VA = dealii::VectorizedArray<Number, simd_length>;
1369 if constexpr (std::is_same_v<VA, Number2>) {
1370 /*
1371 * Vectorized fast access. Indices must be in the range [0,n_internal),
1372 * index must be divisible by simd_length:
1373 */
1374
1375 Assert(row < sparsity_pattern_.n_internal_dofs(),
1376 dealii::ExcMessage(
1377 "Vectorized access only possible in vectorized part"));
1378 Assert(row % simd_length == 0,
1379 dealii::ExcMessage(
1380 "Access only supported for rows at the SIMD granularity"));
1381
1382 Number *store_pos = data_.data();
1383 store_pos +=
1384 sparsity_pattern_.template offset_internal<n_comp>(row, column_index);
1385
1386 if (do_streaming_store)
1387 for (unsigned int d = 0; d < n_comp; ++d)
1388 tensor[d].streaming_store(store_pos + d * warp_size);
1389 else
1390 for (unsigned int d = 0; d < n_comp; ++d)
1391 tensor[d].store(store_pos + d * warp_size);
1392
1393 } else {
1394 /*
1395 * Non-vectorized slow access. Supports all row indices in [0,n_owned):
1396 */
1397
1398 for (unsigned int d = 0; d < n_comp; ++d) {
1399 const auto offset =
1400 sparsity_pattern_.template offset<n_comp>(row, column_index, d);
1401 data_(offset) = tensor[d];
1402 }
1403 }
1404 }
1405
1406
1407 template <typename Number,
1408 int n_comp,
1409 int warp_size,
1410 int simd_length,
1411 typename MemorySpace,
1412 bool writable>
1413 template <typename Number2>
1414 DEAL_II_HOST_DEVICE_ALWAYS_INLINE void
1415 SparseMatrixView<Number,
1416 n_comp,
1417 warp_size,
1418 simd_length,
1419 MemorySpace,
1420 writable>::add_entry(const Number2 entry,
1421 const unsigned int row,
1422 const unsigned int column_index) const
1423 requires(writable)
1424 {
1425 static_assert(
1426 n_comp == 1,
1427 "Attempted to write a scalar value into a tensor-valued matrix entry");
1428
1429 AssertIndexRange(row, sparsity_pattern_.n_rows());
1430 AssertIndexRange(column_index, sparsity_pattern_.row_length(row));
1431
1432 dealii::Tensor<1, n_comp, Number2> tensor;
1433 tensor[0] = entry;
1434
1435 add_tensor<Number2>(tensor, row, column_index);
1436 }
1437
1438
1439 template <typename Number,
1440 int n_comp,
1441 int warp_size,
1442 int simd_length,
1443 typename MemorySpace,
1444 bool writable>
1445 template <typename Number2, typename Tensor>
1446 DEAL_II_HOST_DEVICE_ALWAYS_INLINE void
1447 SparseMatrixView<Number,
1448 n_comp,
1449 warp_size,
1450 simd_length,
1451 MemorySpace,
1452 writable>::add_tensor(const Tensor &tensor,
1453 const unsigned int row,
1454 const unsigned int column_index) const
1455 requires(writable)
1456 {
1457 AssertIndexRange(row, sparsity_pattern_.n_rows());
1458 AssertIndexRange(column_index, sparsity_pattern_.row_length(row));
1459
1460 using VA = dealii::VectorizedArray<Number, simd_length>;
1461 if constexpr (std::is_same_v<VA, Number2>) {
1462 /*
1463 * Vectorized fast access. Indices must be in the range [0,n_internal),
1464 * index must be divisible by simd_length:
1465 */
1466
1467 Assert(row < sparsity_pattern_.n_internal_dofs(),
1468 dealii::ExcMessage(
1469 "Vectorized access only possible in vectorized part"));
1470 Assert(row % simd_length == 0,
1471 dealii::ExcMessage(
1472 "Access only supported for rows at the SIMD granularity"));
1473
1474 Number *store_pos = data_.data();
1475 store_pos +=
1476 sparsity_pattern_.template offset_internal<n_comp>(row, column_index);
1477
1478 for (unsigned int d = 0; d < n_comp; ++d) {
1479 auto temp = tensor[d];
1480 temp.load(store_pos + d * warp_size);
1481 temp += tensor[d];
1482 temp.store(store_pos + d * warp_size);
1483 }
1484
1485 } else {
1486 /*
1487 * Non-vectorized slow access. Supports all row indices in [0,n_owned):
1488 */
1489
1490 for (unsigned int d = 0; d < n_comp; ++d) {
1491 const auto offset =
1492 sparsity_pattern_.template offset<n_comp>(row, column_index, d);
1493 data_(offset) += tensor[d]; /*add*/
1494 ;
1495 }
1496 }
1497 }
1498
1499
1500 template <typename Number,
1501 int n_comp,
1502 int warp_size,
1503 int simd_length,
1504 typename MemorySpace,
1505 bool writable>
1506 void SparseMatrixView<Number,
1507 n_comp,
1508 warp_size,
1509 simd_length,
1510 MemorySpace,
1511 writable>::zero_out_ghost_rows() const
1512 requires(writable)
1513 {
1514 using HostSpace = dealii::MemorySpace::Host;
1515 using DefaultSpace = dealii::MemorySpace::Default;
1516
1517 static_assert(std::is_same_v<MemorySpace, HostSpace> ||
1518 std::is_same_v<MemorySpace, DefaultSpace>,
1519 "Unexpected Kokkos memory space");
1520
1521 Assert(sparse_matrix_->template is_resident<MemorySpace>(),
1522 dealii::ExcMessage("The chosen memory space is not resident."));
1523
1524 sparse_matrix_->template zero_out_ghost_rows_on_memory_space<MemorySpace>();
1525 }
1526
1527
1528 template <typename Number,
1529 int n_comp,
1530 int warp_size,
1531 int simd_length,
1532 typename MemorySpace,
1533 bool writable>
1534 void SparseMatrixView<Number,
1535 n_comp,
1536 warp_size,
1537 simd_length,
1538 MemorySpace,
1539 writable>::update_ghost_rows() const
1540 requires(writable)
1541 {
1542 using HostSpace = dealii::MemorySpace::Host;
1543 using DefaultSpace = dealii::MemorySpace::Default;
1544
1545 static_assert(std::is_same_v<MemorySpace, HostSpace> ||
1546 std::is_same_v<MemorySpace, DefaultSpace>,
1547 "Unexpected Kokkos memory space");
1548
1549 Assert(sparse_matrix_->template is_resident<MemorySpace>(),
1550 dealii::ExcMessage("The chosen memory space is not resident."));
1551
1552 sparse_matrix_->template update_ghost_rows_on_memory_space<MemorySpace>();
1553 }
1554
1555
1556 template <typename Number,
1557 int n_comp,
1558 int warp_size,
1559 int simd_length,
1560 typename MemorySpace,
1561 bool writable>
1562 void SparseMatrixView<Number,
1563 n_comp,
1564 warp_size,
1565 simd_length,
1566 MemorySpace,
1567 writable>::compress(dealii::VectorOperation::values
1568 operation) const
1569 requires(writable)
1570 {
1571 using HostSpace = dealii::MemorySpace::Host;
1572 using DefaultSpace = dealii::MemorySpace::Default;
1573
1574 static_assert(std::is_same_v<MemorySpace, HostSpace> ||
1575 std::is_same_v<MemorySpace, DefaultSpace>,
1576 "Unexpected Kokkos memory space");
1577
1578 Assert(sparse_matrix_->template is_resident<MemorySpace>(),
1579 dealii::ExcMessage("The chosen memory space is not resident."));
1580
1581 sparse_matrix_->template compress_on_memory_space<MemorySpace>(operation);
1582 }
1583
1584
1585 template <typename Number,
1586 int n_comp,
1587 int warp_size,
1588 int simd_length,
1589 typename FM>
1590 void distribute_local_to_global(
1591 const FM &cell_matrix,
1592 const std::vector<dealii::types::global_dof_index> &dof_indices_row,
1593 const std::vector<dealii::types::global_dof_index> &dof_indices_column,
1594 const dealii::AffineConstraints<Number> &affine_constraints
1595 [[maybe_unused]],
1596 SparseMatrix<Number, n_comp, warp_size, simd_length> &sparse_matrix)
1597 {
1598 constexpr bool is_matrix = std::is_same_v<FM, dealii::FullMatrix<Number>>;
1599 constexpr bool is_array =
1600 std::is_same_v<FM, std::array<dealii::FullMatrix<Number>, n_comp>>;
1601 static_assert((n_comp == 1 && is_matrix) || is_array, "not implemented");
1602
1603 if constexpr (is_matrix) {
1604 Assert(cell_matrix.m() == dof_indices_row.size(),
1605 dealii::ExcInternalError());
1606 Assert(cell_matrix.n() == dof_indices_column.size(),
1607 dealii::ExcInternalError());
1608 } else if constexpr (is_array) {
1609 Assert(cell_matrix.size() == n_comp, dealii::ExcInternalError());
1610 for (unsigned int d = 0; d < n_comp; ++d) {
1611 Assert(cell_matrix[d].m() == dof_indices_row.size(),
1612 dealii::ExcInternalError());
1613 Assert(cell_matrix[d].n() == dof_indices_column.size(),
1614 dealii::ExcInternalError());
1615 }
1616 }
1617
1618 const auto sparse_matrix_view = sparse_matrix.view();
1619
1620 const auto &sparsity_pattern = sparse_matrix.sparsity_pattern();
1621 const auto sparsity_pattern_view = sparsity_pattern.view();
1622 const auto &partitioner = sparsity_pattern.partitioner();
1623
1624 /*
1625 * Helper that inserts a single entry into the matrix indexed by (r, c)
1626 * in the cell_matrix and by (i, j) in the sparse matrix and multiplied
1627 * by a weight c_ij.
1628 */
1629
1630 const auto insert_entry =
1631 [&](auto r, auto c, auto i, auto j, auto c_ij) DEAL_II_ALWAYS_INLINE {
1632 if constexpr (is_matrix) {
1633 const Number &entry = cell_matrix(r, c);
1634 if (entry == Number{})
1635 return;
1636 const auto col_idx = sparsity_pattern_view.column_index(i, j);
1637 sparse_matrix_view.add_entry(c_ij * entry, i, col_idx);
1638 } else if constexpr (is_array) {
1639 dealii::Tensor<1, n_comp, Number> entry;
1640 for (unsigned int k = 0; k < n_comp; ++k)
1641 entry[k] = cell_matrix[k](r, c);
1642 if (entry == dealii::Tensor<1, n_comp>{})
1643 return;
1644 const auto col_idx = sparsity_pattern_view.column_index(i, j);
1645 sparse_matrix_view.add_tensor(c_ij * entry, i, col_idx);
1646 }
1647 };
1648
1649 /*
1650 * Helper that iterates over row entries:
1651 */
1652
1653 const auto iterate_over_row_entries =
1654 [&](const auto r, const auto i, const auto c_i) DEAL_II_ALWAYS_INLINE {
1655 /* Iterate over columns: c - column index; j_global, j - dof index */
1656 for (unsigned int c = 0; c < dof_indices_column.size(); ++c) {
1657 const auto j_global = dof_indices_column[c];
1658 if (affine_constraints.is_constrained(j_global)) {
1659 const auto &line =
1660 *affine_constraints.get_constraint_entries(j_global);
1661 for (const auto &[k_global, c_k] : line) {
1662 const auto k = partitioner->global_to_local(k_global);
1663 insert_entry(r, c, i, k, c_i * c_k);
1664 }
1665 } else {
1666 const auto j = partitioner->global_to_local(j_global);
1667 insert_entry(r, c, i, j, c_i);
1668 }
1669 }
1670 };
1671
1672 /* Now, iterate over rows: r - row index; i_global, i - dof index */
1673 for (unsigned int r = 0; r < dof_indices_row.size(); ++r) {
1674 const auto i_global = dof_indices_row[r];
1675 if (affine_constraints.is_constrained(i_global)) {
1676 const auto &line = *affine_constraints.get_constraint_entries(i_global);
1677 for (const auto &[k_global, c_k] : line) {
1678 const auto k = partitioner->global_to_local(k_global);
1679 iterate_over_row_entries(r, k, c_k);
1680 }
1681 } else {
1682 const auto i = partitioner->global_to_local(i_global);
1683 iterate_over_row_entries(r, i, Number(1.));
1684 }
1685 }
1686 }
1687
1688
1689 template <typename Number,
1690 int n_comp,
1691 int warp_size,
1692 int simd_length,
1693 typename FM>
1694 void distribute_local_to_global(
1695 const FM &cell_matrix,
1696 const std::vector<dealii::types::global_dof_index> &dof_indices,
1697 const dealii::AffineConstraints<Number> &affine_constraints,
1698 SparseMatrix<Number, n_comp, warp_size, simd_length> &sparse_matrix)
1699 {
1700 constexpr bool is_matrix = std::is_same_v<FM, dealii::FullMatrix<Number>>;
1701 constexpr bool is_array =
1702 std::is_same_v<FM, std::array<dealii::FullMatrix<Number>, n_comp>>;
1703 static_assert((n_comp == 1 && is_matrix) || is_array, "not implemented");
1704
1705 distribute_local_to_global(cell_matrix,
1706 dof_indices,
1707 dof_indices,
1708 affine_constraints,
1709 sparse_matrix);
1710 }
1711
1712#endif
1713} // namespace ryujin
ryujin::MirroredStorage::transfer_policy
TransferPolicy transfer_policy() const
ryujin::SparseMatrixView::zero_out_ghost_rows
void zero_out_ghost_rows() const
ryujin::SparseMatrixView::read_transposed_tensor
DEAL_II_HOST_DEVICE Tensor read_transposed_tensor(const unsigned int row, const unsigned int column_index) const
ryujin::SparseMatrixView::ACCESSOR_READ_ONLY
ACCESSOR_READ_ONLY(sparsity_pattern)
ryujin::SparseMatrixView::read_transposed_entry
DEAL_II_HOST_DEVICE Number2 read_transposed_entry(const unsigned int row, const unsigned int column_index) const
ryujin::SparseMatrixView::write_tensor
DEAL_II_HOST_DEVICE void write_tensor(const Tensor &tensor, const unsigned int row, const unsigned int column_index, const bool do_streaming_store=false) const
ryujin::SparseMatrixView::reinit
void reinit(SparseMatrix &sparse_matrix)
ryujin::SparseMatrixView::add_tensor
DEAL_II_HOST_DEVICE void add_tensor(const Tensor &tensor, const unsigned int row, const unsigned int column_index) const
ryujin::SparseMatrixView::SparseMatrixView
SparseMatrixView(const SparseMatrix< Number, n_comp, warp_size, simd_length > &sparse_matrix)
ryujin::SparseMatrixView::write_entry
DEAL_II_HOST_DEVICE void write_entry(const Number2 entry, const unsigned int row, const unsigned int column_index, const bool do_streaming_store=false) const
ryujin::SparseMatrixView::SparseMatrixView
DEAL_II_HOST_DEVICE SparseMatrixView(const SparseMatrixView< Number, n_comp, warp_size, simd_length, MemorySpace, other_writable > &other)
ryujin::SparseMatrixView::read_entry
DEAL_II_HOST_DEVICE Number2 read_entry(const unsigned int row, const unsigned int column_index) const
ryujin::SparseMatrixView::add_entry
DEAL_II_HOST_DEVICE void add_entry(const Number2 entry, const unsigned int row, const unsigned int column_index) const
ryujin::SparseMatrixView::update_ghost_rows
void update_ghost_rows() const
ryujin::SparseMatrixView::SparseMatrixView
SparseMatrixView(SparseMatrix< Number, n_comp, warp_size, simd_length > &sparse_matrix)
ryujin::SparseMatrixView::read_tensor
DEAL_II_HOST_DEVICE Tensor read_tensor(const unsigned int row, const unsigned int column_index) const
ryujin::SparseMatrixView::compress
void compress(dealii::VectorOperation::values operation) const
ryujin::SparseMatrix::populate_exchange_buffer_on_memory_space
void populate_exchange_buffer_on_memory_space()
ryujin::SparseMatrix::view
SparseMatrixView< Number, n_comp, warp_size, simd_length, MemorySpace, false > view() const
ryujin::SparseMatrix::compress_on_memory_space
void compress_on_memory_space(dealii::VectorOperation::values operation)
ryujin::SparseMatrix::update_ghost_rows_on_memory_space
void update_ghost_rows_on_memory_space()
ryujin::SparseMatrix::zero_out_ghost_rows_on_memory_space
void zero_out_ghost_rows_on_memory_space()
ryujin::SparseMatrix::view
SparseMatrixView< Number, n_comp, warp_size, simd_length, MemorySpace, true > view()
ryujin::SparseMatrix::ACCESSOR_READ_ONLY
ACCESSOR_READ_ONLY(sparsity_pattern)
ryujin::SparseMatrix::reinit
void reinit(const SparsityPattern< warp_size > &sparsity, const TransferPolicy transfer_policy=TransferPolicy::explicit_transfers)
ryujin::SparseMatrix::SparseMatrix
SparseMatrix(const SparsityPattern< warp_size > &sparsity, const TransferPolicy transfer_policy=TransferPolicy::explicit_transfers)
ryujin::warp_size
constexpr unsigned int warp_size
Definition gpu.h:46
ryujin::TransferPolicy
TransferPolicy
Definition gpu.h:88
ryujin::have_separate_memory_spaces
constexpr bool have_separate_memory_spaces
Definition gpu.h:29
ryujin::TransferPolicy::explicit_transfers
@ explicit_transfers
ryujin::distribute_local_to_global
void distribute_local_to_global(const FullMatrix &cell_matrix, const std::vector< dealii::types::global_dof_index > &dof_indices_row, const std::vector< dealii::types::global_dof_index > &dof_indices_column, const dealii::AffineConstraints< Number > &affine_constraints, SparseMatrix< Number, n_comp, warp_size, simd_length > &sparse_matrix)
sparsity_pattern.h
Generated by doxygen 1.9.8