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