10 #ifndef IFPACK2_BLOCKTRIDICONTAINER_IMPL_HPP
11 #define IFPACK2_BLOCKTRIDICONTAINER_IMPL_HPP
18 #include <Tpetra_Details_extractMpiCommFromTeuchos.hpp>
19 #include <Tpetra_Distributor.hpp>
20 #include <Tpetra_BlockMultiVector.hpp>
22 #include <Kokkos_ArithTraits.hpp>
23 #include <KokkosBatched_Util.hpp>
24 #include <KokkosBatched_Vector.hpp>
25 #include <KokkosBatched_Copy_Decl.hpp>
26 #include <KokkosBatched_Copy_Impl.hpp>
27 #include <KokkosBatched_AddRadial_Decl.hpp>
28 #include <KokkosBatched_AddRadial_Impl.hpp>
29 #include <KokkosBatched_SetIdentity_Decl.hpp>
30 #include <KokkosBatched_SetIdentity_Impl.hpp>
31 #include <KokkosBatched_Gemm_Decl.hpp>
32 #include <KokkosBatched_Gemm_Serial_Impl.hpp>
33 #include <KokkosBatched_Gemm_Team_Impl.hpp>
34 #include <KokkosBatched_Gemv_Decl.hpp>
35 #include <KokkosBatched_Gemv_Team_Impl.hpp>
36 #include <KokkosBatched_Trsm_Decl.hpp>
37 #include <KokkosBatched_Trsm_Serial_Impl.hpp>
38 #include <KokkosBatched_Trsm_Team_Impl.hpp>
39 #include <KokkosBatched_Trsv_Decl.hpp>
40 #include <KokkosBatched_Trsv_Serial_Impl.hpp>
41 #include <KokkosBatched_Trsv_Team_Impl.hpp>
42 #include <KokkosBatched_LU_Decl.hpp>
43 #include <KokkosBatched_LU_Serial_Impl.hpp>
44 #include <KokkosBatched_LU_Team_Impl.hpp>
46 #include <KokkosBlas1_nrm1.hpp>
47 #include <KokkosBlas1_nrm2.hpp>
51 #include "Ifpack2_BlockHelper.hpp"
52 #include "Ifpack2_BlockComputeResidualVector.hpp"
53 #include "Ifpack2_BlockComputeResidualAndSolve.hpp"
60 #if defined(KOKKOS_ENABLE_CUDA) && defined(IFPACK2_BLOCKTRIDICONTAINER_ENABLE_PROFILE)
61 #include "cuda_profiler_api.h"
66 #define IFPACK2_BLOCKTRIDICONTAINER_USE_MPI_3
74 #define IFPACK2_BLOCKTRIDICONTAINER_USE_PINNED_MEMORY_FOR_MPI
78 #define IFPACK2_BLOCKTRIDICONTAINER_USE_CUDA_SPACE
81 #if defined(HAVE_IFPACK2_BLOCKTRIDICONTAINER_SMALL_SCALAR)
82 #define IFPACK2_BLOCKTRIDICONTAINER_USE_SMALL_SCALAR_FOR_BLOCKTRIDIAG
86 #define IFPACK2_BLOCKTRIDICONTAINER_USE_EXEC_SPACE_INSTANCES
90 namespace BlockTriDiContainerDetails {
92 namespace KB = KokkosBatched;
99 template <
typename MemoryTraitsType, Kokkos::MemoryTraitsFlags flag>
100 using MemoryTraits = Kokkos::MemoryTraits<MemoryTraitsType::is_unmanaged |
101 MemoryTraitsType::is_random_access |
104 template <
typename ViewType>
105 using Unmanaged = Kokkos::View<
typename ViewType::data_type,
106 typename ViewType::array_layout,
107 typename ViewType::device_type,
108 MemoryTraits<typename ViewType::memory_traits, Kokkos::Unmanaged>>;
109 template <
typename ViewType>
110 using Atomic = Kokkos::View<
typename ViewType::data_type,
111 typename ViewType::array_layout,
112 typename ViewType::device_type,
113 MemoryTraits<typename ViewType::memory_traits, Kokkos::Atomic>>;
114 template <
typename ViewType>
115 using Const = Kokkos::View<
typename ViewType::const_data_type,
116 typename ViewType::array_layout,
117 typename ViewType::device_type,
118 typename ViewType::memory_traits>;
119 template <
typename ViewType>
120 using ConstUnmanaged = Const<Unmanaged<ViewType>>;
122 template <
typename ViewType>
123 using AtomicUnmanaged = Atomic<Unmanaged<ViewType>>;
125 template <
typename ViewType>
126 using Unmanaged = Kokkos::View<
typename ViewType::data_type,
127 typename ViewType::array_layout,
128 typename ViewType::device_type,
129 MemoryTraits<typename ViewType::memory_traits, Kokkos::Unmanaged>>;
131 template <
typename ViewType>
132 using Scratch = Kokkos::View<
typename ViewType::data_type,
133 typename ViewType::array_layout,
134 typename ViewType::execution_space::scratch_memory_space,
135 MemoryTraits<typename ViewType::memory_traits, Kokkos::Unmanaged>>;
140 template <
typename T>
142 #if defined(IFPACK2_BLOCKTRIDICONTAINER_USE_SMALL_SCALAR_FOR_BLOCKTRIDIAG)
148 #if defined(KOKKOS_ENABLE_CUDA) && defined(IFPACK2_BLOCKTRIDICONTAINER_ENABLE_PROFILE)
149 #define IFPACK2_BLOCKTRIDICONTAINER_PROFILER_REGION_BEGIN \
150 KOKKOS_IMPL_CUDA_SAFE_CALL(cudaProfilerStart());
152 #define IFPACK2_BLOCKTRIDICONTAINER_PROFILER_REGION_END \
153 { KOKKOS_IMPL_CUDA_SAFE_CALL(cudaProfilerStop()); }
155 #define IFPACK2_BLOCKTRIDICONTAINER_PROFILER_REGION_BEGIN
157 #define IFPACK2_BLOCKTRIDICONTAINER_PROFILER_REGION_END
163 template <
typename MatrixType>
166 IFPACK2_BLOCKHELPER_TIMER(
"BlockTriDi::CreateBlockCrsTpetraImporter", CreateBlockCrsTpetraImporter);
168 using tpetra_map_type =
typename impl_type::tpetra_map_type;
169 using tpetra_mv_type =
typename impl_type::tpetra_block_multivector_type;
170 using tpetra_import_type =
typename impl_type::tpetra_import_type;
171 using crs_matrix_type =
typename impl_type::tpetra_crs_matrix_type;
172 using block_crs_matrix_type =
typename impl_type::tpetra_block_crs_matrix_type;
174 auto A_crs = Teuchos::rcp_dynamic_cast<
const crs_matrix_type>(A);
175 auto A_bcrs = Teuchos::rcp_dynamic_cast<
const block_crs_matrix_type>(A);
177 bool hasBlockCrsMatrix = !A_bcrs.is_null();
180 const auto g = hasBlockCrsMatrix ? A_bcrs->getCrsGraph() : *(A_crs->getCrsGraph());
182 const auto blocksize = hasBlockCrsMatrix ? A_bcrs->getBlockSize() : 1;
183 const auto src =
Teuchos::rcp(
new tpetra_map_type(tpetra_mv_type::makePointMap(*g.getDomainMap(), blocksize)));
184 const auto tgt =
Teuchos::rcp(
new tpetra_map_type(tpetra_mv_type::makePointMap(*g.getColMap(), blocksize)));
185 IFPACK2_BLOCKHELPER_TIMER_FENCE(
typename BlockHelperDetails::ImplType<MatrixType>::execution_space)
194 template <
typename MatrixType>
195 struct AsyncableImport {
203 #if !defined(HAVE_IFPACK2_MPI)
204 typedef int MPI_Request;
205 typedef int MPI_Comm;
207 using scalar_type =
typename impl_type::scalar_type;
211 static int isend(
const MPI_Comm comm,
const char *buf,
int count,
int dest,
int tag, MPI_Request *ireq) {
212 #ifdef HAVE_IFPACK2_MPI
214 int ret = MPI_Isend(const_cast<char *>(buf), count, MPI_CHAR, dest, tag, comm, ireq == NULL ? &ureq : ireq);
215 if (ireq == NULL) MPI_Request_free(&ureq);
222 static int irecv(
const MPI_Comm comm,
char *buf,
int count,
int src,
int tag, MPI_Request *ireq) {
223 #ifdef HAVE_IFPACK2_MPI
225 int ret = MPI_Irecv(buf, count, MPI_CHAR, src, tag, comm, ireq == NULL ? &ureq : ireq);
226 if (ireq == NULL) MPI_Request_free(&ureq);
233 static int waitany(
int count, MPI_Request *reqs,
int *index) {
234 #ifdef HAVE_IFPACK2_MPI
235 return MPI_Waitany(count, reqs, index, MPI_STATUS_IGNORE);
241 static int waitall(
int count, MPI_Request *reqs) {
242 #ifdef HAVE_IFPACK2_MPI
243 return MPI_Waitall(count, reqs, MPI_STATUS_IGNORE);
250 using tpetra_map_type =
typename impl_type::tpetra_map_type;
251 using tpetra_import_type =
typename impl_type::tpetra_import_type;
253 using local_ordinal_type =
typename impl_type::local_ordinal_type;
254 using global_ordinal_type =
typename impl_type::global_ordinal_type;
258 using int_1d_view_host = Kokkos::View<int *, Kokkos::HostSpace>;
259 using local_ordinal_type_1d_view_host = Kokkos::View<local_ordinal_type *, Kokkos::HostSpace>;
261 using execution_space =
typename impl_type::execution_space;
262 using memory_space =
typename impl_type::memory_space;
263 using local_ordinal_type_1d_view =
typename impl_type::local_ordinal_type_1d_view;
265 using size_type_1d_view_host = Kokkos::View<size_type *, Kokkos::HostSpace>;
267 #if defined(KOKKOS_ENABLE_CUDA)
268 using impl_scalar_type_1d_view =
269 typename std::conditional<std::is_same<execution_space, Kokkos::Cuda>::value,
270 #if defined(IFPACK2_BLOCKTRIDICONTAINER_USE_PINNED_MEMORY_FOR_MPI)
271 Kokkos::View<impl_scalar_type *, Kokkos::CudaHostPinnedSpace>,
272 #elif defined(IFPACK2_BLOCKTRIDICONTAINER_USE_CUDA_MEMORY_FOR_MPI)
273 Kokkos::View<impl_scalar_type *, Kokkos::CudaSpace>,
274 #else // no experimental macros are defined
275 typename impl_type::impl_scalar_type_1d_view,
277 typename impl_type::impl_scalar_type_1d_view>::type;
279 using impl_scalar_type_1d_view =
typename impl_type::impl_scalar_type_1d_view;
281 using impl_scalar_type_1d_view_host = Kokkos::View<impl_scalar_type *, Kokkos::HostSpace>;
282 using impl_scalar_type_2d_view =
typename impl_type::impl_scalar_type_2d_view;
283 using impl_scalar_type_2d_view_tpetra =
typename impl_type::impl_scalar_type_2d_view_tpetra;
285 #ifdef HAVE_IFPACK2_MPI
289 impl_scalar_type_2d_view_tpetra remote_multivector;
290 local_ordinal_type blocksize;
292 template <
typename T>
293 struct SendRecvPair {
298 SendRecvPair<int_1d_view_host> pids;
299 SendRecvPair<std::vector<MPI_Request>> reqs;
300 SendRecvPair<size_type_1d_view> offset;
301 SendRecvPair<size_type_1d_view_host> offset_host;
302 SendRecvPair<local_ordinal_type_1d_view> lids;
303 SendRecvPair<impl_scalar_type_1d_view> buffer;
304 SendRecvPair<impl_scalar_type_1d_view_host> buffer_host;
306 local_ordinal_type_1d_view dm2cm;
308 #if defined(KOKKOS_ENABLE_CUDA) || defined(KOKKOS_ENABLE_HIP) || defined(KOKKOS_ENABLE_SYCL)
309 using exec_instance_1d_std_vector = std::vector<execution_space>;
310 exec_instance_1d_std_vector exec_instances;
316 const size_type_1d_view &offs) {
318 Kokkos::View<size_t *, Kokkos::HostSpace> lens_host(const_cast<size_t *>(lens.
getRawPtr()), lens.
size());
319 const auto lens_device = Kokkos::create_mirror_view_and_copy(memory_space(), lens_host);
322 const Kokkos::RangePolicy<execution_space> policy(0, offs.extent(0));
323 const local_ordinal_type lens_size = lens_device.extent(0);
324 Kokkos::parallel_scan(
325 "AsyncableImport::RangePolicy::setOffsetValues",
326 policy, KOKKOS_LAMBDA(
const local_ordinal_type &i, size_type &update,
const bool &
final) {
329 update += (i < lens_size ? lens_device[i] : 0);
334 const size_type_1d_view_host &offs) {
336 Kokkos::View<size_t *, Kokkos::HostSpace> lens_host(const_cast<size_t *>(lens.
getRawPtr()), lens.
size());
337 const auto lens_device = Kokkos::create_mirror_view_and_copy(memory_space(), lens_host);
341 for (local_ordinal_type i = 1, iend = offs.extent(0); i < iend; ++i) {
342 offs(i) = offs(i - 1) + lens[i - 1];
347 void createMpiRequests(
const tpetra_import_type &
import) {
348 Tpetra::Distributor &distributor =
import.getDistributor();
351 const auto pids_from = distributor.getProcsFrom();
353 memcpy(pids.recv.data(), pids_from.getRawPtr(),
sizeof(int) * pids.recv.extent(0));
355 const auto pids_to = distributor.getProcsTo();
357 memcpy(pids.send.data(), pids_to.getRawPtr(),
sizeof(int) * pids.send.extent(0));
360 reqs.recv.resize(pids.recv.extent(0));
361 memset(reqs.recv.data(), 0, reqs.recv.size() *
sizeof(MPI_Request));
362 reqs.send.resize(pids.send.extent(0));
363 memset(reqs.send.data(), 0, reqs.send.size() *
sizeof(MPI_Request));
367 const auto lengths_to = distributor.getLengthsTo();
370 const auto lengths_from = distributor.getLengthsFrom();
373 setOffsetValues(lengths_to, offset.send);
374 offset_host.send = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), offset.send);
376 setOffsetValues(lengths_from, offset.recv);
377 offset_host.recv = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), offset.recv);
379 const auto lengths_to = distributor.getLengthsTo();
380 offset_host.send = size_type_1d_view_host(
do_not_initialize_tag(
"offset send"), lengths_to.size() + 1);
382 const auto lengths_from = distributor.getLengthsFrom();
383 offset_host.recv = size_type_1d_view_host(
do_not_initialize_tag(
"offset recv"), lengths_from.size() + 1);
385 setOffsetValuesHost(lengths_to, offset_host.send);
388 setOffsetValuesHost(lengths_from, offset_host.recv);
393 void createSendRecvIDs(
const tpetra_import_type &
import) {
395 const auto remote_lids =
import.getRemoteLIDs();
396 const local_ordinal_type_1d_view_host
397 remote_lids_view_host(const_cast<local_ordinal_type *>(remote_lids.getRawPtr()), remote_lids.size());
399 Kokkos::deep_copy(lids.recv, remote_lids_view_host);
402 auto epids =
import.getExportPIDs();
403 auto elids =
import.getExportLIDs();
406 auto lids_send_host = Kokkos::create_mirror_view(lids.send);
409 for (local_ordinal_type cnt = 0, i = 0, iend = pids.send.extent(0); i < iend; ++i) {
410 const auto pid_send_value = pids.send[i];
411 for (local_ordinal_type j = 0, jend = epids.size(); j < jend; ++j)
412 if (epids[j] == pid_send_value) lids_send_host[cnt++] = elids[j];
413 TEUCHOS_ASSERT(static_cast<size_t>(cnt) == offset_host.send[i + 1]);
415 Kokkos::deep_copy(lids.send, lids_send_host);
418 void createExecutionSpaceInstances() {
419 #if defined(KOKKOS_ENABLE_CUDA) || defined(KOKKOS_ENABLE_HIP) || defined(KOKKOS_ENABLE_SYCL)
421 #if KOKKOS_VERSION >= 40699
423 Kokkos::Experimental::partition_space(execution_space(), std::vector<int>(8, 1));
426 Kokkos::Experimental::partition_space(execution_space(), 1, 1, 1, 1, 1, 1, 1, 1);
434 struct ToMultiVector {};
438 const local_ordinal_type blocksize_,
439 const local_ordinal_type_1d_view dm2cm_) {
440 blocksize = blocksize_;
443 #ifdef HAVE_IFPACK2_MPI
444 comm = Tpetra::Details::extractMpiCommFromTeuchos(*tgt_map->getComm());
446 const tpetra_import_type
import(src_map, tgt_map);
448 createMpiRequests(
import);
449 createSendRecvIDs(
import);
450 createExecutionSpaceInstances();
453 void createDataBuffer(
const local_ordinal_type &num_vectors) {
454 const size_type extent_0 = lids.recv.extent(0) * blocksize;
455 const size_type extent_1 = num_vectors;
456 if (remote_multivector.extent(0) == extent_0 &&
457 remote_multivector.extent(1) == extent_1) {
463 const auto send_buffer_size = offset_host.send[offset_host.send.extent(0) - 1] * blocksize * num_vectors;
464 const auto recv_buffer_size = offset_host.recv[offset_host.recv.extent(0) - 1] * blocksize * num_vectors;
469 if (!Tpetra::Details::Behavior::assumeMpiIsGPUAware()) {
470 buffer_host.send = impl_scalar_type_1d_view_host(
do_not_initialize_tag(
"buffer send"), send_buffer_size);
471 buffer_host.recv = impl_scalar_type_1d_view_host(
do_not_initialize_tag(
"buffer recv"), recv_buffer_size);
477 #ifdef HAVE_IFPACK2_MPI
478 waitall(reqs.recv.size(), reqs.recv.data());
479 waitall(reqs.send.size(), reqs.send.data());
487 #if defined(KOKKOS_ENABLE_CUDA) || defined(KOKKOS_ENABLE_HIP) || defined(KOKKOS_ENABLE_SYCL)
488 template <
typename PackTag>
489 static void copy(
const local_ordinal_type_1d_view &lids_,
490 const impl_scalar_type_1d_view &buffer_,
491 const local_ordinal_type ibeg_,
492 const local_ordinal_type iend_,
493 const impl_scalar_type_2d_view_tpetra &multivector_,
494 const local_ordinal_type blocksize_,
495 const execution_space &exec_instance_) {
496 const local_ordinal_type num_vectors = multivector_.extent(1);
497 const local_ordinal_type mv_blocksize = blocksize_ * num_vectors;
498 const local_ordinal_type idiff = iend_ - ibeg_;
499 const auto abase = buffer_.data() + mv_blocksize * ibeg_;
501 using team_policy_type = Kokkos::TeamPolicy<execution_space>;
502 local_ordinal_type vector_size(0);
505 else if (blocksize_ <= 8)
507 else if (blocksize_ <= 16)
512 const auto work_item_property = Kokkos::Experimental::WorkItemProperty::HintLightWeight;
513 const team_policy_type policy(exec_instance_, idiff, 1, vector_size);
514 Kokkos::parallel_for(
515 Kokkos::Experimental::require(policy, work_item_property),
516 KOKKOS_LAMBDA(
const typename team_policy_type::member_type &member) {
517 const local_ordinal_type i = member.league_rank();
518 Kokkos::parallel_for(Kokkos::TeamThreadRange(member, num_vectors), [&](
const local_ordinal_type &j) {
519 auto aptr = abase + blocksize_ * (i + idiff * j);
520 auto bptr = &multivector_(blocksize_ * lids_(i + ibeg_), j);
521 if (std::is_same<PackTag, ToBuffer>::value)
522 Kokkos::parallel_for(Kokkos::ThreadVectorRange(member, blocksize_), [&](
const local_ordinal_type &k) {
526 Kokkos::parallel_for(Kokkos::ThreadVectorRange(member, blocksize_), [&](
const local_ordinal_type &k) {
533 void asyncSendRecvVar1(
const impl_scalar_type_2d_view_tpetra &mv) {
534 IFPACK2_BLOCKHELPER_TIMER(
"BlockTriDi::AsyncableImport::AsyncSendRecv", AsyncSendRecv);
536 #ifdef HAVE_IFPACK2_MPI
538 const local_ordinal_type num_vectors = mv.extent(1);
539 const local_ordinal_type mv_blocksize = blocksize * num_vectors;
542 for (local_ordinal_type i = 0, iend = pids.recv.extent(0); i < iend; ++i) {
543 if (Tpetra::Details::Behavior::assumeMpiIsGPUAware()) {
545 reinterpret_cast<char *>(buffer.recv.data() + offset_host.recv[i] * mv_blocksize),
546 (offset_host.recv[i + 1] - offset_host.recv[i]) * mv_blocksize *
sizeof(impl_scalar_type),
552 reinterpret_cast<char *>(buffer_host.recv.data() + offset_host.recv[i] * mv_blocksize),
553 (offset_host.recv[i + 1] - offset_host.recv[i]) * mv_blocksize *
sizeof(impl_scalar_type),
561 execution_space().fence();
564 for (local_ordinal_type i = 0; i < static_cast<local_ordinal_type>(pids.send.extent(0)); ++i) {
566 if (i < 8) exec_instances[i % 8].fence();
567 copy<ToBuffer>(lids.send, buffer.send,
568 offset_host.send(i), offset_host.send(i + 1),
571 exec_instances[i % 8]);
572 if (!Tpetra::Details::Behavior::assumeMpiIsGPUAware()) {
574 const local_ordinal_type num_vectors = mv.extent(1);
575 const local_ordinal_type mv_blocksize = blocksize * num_vectors;
577 Kokkos::deep_copy(exec_instances[i % 8],
578 Kokkos::subview(buffer_host.send,
579 Kokkos::pair<local_ordinal_type, local_ordinal_type>(
580 offset_host.send(i) * mv_blocksize,
581 offset_host.send(i + 1) * mv_blocksize)),
582 Kokkos::subview(buffer.send,
583 Kokkos::pair<local_ordinal_type, local_ordinal_type>(
584 offset_host.send(i) * mv_blocksize,
585 offset_host.send(i + 1) * mv_blocksize)));
590 for (local_ordinal_type i = 0; i < static_cast<local_ordinal_type>(pids.send.extent(0)); ++i) {
592 if (i < 8) exec_instances[i % 8].fence();
593 if (Tpetra::Details::Behavior::assumeMpiIsGPUAware()) {
595 reinterpret_cast<const char *>(buffer.send.data() + offset_host.send[i] * mv_blocksize),
596 (offset_host.send[i + 1] - offset_host.send[i]) * mv_blocksize *
sizeof(impl_scalar_type),
602 reinterpret_cast<const char *>(buffer_host.send.data() + offset_host.send[i] * mv_blocksize),
603 (offset_host.send[i + 1] - offset_host.send[i]) * mv_blocksize *
sizeof(impl_scalar_type),
611 for (local_ordinal_type i = 0, iend = pids.recv.extent(0); i < iend; ++i) {
614 MPI_Iprobe(pids.recv[i], 42, comm, &flag, &stat);
616 #endif // HAVE_IFPACK2_MPI
617 IFPACK2_BLOCKHELPER_TIMER_FENCE(execution_space)
620 void syncRecvVar1() {
621 IFPACK2_BLOCKHELPER_TIMER(
"BlockTriDi::AsyncableImport::SyncRecv", SyncRecv);
622 #ifdef HAVE_IFPACK2_MPI
624 for (local_ordinal_type i = 0; i < static_cast<local_ordinal_type>(pids.recv.extent(0)); ++i) {
625 local_ordinal_type idx = i;
628 waitany(pids.recv.extent(0), reqs.recv.data(), &idx);
630 if (!Tpetra::Details::Behavior::assumeMpiIsGPUAware()) {
631 const local_ordinal_type num_vectors = remote_multivector.extent(1);
632 const local_ordinal_type mv_blocksize = blocksize * num_vectors;
635 Kokkos::subview(buffer.recv,
636 Kokkos::pair<local_ordinal_type, local_ordinal_type>(
637 offset_host.recv(idx) * mv_blocksize,
638 offset_host.recv(idx + 1) * mv_blocksize)),
639 Kokkos::subview(buffer_host.recv,
640 Kokkos::pair<local_ordinal_type, local_ordinal_type>(
641 offset_host.recv(idx) * mv_blocksize,
642 offset_host.recv(idx + 1) * mv_blocksize)));
646 copy<ToMultiVector>(lids.recv, buffer.recv,
647 offset_host.recv(idx), offset_host.recv(idx + 1),
648 remote_multivector, blocksize,
649 exec_instances[idx % 8]);
656 waitall(reqs.send.size(), reqs.send.data());
657 #endif // HAVE_IFPACK2_MPI
658 IFPACK2_BLOCKHELPER_TIMER_FENCE(execution_space)
660 #endif // defined(KOKKOS_ENABLE_CUDA|HIP|SYCL)
667 template <
typename PackTag>
668 static void copy(
const local_ordinal_type_1d_view &lids_,
669 const impl_scalar_type_1d_view &buffer_,
670 const local_ordinal_type &ibeg_,
671 const local_ordinal_type &iend_,
672 const impl_scalar_type_2d_view_tpetra &multivector_,
673 const local_ordinal_type blocksize_) {
674 const local_ordinal_type num_vectors = multivector_.extent(1);
675 const local_ordinal_type mv_blocksize = blocksize_ * num_vectors;
676 const local_ordinal_type idiff = iend_ - ibeg_;
677 const auto abase = buffer_.data() + mv_blocksize * ibeg_;
678 if constexpr (BlockHelperDetails::is_device<execution_space>::value) {
679 using team_policy_type = Kokkos::TeamPolicy<execution_space>;
680 local_ordinal_type vector_size(0);
683 else if (blocksize_ <= 8)
685 else if (blocksize_ <= 16)
689 const team_policy_type policy(idiff, 1, vector_size);
690 Kokkos::parallel_for(
691 "AsyncableImport::TeamPolicy::copy",
692 policy, KOKKOS_LAMBDA(
const typename team_policy_type::member_type &member) {
693 const local_ordinal_type i = member.league_rank();
694 Kokkos::parallel_for(Kokkos::TeamThreadRange(member, num_vectors), [&](
const local_ordinal_type &j) {
695 auto aptr = abase + blocksize_ * (i + idiff * j);
696 auto bptr = &multivector_(blocksize_ * lids_(i + ibeg_), j);
697 if (std::is_same<PackTag, ToBuffer>::value)
698 Kokkos::parallel_for(Kokkos::ThreadVectorRange(member, blocksize_), [&](
const local_ordinal_type &k) {
702 Kokkos::parallel_for(Kokkos::ThreadVectorRange(member, blocksize_), [&](
const local_ordinal_type &k) {
708 const Kokkos::RangePolicy<execution_space> policy(0, idiff * num_vectors);
709 Kokkos::parallel_for(
710 "AsyncableImport::RangePolicy::copy",
711 policy, KOKKOS_LAMBDA(
const local_ordinal_type &ij) {
712 const local_ordinal_type i = ij % idiff;
713 const local_ordinal_type j = ij / idiff;
714 auto aptr = abase + blocksize_ * (i + idiff * j);
715 auto bptr = &multivector_(blocksize_ * lids_(i + ibeg_), j);
716 auto from = std::is_same<PackTag, ToBuffer>::value ? bptr : aptr;
717 auto to = std::is_same<PackTag, ToBuffer>::value ? aptr : bptr;
718 memcpy(to, from,
sizeof(impl_scalar_type) * blocksize_);
726 void asyncSendRecvVar0(
const impl_scalar_type_2d_view_tpetra &mv) {
727 IFPACK2_BLOCKHELPER_TIMER(
"BlockTriDi::AsyncableImport::AsyncSendRecv", AsyncSendRecv);
729 #ifdef HAVE_IFPACK2_MPI
731 const local_ordinal_type num_vectors = mv.extent(1);
732 const local_ordinal_type mv_blocksize = blocksize * num_vectors;
735 for (local_ordinal_type i = 0, iend = pids.recv.extent(0); i < iend; ++i) {
736 if (Tpetra::Details::Behavior::assumeMpiIsGPUAware()) {
738 reinterpret_cast<char *>(buffer.recv.data() + offset_host.recv[i] * mv_blocksize),
739 (offset_host.recv[i + 1] - offset_host.recv[i]) * mv_blocksize *
sizeof(impl_scalar_type),
745 reinterpret_cast<char *>(buffer_host.recv.data() + offset_host.recv[i] * mv_blocksize),
746 (offset_host.recv[i + 1] - offset_host.recv[i]) * mv_blocksize *
sizeof(impl_scalar_type),
754 for (local_ordinal_type i = 0, iend = pids.send.extent(0); i < iend; ++i) {
755 copy<ToBuffer>(lids.send, buffer.send, offset_host.send(i), offset_host.send(i + 1),
758 if (Tpetra::Details::Behavior::assumeMpiIsGPUAware()) {
760 reinterpret_cast<const char *>(buffer.send.data() + offset_host.send[i] * mv_blocksize),
761 (offset_host.send[i + 1] - offset_host.send[i]) * mv_blocksize *
sizeof(impl_scalar_type),
767 Kokkos::subview(buffer_host.send,
768 Kokkos::pair<local_ordinal_type, local_ordinal_type>(
769 offset_host.send(i) * mv_blocksize,
770 offset_host.send(i + 1) * mv_blocksize)),
771 Kokkos::subview(buffer.send,
772 Kokkos::pair<local_ordinal_type, local_ordinal_type>(
773 offset_host.send(i) * mv_blocksize,
774 offset_host.send(i + 1) * mv_blocksize)));
776 reinterpret_cast<const char *>(buffer_host.send.data() + offset_host.send[i] * mv_blocksize),
777 (offset_host.send[i + 1] - offset_host.send[i]) * mv_blocksize *
sizeof(impl_scalar_type),
786 for (local_ordinal_type i = 0, iend = pids.recv.extent(0); i < iend; ++i) {
789 MPI_Iprobe(pids.recv[i], 42, comm, &flag, &stat);
792 IFPACK2_BLOCKHELPER_TIMER_FENCE(execution_space)
795 void syncRecvVar0() {
796 IFPACK2_BLOCKHELPER_TIMER(
"BlockTriDi::AsyncableImport::SyncRecv", SyncRecv);
797 #ifdef HAVE_IFPACK2_MPI
799 for (local_ordinal_type i = 0, iend = pids.recv.extent(0); i < iend; ++i) {
800 local_ordinal_type idx = i;
801 waitany(pids.recv.extent(0), reqs.recv.data(), &idx);
802 if (!Tpetra::Details::Behavior::assumeMpiIsGPUAware()) {
803 const local_ordinal_type num_vectors = remote_multivector.extent(1);
804 const local_ordinal_type mv_blocksize = blocksize * num_vectors;
806 Kokkos::subview(buffer.recv,
807 Kokkos::pair<local_ordinal_type, local_ordinal_type>(
808 offset_host.recv(idx) * mv_blocksize,
809 offset_host.recv(idx + 1) * mv_blocksize)),
810 Kokkos::subview(buffer_host.recv,
811 Kokkos::pair<local_ordinal_type, local_ordinal_type>(
812 offset_host.recv(idx) * mv_blocksize,
813 offset_host.recv(idx + 1) * mv_blocksize)));
815 copy<ToMultiVector>(lids.recv, buffer.recv, offset_host.recv(idx), offset_host.recv(idx + 1),
816 remote_multivector, blocksize);
819 waitall(reqs.send.size(), reqs.send.data());
821 IFPACK2_BLOCKHELPER_TIMER_FENCE(execution_space)
827 void asyncSendRecv(
const impl_scalar_type_2d_view_tpetra &mv) {
828 #if defined(KOKKOS_ENABLE_CUDA) || defined(KOKKOS_ENABLE_HIP) || defined(KOKKOS_ENABLE_SYCL)
829 #if defined(IFPACK2_BLOCKTRIDICONTAINER_USE_EXEC_SPACE_INSTANCES)
830 asyncSendRecvVar1(mv);
832 asyncSendRecvVar0(mv);
835 asyncSendRecvVar0(mv);
839 #if defined(KOKKOS_ENABLE_CUDA) || defined(KOKKOS_ENABLE_HIP) || defined(KOKKOS_ENABLE_SYCL)
840 #if defined(IFPACK2_BLOCKTRIDICONTAINER_USE_EXEC_SPACE_INSTANCES)
850 void syncExchange(
const impl_scalar_type_2d_view_tpetra &mv) {
851 IFPACK2_BLOCKHELPER_TIMER(
"BlockTriDi::AsyncableImport::SyncExchange", SyncExchange);
854 IFPACK2_BLOCKHELPER_TIMER_FENCE(execution_space)
857 impl_scalar_type_2d_view_tpetra getRemoteMultiVectorLocalView()
const {
return remote_multivector; }
860 template <
typename ViewType1,
typename ViewType2>
861 struct are_same_struct {
865 are_same_struct(ViewType1 keys1_, ViewType2 keys2_)
868 KOKKOS_INLINE_FUNCTION
869 void operator()(
int i,
unsigned int &count)
const {
870 if (keys1(i) != keys2(i)) count++;
874 template <
typename ViewType1,
typename ViewType2>
875 bool are_same(ViewType1 keys1, ViewType2 keys2) {
876 unsigned int are_same_ = 0;
878 Kokkos::parallel_reduce(Kokkos::RangePolicy<typename ViewType1::execution_space>(0, keys1.extent(0)),
879 are_same_struct(keys1, keys2),
881 return are_same_ == 0;
887 template <
typename MatrixType>
892 using tpetra_map_type =
typename impl_type::tpetra_map_type;
893 using local_ordinal_type =
typename impl_type::local_ordinal_type;
894 using global_ordinal_type =
typename impl_type::global_ordinal_type;
895 using local_ordinal_type_1d_view =
typename impl_type::local_ordinal_type_1d_view;
896 using crs_matrix_type =
typename impl_type::tpetra_crs_matrix_type;
897 using block_crs_matrix_type =
typename impl_type::tpetra_block_crs_matrix_type;
898 using global_indices_array_device_type = Kokkos::View<const global_ordinal_type *, typename tpetra_map_type::device_type>;
900 auto A_crs = Teuchos::rcp_dynamic_cast<
const crs_matrix_type>(A);
901 auto A_bcrs = Teuchos::rcp_dynamic_cast<
const block_crs_matrix_type>(A);
903 bool hasBlockCrsMatrix = !A_bcrs.is_null();
906 const auto g = hasBlockCrsMatrix ? A_bcrs->getCrsGraph() : *(A_crs->getCrsGraph());
908 const auto blocksize = hasBlockCrsMatrix ? A_bcrs->getBlockSize() : 1;
909 const auto domain_map = g.getDomainMap();
910 const auto column_map = g.getColMap();
912 std::vector<global_ordinal_type> gids;
914 Kokkos::Subview<global_indices_array_device_type, std::pair<int, int>> column_map_global_iD_last;
916 bool separate_remotes =
true, found_first =
false, need_owned_permutation =
false;
918 IFPACK2_BLOCKHELPER_TIMER(
"createBlockCrsAsyncImporter::loop_over_local_elements", loop_over_local_elements);
920 global_indices_array_device_type column_map_global_iD = column_map->getMyGlobalIndicesDevice();
921 global_indices_array_device_type domain_map_global_iD = domain_map->getMyGlobalIndicesDevice();
923 if (are_same(domain_map_global_iD, column_map_global_iD)) {
925 separate_remotes =
true;
926 need_owned_permutation =
false;
928 column_map_global_iD_last = Kokkos::subview(column_map_global_iD,
929 std::pair<int, int>(domain_map_global_iD.extent(0), column_map_global_iD.extent(0)));
932 for (
size_t i = 0; i < column_map->getLocalNumElements(); ++i) {
933 const global_ordinal_type gid = column_map->getGlobalElement(i);
934 if (!domain_map->isNodeGlobalElement(gid)) {
937 }
else if (found_first) {
938 separate_remotes =
false;
941 if (!found_first && !need_owned_permutation &&
942 domain_map->getLocalElement(gid) !=
static_cast<local_ordinal_type
>(i)) {
951 need_owned_permutation =
true;
955 IFPACK2_BLOCKHELPER_TIMER_FENCE(
typename BlockHelperDetails::ImplType<MatrixType>::execution_space)
958 if (separate_remotes) {
959 IFPACK2_BLOCKHELPER_TIMER(
"createBlockCrsAsyncImporter::separate_remotes", separate_remotes);
961 const auto parsimonious_col_map = need_owned_permutation ?
Teuchos::rcp(
new tpetra_map_type(invalid, gids.data(), gids.size(), 0, domain_map->getComm())) :
Teuchos::rcp(
new tpetra_map_type(invalid, column_map_global_iD_last, 0, domain_map->getComm()));
962 if (parsimonious_col_map->getGlobalNumElements() > 0) {
964 local_ordinal_type_1d_view dm2cm;
965 if (need_owned_permutation) {
966 dm2cm = local_ordinal_type_1d_view(
do_not_initialize_tag(
"dm2cm"), domain_map->getLocalNumElements());
967 const auto dm2cm_host = Kokkos::create_mirror_view(dm2cm);
968 for (
size_t i = 0; i < domain_map->getLocalNumElements(); ++i)
969 dm2cm_host(i) = domain_map->getLocalElement(column_map->getGlobalElement(i));
970 Kokkos::deep_copy(dm2cm, dm2cm_host);
972 IFPACK2_BLOCKHELPER_TIMER_FENCE(
typename BlockHelperDetails::ImplType<MatrixType>::execution_space)
973 return Teuchos::rcp(
new AsyncableImport<MatrixType>(domain_map, parsimonious_col_map, blocksize, dm2cm));
976 IFPACK2_BLOCKHELPER_TIMER_FENCE(
typename BlockHelperDetails::ImplType<MatrixType>::execution_space)
977 return Teuchos::null;
980 template <
typename local_ordinal_type>
981 local_ordinal_type costTRSM(
const local_ordinal_type block_size) {
982 return block_size * block_size;
985 template <
typename local_ordinal_type>
986 local_ordinal_type costGEMV(
const local_ordinal_type block_size) {
987 return 2 * block_size * block_size;
990 template <
typename local_ordinal_type>
991 local_ordinal_type costTriDiagSolve(
const local_ordinal_type subline_length,
const local_ordinal_type block_size) {
992 return 2 * subline_length * costTRSM(block_size) + 2 * (subline_length - 1) * costGEMV(block_size);
995 template <
typename local_ordinal_type>
996 local_ordinal_type costSolveSchur(
const local_ordinal_type num_parts,
997 const local_ordinal_type num_teams,
998 const local_ordinal_type line_length,
999 const local_ordinal_type block_size,
1000 const local_ordinal_type n_subparts_per_part) {
1001 const local_ordinal_type subline_length = ceil(
double(line_length - (n_subparts_per_part - 1) * 2) / n_subparts_per_part);
1002 if (subline_length < 1) {
1006 const local_ordinal_type p_n_lines = ceil(
double(num_parts) / num_teams);
1007 const local_ordinal_type p_n_sublines = ceil(
double(n_subparts_per_part) * num_parts / num_teams);
1008 const local_ordinal_type p_n_sublines_2 = ceil(
double(n_subparts_per_part - 1) * num_parts / num_teams);
1010 const local_ordinal_type p_costApplyE = p_n_sublines_2 * subline_length * 2 * costGEMV(block_size);
1011 const local_ordinal_type p_costApplyS = p_n_lines * costTriDiagSolve((n_subparts_per_part - 1) * 2, block_size);
1012 const local_ordinal_type p_costApplyAinv = p_n_sublines * costTriDiagSolve(subline_length, block_size);
1013 const local_ordinal_type p_costApplyC = p_n_sublines_2 * 2 * costGEMV(block_size);
1015 if (n_subparts_per_part == 1) {
1016 return p_costApplyAinv;
1018 return p_costApplyE + p_costApplyS + p_costApplyAinv + p_costApplyC;
1021 template <
typename local_ordinal_type>
1022 local_ordinal_type getAutomaticNSubparts(
const local_ordinal_type num_parts,
1023 const local_ordinal_type num_teams,
1024 const local_ordinal_type line_length,
1025 const local_ordinal_type block_size) {
1032 double parallelismSurplus = Kokkos::sqrt((
double)num_teams / num_parts);
1033 double logLineLength = Kokkos::log2((
double)line_length);
1034 (void)logLineLength;
1036 #if defined(KOKKOS_ARCH_AMD_GFX942) || defined(KOKKOS_ARCH_AMD_GFX942_APU)
1038 double modeled = -9.2312 + 4.6946 * parallelismSurplus + 0.4095 * block_size + 0.966 * logLineLength;
1040 if (parallelismSurplus < 0.3)
1042 #elif defined(KOKKOS_ARCH_HOPPER) || defined(KOKKOS_ARCH_BLACKWELL)
1044 double modeled = -9.6053 + 4.7477 * parallelismSurplus + 0.2338 * block_size + 1.0794 * logLineLength;
1046 double maxSplit = (double)line_length / 8;
1047 if (modeled > maxSplit)
1049 #elif defined(KOKKOS_ENABLE_CUDA)
1053 if (parallelismSurplus > 1 && line_length > 64)
1055 #elif defined(KOKKOS_ENABLE_HIP)
1057 double modeled = -8.6214 + 7.3468 * parallelismSurplus + 0.3596 * block_size + 0.6673 * logLineLength;
1061 if (parallelismSurplus > 1 && line_length > 64)
1066 local_ordinal_type n_subparts_per_part = 0.5 + modeled;
1068 if (parallelismSurplus < 0.3)
1069 n_subparts_per_part = 1;
1076 local_ordinal_type min_subparts_per_part = 1;
1077 local_ordinal_type max_subparts_per_part = (line_length + 2) / 3;
1079 if (max_subparts_per_part > 16)
1080 max_subparts_per_part = 16;
1081 if (n_subparts_per_part < min_subparts_per_part)
1082 n_subparts_per_part = min_subparts_per_part;
1083 if (n_subparts_per_part > max_subparts_per_part)
1084 n_subparts_per_part = max_subparts_per_part;
1085 return n_subparts_per_part;
1088 template <
typename ArgActiveExecutionMemorySpace>
1089 struct SolveTridiagsDefaultModeAndAlgo;
1094 template <
typename MatrixType>
1095 BlockHelperDetails::PartInterface<MatrixType>
1097 const Teuchos::RCP<
const typename BlockHelperDetails::ImplType<MatrixType>::tpetra_crs_graph_type> &G,
1099 const typename BlockHelperDetails::ImplType<MatrixType>::local_ordinal_type n_subparts_per_part_in) {
1102 using local_ordinal_type =
typename impl_type::local_ordinal_type;
1103 using local_ordinal_type_1d_view =
typename impl_type::local_ordinal_type_1d_view;
1104 using local_ordinal_type_2d_view =
typename impl_type::local_ordinal_type_2d_view;
1105 using size_type =
typename impl_type::size_type;
1107 auto bA = Teuchos::rcp_dynamic_cast<
const typename BlockHelperDetails::ImplType<MatrixType>::tpetra_block_crs_matrix_type>(A);
1110 const local_ordinal_type blocksize = bA.is_null() ? A->getLocalNumRows() / G->getLocalNumRows() : A->getBlockSize();
1111 constexpr
int vector_length = impl_type::vector_length;
1112 constexpr
int internal_vector_length = impl_type::internal_vector_length;
1114 const auto comm = A->getRowMap()->getComm();
1116 BlockHelperDetails::PartInterface<MatrixType> interf;
1118 const bool jacobi = partitions.size() == 0;
1119 const local_ordinal_type A_n_lclrows = G->getLocalNumRows();
1120 const local_ordinal_type nparts = jacobi ? A_n_lclrows : partitions.size();
1122 typedef std::pair<local_ordinal_type, local_ordinal_type> size_idx_pair_type;
1123 std::vector<size_idx_pair_type> partsz(nparts);
1126 for (local_ordinal_type i = 0; i < nparts; ++i)
1127 partsz[i] = size_idx_pair_type(partitions[i].size(), i);
1128 std::sort(partsz.begin(), partsz.end(),
1129 [](
const size_idx_pair_type &x,
const size_idx_pair_type &y) {
1130 return x.first > y.first;
1134 local_ordinal_type n_subparts_per_part;
1135 if (n_subparts_per_part_in == -1) {
1138 using execution_space =
typename impl_type::execution_space;
1141 if constexpr (impl_type::node_type::is_gpu) {
1142 const int line_length = partsz[0].first;
1144 const local_ordinal_type team_size =
1145 SolveTridiagsDefaultModeAndAlgo<typename execution_space::memory_space>::
1146 recommended_team_size(blocksize, vector_length, internal_vector_length);
1148 const local_ordinal_type num_teams = std::max(1, execution_space().concurrency() / (team_size * vector_length));
1149 n_subparts_per_part = getAutomaticNSubparts(nparts, num_teams, line_length, blocksize);
1150 #ifdef IFPACK2_BLOCKTRIDICONTAINER_USE_PRINTF
1151 printf(
"Automatically chosen n_subparts_per_part = %d for nparts = %d, num_teams = %d, team_size = %d, line_length = %d, and blocksize = %d;\n", n_subparts_per_part, nparts, num_teams, team_size, line_length, blocksize);
1154 n_subparts_per_part = 1;
1155 #ifdef IFPACK2_BLOCKTRIDICONTAINER_USE_PRINTF
1156 printf(
"Automatically chosen n_subparts_per_part = 1 for CPU backend\n");
1160 n_subparts_per_part = n_subparts_per_part_in;
1164 const local_ordinal_type n_sub_parts = nparts * n_subparts_per_part;
1167 const local_ordinal_type n_sub_parts_and_schur = n_sub_parts + nparts * (n_subparts_per_part - 1);
1169 #if defined(BLOCKTRIDICONTAINER_DEBUG)
1170 local_ordinal_type nrows = 0;
1174 for (local_ordinal_type i = 0; i < nparts; ++i) nrows += partitions[i].size();
1176 TEUCHOS_TEST_FOR_EXCEPT_MSG(nrows != A_n_lclrows, BlockHelperDetails::get_msg_prefix(comm) <<
"The #rows implied by the local partition is not "
1177 <<
"the same as getLocalNumRows: " << nrows <<
" vs " << A_n_lclrows);
1181 std::vector<local_ordinal_type> p;
1183 interf.max_partsz = 1;
1184 interf.max_subpartsz = 0;
1185 interf.n_subparts_per_part = 1;
1186 interf.nparts = nparts;
1191 for (local_ordinal_type i = 0; i < nparts; ++i)
1192 p[i] = partsz[i].second;
1194 interf.max_partsz = partsz[0].first;
1196 constexpr local_ordinal_type connection_length = 2;
1197 const local_ordinal_type sub_line_length = (interf.max_partsz - (n_subparts_per_part - 1) * connection_length) / n_subparts_per_part;
1198 const local_ordinal_type last_sub_line_length = interf.max_partsz - (n_subparts_per_part - 1) * (connection_length + sub_line_length);
1200 interf.max_subpartsz = (sub_line_length > last_sub_line_length) ? sub_line_length : last_sub_line_length;
1201 interf.n_subparts_per_part = n_subparts_per_part;
1202 interf.nparts = nparts;
1208 interf.part2rowidx0 = local_ordinal_type_1d_view(
do_not_initialize_tag(
"part2rowidx0"), nparts + 1);
1209 interf.part2packrowidx0 = local_ordinal_type_1d_view(
do_not_initialize_tag(
"part2packrowidx0"), nparts + 1);
1212 interf.part2rowidx0_sub = local_ordinal_type_1d_view(
do_not_initialize_tag(
"part2rowidx0_sub"), n_sub_parts_and_schur + 1);
1213 interf.part2packrowidx0_sub = local_ordinal_type_2d_view(
do_not_initialize_tag(
"part2packrowidx0_sub"), nparts, 2 * n_subparts_per_part);
1214 interf.rowidx2part_sub = local_ordinal_type_1d_view(
do_not_initialize_tag(
"rowidx2part"), A_n_lclrows);
1216 interf.partptr_sub = local_ordinal_type_2d_view(
do_not_initialize_tag(
"partptr_sub"), n_sub_parts_and_schur, 2);
1219 const auto partptr = Kokkos::create_mirror_view(interf.partptr);
1220 const auto partptr_sub = Kokkos::create_mirror_view(interf.partptr_sub);
1222 const auto lclrow = Kokkos::create_mirror_view(interf.lclrow);
1223 const auto part2rowidx0 = Kokkos::create_mirror_view(interf.part2rowidx0);
1224 const auto part2packrowidx0 = Kokkos::create_mirror_view(interf.part2packrowidx0);
1225 const auto rowidx2part = Kokkos::create_mirror_view(interf.rowidx2part);
1227 const auto part2rowidx0_sub = Kokkos::create_mirror_view(interf.part2rowidx0_sub);
1228 const auto part2packrowidx0_sub = Kokkos::create_mirror_view(Kokkos::HostSpace(), interf.part2packrowidx0_sub);
1229 const auto rowidx2part_sub = Kokkos::create_mirror_view(interf.rowidx2part_sub);
1232 interf.row_contiguous =
true;
1234 part2rowidx0(0) = 0;
1235 part2packrowidx0(0) = 0;
1236 local_ordinal_type pack_nrows = 0;
1237 local_ordinal_type pack_nrows_sub = 0;
1239 IFPACK2_BLOCKHELPER_TIMER(
"compute part indices (Jacobi)", Jacobi);
1243 for (local_ordinal_type i = 0; i <= nparts; ++i) {
1244 part2rowidx0(i) = i;
1247 for (local_ordinal_type i = 0; i < nparts; ++i) {
1251 for (local_ordinal_type ip = 0; ip < nparts; ++ip) {
1253 if (ip % vector_length == 0) pack_nrows = 1;
1254 part2packrowidx0(ip + 1) = part2packrowidx0(ip) + ((ip + 1) % vector_length == 0 || ip + 1 == nparts ? pack_nrows : 0);
1256 part2rowidx0_sub(0) = 0;
1257 partptr_sub(0, 0) = 0;
1259 for (local_ordinal_type ip = 0; ip < nparts; ++ip) {
1260 constexpr local_ordinal_type ipnrows = 1;
1261 const local_ordinal_type full_line_length = partptr(ip + 1) - partptr(ip);
1264 "In the part " << ip);
1266 constexpr local_ordinal_type connection_length = 2;
1268 if (full_line_length < n_subparts_per_part + (n_subparts_per_part - 1) * connection_length)
1270 "The part " << ip <<
" is too short to use " << n_subparts_per_part <<
" sub parts.");
1272 const local_ordinal_type sub_line_length = (full_line_length - (n_subparts_per_part - 1) * connection_length) / n_subparts_per_part;
1273 const local_ordinal_type last_sub_line_length = full_line_length - (n_subparts_per_part - 1) * (connection_length + sub_line_length);
1275 if (ip % vector_length == 0) pack_nrows_sub = ipnrows;
1277 for (local_ordinal_type local_sub_ip = 0; local_sub_ip < n_subparts_per_part; ++local_sub_ip) {
1278 const local_ordinal_type sub_ip = nparts * (2 * local_sub_ip) + ip;
1279 const local_ordinal_type schur_ip = nparts * (2 * local_sub_ip + 1) + ip;
1280 if (local_sub_ip != n_subparts_per_part - 1) {
1281 if (local_sub_ip != 0) {
1282 partptr_sub(sub_ip, 0) = partptr_sub(nparts * (2 * local_sub_ip - 1) + ip, 1);
1283 }
else if (ip != 0) {
1284 partptr_sub(sub_ip, 0) = partptr_sub(nparts * 2 * (n_subparts_per_part - 1) + ip - 1, 1);
1286 partptr_sub(sub_ip, 1) = sub_line_length + partptr_sub(sub_ip, 0);
1287 partptr_sub(schur_ip, 0) = partptr_sub(sub_ip, 1);
1288 partptr_sub(schur_ip, 1) = connection_length + partptr_sub(schur_ip, 0);
1290 part2rowidx0_sub(sub_ip + 1) = part2rowidx0_sub(sub_ip) + sub_line_length;
1291 part2rowidx0_sub(sub_ip + 2) = part2rowidx0_sub(sub_ip + 1) + connection_length;
1293 #ifdef IFPACK2_BLOCKTRIDICONTAINER_USE_PRINTF
1294 printf(
"Sub Part index = %d, first LID associated to the sub part = %d, sub part size = %d;\n", sub_ip, partptr_sub(ip, 2 * local_sub_ip), sub_line_length);
1295 printf(
"Sub Part index Schur = %d, first LID associated to the sub part = %d, sub part size = %d;\n", sub_ip + 1, partptr_sub(ip, 2 * local_sub_ip + 1), connection_length);
1298 if (local_sub_ip != 0) {
1299 partptr_sub(sub_ip, 0) = partptr_sub(nparts * (2 * local_sub_ip - 1) + ip, 1);
1300 }
else if (ip != 0) {
1301 partptr_sub(sub_ip, 0) = partptr_sub(nparts * 2 * (n_subparts_per_part - 1) + ip - 1, 1);
1303 partptr_sub(sub_ip, 1) = last_sub_line_length + partptr_sub(sub_ip, 0);
1305 part2rowidx0_sub(sub_ip + 1) = part2rowidx0_sub(sub_ip) + last_sub_line_length;
1307 #ifdef IFPACK2_BLOCKTRIDICONTAINER_USE_PRINTF
1308 printf(
"Sub Part index = %d, first LID associated to the sub part = %d, sub part size = %d;\n", sub_ip, partptr_sub(ip, 2 * local_sub_ip), last_sub_line_length);
1314 #ifdef IFPACK2_BLOCKTRIDICONTAINER_WRITE_MM
1315 std::cout <<
"partptr_sub = " << std::endl;
1316 for (size_type i = 0; i < partptr_sub.extent(0); ++i) {
1317 for (size_type j = 0; j < partptr_sub.extent(1); ++j) {
1318 std::cout << partptr_sub(i, j) <<
" ";
1320 std::cout << std::endl;
1322 std::cout <<
"partptr_sub end" << std::endl;
1326 local_ordinal_type npacks = ceil(
float(nparts) / vector_length);
1328 local_ordinal_type ip_max = nparts > vector_length ? vector_length : nparts;
1329 for (local_ordinal_type ip = 0; ip < ip_max; ++ip) {
1330 part2packrowidx0_sub(ip, 0) = 0;
1332 for (local_ordinal_type ipack = 0; ipack < npacks; ++ipack) {
1334 local_ordinal_type ip_min = ipack * vector_length;
1335 ip_max = nparts > (ipack + 1) * vector_length ? (ipack + 1) * vector_length : nparts;
1336 for (local_ordinal_type ip = ip_min; ip < ip_max; ++ip) {
1337 part2packrowidx0_sub(ip, 0) = part2packrowidx0_sub(ip - vector_length, part2packrowidx0_sub.extent(1) - 1);
1341 for (size_type local_sub_ip = 0; local_sub_ip < part2packrowidx0_sub.extent(1) - 1; ++local_sub_ip) {
1342 local_ordinal_type ip_min = ipack * vector_length;
1343 ip_max = nparts > (ipack + 1) * vector_length ? (ipack + 1) * vector_length : nparts;
1345 const local_ordinal_type full_line_length = partptr(ip_min + 1) - partptr(ip_min);
1347 constexpr local_ordinal_type connection_length = 2;
1349 const local_ordinal_type sub_line_length = (full_line_length - (n_subparts_per_part - 1) * connection_length) / n_subparts_per_part;
1350 const local_ordinal_type last_sub_line_length = full_line_length - (n_subparts_per_part - 1) * (connection_length + sub_line_length);
1352 if (local_sub_ip % 2 == 0) pack_nrows_sub = sub_line_length;
1353 if (local_sub_ip % 2 == 1) pack_nrows_sub = connection_length;
1354 if (local_sub_ip == part2packrowidx0_sub.extent(1) - 2) pack_nrows_sub = last_sub_line_length;
1356 part2packrowidx0_sub(ip_min, local_sub_ip + 1) = part2packrowidx0_sub(ip_min, local_sub_ip) + pack_nrows_sub;
1358 for (local_ordinal_type ip = ip_min + 1; ip < ip_max; ++ip) {
1359 part2packrowidx0_sub(ip, local_sub_ip + 1) = part2packrowidx0_sub(ip_min, local_sub_ip + 1);
1364 Kokkos::deep_copy(interf.part2packrowidx0_sub, part2packrowidx0_sub);
1366 IFPACK2_BLOCKHELPER_TIMER_FENCE(
typename BlockHelperDetails::ImplType<MatrixType>::execution_space)
1368 IFPACK2_BLOCKHELPER_TIMER(
"compute part indices", indices);
1369 for (local_ordinal_type ip = 0; ip < nparts; ++ip) {
1370 const auto *part = &partitions[p[ip]];
1371 const local_ordinal_type ipnrows = part->size();
1372 TEUCHOS_ASSERT(ip == 0 || (ipnrows <= static_cast<local_ordinal_type>(partitions[p[ip - 1]].size())));
1374 BlockHelperDetails::get_msg_prefix(comm)
1375 <<
"partition " << p[ip]
1376 <<
" is empty, which is not allowed.");
1378 part2rowidx0(ip + 1) = part2rowidx0(ip) + ipnrows;
1381 if (ip % vector_length == 0) pack_nrows = ipnrows;
1382 part2packrowidx0(ip + 1) = part2packrowidx0(ip) + ((ip + 1) % vector_length == 0 || ip + 1 == nparts ? pack_nrows : 0);
1383 const local_ordinal_type offset = partptr(ip);
1384 for (local_ordinal_type i = 0; i < ipnrows; ++i) {
1385 const auto lcl_row = (*part)[i];
1387 BlockHelperDetails::get_msg_prefix(comm)
1388 <<
"partitions[" << p[ip] <<
"]["
1389 << i <<
"] = " << lcl_row
1390 <<
" but input matrix implies limits of [0, " << A_n_lclrows - 1
1392 lclrow(offset + i) = lcl_row;
1393 rowidx2part(offset + i) = ip;
1394 if (interf.row_contiguous && offset + i > 0 && lclrow((offset + i) - 1) + 1 != lcl_row)
1395 interf.row_contiguous =
false;
1397 partptr(ip + 1) = offset + ipnrows;
1399 #ifdef IFPACK2_BLOCKTRIDICONTAINER_USE_PRINTF
1400 printf(
"Part index = ip = %d, first LID associated to the part = partptr(ip) = offset = %d, part->size() = ipnrows = %d;\n", ip, offset, ipnrows);
1401 printf(
"partptr(%d+1) = %d\n", ip, partptr(ip + 1));
1405 part2rowidx0_sub(0) = 0;
1406 partptr_sub(0, 0) = 0;
1409 for (local_ordinal_type ip = 0; ip < nparts; ++ip) {
1410 const auto *part = &partitions[p[ip]];
1411 const local_ordinal_type ipnrows = part->size();
1412 const local_ordinal_type full_line_length = partptr(ip + 1) - partptr(ip);
1415 "In the part " << ip);
1417 constexpr local_ordinal_type connection_length = 2;
1419 if (full_line_length < n_subparts_per_part + (n_subparts_per_part - 1) * connection_length)
1421 "The part " << ip <<
" is too short to use " << n_subparts_per_part <<
" sub parts.");
1423 const local_ordinal_type sub_line_length = (full_line_length - (n_subparts_per_part - 1) * connection_length) / n_subparts_per_part;
1424 const local_ordinal_type last_sub_line_length = full_line_length - (n_subparts_per_part - 1) * (connection_length + sub_line_length);
1426 if (ip % vector_length == 0) pack_nrows_sub = ipnrows;
1428 for (local_ordinal_type local_sub_ip = 0; local_sub_ip < n_subparts_per_part; ++local_sub_ip) {
1429 const local_ordinal_type sub_ip = nparts * (2 * local_sub_ip) + ip;
1430 const local_ordinal_type schur_ip = nparts * (2 * local_sub_ip + 1) + ip;
1431 if (local_sub_ip != n_subparts_per_part - 1) {
1432 if (local_sub_ip != 0) {
1433 partptr_sub(sub_ip, 0) = partptr_sub(nparts * (2 * local_sub_ip - 1) + ip, 1);
1434 }
else if (ip != 0) {
1435 partptr_sub(sub_ip, 0) = partptr_sub(nparts * 2 * (n_subparts_per_part - 1) + ip - 1, 1);
1437 partptr_sub(sub_ip, 1) = sub_line_length + partptr_sub(sub_ip, 0);
1438 partptr_sub(schur_ip, 0) = partptr_sub(sub_ip, 1);
1439 partptr_sub(schur_ip, 1) = connection_length + partptr_sub(schur_ip, 0);
1441 part2rowidx0_sub(sub_ip + 1) = part2rowidx0_sub(sub_ip) + sub_line_length;
1442 part2rowidx0_sub(sub_ip + 2) = part2rowidx0_sub(sub_ip + 1) + connection_length;
1444 #ifdef IFPACK2_BLOCKTRIDICONTAINER_USE_PRINTF
1445 printf(
"Sub Part index = %d, first LID associated to the sub part = %d, sub part size = %d;\n", sub_ip, partptr_sub(sub_ip, 0), sub_line_length);
1446 printf(
"Sub Part index Schur = %d, first LID associated to the sub part = %d, sub part size = %d;\n", sub_ip + 1, partptr_sub(ip, 2 * local_sub_ip + 1), connection_length);
1449 if (local_sub_ip != 0) {
1450 partptr_sub(sub_ip, 0) = partptr_sub(nparts * (2 * local_sub_ip - 1) + ip, 1);
1451 }
else if (ip != 0) {
1452 partptr_sub(sub_ip, 0) = partptr_sub(nparts * 2 * (n_subparts_per_part - 1) + ip - 1, 1);
1454 partptr_sub(sub_ip, 1) = last_sub_line_length + partptr_sub(sub_ip, 0);
1456 part2rowidx0_sub(sub_ip + 1) = part2rowidx0_sub(sub_ip) + last_sub_line_length;
1458 #ifdef IFPACK2_BLOCKTRIDICONTAINER_USE_PRINTF
1459 printf(
"Sub Part index = %d, first LID associated to the sub part = %d, sub part size = %d;\n", sub_ip, partptr_sub(sub_ip, 0), last_sub_line_length);
1466 local_ordinal_type npacks = ceil(
float(nparts) / vector_length);
1468 local_ordinal_type ip_max = nparts > vector_length ? vector_length : nparts;
1469 for (local_ordinal_type ip = 0; ip < ip_max; ++ip) {
1470 part2packrowidx0_sub(ip, 0) = 0;
1472 for (local_ordinal_type ipack = 0; ipack < npacks; ++ipack) {
1474 local_ordinal_type ip_min = ipack * vector_length;
1475 ip_max = nparts > (ipack + 1) * vector_length ? (ipack + 1) * vector_length : nparts;
1476 for (local_ordinal_type ip = ip_min; ip < ip_max; ++ip) {
1477 part2packrowidx0_sub(ip, 0) = part2packrowidx0_sub(ip - vector_length, part2packrowidx0_sub.extent(1) - 1);
1481 for (size_type local_sub_ip = 0; local_sub_ip < part2packrowidx0_sub.extent(1) - 1; ++local_sub_ip) {
1482 local_ordinal_type ip_min = ipack * vector_length;
1483 ip_max = nparts > (ipack + 1) * vector_length ? (ipack + 1) * vector_length : nparts;
1485 const local_ordinal_type full_line_length = partptr(ip_min + 1) - partptr(ip_min);
1487 constexpr local_ordinal_type connection_length = 2;
1489 const local_ordinal_type sub_line_length = (full_line_length - (n_subparts_per_part - 1) * connection_length) / n_subparts_per_part;
1490 const local_ordinal_type last_sub_line_length = full_line_length - (n_subparts_per_part - 1) * (connection_length + sub_line_length);
1492 if (local_sub_ip % 2 == 0) pack_nrows_sub = sub_line_length;
1493 if (local_sub_ip % 2 == 1) pack_nrows_sub = connection_length;
1494 if (local_sub_ip == part2packrowidx0_sub.extent(1) - 2) pack_nrows_sub = last_sub_line_length;
1496 part2packrowidx0_sub(ip_min, local_sub_ip + 1) = part2packrowidx0_sub(ip_min, local_sub_ip) + pack_nrows_sub;
1498 for (local_ordinal_type ip = ip_min + 1; ip < ip_max; ++ip) {
1499 part2packrowidx0_sub(ip, local_sub_ip + 1) = part2packrowidx0_sub(ip_min, local_sub_ip + 1);
1504 Kokkos::deep_copy(interf.part2packrowidx0_sub, part2packrowidx0_sub);
1506 IFPACK2_BLOCKHELPER_TIMER_FENCE(
typename BlockHelperDetails::ImplType<MatrixType>::execution_space)
1508 #if defined(BLOCKTRIDICONTAINER_DEBUG)
1511 if (lclrow(0) != 0) interf.row_contiguous =
false;
1513 Kokkos::deep_copy(interf.partptr, partptr);
1514 Kokkos::deep_copy(interf.lclrow, lclrow);
1516 Kokkos::deep_copy(interf.partptr_sub, partptr_sub);
1519 interf.part2rowidx0 = interf.partptr;
1520 Kokkos::deep_copy(interf.part2packrowidx0, part2packrowidx0);
1522 interf.part2packrowidx0_back = part2packrowidx0_sub(part2packrowidx0_sub.extent(0) - 1, part2packrowidx0_sub.extent(1) - 1);
1523 Kokkos::deep_copy(interf.rowidx2part, rowidx2part);
1526 IFPACK2_BLOCKHELPER_TIMER(
"Fill packptr", packptr0);
1527 local_ordinal_type npacks = ceil(
float(nparts) / vector_length) * (part2packrowidx0_sub.extent(1) - 1);
1529 for (local_ordinal_type ip = 1; ip <= nparts; ++ip)
1530 if (part2packrowidx0(ip) != part2packrowidx0(ip - 1))
1534 const auto packptr = Kokkos::create_mirror_view(interf.packptr);
1536 for (local_ordinal_type ip = 1, k = 1; ip <= nparts; ++ip)
1537 if (part2packrowidx0(ip) != part2packrowidx0(ip - 1))
1540 Kokkos::deep_copy(interf.packptr, packptr);
1542 local_ordinal_type npacks_per_subpart = ceil(
float(nparts) / vector_length);
1543 npacks = ceil(
float(nparts) / vector_length) * (part2packrowidx0_sub.extent(1) - 1);
1545 interf.packindices_sub = local_ordinal_type_1d_view(
do_not_initialize_tag(
"packindices_sub"), npacks_per_subpart * n_subparts_per_part);
1546 interf.packindices_schur = local_ordinal_type_2d_view(
do_not_initialize_tag(
"packindices_schur"), npacks_per_subpart, n_subparts_per_part - 1);
1548 const auto packindices_sub = Kokkos::create_mirror_view(interf.packindices_sub);
1549 const auto packindices_schur = Kokkos::create_mirror_view(interf.packindices_schur);
1552 for (local_ordinal_type local_sub_ip = 0; local_sub_ip < n_subparts_per_part - 1; ++local_sub_ip) {
1553 for (local_ordinal_type local_pack_ip = 0; local_pack_ip < npacks_per_subpart; ++local_pack_ip) {
1554 packindices_sub(local_sub_ip * npacks_per_subpart + local_pack_ip) = 2 * local_sub_ip * npacks_per_subpart + local_pack_ip;
1555 packindices_schur(local_pack_ip, local_sub_ip) = 2 * local_sub_ip * npacks_per_subpart + local_pack_ip + npacks_per_subpart;
1559 for (local_ordinal_type local_pack_ip = 0; local_pack_ip < npacks_per_subpart; ++local_pack_ip) {
1560 packindices_sub((n_subparts_per_part - 1) * npacks_per_subpart + local_pack_ip) = 2 * (n_subparts_per_part - 1) * npacks_per_subpart + local_pack_ip;
1563 #ifdef IFPACK2_BLOCKTRIDICONTAINER_WRITE_MM
1564 std::cout <<
"packindices_sub = " << std::endl;
1565 for (size_type i = 0; i < packindices_sub.extent(0); ++i) {
1566 std::cout << packindices_sub(i) <<
" ";
1568 std::cout << std::endl;
1569 std::cout <<
"packindices_sub end" << std::endl;
1571 std::cout <<
"packindices_schur = " << std::endl;
1572 for (size_type i = 0; i < packindices_schur.extent(0); ++i) {
1573 for (size_type j = 0; j < packindices_schur.extent(1); ++j) {
1574 std::cout << packindices_schur(i, j) <<
" ";
1576 std::cout << std::endl;
1579 std::cout <<
"packindices_schur end" << std::endl;
1582 Kokkos::deep_copy(interf.packindices_sub, packindices_sub);
1583 Kokkos::deep_copy(interf.packindices_schur, packindices_schur);
1586 const auto packptr_sub = Kokkos::create_mirror_view(interf.packptr_sub);
1588 for (local_ordinal_type k = 0; k < npacks + 1; ++k)
1589 packptr_sub(k) = packptr(k % npacks_per_subpart) + (k / npacks_per_subpart) * packptr(npacks_per_subpart);
1591 Kokkos::deep_copy(interf.packptr_sub, packptr_sub);
1592 IFPACK2_BLOCKHELPER_TIMER_FENCE(
typename BlockHelperDetails::ImplType<MatrixType>::execution_space)
1594 IFPACK2_BLOCKHELPER_TIMER_FENCE(
typename BlockHelperDetails::ImplType<MatrixType>::execution_space)
1602 template <
typename MatrixType>
1605 using local_ordinal_type_1d_view =
typename impl_type::local_ordinal_type_1d_view;
1607 using size_type_2d_view =
typename impl_type::size_type_2d_view;
1608 using vector_type_3d_view =
typename impl_type::vector_type_3d_view;
1609 using vector_type_4d_view =
typename impl_type::vector_type_4d_view;
1610 using btdm_scalar_type_3d_view =
typename impl_type::btdm_scalar_type_3d_view;
1616 size_type_2d_view flat_td_ptr, pack_td_ptr, pack_td_ptr_schur;
1619 local_ordinal_type_1d_view A_colindsub;
1622 vector_type_3d_view values;
1625 vector_type_3d_view values_schur;
1627 vector_type_4d_view e_values;
1632 size_type_1d_view diag_offsets;
1636 btdm_scalar_type_3d_view d_inv;
1638 bool is_diagonal_only;
1644 template <
typename idx_type>
1645 static KOKKOS_FORCEINLINE_FUNCTION
1647 IndexToRow(
const idx_type &ind) {
return (ind + 1) / 3; }
1650 template <
typename idx_type>
1651 static KOKKOS_FORCEINLINE_FUNCTION
1653 RowToIndex(
const idx_type &row) {
return row > 0 ? 3 * row - 1 : 0; }
1655 template <
typename idx_type>
1656 static KOKKOS_FORCEINLINE_FUNCTION
1658 NumBlocks(
const idx_type &nrows) {
return nrows > 0 ? 3 * nrows - 2 : 0; }
1660 template <
typename idx_type>
1661 static KOKKOS_FORCEINLINE_FUNCTION
1663 NumBlocksSchur(
const idx_type &nrows) {
return nrows > 0 ? 3 * nrows + 2 : 0; }
1669 template <
typename MatrixType>
1672 IFPACK2_BLOCKHELPER_TIMER(
"createBlockTridiags", createBlockTridiags0);
1674 using execution_space =
typename impl_type::execution_space;
1675 using local_ordinal_type =
typename impl_type::local_ordinal_type;
1676 using size_type =
typename impl_type::size_type;
1677 using size_type_2d_view =
typename impl_type::size_type_2d_view;
1679 constexpr
int vector_length = impl_type::vector_length;
1683 const local_ordinal_type ntridiags = interf.partptr_sub.extent(0);
1686 btdm.flat_td_ptr = size_type_2d_view(
do_not_initialize_tag(
"btdm.flat_td_ptr"), interf.nparts, 2 * interf.n_subparts_per_part);
1687 const Kokkos::RangePolicy<execution_space> policy(0, 2 * interf.nparts * interf.n_subparts_per_part);
1688 Kokkos::parallel_scan(
1689 "createBlockTridiags::RangePolicy::flat_td_ptr",
1690 policy, KOKKOS_LAMBDA(
const local_ordinal_type &i, size_type &update,
const bool &
final) {
1691 const local_ordinal_type partidx = i / (2 * interf.n_subparts_per_part);
1692 const local_ordinal_type local_subpartidx = i % (2 * interf.n_subparts_per_part);
1695 btdm.flat_td_ptr(partidx, local_subpartidx) = update;
1697 if (local_subpartidx != (2 * interf.n_subparts_per_part - 1)) {
1698 const local_ordinal_type nrows = interf.partptr_sub(interf.nparts * local_subpartidx + partidx, 1) - interf.partptr_sub(interf.nparts * local_subpartidx + partidx, 0);
1699 if (local_subpartidx % 2 == 0)
1700 update += btdm.NumBlocks(nrows);
1702 update += btdm.NumBlocksSchur(nrows);
1706 const auto nblocks = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), Kokkos::subview(btdm.flat_td_ptr, interf.nparts - 1, 2 * interf.n_subparts_per_part - 1));
1707 btdm.is_diagonal_only = (
static_cast<local_ordinal_type
>(nblocks()) == ntridiags);
1711 if (vector_length == 1) {
1712 btdm.pack_td_ptr = btdm.flat_td_ptr;
1716 local_ordinal_type npacks_per_subpart = 0;
1717 const auto part2packrowidx0 = Kokkos::create_mirror_view(interf.part2packrowidx0);
1718 Kokkos::deep_copy(part2packrowidx0, interf.part2packrowidx0);
1719 for (local_ordinal_type ip = 1; ip <= interf.nparts; ++ip)
1720 if (part2packrowidx0(ip) != part2packrowidx0(ip - 1))
1721 ++npacks_per_subpart;
1723 btdm.pack_td_ptr = size_type_2d_view(
do_not_initialize_tag(
"btdm.pack_td_ptr"), interf.nparts, 2 * interf.n_subparts_per_part);
1724 const Kokkos::RangePolicy<execution_space> policy(0, npacks_per_subpart);
1726 Kokkos::parallel_for(
1727 "createBlockTridiags::RangePolicy::pack_td_ptr",
1728 policy, KOKKOS_LAMBDA(
const local_ordinal_type &i) {
1729 for (local_ordinal_type j = 0; j < 2 * interf.n_subparts_per_part; ++j) {
1730 const local_ordinal_type pack_id = (j == 2 * interf.n_subparts_per_part - 1) ? i + (j - 1) * npacks_per_subpart : i + j * npacks_per_subpart;
1731 const local_ordinal_type nparts_in_pack = interf.packptr_sub(pack_id + 1) - interf.packptr_sub(pack_id);
1733 const local_ordinal_type parti = interf.packptr_sub(pack_id);
1734 const local_ordinal_type partidx = parti % interf.nparts;
1736 for (local_ordinal_type pti = 0; pti < nparts_in_pack; ++pti) {
1737 btdm.pack_td_ptr(partidx + pti, j) = btdm.flat_td_ptr(i, j);
1743 btdm.pack_td_ptr_schur = size_type_2d_view(
do_not_initialize_tag(
"btdm.pack_td_ptr_schur"), interf.nparts, interf.n_subparts_per_part);
1745 const auto host_pack_td_ptr_schur = Kokkos::create_mirror_view(btdm.pack_td_ptr_schur);
1746 constexpr local_ordinal_type connection_length = 2;
1748 host_pack_td_ptr_schur(0, 0) = 0;
1749 for (local_ordinal_type i = 0; i < interf.nparts; ++i) {
1750 if (i % vector_length == 0) {
1752 host_pack_td_ptr_schur(i, 0) = host_pack_td_ptr_schur(i - 1, host_pack_td_ptr_schur.extent(1) - 1);
1753 for (local_ordinal_type j = 0; j < interf.n_subparts_per_part - 1; ++j) {
1754 host_pack_td_ptr_schur(i, j + 1) = host_pack_td_ptr_schur(i, j) + btdm.NumBlocks(connection_length) + (j != 0 ? 1 : 0) + (j != interf.n_subparts_per_part - 2 ? 1 : 0);
1757 for (local_ordinal_type j = 0; j < interf.n_subparts_per_part; ++j) {
1758 host_pack_td_ptr_schur(i, j) = host_pack_td_ptr_schur(i - 1, j);
1763 Kokkos::deep_copy(btdm.pack_td_ptr_schur, host_pack_td_ptr_schur);
1765 #ifdef IFPACK2_BLOCKTRIDICONTAINER_WRITE_MM
1766 const auto host_flat_td_ptr = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), btdm.flat_td_ptr);
1767 std::cout <<
"flat_td_ptr = " << std::endl;
1768 for (size_type i = 0; i < host_flat_td_ptr.extent(0); ++i) {
1769 for (size_type j = 0; j < host_flat_td_ptr.extent(1); ++j) {
1770 std::cout << host_flat_td_ptr(i, j) <<
" ";
1772 std::cout << std::endl;
1774 std::cout <<
"flat_td_ptr end" << std::endl;
1776 const auto host_pack_td_ptr = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), btdm.pack_td_ptr);
1778 std::cout <<
"pack_td_ptr = " << std::endl;
1779 for (size_type i = 0; i < host_pack_td_ptr.extent(0); ++i) {
1780 for (size_type j = 0; j < host_pack_td_ptr.extent(1); ++j) {
1781 std::cout << host_pack_td_ptr(i, j) <<
" ";
1783 std::cout << std::endl;
1785 std::cout <<
"pack_td_ptr end" << std::endl;
1787 std::cout <<
"pack_td_ptr_schur = " << std::endl;
1788 for (size_type i = 0; i < host_pack_td_ptr_schur.extent(0); ++i) {
1789 for (size_type j = 0; j < host_pack_td_ptr_schur.extent(1); ++j) {
1790 std::cout << host_pack_td_ptr_schur(i, j) <<
" ";
1792 std::cout << std::endl;
1794 std::cout <<
"pack_td_ptr_schur end" << std::endl;
1798 IFPACK2_BLOCKHELPER_TIMER_FENCE(
typename BlockHelperDetails::ImplType<MatrixType>::execution_space)
1813 template <
typename MatrixType>
1814 void setTridiagsToIdentity(
const BlockTridiags<MatrixType> &btdm,
1815 const typename BlockHelperDetails::ImplType<MatrixType>::local_ordinal_type_1d_view &packptr) {
1817 using execution_space =
typename impl_type::execution_space;
1818 using local_ordinal_type =
typename impl_type::local_ordinal_type;
1819 using size_type_2d_view =
typename impl_type::size_type_2d_view;
1821 const ConstUnmanaged<size_type_2d_view> pack_td_ptr(btdm.pack_td_ptr);
1822 const local_ordinal_type blocksize = btdm.values.extent(1);
1825 const int vector_length = impl_type::vector_length;
1826 const int internal_vector_length = impl_type::internal_vector_length;
1828 using btdm_scalar_type =
typename impl_type::btdm_scalar_type;
1829 using internal_vector_type =
typename impl_type::internal_vector_type;
1830 using internal_vector_type_4d_view =
1831 typename impl_type::internal_vector_type_4d_view;
1833 using team_policy_type = Kokkos::TeamPolicy<execution_space>;
1834 const internal_vector_type_4d_view values(reinterpret_cast<internal_vector_type *>(btdm.values.data()),
1835 btdm.values.extent(0),
1836 btdm.values.extent(1),
1837 btdm.values.extent(2),
1838 vector_length / internal_vector_length);
1839 const local_ordinal_type vector_loop_size = values.extent(3);
1840 #if defined(KOKKOS_ENABLE_CUDA) && defined(__CUDA_ARCH__)
1841 local_ordinal_type total_team_size(0);
1843 total_team_size = 32;
1844 else if (blocksize <= 9)
1845 total_team_size = 64;
1846 else if (blocksize <= 12)
1847 total_team_size = 96;
1848 else if (blocksize <= 16)
1849 total_team_size = 128;
1850 else if (blocksize <= 20)
1851 total_team_size = 160;
1853 total_team_size = 160;
1854 const local_ordinal_type team_size = total_team_size / vector_loop_size;
1855 const team_policy_type policy(packptr.extent(0) - 1, team_size, vector_loop_size);
1856 #elif defined(KOKKOS_ENABLE_HIP)
1861 local_ordinal_type total_team_size(0);
1863 total_team_size = 32;
1864 else if (blocksize <= 9)
1865 total_team_size = 64;
1866 else if (blocksize <= 12)
1867 total_team_size = 96;
1868 else if (blocksize <= 16)
1869 total_team_size = 128;
1870 else if (blocksize <= 20)
1871 total_team_size = 160;
1873 total_team_size = 160;
1874 const local_ordinal_type team_size = total_team_size / vector_loop_size;
1875 const team_policy_type policy(packptr.extent(0) - 1, team_size, vector_loop_size);
1876 #elif defined(KOKKOS_ENABLE_SYCL)
1878 local_ordinal_type total_team_size(0);
1880 total_team_size = 32;
1881 else if (blocksize <= 9)
1882 total_team_size = 64;
1883 else if (blocksize <= 12)
1884 total_team_size = 96;
1885 else if (blocksize <= 16)
1886 total_team_size = 128;
1887 else if (blocksize <= 20)
1888 total_team_size = 160;
1890 total_team_size = 160;
1891 const local_ordinal_type team_size = total_team_size / vector_loop_size;
1892 const team_policy_type policy(packptr.extent(0) - 1, team_size, vector_loop_size);
1895 const team_policy_type policy(packptr.extent(0) - 1, 1, 1);
1897 Kokkos::parallel_for(
1898 "setTridiagsToIdentity::TeamPolicy",
1899 policy, KOKKOS_LAMBDA(
const typename team_policy_type::member_type &member) {
1900 const local_ordinal_type k = member.league_rank();
1901 const local_ordinal_type ibeg = pack_td_ptr(packptr(k), 0);
1902 const local_ordinal_type iend = pack_td_ptr(packptr(k), pack_td_ptr.extent(1) - 1);
1904 const local_ordinal_type diff = iend - ibeg;
1905 const local_ordinal_type icount = diff / 3 + (diff % 3 > 0);
1906 const btdm_scalar_type one(1);
1907 Kokkos::parallel_for(Kokkos::ThreadVectorRange(member, vector_loop_size), [&](
const int &v) {
1908 Kokkos::parallel_for(Kokkos::TeamThreadRange(member, icount), [&](
const local_ordinal_type &ii) {
1909 const local_ordinal_type i = ibeg + ii * 3;
1910 for (local_ordinal_type j = 0; j < blocksize; ++j) {
1911 values(i, j, j, v) = one;
1922 template <
typename MatrixType>
1924 const Teuchos::RCP<
const typename BlockHelperDetails::ImplType<MatrixType>::tpetra_crs_graph_type> &g,
1925 const BlockHelperDetails::PartInterface<MatrixType> &interf,
1928 const bool overlap_communication_and_computation,
1929 const Teuchos::RCP<AsyncableImport<MatrixType>> &async_importer,
1931 bool use_fused_jacobi) {
1932 IFPACK2_BLOCKHELPER_TIMER(
"BlockTriDi::SymbolicPhase", SymbolicPhase);
1936 using execution_space =
typename impl_type::execution_space;
1938 using local_ordinal_type =
typename impl_type::local_ordinal_type;
1939 using global_ordinal_type =
typename impl_type::global_ordinal_type;
1940 using size_type =
typename impl_type::size_type;
1941 using local_ordinal_type_1d_view =
typename impl_type::local_ordinal_type_1d_view;
1942 using size_type_1d_view =
typename impl_type::size_type_1d_view;
1943 using vector_type_3d_view =
typename impl_type::vector_type_3d_view;
1944 using vector_type_4d_view =
typename impl_type::vector_type_4d_view;
1945 using crs_matrix_type =
typename impl_type::tpetra_crs_matrix_type;
1946 using block_crs_matrix_type =
typename impl_type::tpetra_block_crs_matrix_type;
1947 using btdm_scalar_type_3d_view =
typename impl_type::btdm_scalar_type_3d_view;
1948 using lo_traits = Tpetra::Details::OrdinalTraits<local_ordinal_type>;
1950 constexpr
int vector_length = impl_type::vector_length;
1952 const auto comm = A->getRowMap()->getComm();
1954 auto A_crs = Teuchos::rcp_dynamic_cast<
const crs_matrix_type>(A);
1955 auto A_bcrs = Teuchos::rcp_dynamic_cast<
const block_crs_matrix_type>(A);
1957 bool hasBlockCrsMatrix = !A_bcrs.is_null();
1959 const local_ordinal_type blocksize = hasBlockCrsMatrix ? A->getBlockSize() : A->getLocalNumRows() / g->getLocalNumRows();
1961 const auto partptr = interf.partptr;
1962 const auto lclrow = interf.lclrow;
1963 const auto rowidx2part = interf.rowidx2part;
1964 const auto part2rowidx0 = interf.part2rowidx0;
1965 const auto packptr = interf.packptr;
1968 const local_ordinal_type nrows = Kokkos::create_mirror_view_and_copy(
1969 Kokkos::HostSpace(), Kokkos::subview(partptr, partptr.extent(0) - 1))();
1971 Kokkos::View<local_ordinal_type *, execution_space> col2row(
"col2row", A->getLocalNumCols());
1977 TEUCHOS_ASSERT(!(g->getRowMap().is_null() || g->getColMap().is_null() || g->getDomainMap().is_null()));
1978 #if defined(BLOCKTRIDICONTAINER_DEBUG)
1981 auto rowmapHost = g->getRowMap();
1982 auto colmapHost = g->getColMap();
1983 auto dommapHost = g->getDomainMap();
1984 for (local_ordinal_type lr = 0; lr < nrows; lr++) {
1985 const global_ordinal_type gid = rowmapHost->getGlobalElement(lr);
1987 if (dommapHost->isNodeGlobalElement(gid)) {
1988 const local_ordinal_type lc = colmapHost->getLocalElement(gid);
1990 BlockHelperDetails::get_msg_prefix(comm) <<
"GID " << gid
1991 <<
" gives an invalid local column.");
1996 auto rowmap = g->getRowMap()->getLocalMap();
1997 auto colmap = g->getColMap()->getLocalMap();
1998 auto dommap = g->getDomainMap()->getLocalMap();
2000 const Kokkos::RangePolicy<execution_space> policy(0, nrows);
2001 Kokkos::parallel_for(
2002 "performSymbolicPhase::RangePolicy::col2row",
2003 policy, KOKKOS_LAMBDA(
const local_ordinal_type &lr) {
2004 const global_ordinal_type gid = rowmap.getGlobalElement(lr);
2005 if (dommap.getLocalElement(gid) != lo_traits::invalid()) {
2006 const local_ordinal_type lc = colmap.getLocalElement(gid);
2014 const auto local_graph = g->getLocalGraphDevice();
2015 const auto local_graph_rowptr = local_graph.row_map;
2016 TEUCHOS_ASSERT(local_graph_rowptr.size() ==
static_cast<size_t>(nrows + 1));
2017 const auto local_graph_colidx = local_graph.entries;
2021 Kokkos::View<local_ordinal_type *, execution_space> lclrow2idx(
"lclrow2idx", nrows);
2023 const Kokkos::RangePolicy<execution_space> policy(0, nrows);
2024 Kokkos::parallel_for(
2025 "performSymbolicPhase::RangePolicy::lclrow2idx",
2026 policy, KOKKOS_LAMBDA(
const local_ordinal_type &i) {
2027 lclrow2idx(lclrow(i)) = i;
2032 size_type D_nnz, R_nnz_owned, R_nnz_remote;
2034 const Kokkos::RangePolicy<execution_space> policy(0, nrows);
2035 Kokkos::parallel_reduce
2038 policy, KOKKOS_LAMBDA(
const local_ordinal_type &lr, size_type &update_D_nnz, size_type &update_R_nnz_owned, size_type &update_R_nnz_remote) {
2040 const local_ordinal_type ri0 = lclrow2idx(lr);
2041 const local_ordinal_type pi0 = rowidx2part(ri0);
2042 for (size_type j = local_graph_rowptr(lr); j < local_graph_rowptr(lr + 1); ++j) {
2043 const local_ordinal_type lc = local_graph_colidx(j);
2044 const local_ordinal_type lc2r = col2row(lc);
2045 bool incr_R =
false;
2047 if (lc2r == (local_ordinal_type)-1) {
2051 const local_ordinal_type ri = lclrow2idx(lc2r);
2052 const local_ordinal_type pi = rowidx2part(ri);
2060 if (ri0 + 1 >= ri && ri0 <= ri + 1)
2067 ++update_R_nnz_owned;
2069 ++update_R_nnz_remote;
2073 D_nnz, R_nnz_owned, R_nnz_remote);
2076 if (!overlap_communication_and_computation) {
2077 R_nnz_owned += R_nnz_remote;
2083 const auto flat_td_ptr = btdm.flat_td_ptr;
2085 btdm.A_colindsub = local_ordinal_type_1d_view(
"btdm.A_colindsub", D_nnz);
2086 const auto D_A_colindsub = btdm.A_colindsub;
2088 #if defined(BLOCKTRIDICONTAINER_DEBUG)
2092 const local_ordinal_type nparts = partptr.extent(0) - 1;
2095 const Kokkos::RangePolicy<execution_space> policy(0, nparts);
2096 Kokkos::parallel_for(
2097 "performSymbolicPhase::RangePolicy<execution_space>::D_graph",
2098 policy, KOKKOS_LAMBDA(
const local_ordinal_type &pi0) {
2099 const local_ordinal_type part_ri0 = part2rowidx0(pi0);
2100 local_ordinal_type offset = 0;
2101 for (local_ordinal_type ri0 = partptr(pi0); ri0 < partptr(pi0 + 1); ++ri0) {
2102 const local_ordinal_type td_row_os = btdm.RowToIndex(ri0 - part_ri0) + offset;
2104 const local_ordinal_type lr0 = lclrow(ri0);
2105 const size_type j0 = local_graph_rowptr(lr0);
2106 for (size_type j = j0; j < local_graph_rowptr(lr0 + 1); ++j) {
2107 const local_ordinal_type lc = local_graph_colidx(j);
2108 const local_ordinal_type lc2r = col2row[lc];
2109 if (lc2r == (local_ordinal_type)-1)
continue;
2110 const local_ordinal_type ri = lclrow2idx[lc2r];
2111 const local_ordinal_type pi = rowidx2part(ri);
2112 if (pi != pi0)
continue;
2113 if (ri + 1 < ri0 || ri > ri0 + 1)
continue;
2114 const local_ordinal_type row_entry = j - j0;
2115 D_A_colindsub(flat_td_ptr(pi0, 0) + ((td_row_os + ri) - ri0)) = row_entry;
2120 #if defined(BLOCKTRIDICONTAINER_DEBUG)
2122 auto D_A_colindsub_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), D_A_colindsub);
2123 for (
size_t i = 0; i < D_A_colindsub_host.extent(0); ++i)
2130 const auto pack_td_ptr_last = Kokkos::subview(btdm.pack_td_ptr, btdm.pack_td_ptr.extent(0) - 1, btdm.pack_td_ptr.extent(1) - 1);
2131 const auto num_packed_blocks = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), pack_td_ptr_last);
2132 btdm.values = vector_type_3d_view(
"btdm.values", num_packed_blocks(), blocksize, blocksize);
2134 if (interf.n_subparts_per_part > 1) {
2135 const auto pack_td_ptr_schur_last = Kokkos::subview(btdm.pack_td_ptr_schur, btdm.pack_td_ptr_schur.extent(0) - 1, btdm.pack_td_ptr_schur.extent(1) - 1);
2136 const auto num_packed_blocks_schur = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), pack_td_ptr_schur_last);
2137 btdm.values_schur = vector_type_3d_view(
"btdm.values_schur", num_packed_blocks_schur(), blocksize, blocksize);
2140 if (vector_length > 1) setTridiagsToIdentity(btdm, interf.packptr);
2146 amd.rowptr = size_type_1d_view(
"amd.rowptr", nrows + 1);
2147 amd.A_colindsub = local_ordinal_type_1d_view(
do_not_initialize_tag(
"amd.A_colindsub"), R_nnz_owned);
2149 const auto R_rowptr = amd.rowptr;
2150 const auto R_A_colindsub = amd.A_colindsub;
2152 amd.rowptr_remote = size_type_1d_view(
"amd.rowptr_remote", overlap_communication_and_computation ? nrows + 1 : 0);
2153 amd.A_colindsub_remote = local_ordinal_type_1d_view(
do_not_initialize_tag(
"amd.A_colindsub_remote"), R_nnz_remote);
2155 const auto R_rowptr_remote = amd.rowptr_remote;
2156 const auto R_A_colindsub_remote = amd.A_colindsub_remote;
2159 const Kokkos::RangePolicy<execution_space> policy(0, nrows);
2160 Kokkos::parallel_for(
2161 "performSymbolicPhase::RangePolicy<execution_space>::R_graph_count",
2162 policy, KOKKOS_LAMBDA(
const local_ordinal_type &lr) {
2163 const local_ordinal_type ri0 = lclrow2idx[lr];
2164 const local_ordinal_type pi0 = rowidx2part(ri0);
2165 const size_type j0 = local_graph_rowptr(lr);
2166 for (size_type j = j0; j < local_graph_rowptr(lr + 1); ++j) {
2167 const local_ordinal_type lc = local_graph_colidx(j);
2168 const local_ordinal_type lc2r = col2row[lc];
2169 if (lc2r != (local_ordinal_type)-1) {
2170 const local_ordinal_type ri = lclrow2idx[lc2r];
2171 const local_ordinal_type pi = rowidx2part(ri);
2172 if (pi == pi0 && ri + 1 >= ri0 && ri <= ri0 + 1) {
2177 if (!overlap_communication_and_computation || lc < nrows) {
2180 ++R_rowptr_remote(lr);
2189 Kokkos::RangePolicy<execution_space> policy(0, nrows + 1);
2190 Kokkos::parallel_scan(
2191 "performSymbolicPhase::RangePolicy<execution_space>::R_graph_fill",
2192 policy, KOKKOS_LAMBDA(
const local_ordinal_type &lr, update_type &update,
const bool &
final) {
2194 val.v[0] = R_rowptr(lr);
2195 if (overlap_communication_and_computation)
2196 val.v[1] = R_rowptr_remote(lr);
2199 R_rowptr(lr) = update.v[0];
2200 if (overlap_communication_and_computation)
2201 R_rowptr_remote(lr) = update.v[1];
2204 const local_ordinal_type ri0 = lclrow2idx[lr];
2205 const local_ordinal_type pi0 = rowidx2part(ri0);
2207 size_type cnt_rowptr = R_rowptr(lr);
2208 size_type cnt_rowptr_remote = overlap_communication_and_computation ? R_rowptr_remote(lr) : 0;
2210 const size_type j0 = local_graph_rowptr(lr);
2211 for (size_type j = j0; j < local_graph_rowptr(lr + 1); ++j) {
2212 const local_ordinal_type lc = local_graph_colidx(j);
2213 const local_ordinal_type lc2r = col2row[lc];
2214 if (lc2r != (local_ordinal_type)-1) {
2215 const local_ordinal_type ri = lclrow2idx[lc2r];
2216 const local_ordinal_type pi = rowidx2part(ri);
2217 if (pi == pi0 && ri + 1 >= ri0 && ri <= ri0 + 1)
2220 const local_ordinal_type row_entry = j - j0;
2221 if (!overlap_communication_and_computation || lc < nrows)
2222 R_A_colindsub(cnt_rowptr++) = row_entry;
2224 R_A_colindsub_remote(cnt_rowptr_remote++) = row_entry;
2234 auto r_rowptr_end = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), Kokkos::subview(R_rowptr, nrows));
2236 if (overlap_communication_and_computation) {
2237 auto r_rowptr_remote_end = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), Kokkos::subview(R_rowptr_remote, nrows));
2243 if (hasBlockCrsMatrix)
2244 amd.tpetra_values = (
const_cast<block_crs_matrix_type *
>(A_bcrs.get())->getValuesDeviceNonConst());
2246 amd.tpetra_values = (
const_cast<crs_matrix_type *
>(A_crs.get()))->getLocalValuesDevice(Tpetra::Access::ReadWrite);
2252 if (interf.n_subparts_per_part > 1)
2253 btdm.e_values = vector_type_4d_view(
"btdm.e_values", 2, interf.part2packrowidx0_back, blocksize, blocksize);
2263 if (BlockHelperDetails::is_device<execution_space>::value && !useSeqMethod && hasBlockCrsMatrix) {
2264 bool is_async_importer_active = !async_importer.is_null();
2265 local_ordinal_type_1d_view dm2cm = is_async_importer_active ? async_importer->dm2cm : local_ordinal_type_1d_view();
2266 bool ownedRemoteSeparate = overlap_communication_and_computation || !is_async_importer_active;
2267 BlockHelperDetails::precompute_A_x_offsets<MatrixType>(amd, interf, g, dm2cm, blocksize, ownedRemoteSeparate);
2271 if (use_fused_jacobi) {
2272 btdm.d_inv = btdm_scalar_type_3d_view(
do_not_initialize_tag(
"btdm.d_inv"), interf.nparts, blocksize, blocksize);
2273 auto rowptrs = A_bcrs->getCrsGraph().getLocalRowPtrsDevice();
2274 auto entries = A_bcrs->getCrsGraph().getLocalIndicesDevice();
2275 btdm.diag_offsets = BlockHelperDetails::findDiagOffsets<execution_space, size_type_1d_view>(rowptrs, entries, interf.nparts, blocksize);
2277 IFPACK2_BLOCKHELPER_TIMER_FENCE(
typename BlockHelperDetails::ImplType<MatrixType>::execution_space)
2283 template <
typename ArgActiveExecutionMemorySpace>
2288 typedef KB::Mode::Serial mode_type;
2289 #if defined(__KOKKOSBATCHED_INTEL_MKL_COMPACT_BATCHED__)
2290 typedef KB::Algo::Level3::CompactMKL algo_type;
2292 typedef KB::Algo::Level3::Blocked algo_type;
2294 static int recommended_team_size(
const int ,
2301 #if defined(KOKKOS_ENABLE_CUDA)
2302 static inline int ExtractAndFactorizeRecommendedCudaTeamSize(
const int blksize,
2303 const int vector_length,
2304 const int internal_vector_length) {
2305 const int vector_size = vector_length / internal_vector_length;
2306 int total_team_size(0);
2308 total_team_size = 32;
2309 else if (blksize <= 9)
2310 total_team_size = 32;
2311 else if (blksize <= 12)
2312 total_team_size = 96;
2313 else if (blksize <= 16)
2314 total_team_size = 128;
2315 else if (blksize <= 20)
2316 total_team_size = 160;
2318 total_team_size = 160;
2319 return 2 * total_team_size / vector_size;
2322 struct ExtractAndFactorizeTridiagsDefaultModeAndAlgo<Kokkos::CudaSpace> {
2323 typedef KB::Mode::Team mode_type;
2324 typedef KB::Algo::Level3::Unblocked algo_type;
2325 static int recommended_team_size(
const int blksize,
2326 const int vector_length,
2327 const int internal_vector_length) {
2328 return ExtractAndFactorizeRecommendedCudaTeamSize(blksize, vector_length, internal_vector_length);
2332 struct ExtractAndFactorizeTridiagsDefaultModeAndAlgo<Kokkos::CudaUVMSpace> {
2333 typedef KB::Mode::Team mode_type;
2334 typedef KB::Algo::Level3::Unblocked algo_type;
2335 static int recommended_team_size(
const int blksize,
2336 const int vector_length,
2337 const int internal_vector_length) {
2338 return ExtractAndFactorizeRecommendedCudaTeamSize(blksize, vector_length, internal_vector_length);
2343 #if defined(KOKKOS_ENABLE_HIP)
2344 static inline int ExtractAndFactorizeRecommendedHIPTeamSize(
const int blksize,
2345 const int vector_length,
2346 const int internal_vector_length) {
2347 const int vector_size = vector_length / internal_vector_length;
2348 int total_team_size(0);
2350 total_team_size = 32;
2351 else if (blksize <= 9)
2352 total_team_size = 32;
2353 else if (blksize <= 12)
2354 total_team_size = 96;
2355 else if (blksize <= 16)
2356 total_team_size = 128;
2357 else if (blksize <= 20)
2358 total_team_size = 160;
2360 total_team_size = 160;
2361 return 2 * total_team_size / vector_size;
2364 struct ExtractAndFactorizeTridiagsDefaultModeAndAlgo<Kokkos::HIPSpace> {
2365 typedef KB::Mode::Team mode_type;
2366 typedef KB::Algo::Level3::Unblocked algo_type;
2367 static int recommended_team_size(
const int blksize,
2368 const int vector_length,
2369 const int internal_vector_length) {
2370 return ExtractAndFactorizeRecommendedHIPTeamSize(blksize, vector_length, internal_vector_length);
2374 struct ExtractAndFactorizeTridiagsDefaultModeAndAlgo<Kokkos::HIPHostPinnedSpace> {
2375 typedef KB::Mode::Team mode_type;
2376 typedef KB::Algo::Level3::Unblocked algo_type;
2377 static int recommended_team_size(
const int blksize,
2378 const int vector_length,
2379 const int internal_vector_length) {
2380 return ExtractAndFactorizeRecommendedHIPTeamSize(blksize, vector_length, internal_vector_length);
2385 #if defined(KOKKOS_ENABLE_SYCL)
2386 static inline int ExtractAndFactorizeRecommendedSYCLTeamSize(
const int blksize,
2387 const int vector_length,
2388 const int internal_vector_length) {
2389 const int vector_size = vector_length / internal_vector_length;
2390 int total_team_size(0);
2392 total_team_size = 32;
2393 else if (blksize <= 9)
2394 total_team_size = 32;
2395 else if (blksize <= 12)
2396 total_team_size = 96;
2397 else if (blksize <= 16)
2398 total_team_size = 128;
2399 else if (blksize <= 20)
2400 total_team_size = 160;
2402 total_team_size = 160;
2403 return 2 * total_team_size / vector_size;
2406 struct ExtractAndFactorizeTridiagsDefaultModeAndAlgo<Kokkos::Experimental::SYCLDeviceUSMSpace> {
2407 typedef KB::Mode::Team mode_type;
2408 typedef KB::Algo::Level3::Unblocked algo_type;
2409 static int recommended_team_size(
const int blksize,
2410 const int vector_length,
2411 const int internal_vector_length) {
2412 return ExtractAndFactorizeRecommendedSYCLTeamSize(blksize, vector_length, internal_vector_length);
2416 struct ExtractAndFactorizeTridiagsDefaultModeAndAlgo<Kokkos::Experimental::SYCLSharedUSMSpace> {
2417 typedef KB::Mode::Team mode_type;
2418 typedef KB::Algo::Level3::Unblocked algo_type;
2419 static int recommended_team_size(
const int blksize,
2420 const int vector_length,
2421 const int internal_vector_length) {
2422 return ExtractAndFactorizeRecommendedSYCLTeamSize(blksize, vector_length, internal_vector_length);
2427 template <
typename impl_type,
typename WWViewType>
2428 KOKKOS_INLINE_FUNCTION
void
2429 solveMultiVector(
const typename Kokkos::TeamPolicy<typename impl_type::execution_space>::member_type &member,
2430 const typename impl_type::local_ordinal_type & ,
2431 const typename impl_type::local_ordinal_type &i0,
2432 const typename impl_type::local_ordinal_type &r0,
2433 const typename impl_type::local_ordinal_type &nrows,
2434 const typename impl_type::local_ordinal_type &v,
2435 const ConstUnmanaged<typename impl_type::internal_vector_type_4d_view> D_internal_vector_values,
2436 const Unmanaged<typename impl_type::internal_vector_type_4d_view> X_internal_vector_values,
2437 const WWViewType &WW,
2438 const bool skip_first_pass =
false) {
2439 using execution_space =
typename impl_type::execution_space;
2440 using team_policy_type = Kokkos::TeamPolicy<execution_space>;
2441 using member_type =
typename team_policy_type::member_type;
2442 using local_ordinal_type =
typename impl_type::local_ordinal_type;
2444 typedef SolveTridiagsDefaultModeAndAlgo<typename execution_space::memory_space> default_mode_and_algo_type;
2446 typedef typename default_mode_and_algo_type::mode_type default_mode_type;
2447 typedef typename default_mode_and_algo_type::multi_vector_algo_type default_algo_type;
2449 using btdm_magnitude_type =
typename impl_type::btdm_magnitude_type;
2452 const auto one = Kokkos::ArithTraits<btdm_magnitude_type>::one();
2453 const auto zero = Kokkos::ArithTraits<btdm_magnitude_type>::zero();
2456 auto A = Kokkos::subview(D_internal_vector_values, i0, Kokkos::ALL(), Kokkos::ALL(), v);
2457 auto X1 = Kokkos::subview(X_internal_vector_values, r0, Kokkos::ALL(), Kokkos::ALL(), v);
2460 local_ordinal_type i = i0, r = r0;
2464 if (skip_first_pass) {
2465 i += (nrows - 2) * 3;
2467 A.assign_data(&D_internal_vector_values(i + 2, 0, 0, v));
2468 X2.assign_data(&X_internal_vector_values(++r, 0, 0, v));
2469 A.assign_data(&D_internal_vector_values(i + 3, 0, 0, v));
2470 KB::Trsm<member_type,
2471 KB::Side::Left, KB::Uplo::Lower, KB::Trans::NoTranspose, KB::Diag::Unit,
2472 default_mode_type, default_algo_type>::invoke(member, one, A, X2);
2473 X1.assign_data(X2.data());
2476 KB::Trsm<member_type,
2477 KB::Side::Left, KB::Uplo::Lower, KB::Trans::NoTranspose, KB::Diag::Unit,
2478 default_mode_type, default_algo_type>::invoke(member, one, A, X1);
2479 for (local_ordinal_type tr = 1; tr < nrows; ++tr, i += 3) {
2480 A.assign_data(&D_internal_vector_values(i + 2, 0, 0, v));
2481 X2.assign_data(&X_internal_vector_values(++r, 0, 0, v));
2482 member.team_barrier();
2483 KB::Gemm<member_type,
2484 KB::Trans::NoTranspose, KB::Trans::NoTranspose,
2485 default_mode_type, default_algo_type>::invoke(member, -one, A, X1, one, X2);
2486 A.assign_data(&D_internal_vector_values(i + 3, 0, 0, v));
2487 KB::Trsm<member_type,
2488 KB::Side::Left, KB::Uplo::Lower, KB::Trans::NoTranspose, KB::Diag::Unit,
2489 default_mode_type, default_algo_type>::invoke(member, one, A, X2);
2490 X1.assign_data(X2.data());
2495 KB::Trsm<member_type,
2496 KB::Side::Left, KB::Uplo::Upper, KB::Trans::NoTranspose, KB::Diag::NonUnit,
2497 default_mode_type, default_algo_type>::invoke(member, one, A, X1);
2498 for (local_ordinal_type tr = nrows; tr > 1; --tr) {
2500 A.assign_data(&D_internal_vector_values(i + 1, 0, 0, v));
2501 X2.assign_data(&X_internal_vector_values(--r, 0, 0, v));
2502 member.team_barrier();
2503 KB::Gemm<member_type,
2504 KB::Trans::NoTranspose, KB::Trans::NoTranspose,
2505 default_mode_type, default_algo_type>::invoke(member, -one, A, X1, one, X2);
2507 A.assign_data(&D_internal_vector_values(i, 0, 0, v));
2508 KB::Trsm<member_type,
2509 KB::Side::Left, KB::Uplo::Upper, KB::Trans::NoTranspose, KB::Diag::NonUnit,
2510 default_mode_type, default_algo_type>::invoke(member, one, A, X2);
2511 X1.assign_data(X2.data());
2515 auto W = Kokkos::subview(WW, Kokkos::ALL(), Kokkos::ALL(), v);
2516 KB::Copy<member_type, KB::Trans::NoTranspose, default_mode_type>::invoke(member, X1, W);
2517 member.team_barrier();
2518 KB::Gemm<member_type,
2519 KB::Trans::NoTranspose, KB::Trans::NoTranspose,
2520 default_mode_type, default_algo_type>::invoke(member, one, A, W, zero, X1);
2524 template <
typename impl_type,
typename WWViewType,
typename XViewType>
2525 KOKKOS_INLINE_FUNCTION
void
2526 solveSingleVectorNew(
const typename Kokkos::TeamPolicy<typename impl_type::execution_space>::member_type &member,
2527 const typename impl_type::local_ordinal_type &blocksize,
2528 const typename impl_type::local_ordinal_type &i0,
2529 const typename impl_type::local_ordinal_type &r0,
2530 const typename impl_type::local_ordinal_type &nrows,
2531 const typename impl_type::local_ordinal_type &v,
2532 const ConstUnmanaged<typename impl_type::internal_vector_type_4d_view> D_internal_vector_values,
2533 const XViewType &X_internal_vector_values,
2534 const WWViewType &WW) {
2535 using execution_space =
typename impl_type::execution_space;
2538 using local_ordinal_type =
typename impl_type::local_ordinal_type;
2540 typedef SolveTridiagsDefaultModeAndAlgo<typename execution_space::memory_space> default_mode_and_algo_type;
2542 typedef typename default_mode_and_algo_type::mode_type default_mode_type;
2543 typedef typename default_mode_and_algo_type::single_vector_algo_type default_algo_type;
2545 using btdm_magnitude_type =
typename impl_type::btdm_magnitude_type;
2548 auto A = D_internal_vector_values.data();
2549 auto X = X_internal_vector_values.data();
2552 const auto one = Kokkos::ArithTraits<btdm_magnitude_type>::one();
2553 const auto zero = Kokkos::ArithTraits<btdm_magnitude_type>::zero();
2557 const local_ordinal_type astep = D_internal_vector_values.stride(0);
2558 const local_ordinal_type as0 = D_internal_vector_values.stride(1);
2559 const local_ordinal_type as1 = D_internal_vector_values.stride(2);
2560 const local_ordinal_type xstep = X_internal_vector_values.stride(0);
2561 const local_ordinal_type xs0 = X_internal_vector_values.stride(1);
2564 A += i0 * astep + v;
2565 X += r0 * xstep + v;
2570 KOKKOSBATCHED_TRSV_LOWER_NO_TRANSPOSE_INTERNAL_INVOKE(default_mode_type, default_algo_type,
2573 blocksize, blocksize,
2578 for (local_ordinal_type tr = 1; tr < nrows; ++tr) {
2579 member.team_barrier();
2580 KOKKOSBATCHED_GEMV_NO_TRANSPOSE_INTERNAL_INVOKE(default_mode_type, default_algo_type,
2582 blocksize, blocksize,
2584 A + 2 * astep, as0, as1,
2587 X + 1 * xstep, xs0);
2588 KOKKOSBATCHED_TRSV_LOWER_NO_TRANSPOSE_INTERNAL_INVOKE(default_mode_type, default_algo_type,
2591 blocksize, blocksize,
2593 A + 3 * astep, as0, as1,
2594 X + 1 * xstep, xs0);
2601 KOKKOSBATCHED_TRSV_UPPER_NO_TRANSPOSE_INTERNAL_INVOKE(default_mode_type, default_algo_type,
2604 blocksize, blocksize,
2609 for (local_ordinal_type tr = nrows; tr > 1; --tr) {
2611 member.team_barrier();
2612 KOKKOSBATCHED_GEMV_NO_TRANSPOSE_INTERNAL_INVOKE(default_mode_type, default_algo_type,
2614 blocksize, blocksize,
2616 A + 1 * astep, as0, as1,
2619 X - 1 * xstep, xs0);
2620 KOKKOSBATCHED_TRSV_UPPER_NO_TRANSPOSE_INTERNAL_INVOKE(default_mode_type, default_algo_type,
2623 blocksize, blocksize,
2626 X - 1 * xstep, xs0);
2632 const local_ordinal_type ws0 = WW.stride(0);
2633 auto W = WW.data() + v;
2634 KOKKOSBATCHED_COPY_VECTOR_NO_TRANSPOSE_INTERNAL_INVOKE(default_mode_type,
2635 member, blocksize, X, xs0, W, ws0);
2636 member.team_barrier();
2637 KOKKOSBATCHED_GEMV_NO_TRANSPOSE_INTERNAL_INVOKE(default_mode_type, default_algo_type,
2639 blocksize, blocksize,
2648 template <
typename local_ordinal_type,
typename ViewType>
2649 void writeBTDValuesToFile(
const local_ordinal_type &n_parts,
const ViewType &scalar_values_device, std::string fileName) {
2650 #ifdef IFPACK2_BLOCKTRIDICONTAINER_WRITE_MM
2651 auto scalar_values = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), scalar_values_device);
2652 std::ofstream myfile;
2653 myfile.open(fileName);
2655 const local_ordinal_type n_parts_per_pack = n_parts < (local_ordinal_type)scalar_values.extent(3) ? n_parts : scalar_values.extent(3);
2656 local_ordinal_type nnz = scalar_values.extent(0) * scalar_values.extent(1) * scalar_values.extent(2) * n_parts_per_pack;
2657 const local_ordinal_type n_blocks = scalar_values.extent(0) * n_parts_per_pack;
2658 const local_ordinal_type n_blocks_per_part = n_blocks / n_parts;
2660 const local_ordinal_type block_size = scalar_values.extent(1);
2662 const local_ordinal_type n_rows_per_part = (n_blocks_per_part + 2) / 3 * block_size;
2663 const local_ordinal_type n_rows = n_rows_per_part * n_parts;
2665 const local_ordinal_type n_packs = ceil(
float(n_parts) / n_parts_per_pack);
2667 myfile <<
"%%MatrixMarket matrix coordinate real general" << std::endl;
2668 myfile <<
"%%nnz = " << nnz;
2669 myfile <<
" block size = " << block_size;
2670 myfile <<
" number of blocks = " << n_blocks;
2671 myfile <<
" number of parts = " << n_parts;
2672 myfile <<
" number of blocks per part = " << n_blocks_per_part;
2673 myfile <<
" number of rows = " << n_rows;
2674 myfile <<
" number of cols = " << n_rows;
2675 myfile <<
" number of packs = " << n_packs << std::endl;
2677 myfile << n_rows <<
" " << n_rows <<
" " << nnz << std::setprecision(9) << std::endl;
2679 local_ordinal_type current_part_idx, current_block_idx, current_row_offset, current_col_offset, current_row, current_col;
2680 for (local_ordinal_type i_pack = 0; i_pack < n_packs; ++i_pack) {
2681 for (local_ordinal_type i_part_in_pack = 0; i_part_in_pack < n_parts_per_pack; ++i_part_in_pack) {
2682 current_part_idx = i_part_in_pack + i_pack * n_parts_per_pack;
2683 for (local_ordinal_type i_block_in_part = 0; i_block_in_part < n_blocks_per_part; ++i_block_in_part) {
2684 current_block_idx = i_block_in_part + i_pack * n_blocks_per_part;
2685 if (current_block_idx >= (local_ordinal_type)scalar_values.extent(0))
2687 if (i_block_in_part % 3 == 0) {
2688 current_row_offset = i_block_in_part / 3 * block_size;
2689 current_col_offset = i_block_in_part / 3 * block_size;
2690 }
else if (i_block_in_part % 3 == 1) {
2691 current_row_offset = (i_block_in_part - 1) / 3 * block_size;
2692 current_col_offset = ((i_block_in_part - 1) / 3 + 1) * block_size;
2693 }
else if (i_block_in_part % 3 == 2) {
2694 current_row_offset = ((i_block_in_part - 2) / 3 + 1) * block_size;
2695 current_col_offset = (i_block_in_part - 2) / 3 * block_size;
2697 current_row_offset += current_part_idx * n_rows_per_part;
2698 current_col_offset += current_part_idx * n_rows_per_part;
2699 for (local_ordinal_type i_in_block = 0; i_in_block < block_size; ++i_in_block) {
2700 for (local_ordinal_type j_in_block = 0; j_in_block < block_size; ++j_in_block) {
2701 current_row = current_row_offset + i_in_block + 1;
2702 current_col = current_col_offset + j_in_block + 1;
2703 myfile << current_row <<
" " << current_col <<
" " << scalar_values(current_block_idx, i_in_block, j_in_block, i_part_in_pack) << std::endl;
2714 template <
typename local_ordinal_type,
typename ViewType>
2715 void write4DMultiVectorValuesToFile(
const local_ordinal_type &n_parts,
const ViewType &scalar_values_device, std::string fileName) {
2716 #ifdef IFPACK2_BLOCKTRIDICONTAINER_WRITE_MM
2717 auto scalar_values = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), scalar_values_device);
2718 std::ofstream myfile;
2719 myfile.open(fileName);
2721 const local_ordinal_type n_parts_per_pack = n_parts < scalar_values.extent(3) ? n_parts : scalar_values.extent(3);
2722 const local_ordinal_type n_blocks = scalar_values.extent(0) * n_parts_per_pack;
2723 const local_ordinal_type n_blocks_per_part = n_blocks / n_parts;
2725 const local_ordinal_type block_size = scalar_values.extent(1);
2726 const local_ordinal_type n_cols = scalar_values.extent(2);
2728 const local_ordinal_type n_rows_per_part = n_blocks_per_part * block_size;
2729 const local_ordinal_type n_rows = n_rows_per_part * n_parts;
2731 const local_ordinal_type n_packs = ceil(
float(n_parts) / n_parts_per_pack);
2733 myfile <<
"%%MatrixMarket matrix array real general" << std::endl;
2734 myfile <<
"%%block size = " << block_size;
2735 myfile <<
" number of blocks = " << n_blocks;
2736 myfile <<
" number of parts = " << n_parts;
2737 myfile <<
" number of blocks per part = " << n_blocks_per_part;
2738 myfile <<
" number of rows = " << n_rows;
2739 myfile <<
" number of cols = " << n_cols;
2740 myfile <<
" number of packs = " << n_packs << std::endl;
2742 myfile << n_rows <<
" " << n_cols << std::setprecision(9) << std::endl;
2744 local_ordinal_type current_part_idx, current_block_idx, current_row_offset;
2745 (void)current_row_offset;
2746 (void)current_part_idx;
2747 for (local_ordinal_type j_in_block = 0; j_in_block < n_cols; ++j_in_block) {
2748 for (local_ordinal_type i_pack = 0; i_pack < n_packs; ++i_pack) {
2749 for (local_ordinal_type i_part_in_pack = 0; i_part_in_pack < n_parts_per_pack; ++i_part_in_pack) {
2750 current_part_idx = i_part_in_pack + i_pack * n_parts_per_pack;
2751 for (local_ordinal_type i_block_in_part = 0; i_block_in_part < n_blocks_per_part; ++i_block_in_part) {
2752 current_block_idx = i_block_in_part + i_pack * n_blocks_per_part;
2754 if (current_block_idx >= (local_ordinal_type)scalar_values.extent(0))
2756 for (local_ordinal_type i_in_block = 0; i_in_block < block_size; ++i_in_block) {
2757 myfile << scalar_values(current_block_idx, i_in_block, j_in_block, i_part_in_pack) << std::endl;
2767 template <
typename local_ordinal_type,
typename ViewType>
2768 void write5DMultiVectorValuesToFile(
const local_ordinal_type &n_parts,
const ViewType &scalar_values_device, std::string fileName) {
2769 #ifdef IFPACK2_BLOCKTRIDICONTAINER_WRITE_MM
2770 auto scalar_values = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), scalar_values_device);
2771 std::ofstream myfile;
2772 myfile.open(fileName);
2774 const local_ordinal_type n_parts_per_pack = n_parts < scalar_values.extent(4) ? n_parts : scalar_values.extent(4);
2775 const local_ordinal_type n_blocks = scalar_values.extent(1) * n_parts_per_pack;
2776 const local_ordinal_type n_blocks_per_part = n_blocks / n_parts;
2778 const local_ordinal_type block_size = scalar_values.extent(2);
2779 const local_ordinal_type n_blocks_cols = scalar_values.extent(0);
2780 const local_ordinal_type n_cols = n_blocks_cols * block_size;
2782 const local_ordinal_type n_rows_per_part = n_blocks_per_part * block_size;
2783 const local_ordinal_type n_rows = n_rows_per_part * n_parts;
2785 const local_ordinal_type n_packs = ceil(
float(n_parts) / n_parts_per_pack);
2787 myfile <<
"%%MatrixMarket matrix array real general" << std::endl;
2788 myfile <<
"%%block size = " << block_size;
2789 myfile <<
" number of blocks = " << n_blocks;
2790 myfile <<
" number of parts = " << n_parts;
2791 myfile <<
" number of blocks per part = " << n_blocks_per_part;
2792 myfile <<
" number of rows = " << n_rows;
2793 myfile <<
" number of cols = " << n_cols;
2794 myfile <<
" number of packs = " << n_packs << std::endl;
2796 myfile << n_rows <<
" " << n_cols << std::setprecision(9) << std::endl;
2798 local_ordinal_type current_part_idx, current_block_idx, current_row_offset;
2799 (void)current_row_offset;
2800 (void)current_part_idx;
2801 for (local_ordinal_type i_block_col = 0; i_block_col < n_blocks_cols; ++i_block_col) {
2802 for (local_ordinal_type j_in_block = 0; j_in_block < block_size; ++j_in_block) {
2803 for (local_ordinal_type i_pack = 0; i_pack < n_packs; ++i_pack) {
2804 for (local_ordinal_type i_part_in_pack = 0; i_part_in_pack < n_parts_per_pack; ++i_part_in_pack) {
2805 current_part_idx = i_part_in_pack + i_pack * n_parts_per_pack;
2806 for (local_ordinal_type i_block_in_part = 0; i_block_in_part < n_blocks_per_part; ++i_block_in_part) {
2807 current_block_idx = i_block_in_part + i_pack * n_blocks_per_part;
2809 if (current_block_idx >= (local_ordinal_type)scalar_values.extent(1))
2811 for (local_ordinal_type i_in_block = 0; i_in_block < block_size; ++i_in_block) {
2812 myfile << scalar_values(i_block_col, current_block_idx, i_in_block, j_in_block, i_part_in_pack) << std::endl;
2823 template <
typename local_ordinal_type,
typename member_type,
typename ViewType1,
typename ViewType2>
2824 KOKKOS_INLINE_FUNCTION
void
2825 copy3DView(
const member_type &member,
const ViewType1 &view1,
const ViewType2 &view2) {
2838 Kokkos::Experimental::local_deep_copy(member, view1, view2);
2840 template <
typename MatrixType,
int ScratchLevel>
2841 struct ExtractAndFactorizeTridiags {
2843 using impl_type = BlockHelperDetails::ImplType<MatrixType>;
2845 using execution_space =
typename impl_type::execution_space;
2846 using memory_space =
typename impl_type::memory_space;
2848 using local_ordinal_type =
typename impl_type::local_ordinal_type;
2851 using magnitude_type =
typename impl_type::magnitude_type;
2853 using row_matrix_type =
typename impl_type::tpetra_row_matrix_type;
2854 using crs_graph_type =
typename impl_type::tpetra_crs_graph_type;
2856 using local_ordinal_type_1d_view =
typename impl_type::local_ordinal_type_1d_view;
2857 using local_ordinal_type_2d_view =
typename impl_type::local_ordinal_type_2d_view;
2859 using size_type_2d_view =
typename impl_type::size_type_2d_view;
2860 using impl_scalar_type_1d_view_tpetra =
typename impl_type::impl_scalar_type_1d_view_tpetra;
2862 using btdm_scalar_type =
typename impl_type::btdm_scalar_type;
2863 using btdm_magnitude_type =
typename impl_type::btdm_magnitude_type;
2864 using vector_type_3d_view =
typename impl_type::vector_type_3d_view;
2865 using vector_type_4d_view =
typename impl_type::vector_type_4d_view;
2866 using internal_vector_type_4d_view =
typename impl_type::internal_vector_type_4d_view;
2867 using internal_vector_type_5d_view =
typename impl_type::internal_vector_type_5d_view;
2868 using btdm_scalar_type_2d_view =
typename impl_type::btdm_scalar_type_2d_view;
2869 using btdm_scalar_type_3d_view =
typename impl_type::btdm_scalar_type_3d_view;
2870 using btdm_scalar_type_4d_view =
typename impl_type::btdm_scalar_type_4d_view;
2871 using btdm_scalar_type_5d_view =
typename impl_type::btdm_scalar_type_5d_view;
2872 using internal_vector_scratch_type_3d_view = Scratch<typename impl_type::internal_vector_type_3d_view>;
2873 using btdm_scalar_scratch_type_3d_view = Scratch<typename impl_type::btdm_scalar_type_3d_view>;
2874 using tpetra_block_access_view_type =
typename impl_type::tpetra_block_access_view_type;
2875 using local_crs_graph_type =
typename impl_type::local_crs_graph_type;
2876 using colinds_view =
typename local_crs_graph_type::entries_type;
2878 using internal_vector_type =
typename impl_type::internal_vector_type;
2879 static constexpr
int vector_length = impl_type::vector_length;
2880 static constexpr
int internal_vector_length = impl_type::internal_vector_length;
2881 static_assert(vector_length >= internal_vector_length,
"Ifpack2 BlockTriDi Numeric: vector_length must be at least as large as internal_vector_length");
2882 static_assert(vector_length % internal_vector_length == 0,
"Ifpack2 BlockTriDi Numeric: vector_length must be divisible by internal_vector_length");
2887 static constexpr
int half_vector_length = impl_type::half_vector_length;
2890 using team_policy_type = Kokkos::TeamPolicy<execution_space>;
2891 using member_type =
typename team_policy_type::member_type;
2895 const ConstUnmanaged<local_ordinal_type_1d_view> partptr, lclrow, packptr, packindices_sub, packptr_sub;
2896 const ConstUnmanaged<local_ordinal_type_2d_view> partptr_sub, part2packrowidx0_sub, packindices_schur;
2897 const local_ordinal_type max_partsz;
2899 using size_type_1d_view_tpetra = Kokkos::View<size_t *, typename impl_type::node_device_type>;
2900 ConstUnmanaged<size_type_1d_view_tpetra> A_block_rowptr;
2901 ConstUnmanaged<size_type_1d_view_tpetra> A_point_rowptr;
2902 ConstUnmanaged<impl_scalar_type_1d_view_tpetra> A_values;
2904 const ConstUnmanaged<size_type_2d_view> pack_td_ptr, flat_td_ptr, pack_td_ptr_schur;
2905 const ConstUnmanaged<local_ordinal_type_1d_view> A_colindsub;
2906 const Unmanaged<internal_vector_type_4d_view> internal_vector_values, internal_vector_values_schur;
2907 const Unmanaged<internal_vector_type_5d_view> e_internal_vector_values;
2908 const Unmanaged<btdm_scalar_type_4d_view> scalar_values, scalar_values_schur;
2909 const Unmanaged<btdm_scalar_type_5d_view> e_scalar_values;
2910 const Unmanaged<btdm_scalar_type_3d_view> d_inv;
2911 const Unmanaged<size_type_1d_view> diag_offsets;
2913 const local_ordinal_type blocksize, blocksize_square;
2915 const magnitude_type tiny;
2916 const local_ordinal_type vector_loop_size;
2918 bool hasBlockCrsMatrix;
2921 ExtractAndFactorizeTridiags(
const BlockTridiags<MatrixType> &btdm_,
2922 const BlockHelperDetails::PartInterface<MatrixType> &interf_,
2925 const magnitude_type &tiny_)
2927 partptr(interf_.partptr)
2928 , lclrow(interf_.lclrow)
2929 , packptr(interf_.packptr)
2930 , packindices_sub(interf_.packindices_sub)
2931 , packptr_sub(interf_.packptr_sub)
2932 , partptr_sub(interf_.partptr_sub)
2933 , part2packrowidx0_sub(interf_.part2packrowidx0_sub)
2934 , packindices_schur(interf_.packindices_schur)
2935 , max_partsz(interf_.max_partsz)
2938 pack_td_ptr(btdm_.pack_td_ptr)
2939 , flat_td_ptr(btdm_.flat_td_ptr)
2940 , pack_td_ptr_schur(btdm_.pack_td_ptr_schur)
2941 , A_colindsub(btdm_.A_colindsub)
2942 , internal_vector_values((internal_vector_type *)btdm_.values.data(),
2943 btdm_.values.extent(0),
2944 btdm_.values.extent(1),
2945 btdm_.values.extent(2),
2946 vector_length / internal_vector_length)
2947 , internal_vector_values_schur((internal_vector_type *)btdm_.values_schur.data(),
2948 btdm_.values_schur.extent(0),
2949 btdm_.values_schur.extent(1),
2950 btdm_.values_schur.extent(2),
2951 vector_length / internal_vector_length)
2952 , e_internal_vector_values((internal_vector_type *)btdm_.e_values.data(),
2953 btdm_.e_values.extent(0),
2954 btdm_.e_values.extent(1),
2955 btdm_.e_values.extent(2),
2956 btdm_.e_values.extent(3),
2957 vector_length / internal_vector_length)
2958 , scalar_values((btdm_scalar_type *)btdm_.values.data(),
2959 btdm_.values.extent(0),
2960 btdm_.values.extent(1),
2961 btdm_.values.extent(2),
2963 , scalar_values_schur((btdm_scalar_type *)btdm_.values_schur.data(),
2964 btdm_.values_schur.extent(0),
2965 btdm_.values_schur.extent(1),
2966 btdm_.values_schur.extent(2),
2968 , e_scalar_values((btdm_scalar_type *)btdm_.e_values.data(),
2969 btdm_.e_values.extent(0),
2970 btdm_.e_values.extent(1),
2971 btdm_.e_values.extent(2),
2972 btdm_.e_values.extent(3),
2974 , d_inv(btdm_.d_inv)
2975 , diag_offsets(btdm_.diag_offsets)
2976 , blocksize(btdm_.values.extent(1))
2977 , blocksize_square(blocksize * blocksize)
2981 , vector_loop_size(vector_length / internal_vector_length) {
2982 using crs_matrix_type =
typename impl_type::tpetra_crs_matrix_type;
2983 using block_crs_matrix_type =
typename impl_type::tpetra_block_crs_matrix_type;
2985 auto A_crs = Teuchos::rcp_dynamic_cast<
const crs_matrix_type>(A_);
2986 auto A_bcrs = Teuchos::rcp_dynamic_cast<
const block_crs_matrix_type>(A_);
2988 hasBlockCrsMatrix = !A_bcrs.is_null();
2990 A_block_rowptr = G_->getLocalGraphDevice().row_map;
2991 if (hasBlockCrsMatrix) {
2992 A_values =
const_cast<block_crs_matrix_type *
>(A_bcrs.get())->getValuesDeviceNonConst();
2994 A_point_rowptr = A_crs->getCrsGraph()->getLocalGraphDevice().row_map;
2995 A_values = A_crs->getLocalValuesDevice(Tpetra::Access::ReadOnly);
3000 KOKKOS_INLINE_FUNCTION
3002 extract(local_ordinal_type partidx,
3003 local_ordinal_type local_subpartidx,
3004 local_ordinal_type npacks)
const {
3005 #ifdef IFPACK2_BLOCKTRIDICONTAINER_USE_PRINTF
3006 printf(
"extract partidx = %d, local_subpartidx = %d, npacks = %d;\n", partidx, local_subpartidx, npacks);
3008 using tlb = BlockHelperDetails::TpetraLittleBlock<Tpetra::Impl::BlockCrsMatrixLittleBlockArrayLayout>;
3009 const size_type kps = pack_td_ptr(partidx, local_subpartidx);
3010 local_ordinal_type kfs[vector_length] = {};
3011 local_ordinal_type ri0[vector_length] = {};
3012 local_ordinal_type nrows[vector_length] = {};
3014 for (local_ordinal_type vi = 0; vi < npacks; ++vi, ++partidx) {
3015 kfs[vi] = flat_td_ptr(partidx, local_subpartidx);
3016 ri0[vi] = partptr_sub(pack_td_ptr.extent(0) * local_subpartidx + partidx, 0);
3017 nrows[vi] = partptr_sub(pack_td_ptr.extent(0) * local_subpartidx + partidx, 1) - ri0[vi];
3018 #ifdef IFPACK2_BLOCKTRIDICONTAINER_USE_PRINTF
3019 printf(
"kfs[%d] = %d;\n", vi, kfs[vi]);
3020 printf(
"ri0[%d] = %d;\n", vi, ri0[vi]);
3021 printf(
"nrows[%d] = %d;\n", vi, nrows[vi]);
3024 local_ordinal_type tr_min = 0;
3025 local_ordinal_type tr_max = nrows[0];
3026 if (local_subpartidx % 2 == 1) {
3030 #ifdef IFPACK2_BLOCKTRIDICONTAINER_USE_PRINTF
3031 printf(
"tr_min = %d and tr_max = %d;\n", tr_min, tr_max);
3033 for (local_ordinal_type tr = tr_min, j = 0; tr < tr_max; ++tr) {
3034 for (local_ordinal_type e = 0; e < 3; ++e) {
3035 if (hasBlockCrsMatrix) {
3036 const impl_scalar_type *block[vector_length] = {};
3037 for (local_ordinal_type vi = 0; vi < npacks; ++vi) {
3038 const size_type Aj = A_block_rowptr(lclrow(ri0[vi] + tr)) + A_colindsub(kfs[vi] + j);
3040 block[vi] = &A_values(Aj * blocksize_square);
3042 const size_type pi = kps + j;
3043 #ifdef IFPACK2_BLOCKTRIDICONTAINER_USE_PRINTF
3044 printf(
"Extract pi = %ld, ri0 + tr = %d, kfs + j = %d\n", pi, ri0[0] + tr, kfs[0] + j);
3047 for (local_ordinal_type ii = 0; ii < blocksize; ++ii) {
3048 for (local_ordinal_type jj = 0; jj < blocksize; ++jj) {
3049 const auto idx = tlb::getFlatIndex(ii, jj, blocksize);
3050 auto &v = internal_vector_values(pi, ii, jj, 0);
3051 for (local_ordinal_type vi = 0; vi < npacks; ++vi) {
3052 v[vi] =
static_cast<btdm_scalar_type
>(block[vi][idx]);
3057 const size_type pi = kps + j;
3059 for (local_ordinal_type vi = 0; vi < npacks; ++vi) {
3060 const size_type Aj_c = A_colindsub(kfs[vi] + j);
3062 for (local_ordinal_type ii = 0; ii < blocksize; ++ii) {
3063 auto point_row_offset = A_point_rowptr(lclrow(ri0[vi] + tr) * blocksize + ii);
3065 for (local_ordinal_type jj = 0; jj < blocksize; ++jj) {
3066 scalar_values(pi, ii, jj, vi) = A_values(point_row_offset + Aj_c * blocksize + jj);
3072 if (nrows[0] == 1)
break;
3073 if (local_subpartidx % 2 == 0) {
3074 if (e == 1 && (tr == 0 || tr + 1 == nrows[0]))
break;
3075 for (local_ordinal_type vi = 1; vi < npacks; ++vi) {
3076 if ((e == 0 && nrows[vi] == 1) || (e == 1 && tr + 1 == nrows[vi])) {
3082 if (e == 0 && (tr == -1 || tr == nrows[0]))
break;
3083 for (local_ordinal_type vi = 1; vi < npacks; ++vi) {
3084 if ((e == 0 && nrows[vi] == 1) || (e == 0 && tr == nrows[vi])) {
3094 KOKKOS_INLINE_FUNCTION
3096 extract(
const member_type &member,
3097 const local_ordinal_type &partidxbeg,
3098 local_ordinal_type local_subpartidx,
3099 const local_ordinal_type &npacks,
3100 const local_ordinal_type &vbeg)
const {
3101 #ifdef IFPACK2_BLOCKTRIDICONTAINER_USE_PRINTF
3102 printf(
"extract partidxbeg = %d, local_subpartidx = %d, npacks = %d, vbeg = %d;\n", partidxbeg, local_subpartidx, npacks, vbeg);
3104 using tlb = BlockHelperDetails::TpetraLittleBlock<Tpetra::Impl::BlockCrsMatrixLittleBlockArrayLayout>;
3105 local_ordinal_type kfs_vals[internal_vector_length] = {};
3106 local_ordinal_type ri0_vals[internal_vector_length] = {};
3107 local_ordinal_type nrows_vals[internal_vector_length] = {};
3109 const size_type kps = pack_td_ptr(partidxbeg, local_subpartidx);
3110 for (local_ordinal_type v = vbeg, vi = 0; v < npacks && vi < internal_vector_length; ++v, ++vi) {
3111 kfs_vals[vi] = flat_td_ptr(partidxbeg + vi, local_subpartidx);
3112 ri0_vals[vi] = partptr_sub(pack_td_ptr.extent(0) * local_subpartidx + partidxbeg + vi, 0);
3113 nrows_vals[vi] = partptr_sub(pack_td_ptr.extent(0) * local_subpartidx + partidxbeg + vi, 1) - ri0_vals[vi];
3114 #ifdef IFPACK2_BLOCKTRIDICONTAINER_USE_PRINTF
3115 printf(
"kfs_vals[%d] = %d;\n", vi, kfs_vals[vi]);
3116 printf(
"ri0_vals[%d] = %d;\n", vi, ri0_vals[vi]);
3117 printf(
"nrows_vals[%d] = %d;\n", vi, nrows_vals[vi]);
3121 local_ordinal_type j_vals[internal_vector_length] = {};
3123 local_ordinal_type tr_min = 0;
3124 local_ordinal_type tr_max = nrows_vals[0];
3125 if (local_subpartidx % 2 == 1) {
3129 #ifdef IFPACK2_BLOCKTRIDICONTAINER_USE_PRINTF
3130 printf(
"tr_min = %d and tr_max = %d;\n", tr_min, tr_max);
3132 for (local_ordinal_type tr = tr_min; tr < tr_max; ++tr) {
3133 for (local_ordinal_type v = vbeg, vi = 0; v < npacks && vi < internal_vector_length; ++v, ++vi) {
3134 const local_ordinal_type nrows = (local_subpartidx % 2 == 0 ? nrows_vals[vi] : nrows_vals[vi]);
3135 if ((local_subpartidx % 2 == 0 && tr < nrows) || (local_subpartidx % 2 == 1 && tr < nrows + 1)) {
3136 auto &j = j_vals[vi];
3137 const local_ordinal_type kfs = kfs_vals[vi];
3138 const local_ordinal_type ri0 = ri0_vals[vi];
3139 local_ordinal_type lbeg, lend;
3140 if (local_subpartidx % 2 == 0) {
3141 lbeg = (tr == tr_min ? 1 : 0);
3142 lend = (tr == nrows - 1 ? 2 : 3);
3149 }
else if (tr == nrows) {
3154 if (hasBlockCrsMatrix) {
3155 for (local_ordinal_type l = lbeg; l < lend; ++l, ++j) {
3156 const size_type Aj = A_block_rowptr(lclrow(ri0 + tr)) + A_colindsub(kfs + j);
3157 const impl_scalar_type *block = &A_values(Aj * blocksize_square);
3158 const size_type pi = kps + j;
3159 #ifdef IFPACK2_BLOCKTRIDICONTAINER_USE_PRINTF
3160 printf(
"Extract pi = %ld, ri0 + tr = %d, kfs + j = %d, tr = %d, lbeg = %d, lend = %d, l = %d\n", pi, ri0 + tr, kfs + j, tr, lbeg, lend, l);
3162 Kokkos::parallel_for(Kokkos::TeamThreadRange(member, blocksize),
3163 [&](
const local_ordinal_type &ii) {
3164 for (local_ordinal_type jj = 0; jj < blocksize; ++jj) {
3165 scalar_values(pi, ii, jj, v) =
static_cast<btdm_scalar_type
>(block[tlb::getFlatIndex(ii, jj, blocksize)]);
3170 for (local_ordinal_type l = lbeg; l < lend; ++l, ++j) {
3171 const size_type Aj_c = A_colindsub(kfs + j);
3172 const size_type pi = kps + j;
3173 Kokkos::parallel_for(Kokkos::TeamThreadRange(member, blocksize),
3174 [&](
const local_ordinal_type &ii) {
3175 auto point_row_offset = A_point_rowptr(lclrow(ri0 + tr) * blocksize + ii);
3176 for (local_ordinal_type jj = 0; jj < blocksize; ++jj) {
3177 scalar_values(pi, ii, jj, v) = A_values(point_row_offset + Aj_c * blocksize + jj);
3187 template <
typename AAViewType,
3188 typename WWViewType>
3189 KOKKOS_INLINE_FUNCTION
void
3190 factorize_subline(
const member_type &member,
3191 const local_ordinal_type &i0,
3192 const local_ordinal_type &nrows,
3193 const local_ordinal_type &v,
3194 const AAViewType &AA,
3195 const WWViewType &WW)
const {
3196 typedef ExtractAndFactorizeTridiagsDefaultModeAndAlgo<typename execution_space::memory_space> default_mode_and_algo_type;
3198 typedef typename default_mode_and_algo_type::mode_type default_mode_type;
3199 typedef typename default_mode_and_algo_type::algo_type default_algo_type;
3202 const auto one = Kokkos::ArithTraits<btdm_magnitude_type>::one();
3204 #ifdef IFPACK2_BLOCKTRIDICONTAINER_USE_PRINTF
3205 printf(
"i0 = %d, nrows = %d, v = %d, AA.extent(0) = %ld;\n", i0, nrows, v, AA.extent(0));
3209 auto A = Kokkos::subview(AA, i0, Kokkos::ALL(), Kokkos::ALL(), v);
3211 default_mode_type, KB::Algo::LU::Unblocked>::invoke(member, A, tiny);
3216 local_ordinal_type i = i0;
3217 for (local_ordinal_type tr = 1; tr < nrows; ++tr, i += 3) {
3218 #ifdef IFPACK2_BLOCKTRIDICONTAINER_USE_PRINTF
3219 printf(
"tr = %d, i = %d;\n", tr, i);
3221 B.assign_data(&AA(i + 1, 0, 0, v));
3222 KB::Trsm<member_type,
3223 KB::Side::Left, KB::Uplo::Lower, KB::Trans::NoTranspose, KB::Diag::Unit,
3224 default_mode_type, default_algo_type>::invoke(member, one, A, B);
3225 C.assign_data(&AA(i + 2, 0, 0, v));
3226 KB::Trsm<member_type,
3227 KB::Side::Right, KB::Uplo::Upper, KB::Trans::NoTranspose, KB::Diag::NonUnit,
3228 default_mode_type, default_algo_type>::invoke(member, one, A, C);
3229 A.assign_data(&AA(i + 3, 0, 0, v));
3231 member.team_barrier();
3232 KB::Gemm<member_type,
3233 KB::Trans::NoTranspose, KB::Trans::NoTranspose,
3234 default_mode_type, default_algo_type>::invoke(member, -one, C, B, one, A);
3236 default_mode_type, KB::Algo::LU::Unblocked>::invoke(member, A, tiny);
3240 auto W = Kokkos::subview(WW, Kokkos::ALL(), Kokkos::ALL(), v);
3241 KB::Copy<member_type, KB::Trans::NoTranspose, default_mode_type>::invoke(member, A, W);
3242 KB::SetIdentity<member_type, default_mode_type>::invoke(member, A);
3243 member.team_barrier();
3244 KB::Trsm<member_type,
3245 KB::Side::Left, KB::Uplo::Lower, KB::Trans::NoTranspose, KB::Diag::Unit,
3246 default_mode_type, default_algo_type>::invoke(member, one, W, A);
3247 KB::Trsm<member_type,
3248 KB::Side::Left, KB::Uplo::Upper, KB::Trans::NoTranspose, KB::Diag::NonUnit,
3249 default_mode_type, default_algo_type>::invoke(member, one, W, A);
3254 struct ExtractAndFactorizeSubLineTag {};
3255 struct ExtractAndFactorizeFusedJacobiTag {};
3256 struct ExtractBCDTag {};
3257 struct ComputeETag {};
3258 struct ComputeSchurTag {};
3259 struct FactorizeSchurTag {};
3261 KOKKOS_INLINE_FUNCTION
3263 operator()(
const ExtractAndFactorizeSubLineTag &,
const member_type &member)
const {
3265 const local_ordinal_type packidx = packindices_sub(member.league_rank());
3267 const local_ordinal_type subpartidx = packptr_sub(packidx);
3268 const local_ordinal_type n_parts = part2packrowidx0_sub.extent(0);
3269 const local_ordinal_type local_subpartidx = subpartidx / n_parts;
3270 const local_ordinal_type partidx = subpartidx % n_parts;
3272 const local_ordinal_type npacks = packptr_sub(packidx + 1) - subpartidx;
3273 const local_ordinal_type i0 = pack_td_ptr(partidx, local_subpartidx);
3274 const local_ordinal_type nrows = partptr_sub(subpartidx, 1) - partptr_sub(subpartidx, 0);
3276 internal_vector_scratch_type_3d_view
3277 WW(member.team_scratch(ScratchLevel), blocksize, blocksize, vector_loop_size);
3279 #ifdef IFPACK2_BLOCKTRIDICONTAINER_USE_PRINTF
3280 printf(
"rank = %d, i0 = %d, npacks = %d, nrows = %d, packidx = %d, subpartidx = %d, partidx = %d, local_subpartidx = %d;\n", member.league_rank(), i0, npacks, nrows, packidx, subpartidx, partidx, local_subpartidx);
3281 printf(
"vector_loop_size = %d\n", vector_loop_size);
3284 if (vector_loop_size == 1) {
3285 extract(partidx, local_subpartidx, npacks);
3286 factorize_subline(member, i0, nrows, 0, internal_vector_values, WW);
3288 Kokkos::parallel_for(Kokkos::ThreadVectorRange(member, vector_loop_size),
3289 [&](
const local_ordinal_type &v) {
3290 const local_ordinal_type vbeg = v * internal_vector_length;
3291 #ifdef IFPACK2_BLOCKTRIDICONTAINER_USE_PRINTF
3292 printf(
"i0 = %d, npacks = %d, vbeg = %d;\n", i0, npacks, vbeg);
3295 extract(member, partidx + vbeg, local_subpartidx, npacks, vbeg);
3298 member.team_barrier();
3299 factorize_subline(member, i0, nrows, v, internal_vector_values, WW);
3304 KOKKOS_INLINE_FUNCTION
3306 operator()(
const ExtractAndFactorizeFusedJacobiTag &,
const member_type &member)
const {
3307 using default_mode_and_algo_type = ExtractAndFactorizeTridiagsDefaultModeAndAlgo<typename execution_space::memory_space>;
3308 using default_mode_type =
typename default_mode_and_algo_type::mode_type;
3309 using default_algo_type =
typename default_mode_and_algo_type::algo_type;
3312 btdm_scalar_scratch_type_3d_view WW1(member.team_scratch(ScratchLevel), half_vector_length, blocksize, blocksize);
3313 btdm_scalar_scratch_type_3d_view WW2(member.team_scratch(ScratchLevel), half_vector_length, blocksize, blocksize);
3314 const auto one = Kokkos::ArithTraits<btdm_magnitude_type>::one();
3315 const local_ordinal_type nrows = lclrow.extent(0);
3316 Kokkos::parallel_for(Kokkos::ThreadVectorRange(member, half_vector_length),
3317 [&](
const local_ordinal_type &v) {
3318 local_ordinal_type row = member.league_rank() * half_vector_length + v;
3320 auto W1 = Kokkos::subview(WW1, v, Kokkos::ALL(), Kokkos::ALL());
3321 auto W2 = Kokkos::subview(WW2, v, Kokkos::ALL(), Kokkos::ALL());
3324 const impl_scalar_type *A_diag = A_values.data() + diag_offsets(row);
3327 Kokkos::parallel_for(Kokkos::TeamThreadRange(member, blocksize * blocksize),
3329 W1.data()[i] = A_diag[i];
3332 KB::SetIdentity<member_type, default_mode_type>::invoke(member, W2);
3337 KB::SetIdentity<member_type, default_mode_type>::invoke(member, W1);
3339 member.team_barrier();
3341 KB::LU<member_type, default_mode_type, KB::Algo::LU::Unblocked>::invoke(member, W1, tiny);
3342 member.team_barrier();
3343 KB::Trsm<member_type,
3344 KB::Side::Left, KB::Uplo::Lower, KB::Trans::NoTranspose, KB::Diag::Unit,
3345 default_mode_type, default_algo_type>::invoke(member, one, W1, W2);
3346 KB::Trsm<member_type,
3347 KB::Side::Left, KB::Uplo::Upper, KB::Trans::NoTranspose, KB::Diag::NonUnit,
3348 default_mode_type, default_algo_type>::invoke(member, one, W1, W2);
3349 member.team_barrier();
3351 Kokkos::parallel_for(Kokkos::TeamThreadRange(member, blocksize * blocksize),
3353 auto d_inv_block = &d_inv(row, 0, 0);
3354 d_inv_block[i] = W2.data()[i];
3360 KOKKOS_INLINE_FUNCTION
3362 operator()(
const ExtractBCDTag &,
const member_type &member)
const {
3364 const local_ordinal_type packindices_schur_i = member.league_rank() % packindices_schur.extent(0);
3365 const local_ordinal_type packindices_schur_j = member.league_rank() / packindices_schur.extent(0);
3366 const local_ordinal_type packidx = packindices_schur(packindices_schur_i, packindices_schur_j);
3368 const local_ordinal_type subpartidx = packptr_sub(packidx);
3369 const local_ordinal_type n_parts = part2packrowidx0_sub.extent(0);
3370 const local_ordinal_type local_subpartidx = subpartidx / n_parts;
3371 const local_ordinal_type partidx = subpartidx % n_parts;
3373 const local_ordinal_type npacks = packptr_sub(packidx + 1) - subpartidx;
3377 if (vector_loop_size == 1) {
3378 extract(partidx, local_subpartidx, npacks);
3380 Kokkos::parallel_for(Kokkos::ThreadVectorRange(member, vector_loop_size),
3381 [&](
const local_ordinal_type &v) {
3382 const local_ordinal_type vbeg = v * internal_vector_length;
3383 #ifdef IFPACK2_BLOCKTRIDICONTAINER_USE_PRINTF
3384 const local_ordinal_type i0 = pack_td_ptr(partidx, local_subpartidx);
3385 printf(
"i0 = %d, npacks = %d, vbeg = %d;\n", i0, npacks, vbeg);
3388 extract(member, partidx + vbeg, local_subpartidx, npacks, vbeg);
3392 member.team_barrier();
3394 const size_type kps1 = pack_td_ptr(partidx, local_subpartidx);
3395 const size_type kps2 = pack_td_ptr(partidx, local_subpartidx + 1) - 1;
3397 const local_ordinal_type r1 = part2packrowidx0_sub(partidx, local_subpartidx) - 1;
3398 const local_ordinal_type r2 = part2packrowidx0_sub(partidx, local_subpartidx) + 2;
3400 #ifdef IFPACK2_BLOCKTRIDICONTAINER_USE_PRINTF
3401 printf(
"Copy for Schur complement part id = %d from kps1 = %ld to r1 = %d and from kps2 = %ld to r2 = %d partidx = %d local_subpartidx = %d;\n", packidx, kps1, r1, kps2, r2, partidx, local_subpartidx);
3405 copy3DView<local_ordinal_type>(member, Kokkos::subview(e_internal_vector_values, 0, r1, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()),
3406 Kokkos::subview(internal_vector_values, kps1, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()));
3408 copy3DView<local_ordinal_type>(member, Kokkos::subview(e_internal_vector_values, 1, r2, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()),
3409 Kokkos::subview(internal_vector_values, kps2, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()));
3412 KOKKOS_INLINE_FUNCTION
3414 operator()(
const ComputeETag &,
const member_type &member)
const {
3416 const local_ordinal_type packidx = packindices_sub(member.league_rank());
3418 const local_ordinal_type subpartidx = packptr_sub(packidx);
3419 const local_ordinal_type n_parts = part2packrowidx0_sub.extent(0);
3420 const local_ordinal_type local_subpartidx = subpartidx / n_parts;
3421 const local_ordinal_type partidx = subpartidx % n_parts;
3423 const local_ordinal_type npacks = packptr_sub(packidx + 1) - subpartidx;
3424 const local_ordinal_type i0 = pack_td_ptr(partidx, local_subpartidx);
3425 const local_ordinal_type r0 = part2packrowidx0_sub(partidx, local_subpartidx);
3426 const local_ordinal_type nrows = partptr_sub(subpartidx, 1) - partptr_sub(subpartidx, 0);
3427 const local_ordinal_type num_vectors = blocksize;
3431 internal_vector_scratch_type_3d_view
3432 WW(member.team_scratch(ScratchLevel), blocksize, num_vectors, vector_loop_size);
3433 if (local_subpartidx == 0) {
3434 Kokkos::parallel_for(Kokkos::ThreadVectorRange(member, vector_loop_size), [&](
const int &v) {
3435 solveMultiVector<impl_type, internal_vector_scratch_type_3d_view>(member, blocksize, i0, r0, nrows, v, internal_vector_values, Kokkos::subview(e_internal_vector_values, 0, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()), WW,
true);
3437 }
else if (local_subpartidx == (local_ordinal_type)part2packrowidx0_sub.extent(1) - 2) {
3438 Kokkos::parallel_for(Kokkos::ThreadVectorRange(member, vector_loop_size), [&](
const int &v) {
3439 solveMultiVector<impl_type, internal_vector_scratch_type_3d_view>(member, blocksize, i0, r0, nrows, v, internal_vector_values, Kokkos::subview(e_internal_vector_values, 1, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()), WW);
3442 Kokkos::parallel_for(Kokkos::ThreadVectorRange(member, vector_loop_size), [&](
const int &v) {
3443 solveMultiVector<impl_type, internal_vector_scratch_type_3d_view>(member, blocksize, i0, r0, nrows, v, internal_vector_values, Kokkos::subview(e_internal_vector_values, 0, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()), WW,
true);
3444 solveMultiVector<impl_type, internal_vector_scratch_type_3d_view>(member, blocksize, i0, r0, nrows, v, internal_vector_values, Kokkos::subview(e_internal_vector_values, 1, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()), WW);
3449 KOKKOS_INLINE_FUNCTION
3451 operator()(
const ComputeSchurTag &,
const member_type &member)
const {
3453 const local_ordinal_type packindices_schur_i = member.league_rank() % packindices_schur.extent(0);
3454 const local_ordinal_type packindices_schur_j = member.league_rank() / packindices_schur.extent(0);
3455 const local_ordinal_type packidx = packindices_schur(packindices_schur_i, packindices_schur_j);
3457 const local_ordinal_type subpartidx = packptr_sub(packidx);
3458 const local_ordinal_type n_parts = part2packrowidx0_sub.extent(0);
3459 const local_ordinal_type local_subpartidx = subpartidx / n_parts;
3460 const local_ordinal_type partidx = subpartidx % n_parts;
3463 const local_ordinal_type i0 = pack_td_ptr(partidx, local_subpartidx);
3469 const local_ordinal_type local_subpartidx_schur = (local_subpartidx - 1) / 2;
3470 const local_ordinal_type i0_schur = local_subpartidx_schur == 0 ? pack_td_ptr_schur(partidx, local_subpartidx_schur) : pack_td_ptr_schur(partidx, local_subpartidx_schur) + 1;
3471 const local_ordinal_type i0_offset = local_subpartidx_schur == 0 ? i0 + 2 : i0 + 2;
3473 for (local_ordinal_type i = 0; i < 4; ++i) {
3474 copy3DView<local_ordinal_type>(member, Kokkos::subview(internal_vector_values_schur, i0_schur + i, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()),
3475 Kokkos::subview(internal_vector_values, i0_offset + i, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()));
3478 member.team_barrier();
3480 const auto one = Kokkos::ArithTraits<btdm_magnitude_type>::one();
3482 const size_type c_kps1 = pack_td_ptr(partidx, local_subpartidx) + 1;
3483 const size_type c_kps2 = pack_td_ptr(partidx, local_subpartidx + 1) - 2;
3485 const local_ordinal_type e_r1 = part2packrowidx0_sub(partidx, local_subpartidx) - 1;
3486 const local_ordinal_type e_r2 = part2packrowidx0_sub(partidx, local_subpartidx) + 2;
3488 typedef ExtractAndFactorizeTridiagsDefaultModeAndAlgo<typename execution_space::memory_space> default_mode_and_algo_type;
3490 typedef typename default_mode_and_algo_type::mode_type default_mode_type;
3491 typedef typename default_mode_and_algo_type::algo_type default_algo_type;
3493 Kokkos::parallel_for(Kokkos::ThreadVectorRange(member, vector_loop_size), [&](
const int &v) {
3494 for (size_type i = 0; i < pack_td_ptr_schur(partidx, local_subpartidx_schur + 1) - pack_td_ptr_schur(partidx, local_subpartidx_schur); ++i) {
3495 local_ordinal_type e_r, e_c, c_kps;
3497 if (local_subpartidx_schur == 0) {
3502 }
else if (i == 3) {
3506 }
else if (i == 4) {
3518 }
else if (i == 1) {
3522 }
else if (i == 4) {
3526 }
else if (i == 5) {
3535 auto S = Kokkos::subview(internal_vector_values_schur, pack_td_ptr_schur(partidx, local_subpartidx_schur) + i, Kokkos::ALL(), Kokkos::ALL(), v);
3536 auto C = Kokkos::subview(internal_vector_values, c_kps, Kokkos::ALL(), Kokkos::ALL(), v);
3537 auto E = Kokkos::subview(e_internal_vector_values, e_c, e_r, Kokkos::ALL(), Kokkos::ALL(), v);
3538 KB::Gemm<member_type,
3539 KB::Trans::NoTranspose, KB::Trans::NoTranspose,
3540 default_mode_type, default_algo_type>::invoke(member, -one, C, E, one, S);
3545 KOKKOS_INLINE_FUNCTION
3547 operator()(
const FactorizeSchurTag &,
const member_type &member)
const {
3548 const local_ordinal_type packidx = packindices_schur(member.league_rank(), 0);
3550 const local_ordinal_type subpartidx = packptr_sub(packidx);
3552 const local_ordinal_type n_parts = part2packrowidx0_sub.extent(0);
3553 const local_ordinal_type partidx = subpartidx % n_parts;
3555 const local_ordinal_type i0 = pack_td_ptr_schur(partidx, 0);
3556 const local_ordinal_type nrows = 2 * (pack_td_ptr_schur.extent(1) - 1);
3558 internal_vector_scratch_type_3d_view
3559 WW(member.team_scratch(ScratchLevel), blocksize, blocksize, vector_loop_size);
3561 #ifdef IFPACK2_BLOCKTRIDICONTAINER_USE_PRINTF
3562 printf(
"FactorizeSchurTag rank = %d, i0 = %d, nrows = %d, vector_loop_size = %d;\n", member.league_rank(), i0, nrows, vector_loop_size);
3565 if (vector_loop_size == 1) {
3566 factorize_subline(member, i0, nrows, 0, internal_vector_values_schur, WW);
3568 Kokkos::parallel_for(Kokkos::ThreadVectorRange(member, vector_loop_size),
3569 [&](
const local_ordinal_type &v) {
3570 factorize_subline(member, i0, nrows, v, internal_vector_values_schur, WW);
3576 IFPACK2_BLOCKTRIDICONTAINER_PROFILER_REGION_BEGIN;
3577 const local_ordinal_type team_size =
3578 ExtractAndFactorizeTridiagsDefaultModeAndAlgo<typename execution_space::memory_space>::
3579 recommended_team_size(blocksize, vector_length, internal_vector_length);
3580 const local_ordinal_type per_team_scratch = internal_vector_scratch_type_3d_view::
3581 shmem_size(blocksize, blocksize, vector_loop_size);
3584 #ifdef IFPACK2_BLOCKTRIDICONTAINER_USE_PRINTF
3585 printf(
"Start ExtractAndFactorizeSubLineTag\n");
3587 IFPACK2_BLOCKHELPER_TIMER(
"BlockTriDi::NumericPhase::ExtractAndFactorizeSubLineTag", ExtractAndFactorizeSubLineTag0);
3588 Kokkos::TeamPolicy<execution_space, ExtractAndFactorizeSubLineTag>
3589 policy(packindices_sub.extent(0), team_size, vector_loop_size);
3591 const local_ordinal_type n_parts = part2packrowidx0_sub.extent(0);
3592 writeBTDValuesToFile(n_parts, scalar_values,
"before.mm");
3594 policy.set_scratch_size(ScratchLevel, Kokkos::PerTeam(per_team_scratch));
3595 Kokkos::parallel_for(
"ExtractAndFactorize::TeamPolicy::run<ExtractAndFactorizeSubLineTag>",
3597 execution_space().fence();
3599 writeBTDValuesToFile(n_parts, scalar_values,
"after.mm");
3600 #ifdef IFPACK2_BLOCKTRIDICONTAINER_USE_PRINTF
3601 printf(
"End ExtractAndFactorizeSubLineTag\n");
3605 if (packindices_schur.extent(1) > 0) {
3607 #ifdef IFPACK2_BLOCKTRIDICONTAINER_USE_PRINTF
3608 printf(
"Start ExtractBCDTag\n");
3610 Kokkos::deep_copy(e_scalar_values, Kokkos::ArithTraits<btdm_magnitude_type>::zero());
3611 Kokkos::deep_copy(scalar_values_schur, Kokkos::ArithTraits<btdm_magnitude_type>::zero());
3613 write5DMultiVectorValuesToFile(part2packrowidx0_sub.extent(0), e_scalar_values,
"e_scalar_values_before_extract.mm");
3616 IFPACK2_BLOCKHELPER_TIMER(
"BlockTriDi::NumericPhase::ExtractBCDTag", ExtractBCDTag0);
3617 Kokkos::TeamPolicy<execution_space, ExtractBCDTag>
3618 policy(packindices_schur.extent(0) * packindices_schur.extent(1), team_size, vector_loop_size);
3620 policy.set_scratch_size(ScratchLevel, Kokkos::PerTeam(per_team_scratch));
3621 Kokkos::parallel_for(
"ExtractAndFactorize::TeamPolicy::run<ExtractBCDTag>",
3623 execution_space().fence();
3626 #ifdef IFPACK2_BLOCKTRIDICONTAINER_USE_PRINTF
3627 printf(
"End ExtractBCDTag\n");
3629 writeBTDValuesToFile(part2packrowidx0_sub.extent(0), scalar_values,
"after_extraction_of_BCD.mm");
3630 #ifdef IFPACK2_BLOCKTRIDICONTAINER_USE_PRINTF
3631 printf(
"Start ComputeETag\n");
3633 write5DMultiVectorValuesToFile(part2packrowidx0_sub.extent(0), e_scalar_values,
"e_scalar_values_after_extract.mm");
3635 IFPACK2_BLOCKHELPER_TIMER(
"BlockTriDi::NumericPhase::ComputeETag", ComputeETag0);
3636 Kokkos::TeamPolicy<execution_space, ComputeETag>
3637 policy(packindices_sub.extent(0), team_size, vector_loop_size);
3639 policy.set_scratch_size(ScratchLevel, Kokkos::PerTeam(per_team_scratch));
3640 Kokkos::parallel_for(
"ExtractAndFactorize::TeamPolicy::run<ComputeETag>",
3642 execution_space().fence();
3644 write5DMultiVectorValuesToFile(part2packrowidx0_sub.extent(0), e_scalar_values,
"e_scalar_values_after_compute.mm");
3646 #ifdef IFPACK2_BLOCKTRIDICONTAINER_USE_PRINTF
3647 printf(
"End ComputeETag\n");
3652 #ifdef IFPACK2_BLOCKTRIDICONTAINER_USE_PRINTF
3653 printf(
"Start ComputeSchurTag\n");
3655 IFPACK2_BLOCKHELPER_TIMER(
"BlockTriDi::NumericPhase::ComputeSchurTag", ComputeSchurTag0);
3656 writeBTDValuesToFile(part2packrowidx0_sub.extent(0), scalar_values_schur,
"before_schur.mm");
3657 Kokkos::TeamPolicy<execution_space, ComputeSchurTag>
3658 policy(packindices_schur.extent(0) * packindices_schur.extent(1), team_size, vector_loop_size);
3660 Kokkos::parallel_for(
"ExtractAndFactorize::TeamPolicy::run<ComputeSchurTag>",
3662 writeBTDValuesToFile(part2packrowidx0_sub.extent(0), scalar_values_schur,
"after_schur.mm");
3663 execution_space().fence();
3664 #ifdef IFPACK2_BLOCKTRIDICONTAINER_USE_PRINTF
3665 printf(
"End ComputeSchurTag\n");
3670 #ifdef IFPACK2_BLOCKTRIDICONTAINER_USE_PRINTF
3671 printf(
"Start FactorizeSchurTag\n");
3673 IFPACK2_BLOCKHELPER_TIMER(
"BlockTriDi::NumericPhase::FactorizeSchurTag", FactorizeSchurTag0);
3674 Kokkos::TeamPolicy<execution_space, FactorizeSchurTag>
3675 policy(packindices_schur.extent(0), team_size, vector_loop_size);
3676 policy.set_scratch_size(ScratchLevel, Kokkos::PerTeam(per_team_scratch));
3677 Kokkos::parallel_for(
"ExtractAndFactorize::TeamPolicy::run<FactorizeSchurTag>",
3679 execution_space().fence();
3680 writeBTDValuesToFile(part2packrowidx0_sub.extent(0), scalar_values_schur,
"after_factor_schur.mm");
3681 #ifdef IFPACK2_BLOCKTRIDICONTAINER_USE_PRINTF
3682 printf(
"End FactorizeSchurTag\n");
3687 IFPACK2_BLOCKTRIDICONTAINER_PROFILER_REGION_END;
3690 void run_fused_jacobi() {
3691 IFPACK2_BLOCKTRIDICONTAINER_PROFILER_REGION_BEGIN;
3692 const local_ordinal_type team_size =
3693 ExtractAndFactorizeTridiagsDefaultModeAndAlgo<typename execution_space::memory_space>::
3694 recommended_team_size(blocksize, half_vector_length, 1);
3695 const local_ordinal_type per_team_scratch =
3696 btdm_scalar_scratch_type_3d_view::shmem_size(blocksize, blocksize, 2 * half_vector_length);
3698 IFPACK2_BLOCKHELPER_TIMER(
"BlockTriDi::NumericPhase::ExtractAndFactorizeFusedJacobi", ExtractAndFactorizeFusedJacobiTag);
3699 Kokkos::TeamPolicy<execution_space, ExtractAndFactorizeFusedJacobiTag>
3700 policy((lclrow.extent(0) + half_vector_length - 1) / half_vector_length, team_size, half_vector_length);
3702 policy.set_scratch_size(ScratchLevel, Kokkos::PerTeam(per_team_scratch));
3703 Kokkos::parallel_for(
"ExtractAndFactorize::TeamPolicy::run<ExtractAndFactorizeFusedJacobiTag>",
3706 IFPACK2_BLOCKTRIDICONTAINER_PROFILER_REGION_END;
3713 template <
typename MatrixType>
3715 const Teuchos::RCP<
const typename BlockHelperDetails::ImplType<MatrixType>::tpetra_crs_graph_type> &G,
3716 const BlockHelperDetails::PartInterface<MatrixType> &interf,
3718 const typename BlockHelperDetails::ImplType<MatrixType>::magnitude_type tiny,
3719 bool use_fused_jacobi) {
3721 using execution_space =
typename impl_type::execution_space;
3722 using team_policy_type = Kokkos::TeamPolicy<execution_space>;
3723 using internal_vector_scratch_type_3d_view = Scratch<typename impl_type::internal_vector_type_3d_view>;
3724 using btdm_scalar_scratch_type_3d_view = Scratch<typename impl_type::btdm_scalar_type_3d_view>;
3726 IFPACK2_BLOCKHELPER_TIMER(
"BlockTriDi::NumericPhase", NumericPhase);
3728 int blocksize = btdm.values.extent(1);
3731 int scratch_required;
3732 if (!use_fused_jacobi) {
3734 scratch_required = internal_vector_scratch_type_3d_view::shmem_size(blocksize, blocksize, impl_type::vector_length / impl_type::internal_vector_length);
3737 scratch_required = btdm_scalar_scratch_type_3d_view::shmem_size(blocksize, blocksize, 2 * impl_type::half_vector_length);
3740 int max_scratch = team_policy_type::scratch_size_max(0);
3742 if (scratch_required < max_scratch) {
3744 ExtractAndFactorizeTridiags<MatrixType, 0>
function(btdm, interf, A, G, tiny);
3745 if (!use_fused_jacobi)
3748 function.run_fused_jacobi();
3751 ExtractAndFactorizeTridiags<MatrixType, 1>
function(btdm, interf, A, G, tiny);
3752 if (!use_fused_jacobi)
3755 function.run_fused_jacobi();
3757 IFPACK2_BLOCKHELPER_TIMER_FENCE(
typename BlockHelperDetails::ImplType<MatrixType>::execution_space)
3763 template <
typename MatrixType>
3767 using execution_space =
typename impl_type::execution_space;
3768 using memory_space =
typename impl_type::memory_space;
3770 using local_ordinal_type =
typename impl_type::local_ordinal_type;
3772 using btdm_scalar_type =
typename impl_type::btdm_scalar_type;
3773 using tpetra_multivector_type =
typename impl_type::tpetra_multivector_type;
3774 using local_ordinal_type_1d_view =
typename impl_type::local_ordinal_type_1d_view;
3775 using vector_type_3d_view =
typename impl_type::vector_type_3d_view;
3776 using impl_scalar_type_2d_view_tpetra =
typename impl_type::impl_scalar_type_2d_view_tpetra;
3777 using const_impl_scalar_type_2d_view_tpetra =
typename impl_scalar_type_2d_view_tpetra::const_type;
3778 static constexpr
int vector_length = impl_type::vector_length;
3780 using member_type =
typename Kokkos::TeamPolicy<execution_space>::member_type;
3784 const ConstUnmanaged<local_ordinal_type_1d_view> partptr;
3785 const ConstUnmanaged<local_ordinal_type_1d_view> packptr;
3786 const ConstUnmanaged<local_ordinal_type_1d_view> part2packrowidx0;
3787 const ConstUnmanaged<local_ordinal_type_1d_view> part2rowidx0;
3788 const ConstUnmanaged<local_ordinal_type_1d_view> lclrow;
3789 const local_ordinal_type blocksize;
3790 const local_ordinal_type num_vectors;
3793 vector_type_3d_view packed_multivector;
3794 const_impl_scalar_type_2d_view_tpetra scalar_multivector;
3796 template <
typename TagType>
3797 KOKKOS_INLINE_FUNCTION
void copy_multivectors(
const local_ordinal_type &j,
3798 const local_ordinal_type &vi,
3799 const local_ordinal_type &pri,
3800 const local_ordinal_type &ri0)
const {
3801 for (local_ordinal_type col = 0; col < num_vectors; ++col)
3802 for (local_ordinal_type i = 0; i < blocksize; ++i)
3803 packed_multivector(pri, i, col)[vi] =
static_cast<btdm_scalar_type
>(scalar_multivector(blocksize * lclrow(ri0 + j) + i, col));
3808 const vector_type_3d_view &pmv)
3809 : partptr(interf.partptr)
3810 , packptr(interf.packptr)
3811 , part2packrowidx0(interf.part2packrowidx0)
3812 , part2rowidx0(interf.part2rowidx0)
3813 , lclrow(interf.lclrow)
3814 , blocksize(pmv.extent(1))
3815 , num_vectors(pmv.extent(2))
3816 , packed_multivector(pmv) {}
3819 KOKKOS_INLINE_FUNCTION
3821 operator()(
const local_ordinal_type &packidx)
const {
3822 local_ordinal_type partidx = packptr(packidx);
3823 local_ordinal_type npacks = packptr(packidx + 1) - partidx;
3824 const local_ordinal_type pri0 = part2packrowidx0(partidx);
3826 local_ordinal_type ri0[vector_length] = {};
3827 local_ordinal_type nrows[vector_length] = {};
3828 for (local_ordinal_type v = 0; v < npacks; ++v, ++partidx) {
3829 ri0[v] = part2rowidx0(partidx);
3830 nrows[v] = part2rowidx0(partidx + 1) - ri0[v];
3832 for (local_ordinal_type j = 0; j < nrows[0]; ++j) {
3833 local_ordinal_type cnt = 1;
3834 for (; cnt < npacks && j != nrows[cnt]; ++cnt)
3837 const local_ordinal_type pri = pri0 + j;
3838 for (local_ordinal_type col = 0; col < num_vectors; ++col)
3839 for (local_ordinal_type i = 0; i < blocksize; ++i)
3840 for (local_ordinal_type v = 0; v < npacks; ++v)
3841 packed_multivector(pri, i, col)[v] =
static_cast<btdm_scalar_type
>(scalar_multivector(blocksize * lclrow(ri0[v] + j) + i, col));
3845 KOKKOS_INLINE_FUNCTION
3847 operator()(
const member_type &member)
const {
3848 const local_ordinal_type packidx = member.league_rank();
3849 const local_ordinal_type partidx_begin = packptr(packidx);
3850 const local_ordinal_type npacks = packptr(packidx + 1) - partidx_begin;
3851 const local_ordinal_type pri0 = part2packrowidx0(partidx_begin);
3852 Kokkos::parallel_for(Kokkos::ThreadVectorRange(member, npacks), [&](
const local_ordinal_type &v) {
3853 const local_ordinal_type partidx = partidx_begin + v;
3854 const local_ordinal_type ri0 = part2rowidx0(partidx);
3855 const local_ordinal_type nrows = part2rowidx0(partidx + 1) - ri0;
3858 const local_ordinal_type pri = pri0;
3859 for (local_ordinal_type col = 0; col < num_vectors; ++col) {
3860 Kokkos::parallel_for(Kokkos::TeamThreadRange(member, blocksize), [&](
const local_ordinal_type &i) {
3861 packed_multivector(pri, i, col)[v] =
static_cast<btdm_scalar_type
>(scalar_multivector(blocksize * lclrow(ri0) + i, col));
3865 Kokkos::parallel_for(Kokkos::TeamThreadRange(member, nrows), [&](
const local_ordinal_type &j) {
3866 const local_ordinal_type pri = pri0 + j;
3867 for (local_ordinal_type col = 0; col < num_vectors; ++col)
3868 for (local_ordinal_type i = 0; i < blocksize; ++i)
3869 packed_multivector(pri, i, col)[v] = static_cast<btdm_scalar_type>(scalar_multivector(blocksize * lclrow(ri0 + j) + i, col));
3875 void run(
const const_impl_scalar_type_2d_view_tpetra &scalar_multivector_) {
3876 IFPACK2_BLOCKTRIDICONTAINER_PROFILER_REGION_BEGIN;
3877 IFPACK2_BLOCKHELPER_TIMER(
"BlockTriDi::MultiVectorConverter", MultiVectorConverter0);
3879 scalar_multivector = scalar_multivector_;
3880 if constexpr (BlockHelperDetails::is_device<execution_space>::value) {
3881 const local_ordinal_type vl = vector_length;
3882 const Kokkos::TeamPolicy<execution_space> policy(packptr.extent(0) - 1, Kokkos::AUTO(), vl);
3883 Kokkos::parallel_for(
"MultiVectorConverter::TeamPolicy", policy, *
this);
3885 const Kokkos::RangePolicy<execution_space> policy(0, packptr.extent(0) - 1);
3886 Kokkos::parallel_for(
"MultiVectorConverter::RangePolicy", policy, *
this);
3888 IFPACK2_BLOCKTRIDICONTAINER_PROFILER_REGION_END;
3889 IFPACK2_BLOCKHELPER_TIMER_FENCE(execution_space)
3898 struct SolveTridiagsDefaultModeAndAlgo<Kokkos::HostSpace> {
3899 typedef KB::Mode::Serial mode_type;
3900 typedef KB::Algo::Level2::Unblocked single_vector_algo_type;
3901 #if defined(__KOKKOSBATCHED_INTEL_MKL_COMPACT_BATCHED__)
3902 typedef KB::Algo::Level3::CompactMKL multi_vector_algo_type;
3904 typedef KB::Algo::Level3::Blocked multi_vector_algo_type;
3906 static int recommended_team_size(
const int ,
3913 #if defined(KOKKOS_ENABLE_CUDA)
3914 static inline int SolveTridiagsRecommendedCudaTeamSize(
const int blksize,
3915 const int vector_length,
3916 const int internal_vector_length) {
3917 const int vector_size = vector_length / internal_vector_length;
3918 int total_team_size(0);
3920 total_team_size = 32;
3921 else if (blksize <= 9)
3922 total_team_size = 32;
3923 else if (blksize <= 12)
3924 total_team_size = 96;
3925 else if (blksize <= 16)
3926 total_team_size = 128;
3927 else if (blksize <= 20)
3928 total_team_size = 160;
3930 total_team_size = 160;
3931 return total_team_size / vector_size;
3935 struct SolveTridiagsDefaultModeAndAlgo<Kokkos::CudaSpace> {
3936 typedef KB::Mode::Team mode_type;
3937 typedef KB::Algo::Level2::Unblocked single_vector_algo_type;
3938 typedef KB::Algo::Level3::Unblocked multi_vector_algo_type;
3939 static int recommended_team_size(
const int blksize,
3940 const int vector_length,
3941 const int internal_vector_length) {
3942 return SolveTridiagsRecommendedCudaTeamSize(blksize, vector_length, internal_vector_length);
3946 struct SolveTridiagsDefaultModeAndAlgo<Kokkos::CudaUVMSpace> {
3947 typedef KB::Mode::Team mode_type;
3948 typedef KB::Algo::Level2::Unblocked single_vector_algo_type;
3949 typedef KB::Algo::Level3::Unblocked multi_vector_algo_type;
3950 static int recommended_team_size(
const int blksize,
3951 const int vector_length,
3952 const int internal_vector_length) {
3953 return SolveTridiagsRecommendedCudaTeamSize(blksize, vector_length, internal_vector_length);
3958 #if defined(KOKKOS_ENABLE_HIP)
3959 static inline int SolveTridiagsRecommendedHIPTeamSize(
const int blksize,
3960 const int vector_length,
3961 const int internal_vector_length) {
3962 const int vector_size = vector_length / internal_vector_length;
3963 int total_team_size(0);
3965 total_team_size = 32;
3966 else if (blksize <= 9)
3967 total_team_size = 32;
3968 else if (blksize <= 12)
3969 total_team_size = 96;
3970 else if (blksize <= 16)
3971 total_team_size = 128;
3972 else if (blksize <= 20)
3973 total_team_size = 160;
3975 total_team_size = 160;
3976 return total_team_size / vector_size;
3980 struct SolveTridiagsDefaultModeAndAlgo<Kokkos::HIPSpace> {
3981 typedef KB::Mode::Team mode_type;
3982 typedef KB::Algo::Level2::Unblocked single_vector_algo_type;
3983 typedef KB::Algo::Level3::Unblocked multi_vector_algo_type;
3984 static int recommended_team_size(
const int blksize,
3985 const int vector_length,
3986 const int internal_vector_length) {
3987 return SolveTridiagsRecommendedHIPTeamSize(blksize, vector_length, internal_vector_length);
3991 struct SolveTridiagsDefaultModeAndAlgo<Kokkos::HIPHostPinnedSpace> {
3992 typedef KB::Mode::Team mode_type;
3993 typedef KB::Algo::Level2::Unblocked single_vector_algo_type;
3994 typedef KB::Algo::Level3::Unblocked multi_vector_algo_type;
3995 static int recommended_team_size(
const int blksize,
3996 const int vector_length,
3997 const int internal_vector_length) {
3998 return SolveTridiagsRecommendedHIPTeamSize(blksize, vector_length, internal_vector_length);
4003 #if defined(KOKKOS_ENABLE_SYCL)
4004 static inline int SolveTridiagsRecommendedSYCLTeamSize(
const int blksize,
4005 const int vector_length,
4006 const int internal_vector_length) {
4007 const int vector_size = vector_length / internal_vector_length;
4008 int total_team_size(0);
4010 total_team_size = 32;
4011 else if (blksize <= 9)
4012 total_team_size = 32;
4013 else if (blksize <= 12)
4014 total_team_size = 96;
4015 else if (blksize <= 16)
4016 total_team_size = 128;
4017 else if (blksize <= 20)
4018 total_team_size = 160;
4020 total_team_size = 160;
4021 return total_team_size / vector_size;
4025 struct SolveTridiagsDefaultModeAndAlgo<Kokkos::Experimental::SYCLSharedUSMSpace> {
4026 typedef KB::Mode::Team mode_type;
4027 typedef KB::Algo::Level2::Unblocked single_vector_algo_type;
4028 typedef KB::Algo::Level3::Unblocked multi_vector_algo_type;
4029 static int recommended_team_size(
const int blksize,
4030 const int vector_length,
4031 const int internal_vector_length) {
4032 return SolveTridiagsRecommendedSYCLTeamSize(blksize, vector_length, internal_vector_length);
4036 struct SolveTridiagsDefaultModeAndAlgo<Kokkos::Experimental::SYCLDeviceUSMSpace> {
4037 typedef KB::Mode::Team mode_type;
4038 typedef KB::Algo::Level2::Unblocked single_vector_algo_type;
4039 typedef KB::Algo::Level3::Unblocked multi_vector_algo_type;
4040 static int recommended_team_size(
const int blksize,
4041 const int vector_length,
4042 const int internal_vector_length) {
4043 return SolveTridiagsRecommendedSYCLTeamSize(blksize, vector_length, internal_vector_length);
4048 template <
typename MatrixType>
4049 struct SolveTridiags {
4051 using impl_type = BlockHelperDetails::ImplType<MatrixType>;
4052 using execution_space =
typename impl_type::execution_space;
4054 using local_ordinal_type =
typename impl_type::local_ordinal_type;
4057 using magnitude_type =
typename impl_type::magnitude_type;
4058 using btdm_scalar_type =
typename impl_type::btdm_scalar_type;
4059 using btdm_magnitude_type =
typename impl_type::btdm_magnitude_type;
4061 using local_ordinal_type_1d_view =
typename impl_type::local_ordinal_type_1d_view;
4062 using local_ordinal_type_2d_view =
typename impl_type::local_ordinal_type_2d_view;
4063 using size_type_2d_view =
typename impl_type::size_type_2d_view;
4065 using vector_type_3d_view =
typename impl_type::vector_type_3d_view;
4066 using internal_vector_type_4d_view =
typename impl_type::internal_vector_type_4d_view;
4067 using internal_vector_type_5d_view =
typename impl_type::internal_vector_type_5d_view;
4068 using btdm_scalar_type_4d_view =
typename impl_type::btdm_scalar_type_4d_view;
4070 using internal_vector_scratch_type_3d_view = Scratch<typename impl_type::internal_vector_type_3d_view>;
4072 using internal_vector_type =
typename impl_type::internal_vector_type;
4073 static constexpr
int vector_length = impl_type::vector_length;
4074 static constexpr
int internal_vector_length = impl_type::internal_vector_length;
4077 using impl_scalar_type_1d_view =
typename impl_type::impl_scalar_type_1d_view;
4078 using impl_scalar_type_2d_view_tpetra =
typename impl_type::impl_scalar_type_2d_view_tpetra;
4081 using team_policy_type = Kokkos::TeamPolicy<execution_space>;
4082 using member_type =
typename team_policy_type::member_type;
4086 local_ordinal_type n_subparts_per_part;
4087 const ConstUnmanaged<local_ordinal_type_1d_view> partptr;
4088 const ConstUnmanaged<local_ordinal_type_1d_view> packptr;
4089 const ConstUnmanaged<local_ordinal_type_1d_view> packindices_sub;
4090 const ConstUnmanaged<local_ordinal_type_2d_view> packindices_schur;
4091 const ConstUnmanaged<local_ordinal_type_1d_view> part2packrowidx0;
4092 const ConstUnmanaged<local_ordinal_type_2d_view> part2packrowidx0_sub;
4093 const ConstUnmanaged<local_ordinal_type_1d_view> lclrow;
4094 const ConstUnmanaged<local_ordinal_type_1d_view> packptr_sub;
4096 const ConstUnmanaged<local_ordinal_type_2d_view> partptr_sub;
4097 const ConstUnmanaged<size_type_2d_view> pack_td_ptr_schur;
4100 const ConstUnmanaged<size_type_2d_view> pack_td_ptr;
4103 const ConstUnmanaged<internal_vector_type_4d_view> D_internal_vector_values;
4104 const Unmanaged<internal_vector_type_4d_view> X_internal_vector_values;
4105 const Unmanaged<btdm_scalar_type_4d_view> X_internal_scalar_values;
4107 internal_vector_type_4d_view X_internal_vector_values_schur;
4109 const ConstUnmanaged<internal_vector_type_4d_view> D_internal_vector_values_schur;
4110 const ConstUnmanaged<internal_vector_type_5d_view> e_internal_vector_values;
4112 const local_ordinal_type vector_loop_size;
4115 Unmanaged<impl_scalar_type_2d_view_tpetra> Y_scalar_multivector;
4116 #if defined(__CUDA_ARCH__) || defined(__HIP_DEVICE_COMPILE__) || defined(__SYCL_DEVICE_ONLY__)
4117 AtomicUnmanaged<impl_scalar_type_1d_view> Z_scalar_vector;
4119 Unmanaged<impl_scalar_type_1d_view> Z_scalar_vector;
4121 const impl_scalar_type df;
4122 const bool compute_diff;
4125 SolveTridiags(
const BlockHelperDetails::PartInterface<MatrixType> &interf,
4126 const BlockTridiags<MatrixType> &btdm,
4127 const vector_type_3d_view &pmv,
4128 const impl_scalar_type damping_factor,
4129 const bool is_norm_manager_active)
4131 n_subparts_per_part(interf.n_subparts_per_part)
4132 , partptr(interf.partptr)
4133 , packptr(interf.packptr)
4134 , packindices_sub(interf.packindices_sub)
4135 , packindices_schur(interf.packindices_schur)
4136 , part2packrowidx0(interf.part2packrowidx0)
4137 , part2packrowidx0_sub(interf.part2packrowidx0_sub)
4138 , lclrow(interf.lclrow)
4139 , packptr_sub(interf.packptr_sub)
4140 , partptr_sub(interf.partptr_sub)
4141 , pack_td_ptr_schur(btdm.pack_td_ptr_schur)
4144 pack_td_ptr(btdm.pack_td_ptr)
4145 , D_internal_vector_values((internal_vector_type *)btdm.values.data(),
4146 btdm.values.extent(0),
4147 btdm.values.extent(1),
4148 btdm.values.extent(2),
4149 vector_length / internal_vector_length)
4150 , X_internal_vector_values((internal_vector_type *)pmv.data(),
4154 vector_length / internal_vector_length)
4155 , X_internal_scalar_values((btdm_scalar_type *)pmv.data(),
4161 2 * (n_subparts_per_part - 1) * part2packrowidx0_sub.extent(0),
4164 vector_length / internal_vector_length)
4165 , D_internal_vector_values_schur((internal_vector_type *)btdm.values_schur.data(),
4166 btdm.values_schur.extent(0),
4167 btdm.values_schur.extent(1),
4168 btdm.values_schur.extent(2),
4169 vector_length / internal_vector_length)
4170 , e_internal_vector_values((internal_vector_type *)btdm.e_values.data(),
4171 btdm.e_values.extent(0),
4172 btdm.e_values.extent(1),
4173 btdm.e_values.extent(2),
4174 btdm.e_values.extent(3),
4175 vector_length / internal_vector_length)
4176 , vector_loop_size(vector_length / internal_vector_length)
4177 , Y_scalar_multivector()
4179 , df(damping_factor)
4180 , compute_diff(is_norm_manager_active) {}
4184 KOKKOS_INLINE_FUNCTION
4186 copyToFlatMultiVector(
const member_type &member,
4187 const local_ordinal_type partidxbeg,
4188 const local_ordinal_type npacks,
4189 const local_ordinal_type pri0,
4190 const local_ordinal_type v,
4191 const local_ordinal_type blocksize,
4192 const local_ordinal_type num_vectors)
const {
4193 const local_ordinal_type vbeg = v * internal_vector_length;
4194 if (vbeg < npacks) {
4195 local_ordinal_type ri0_vals[internal_vector_length] = {};
4196 local_ordinal_type nrows_vals[internal_vector_length] = {};
4197 for (local_ordinal_type vv = vbeg, vi = 0; vv < npacks && vi < internal_vector_length; ++vv, ++vi) {
4198 const local_ordinal_type partidx = partidxbeg + vv;
4199 ri0_vals[vi] = partptr(partidx);
4200 nrows_vals[vi] = partptr(partidx + 1) - ri0_vals[vi];
4203 impl_scalar_type z_partial_sum(0);
4204 if (nrows_vals[0] == 1) {
4205 const local_ordinal_type j = 0, pri = pri0;
4207 for (local_ordinal_type vv = vbeg, vi = 0; vv < npacks && vi < internal_vector_length; ++vv, ++vi) {
4208 const local_ordinal_type ri0 = ri0_vals[vi];
4209 const local_ordinal_type nrows = nrows_vals[vi];
4211 Kokkos::parallel_for(Kokkos::TeamThreadRange(member, blocksize),
4212 [&](
const local_ordinal_type &i) {
4213 const local_ordinal_type row = blocksize * lclrow(ri0 + j) + i;
4214 for (local_ordinal_type col = 0; col < num_vectors; ++col) {
4215 impl_scalar_type &y = Y_scalar_multivector(row, col);
4216 const impl_scalar_type yd = X_internal_vector_values(pri, i, col, v)[vi] - y;
4220 const auto yd_abs = Kokkos::ArithTraits<impl_scalar_type>::abs(yd);
4221 z_partial_sum += yd_abs * yd_abs;
4229 Kokkos::parallel_for(Kokkos::TeamThreadRange(member, nrows_vals[0]),
4230 [&](
const local_ordinal_type &j) {
4231 const local_ordinal_type pri = pri0 + j;
4232 for (local_ordinal_type vv = vbeg, vi = 0; vv < npacks && vi < internal_vector_length; ++vv, ++vi) {
4233 const local_ordinal_type ri0 = ri0_vals[vi];
4234 const local_ordinal_type nrows = nrows_vals[vi];
4236 for (local_ordinal_type col = 0; col < num_vectors; ++col) {
4237 for (local_ordinal_type i = 0; i < blocksize; ++i) {
4238 const local_ordinal_type row = blocksize * lclrow(ri0 + j) + i;
4239 impl_scalar_type &y = Y_scalar_multivector(row, col);
4240 const impl_scalar_type yd = X_internal_vector_values(pri, i, col, v)[vi] - y;
4244 const auto yd_abs = Kokkos::ArithTraits<impl_scalar_type>::abs(yd);
4245 z_partial_sum += yd_abs * yd_abs;
4254 Z_scalar_vector(member.league_rank()) += z_partial_sum;
4261 template <
typename WWViewType>
4262 KOKKOS_INLINE_FUNCTION
void
4263 solveSingleVector(
const member_type &member,
4264 const local_ordinal_type &blocksize,
4265 const local_ordinal_type &i0,
4266 const local_ordinal_type &r0,
4267 const local_ordinal_type &nrows,
4268 const local_ordinal_type &v,
4269 const WWViewType &WW)
const {
4270 typedef SolveTridiagsDefaultModeAndAlgo<typename execution_space::memory_space> default_mode_and_algo_type;
4272 typedef typename default_mode_and_algo_type::mode_type default_mode_type;
4273 typedef typename default_mode_and_algo_type::single_vector_algo_type default_algo_type;
4276 auto A = D_internal_vector_values.data();
4277 auto X = X_internal_vector_values.data();
4280 const auto one = Kokkos::ArithTraits<btdm_magnitude_type>::one();
4281 const auto zero = Kokkos::ArithTraits<btdm_magnitude_type>::zero();
4285 const local_ordinal_type astep = D_internal_vector_values.stride(0);
4286 const local_ordinal_type as0 = D_internal_vector_values.stride(1);
4287 const local_ordinal_type as1 = D_internal_vector_values.stride(2);
4288 const local_ordinal_type xstep = X_internal_vector_values.stride(0);
4289 const local_ordinal_type xs0 = X_internal_vector_values.stride(1);
4292 A += i0 * astep + v;
4293 X += r0 * xstep + v;
4298 KOKKOSBATCHED_TRSV_LOWER_NO_TRANSPOSE_INTERNAL_INVOKE(default_mode_type, default_algo_type,
4301 blocksize, blocksize,
4306 for (local_ordinal_type tr = 1; tr < nrows; ++tr) {
4307 member.team_barrier();
4308 KOKKOSBATCHED_GEMV_NO_TRANSPOSE_INTERNAL_INVOKE(default_mode_type, default_algo_type,
4310 blocksize, blocksize,
4312 A + 2 * astep, as0, as1,
4315 X + 1 * xstep, xs0);
4316 KOKKOSBATCHED_TRSV_LOWER_NO_TRANSPOSE_INTERNAL_INVOKE(default_mode_type, default_algo_type,
4319 blocksize, blocksize,
4321 A + 3 * astep, as0, as1,
4322 X + 1 * xstep, xs0);
4329 KOKKOSBATCHED_TRSV_UPPER_NO_TRANSPOSE_INTERNAL_INVOKE(default_mode_type, default_algo_type,
4332 blocksize, blocksize,
4337 for (local_ordinal_type tr = nrows; tr > 1; --tr) {
4339 member.team_barrier();
4340 KOKKOSBATCHED_GEMV_NO_TRANSPOSE_INTERNAL_INVOKE(default_mode_type, default_algo_type,
4342 blocksize, blocksize,
4344 A + 1 * astep, as0, as1,
4347 X - 1 * xstep, xs0);
4348 KOKKOSBATCHED_TRSV_UPPER_NO_TRANSPOSE_INTERNAL_INVOKE(default_mode_type, default_algo_type,
4351 blocksize, blocksize,
4354 X - 1 * xstep, xs0);
4360 const local_ordinal_type ws0 = WW.stride(0);
4361 auto W = WW.data() + v;
4362 KOKKOSBATCHED_COPY_VECTOR_NO_TRANSPOSE_INTERNAL_INVOKE(default_mode_type,
4363 member, blocksize, X, xs0, W, ws0);
4364 member.team_barrier();
4365 KOKKOSBATCHED_GEMV_NO_TRANSPOSE_INTERNAL_INVOKE(default_mode_type, default_algo_type,
4367 blocksize, blocksize,
4376 template <
typename WWViewType>
4377 KOKKOS_INLINE_FUNCTION
void
4378 solveMultiVector(
const member_type &member,
4379 const local_ordinal_type & ,
4380 const local_ordinal_type &i0,
4381 const local_ordinal_type &r0,
4382 const local_ordinal_type &nrows,
4383 const local_ordinal_type &v,
4384 const WWViewType &WW)
const {
4385 typedef SolveTridiagsDefaultModeAndAlgo<typename execution_space::memory_space> default_mode_and_algo_type;
4387 typedef typename default_mode_and_algo_type::mode_type default_mode_type;
4388 typedef typename default_mode_and_algo_type::multi_vector_algo_type default_algo_type;
4391 const auto one = Kokkos::ArithTraits<btdm_magnitude_type>::one();
4392 const auto zero = Kokkos::ArithTraits<btdm_magnitude_type>::zero();
4395 auto A = Kokkos::subview(D_internal_vector_values, i0, Kokkos::ALL(), Kokkos::ALL(), v);
4396 auto X1 = Kokkos::subview(X_internal_vector_values, r0, Kokkos::ALL(), Kokkos::ALL(), v);
4399 local_ordinal_type i = i0, r = r0;
4403 KB::Trsm<member_type,
4404 KB::Side::Left, KB::Uplo::Lower, KB::Trans::NoTranspose, KB::Diag::Unit,
4405 default_mode_type, default_algo_type>::invoke(member, one, A, X1);
4406 for (local_ordinal_type tr = 1; tr < nrows; ++tr, i += 3) {
4407 A.assign_data(&D_internal_vector_values(i + 2, 0, 0, v));
4408 X2.assign_data(&X_internal_vector_values(++r, 0, 0, v));
4409 member.team_barrier();
4410 KB::Gemm<member_type,
4411 KB::Trans::NoTranspose, KB::Trans::NoTranspose,
4412 default_mode_type, default_algo_type>::invoke(member, -one, A, X1, one, X2);
4413 A.assign_data(&D_internal_vector_values(i + 3, 0, 0, v));
4414 KB::Trsm<member_type,
4415 KB::Side::Left, KB::Uplo::Lower, KB::Trans::NoTranspose, KB::Diag::Unit,
4416 default_mode_type, default_algo_type>::invoke(member, one, A, X2);
4417 X1.assign_data(X2.data());
4421 KB::Trsm<member_type,
4422 KB::Side::Left, KB::Uplo::Upper, KB::Trans::NoTranspose, KB::Diag::NonUnit,
4423 default_mode_type, default_algo_type>::invoke(member, one, A, X1);
4424 for (local_ordinal_type tr = nrows; tr > 1; --tr) {
4426 A.assign_data(&D_internal_vector_values(i + 1, 0, 0, v));
4427 X2.assign_data(&X_internal_vector_values(--r, 0, 0, v));
4428 member.team_barrier();
4429 KB::Gemm<member_type,
4430 KB::Trans::NoTranspose, KB::Trans::NoTranspose,
4431 default_mode_type, default_algo_type>::invoke(member, -one, A, X1, one, X2);
4433 A.assign_data(&D_internal_vector_values(i, 0, 0, v));
4434 KB::Trsm<member_type,
4435 KB::Side::Left, KB::Uplo::Upper, KB::Trans::NoTranspose, KB::Diag::NonUnit,
4436 default_mode_type, default_algo_type>::invoke(member, one, A, X2);
4437 X1.assign_data(X2.data());
4441 auto W = Kokkos::subview(WW, Kokkos::ALL(), Kokkos::ALL(), v);
4442 KB::Copy<member_type, KB::Trans::NoTranspose, default_mode_type>::invoke(member, X1, W);
4443 member.team_barrier();
4444 KB::Gemm<member_type,
4445 KB::Trans::NoTranspose, KB::Trans::NoTranspose,
4446 default_mode_type, default_algo_type>::invoke(member, one, A, W, zero, X1);
4451 struct SingleVectorTag {};
4453 struct MultiVectorTag {};
4456 struct SingleVectorSubLineTag {};
4458 struct MultiVectorSubLineTag {};
4460 struct SingleVectorApplyCTag {};
4462 struct MultiVectorApplyCTag {};
4464 struct SingleVectorSchurTag {};
4466 struct MultiVectorSchurTag {};
4468 struct SingleVectorApplyETag {};
4470 struct MultiVectorApplyETag {};
4472 struct SingleVectorCopyToFlatTag {};
4474 struct SingleZeroingTag {};
4477 KOKKOS_INLINE_FUNCTION
void
4478 operator()(
const SingleVectorTag<B> &,
const member_type &member)
const {
4479 const local_ordinal_type packidx = member.league_rank();
4480 const local_ordinal_type partidx = packptr(packidx);
4481 const local_ordinal_type npacks = packptr(packidx + 1) - partidx;
4482 const local_ordinal_type pri0 = part2packrowidx0(partidx);
4483 const local_ordinal_type i0 = pack_td_ptr(partidx, 0);
4484 const local_ordinal_type r0 = part2packrowidx0(partidx);
4485 const local_ordinal_type nrows = partptr(partidx + 1) - partptr(partidx);
4486 const local_ordinal_type blocksize = (B == 0 ? D_internal_vector_values.extent(1) : B);
4487 const local_ordinal_type num_vectors = 1;
4488 internal_vector_scratch_type_3d_view
4489 WW(member.team_scratch(0), blocksize, 1, vector_loop_size);
4490 Kokkos::single(Kokkos::PerTeam(member), [&]() {
4491 Z_scalar_vector(member.league_rank()) = impl_scalar_type(0);
4493 Kokkos::parallel_for(Kokkos::ThreadVectorRange(member, vector_loop_size), [&](
const int &v) {
4494 solveSingleVector(member, blocksize, i0, r0, nrows, v, WW);
4495 copyToFlatMultiVector(member, partidx, npacks, pri0, v, blocksize, num_vectors);
4500 KOKKOS_INLINE_FUNCTION
void
4501 operator()(
const MultiVectorTag<B> &,
const member_type &member)
const {
4502 const local_ordinal_type packidx = member.league_rank();
4503 const local_ordinal_type partidx = packptr(packidx);
4504 const local_ordinal_type npacks = packptr(packidx + 1) - partidx;
4505 const local_ordinal_type pri0 = part2packrowidx0(partidx);
4506 const local_ordinal_type i0 = pack_td_ptr(partidx, 0);
4507 const local_ordinal_type r0 = part2packrowidx0(partidx);
4508 const local_ordinal_type nrows = partptr(partidx + 1) - partptr(partidx);
4509 const local_ordinal_type blocksize = (B == 0 ? D_internal_vector_values.extent(1) : B);
4510 const local_ordinal_type num_vectors = X_internal_vector_values.extent(2);
4512 internal_vector_scratch_type_3d_view
4513 WW(member.team_scratch(0), blocksize, num_vectors, vector_loop_size);
4514 Kokkos::single(Kokkos::PerTeam(member), [&]() {
4515 Z_scalar_vector(member.league_rank()) = impl_scalar_type(0);
4517 Kokkos::parallel_for(Kokkos::ThreadVectorRange(member, vector_loop_size), [&](
const int &v) {
4518 solveMultiVector(member, blocksize, i0, r0, nrows, v, WW);
4519 copyToFlatMultiVector(member, partidx, npacks, pri0, v, blocksize, num_vectors);
4524 KOKKOS_INLINE_FUNCTION
void
4525 operator()(
const SingleVectorSubLineTag<B> &,
const member_type &member)
const {
4527 const local_ordinal_type packidx = packindices_sub(member.league_rank());
4529 const local_ordinal_type subpartidx = packptr_sub(packidx);
4530 const local_ordinal_type n_parts = part2packrowidx0_sub.extent(0);
4531 const local_ordinal_type local_subpartidx = subpartidx / n_parts;
4532 const local_ordinal_type partidx = subpartidx % n_parts;
4534 const local_ordinal_type npacks = packptr_sub(packidx + 1) - subpartidx;
4535 const local_ordinal_type i0 = pack_td_ptr(partidx, local_subpartidx);
4536 const local_ordinal_type r0 = part2packrowidx0_sub(partidx, local_subpartidx);
4537 const local_ordinal_type nrows = partptr_sub(subpartidx, 1) - partptr_sub(subpartidx, 0);
4538 const local_ordinal_type blocksize = e_internal_vector_values.extent(2);
4544 internal_vector_scratch_type_3d_view
4545 WW(member.team_scratch(0), blocksize, 1, vector_loop_size);
4547 Kokkos::parallel_for(Kokkos::ThreadVectorRange(member, vector_loop_size), [&](
const int &v) {
4548 solveSingleVectorNew<impl_type, internal_vector_scratch_type_3d_view>(member, blocksize, i0, r0, nrows, v, D_internal_vector_values, X_internal_vector_values, WW);
4553 KOKKOS_INLINE_FUNCTION
void
4554 operator()(
const SingleVectorApplyCTag<B> &,
const member_type &member)
const {
4557 const local_ordinal_type packidx = packindices_sub(member.league_rank());
4559 const local_ordinal_type subpartidx = packptr_sub(packidx);
4560 const local_ordinal_type n_parts = part2packrowidx0_sub.extent(0);
4561 const local_ordinal_type local_subpartidx = subpartidx / n_parts;
4562 const local_ordinal_type partidx = subpartidx % n_parts;
4563 const local_ordinal_type blocksize = e_internal_vector_values.extent(2);
4566 const local_ordinal_type i0 = pack_td_ptr(partidx, local_subpartidx);
4567 const local_ordinal_type r0 = part2packrowidx0_sub(partidx, local_subpartidx);
4568 const local_ordinal_type nrows = partptr_sub(subpartidx, 1) - partptr_sub(subpartidx, 0);
4570 internal_vector_scratch_type_3d_view
4571 WW(member.team_scratch(0), blocksize, blocksize, vector_loop_size);
4575 const local_ordinal_type local_subpartidx_schur = (local_subpartidx - 1) / 2;
4576 const local_ordinal_type i0_schur = local_subpartidx_schur == 0 ? pack_td_ptr_schur(partidx, local_subpartidx_schur) : pack_td_ptr_schur(partidx, local_subpartidx_schur) + 1;
4577 const local_ordinal_type i0_offset = local_subpartidx_schur == 0 ? i0 + 2 : i0 + 2;
4582 const auto one = Kokkos::ArithTraits<btdm_magnitude_type>::one();
4584 const size_type c_kps2 = local_subpartidx > 0 ? pack_td_ptr(partidx, local_subpartidx) - 2 : 0;
4585 const size_type c_kps1 = pack_td_ptr(partidx, local_subpartidx + 1) + 1;
4587 typedef SolveTridiagsDefaultModeAndAlgo<typename execution_space::memory_space> default_mode_and_algo_type;
4589 typedef typename default_mode_and_algo_type::mode_type default_mode_type;
4590 typedef typename default_mode_and_algo_type::single_vector_algo_type default_algo_type;
4592 if (local_subpartidx == 0) {
4593 Kokkos::parallel_for(Kokkos::ThreadVectorRange(member, vector_loop_size), [&](
const int &v) {
4594 auto v_1 = Kokkos::subview(X_internal_vector_values, r0 + nrows - 1, Kokkos::ALL(), 0, v);
4595 auto v_2 = Kokkos::subview(X_internal_vector_values, r0 + nrows, Kokkos::ALL(), 0, v);
4596 auto C = Kokkos::subview(D_internal_vector_values, c_kps1, Kokkos::ALL(), Kokkos::ALL(), v);
4598 KOKKOSBATCHED_GEMV_NO_TRANSPOSE_INTERNAL_INVOKE(default_mode_type, default_algo_type,
4600 blocksize, blocksize,
4602 C.data(), C.stride(0), C.stride(1),
4603 v_1.data(), v_1.stride(0),
4605 v_2.data(), v_2.stride(0));
4607 }
else if (local_subpartidx == (local_ordinal_type)part2packrowidx0_sub.extent(1) - 2) {
4608 Kokkos::parallel_for(Kokkos::ThreadVectorRange(member, vector_loop_size), [&](
const int &v) {
4609 auto v_1 = Kokkos::subview(X_internal_vector_values, r0, Kokkos::ALL(), 0, v);
4610 auto v_2 = Kokkos::subview(X_internal_vector_values, r0 - 1, Kokkos::ALL(), 0, v);
4611 auto C = Kokkos::subview(D_internal_vector_values, c_kps2, Kokkos::ALL(), Kokkos::ALL(), v);
4613 KOKKOSBATCHED_GEMV_NO_TRANSPOSE_INTERNAL_INVOKE(default_mode_type, default_algo_type,
4615 blocksize, blocksize,
4617 C.data(), C.stride(0), C.stride(1),
4618 v_1.data(), v_1.stride(0),
4620 v_2.data(), v_2.stride(0));
4623 Kokkos::parallel_for(Kokkos::ThreadVectorRange(member, vector_loop_size), [&](
const int &v) {
4625 auto v_1 = Kokkos::subview(X_internal_vector_values, r0 + nrows - 1, Kokkos::ALL(), 0, v);
4626 auto v_2 = Kokkos::subview(X_internal_vector_values, r0 + nrows, Kokkos::ALL(), 0, v);
4627 auto C = Kokkos::subview(D_internal_vector_values, c_kps1, Kokkos::ALL(), Kokkos::ALL(), v);
4629 KOKKOSBATCHED_GEMV_NO_TRANSPOSE_INTERNAL_INVOKE(default_mode_type, default_algo_type,
4631 blocksize, blocksize,
4633 C.data(), C.stride(0), C.stride(1),
4634 v_1.data(), v_1.stride(0),
4636 v_2.data(), v_2.stride(0));
4639 auto v_1 = Kokkos::subview(X_internal_vector_values, r0, Kokkos::ALL(), 0, v);
4640 auto v_2 = Kokkos::subview(X_internal_vector_values, r0 - 1, Kokkos::ALL(), 0, v);
4641 auto C = Kokkos::subview(D_internal_vector_values, c_kps2, Kokkos::ALL(), Kokkos::ALL(), v);
4643 KOKKOSBATCHED_GEMV_NO_TRANSPOSE_INTERNAL_INVOKE(default_mode_type, default_algo_type,
4645 blocksize, blocksize,
4647 C.data(), C.stride(0), C.stride(1),
4648 v_1.data(), v_1.stride(0),
4650 v_2.data(), v_2.stride(0));
4657 KOKKOS_INLINE_FUNCTION
void
4658 operator()(
const SingleVectorSchurTag<B> &,
const member_type &member)
const {
4659 const local_ordinal_type packidx = packindices_sub(member.league_rank());
4661 const local_ordinal_type partidx = packptr_sub(packidx);
4663 const local_ordinal_type blocksize = e_internal_vector_values.extent(2);
4665 const local_ordinal_type i0_schur = pack_td_ptr_schur(partidx, 0);
4666 const local_ordinal_type nrows = 2 * (n_subparts_per_part - 1);
4668 const local_ordinal_type r0_schur = nrows * member.league_rank();
4670 internal_vector_scratch_type_3d_view
4671 WW(member.team_scratch(0), blocksize, blocksize, vector_loop_size);
4673 for (local_ordinal_type schur_sub_part = 0; schur_sub_part < n_subparts_per_part - 1; ++schur_sub_part) {
4674 const local_ordinal_type r0 = part2packrowidx0_sub(partidx, 2 * schur_sub_part + 1);
4675 for (local_ordinal_type i = 0; i < 2; ++i) {
4676 copy3DView<local_ordinal_type>(member,
4677 Kokkos::subview(X_internal_vector_values_schur, r0_schur + 2 * schur_sub_part + i, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()),
4678 Kokkos::subview(X_internal_vector_values, r0 + i, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()));
4682 Kokkos::parallel_for(Kokkos::ThreadVectorRange(member, vector_loop_size), [&](
const int &v) {
4683 solveSingleVectorNew<impl_type, internal_vector_scratch_type_3d_view>(member, blocksize, i0_schur, r0_schur, nrows, v, D_internal_vector_values_schur, X_internal_vector_values_schur, WW);
4686 for (local_ordinal_type schur_sub_part = 0; schur_sub_part < n_subparts_per_part - 1; ++schur_sub_part) {
4687 const local_ordinal_type r0 = part2packrowidx0_sub(partidx, 2 * schur_sub_part + 1);
4688 for (local_ordinal_type i = 0; i < 2; ++i) {
4689 copy3DView<local_ordinal_type>(member,
4690 Kokkos::subview(X_internal_vector_values, r0 + i, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()),
4691 Kokkos::subview(X_internal_vector_values_schur, r0_schur + 2 * schur_sub_part + i, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()));
4697 KOKKOS_INLINE_FUNCTION
void
4698 operator()(
const SingleVectorApplyETag<B> &,
const member_type &member)
const {
4699 const local_ordinal_type packidx = packindices_sub(member.league_rank());
4701 const local_ordinal_type subpartidx = packptr_sub(packidx);
4702 const local_ordinal_type n_parts = part2packrowidx0_sub.extent(0);
4703 const local_ordinal_type local_subpartidx = subpartidx / n_parts;
4704 const local_ordinal_type partidx = subpartidx % n_parts;
4705 const local_ordinal_type blocksize = e_internal_vector_values.extent(2);
4707 const local_ordinal_type r0 = part2packrowidx0_sub(partidx, local_subpartidx);
4708 const local_ordinal_type nrows = partptr_sub(subpartidx, 1) - partptr_sub(subpartidx, 0);
4710 internal_vector_scratch_type_3d_view
4711 WW(member.team_scratch(0), blocksize, blocksize, vector_loop_size);
4715 const auto one = Kokkos::ArithTraits<btdm_magnitude_type>::one();
4717 typedef SolveTridiagsDefaultModeAndAlgo<typename execution_space::memory_space> default_mode_and_algo_type;
4719 typedef typename default_mode_and_algo_type::mode_type default_mode_type;
4720 typedef typename default_mode_and_algo_type::single_vector_algo_type default_algo_type;
4722 if (local_subpartidx == 0) {
4723 Kokkos::parallel_for(Kokkos::ThreadVectorRange(member, vector_loop_size), [&](
const int &v) {
4724 auto v_2 = Kokkos::subview(X_internal_vector_values, r0 + nrows, Kokkos::ALL(), 0, v);
4726 for (local_ordinal_type row = 0; row < nrows; ++row) {
4727 auto v_1 = Kokkos::subview(X_internal_vector_values, r0 + row, Kokkos::ALL(), 0, v);
4728 auto E = Kokkos::subview(e_internal_vector_values, 0, r0 + row, Kokkos::ALL(), Kokkos::ALL(), v);
4730 KOKKOSBATCHED_GEMV_NO_TRANSPOSE_INTERNAL_INVOKE(default_mode_type, default_algo_type,
4732 blocksize, blocksize,
4734 E.data(), E.stride(0), E.stride(1),
4735 v_2.data(), v_2.stride(0),
4737 v_1.data(), v_1.stride(0));
4740 }
else if (local_subpartidx == (local_ordinal_type)part2packrowidx0_sub.extent(1) - 2) {
4741 Kokkos::parallel_for(Kokkos::ThreadVectorRange(member, vector_loop_size), [&](
const int &v) {
4742 auto v_2 = Kokkos::subview(X_internal_vector_values, r0 - 1, Kokkos::ALL(), 0, v);
4744 for (local_ordinal_type row = 0; row < nrows; ++row) {
4745 auto v_1 = Kokkos::subview(X_internal_vector_values, r0 + row, Kokkos::ALL(), 0, v);
4746 auto E = Kokkos::subview(e_internal_vector_values, 1, r0 + row, Kokkos::ALL(), Kokkos::ALL(), v);
4748 KOKKOSBATCHED_GEMV_NO_TRANSPOSE_INTERNAL_INVOKE(default_mode_type, default_algo_type,
4750 blocksize, blocksize,
4752 E.data(), E.stride(0), E.stride(1),
4753 v_2.data(), v_2.stride(0),
4755 v_1.data(), v_1.stride(0));
4759 Kokkos::parallel_for(Kokkos::ThreadVectorRange(member, vector_loop_size), [&](
const int &v) {
4761 auto v_2 = Kokkos::subview(X_internal_vector_values, r0 + nrows, Kokkos::ALL(), 0, v);
4763 for (local_ordinal_type row = 0; row < nrows; ++row) {
4764 auto v_1 = Kokkos::subview(X_internal_vector_values, r0 + row, Kokkos::ALL(), 0, v);
4765 auto E = Kokkos::subview(e_internal_vector_values, 0, r0 + row, Kokkos::ALL(), Kokkos::ALL(), v);
4767 KOKKOSBATCHED_GEMV_NO_TRANSPOSE_INTERNAL_INVOKE(default_mode_type, default_algo_type,
4769 blocksize, blocksize,
4771 E.data(), E.stride(0), E.stride(1),
4772 v_2.data(), v_2.stride(0),
4774 v_1.data(), v_1.stride(0));
4778 auto v_2 = Kokkos::subview(X_internal_vector_values, r0 - 1, Kokkos::ALL(), 0, v);
4780 for (local_ordinal_type row = 0; row < nrows; ++row) {
4781 auto v_1 = Kokkos::subview(X_internal_vector_values, r0 + row, Kokkos::ALL(), 0, v);
4782 auto E = Kokkos::subview(e_internal_vector_values, 1, r0 + row, Kokkos::ALL(), Kokkos::ALL(), v);
4784 KOKKOSBATCHED_GEMV_NO_TRANSPOSE_INTERNAL_INVOKE(default_mode_type, default_algo_type,
4786 blocksize, blocksize,
4788 E.data(), E.stride(0), E.stride(1),
4789 v_2.data(), v_2.stride(0),
4791 v_1.data(), v_1.stride(0));
4799 KOKKOS_INLINE_FUNCTION
void
4800 operator()(
const SingleVectorCopyToFlatTag<B> &,
const member_type &member)
const {
4801 const local_ordinal_type packidx = member.league_rank();
4802 const local_ordinal_type partidx = packptr(packidx);
4803 const local_ordinal_type npacks = packptr(packidx + 1) - partidx;
4804 const local_ordinal_type pri0 = part2packrowidx0(partidx);
4805 const local_ordinal_type blocksize = (B == 0 ? D_internal_vector_values.extent(1) : B);
4806 const local_ordinal_type num_vectors = 1;
4808 Kokkos::parallel_for(Kokkos::ThreadVectorRange(member, vector_loop_size), [&](
const int &v) {
4809 copyToFlatMultiVector(member, partidx, npacks, pri0, v, blocksize, num_vectors);
4814 KOKKOS_INLINE_FUNCTION
void
4815 operator()(
const SingleZeroingTag<B> &,
const member_type &member)
const {
4816 Kokkos::single(Kokkos::PerTeam(member), [&]() {
4817 Z_scalar_vector(member.league_rank()) = impl_scalar_type(0);
4821 void run(
const impl_scalar_type_2d_view_tpetra &Y,
4822 const impl_scalar_type_1d_view &Z) {
4823 IFPACK2_BLOCKTRIDICONTAINER_PROFILER_REGION_BEGIN;
4824 IFPACK2_BLOCKHELPER_TIMER(
"BlockTriDi::SolveTridiags", SolveTridiags);
4827 this->Y_scalar_multivector = Y;
4828 this->Z_scalar_vector = Z;
4830 const local_ordinal_type num_vectors = X_internal_vector_values.extent(2);
4831 const local_ordinal_type blocksize = D_internal_vector_values.extent(1);
4833 const local_ordinal_type team_size =
4834 SolveTridiagsDefaultModeAndAlgo<typename execution_space::memory_space>::
4835 recommended_team_size(blocksize, vector_length, internal_vector_length);
4836 const int per_team_scratch = internal_vector_scratch_type_3d_view ::shmem_size(blocksize, num_vectors, vector_loop_size);
4838 #if defined(KOKKOS_ENABLE_DEPRECATED_CODE)
4839 #define BLOCKTRIDICONTAINER_DETAILS_SOLVETRIDIAGS(B) \
4840 if (num_vectors == 1) { \
4841 const Kokkos::TeamPolicy<execution_space, SingleVectorTag<B>> \
4842 policy(packptr.extent(0) - 1, team_size, vector_loop_size); \
4843 Kokkos::parallel_for("SolveTridiags::TeamPolicy::run<SingleVector>", \
4844 policy.set_scratch_size(0, Kokkos::PerTeam(per_team_scratch)), *this); \
4846 const Kokkos::TeamPolicy<execution_space, MultiVectorTag<B>> \
4847 policy(packptr.extent(0) - 1, team_size, vector_loop_size); \
4848 Kokkos::parallel_for("SolveTridiags::TeamPolicy::run<MultiVector>", \
4849 policy.set_scratch_size(0, Kokkos::PerTeam(per_team_scratch)), *this); \
4853 #define BLOCKTRIDICONTAINER_DETAILS_SOLVETRIDIAGS(B) \
4854 if (num_vectors == 1) { \
4855 if (packindices_schur.extent(1) <= 0) { \
4856 Kokkos::TeamPolicy<execution_space, SingleVectorTag<B>> \
4857 policy(packptr.extent(0) - 1, team_size, vector_loop_size); \
4858 policy.set_scratch_size(0, Kokkos::PerTeam(per_team_scratch)); \
4859 Kokkos::parallel_for("SolveTridiags::TeamPolicy::run<SingleVector>", \
4863 Kokkos::TeamPolicy<execution_space, SingleZeroingTag<B>> \
4864 policy(packptr.extent(0) - 1, team_size, vector_loop_size); \
4865 Kokkos::parallel_for("SolveTridiags::TeamPolicy::run<SingleZeroingTag>", \
4869 IFPACK2_BLOCKHELPER_TIMER("BlockTriDi::ApplyInverseJacobi::SingleVectorSubLineTag", SingleVectorSubLineTag0); \
4870 write4DMultiVectorValuesToFile(part2packrowidx0_sub.extent(0), X_internal_scalar_values, "x_scalar_values_before_SingleVectorSubLineTag.mm"); \
4871 Kokkos::TeamPolicy<execution_space, SingleVectorSubLineTag<B>> \
4872 policy(packindices_sub.extent(0), team_size, vector_loop_size); \
4873 policy.set_scratch_size(0, Kokkos::PerTeam(per_team_scratch)); \
4874 Kokkos::parallel_for("SolveTridiags::TeamPolicy::run<SingleVector>", \
4876 write4DMultiVectorValuesToFile(part2packrowidx0_sub.extent(0), X_internal_scalar_values, "x_scalar_values_after_SingleVectorSubLineTag.mm"); \
4877 IFPACK2_BLOCKHELPER_TIMER_FENCE(execution_space) \
4880 IFPACK2_BLOCKHELPER_TIMER("BlockTriDi::ApplyInverseJacobi::SingleVectorApplyCTag", SingleVectorApplyCTag0); \
4881 write4DMultiVectorValuesToFile(part2packrowidx0_sub.extent(0), X_internal_scalar_values, "x_scalar_values_before_SingleVectorApplyCTag.mm"); \
4882 Kokkos::TeamPolicy<execution_space, SingleVectorApplyCTag<B>> \
4883 policy(packindices_sub.extent(0), team_size, vector_loop_size); \
4884 policy.set_scratch_size(0, Kokkos::PerTeam(per_team_scratch)); \
4885 Kokkos::parallel_for("SolveTridiags::TeamPolicy::run<SingleVector>", \
4887 write4DMultiVectorValuesToFile(part2packrowidx0_sub.extent(0), X_internal_scalar_values, "x_scalar_values_after_SingleVectorApplyCTag.mm"); \
4888 IFPACK2_BLOCKHELPER_TIMER_FENCE(execution_space) \
4891 IFPACK2_BLOCKHELPER_TIMER("BlockTriDi::ApplyInverseJacobi::SingleVectorSchurTag", SingleVectorSchurTag0); \
4892 write4DMultiVectorValuesToFile(part2packrowidx0_sub.extent(0), X_internal_scalar_values, "x_scalar_values_before_SingleVectorSchurTag.mm"); \
4893 Kokkos::TeamPolicy<execution_space, SingleVectorSchurTag<B>> \
4894 policy(packindices_schur.extent(0), team_size, vector_loop_size); \
4895 policy.set_scratch_size(0, Kokkos::PerTeam(per_team_scratch)); \
4896 Kokkos::parallel_for("SolveTridiags::TeamPolicy::run<SingleVector>", \
4898 write4DMultiVectorValuesToFile(part2packrowidx0_sub.extent(0), X_internal_scalar_values, "x_scalar_values_after_SingleVectorSchurTag.mm"); \
4899 IFPACK2_BLOCKHELPER_TIMER_FENCE(execution_space) \
4902 IFPACK2_BLOCKHELPER_TIMER("BlockTriDi::ApplyInverseJacobi::SingleVectorApplyETag", SingleVectorApplyETag0); \
4903 write4DMultiVectorValuesToFile(part2packrowidx0_sub.extent(0), X_internal_scalar_values, "x_scalar_values_before_SingleVectorApplyETag.mm"); \
4904 Kokkos::TeamPolicy<execution_space, SingleVectorApplyETag<B>> \
4905 policy(packindices_sub.extent(0), team_size, vector_loop_size); \
4906 policy.set_scratch_size(0, Kokkos::PerTeam(per_team_scratch)); \
4907 Kokkos::parallel_for("SolveTridiags::TeamPolicy::run<SingleVector>", \
4909 write4DMultiVectorValuesToFile(part2packrowidx0_sub.extent(0), X_internal_scalar_values, "x_scalar_values_after_SingleVectorApplyETag.mm"); \
4910 IFPACK2_BLOCKHELPER_TIMER_FENCE(execution_space) \
4913 Kokkos::TeamPolicy<execution_space, SingleVectorCopyToFlatTag<B>> \
4914 policy(packptr.extent(0) - 1, team_size, vector_loop_size); \
4915 Kokkos::parallel_for("SolveTridiags::TeamPolicy::run<SingleVectorCopyToFlatTag>", \
4920 Kokkos::TeamPolicy<execution_space, MultiVectorTag<B>> \
4921 policy(packptr.extent(0) - 1, team_size, vector_loop_size); \
4922 policy.set_scratch_size(0, Kokkos::PerTeam(per_team_scratch)); \
4923 Kokkos::parallel_for("SolveTridiags::TeamPolicy::run<MultiVector>", \
4928 switch (blocksize) {
4929 case 3: BLOCKTRIDICONTAINER_DETAILS_SOLVETRIDIAGS(3);
4930 case 5: BLOCKTRIDICONTAINER_DETAILS_SOLVETRIDIAGS(5);
4931 case 6: BLOCKTRIDICONTAINER_DETAILS_SOLVETRIDIAGS(6);
4932 case 7: BLOCKTRIDICONTAINER_DETAILS_SOLVETRIDIAGS(7);
4933 case 10: BLOCKTRIDICONTAINER_DETAILS_SOLVETRIDIAGS(10);
4934 case 11: BLOCKTRIDICONTAINER_DETAILS_SOLVETRIDIAGS(11);
4935 case 12: BLOCKTRIDICONTAINER_DETAILS_SOLVETRIDIAGS(12);
4936 case 13: BLOCKTRIDICONTAINER_DETAILS_SOLVETRIDIAGS(13);
4937 case 16: BLOCKTRIDICONTAINER_DETAILS_SOLVETRIDIAGS(16);
4938 case 17: BLOCKTRIDICONTAINER_DETAILS_SOLVETRIDIAGS(17);
4939 case 18: BLOCKTRIDICONTAINER_DETAILS_SOLVETRIDIAGS(18);
4940 case 19: BLOCKTRIDICONTAINER_DETAILS_SOLVETRIDIAGS(19);
4941 default: BLOCKTRIDICONTAINER_DETAILS_SOLVETRIDIAGS(0);
4943 #undef BLOCKTRIDICONTAINER_DETAILS_SOLVETRIDIAGS
4945 IFPACK2_BLOCKTRIDICONTAINER_PROFILER_REGION_END;
4946 IFPACK2_BLOCKHELPER_TIMER_FENCE(execution_space)
4953 template <
typename MatrixType>
4955 const Teuchos::RCP<
const typename BlockHelperDetails::ImplType<MatrixType>::tpetra_row_matrix_type> &A,
4956 const Teuchos::RCP<
const typename BlockHelperDetails::ImplType<MatrixType>::tpetra_crs_graph_type> &G,
4957 const Teuchos::RCP<
const typename BlockHelperDetails::ImplType<MatrixType>::tpetra_import_type> &tpetra_importer,
4958 const Teuchos::RCP<AsyncableImport<MatrixType>> &async_importer,
4959 const bool overlap_communication_and_computation,
4961 const typename BlockHelperDetails::ImplType<MatrixType>::tpetra_multivector_type &X,
4962 typename BlockHelperDetails::ImplType<MatrixType>::tpetra_multivector_type &Y,
4963 typename BlockHelperDetails::ImplType<MatrixType>::tpetra_multivector_type &Z,
4964 typename BlockHelperDetails::ImplType<MatrixType>::impl_scalar_type_1d_view &W,
4966 const BlockHelperDetails::PartInterface<MatrixType> &interf,
4969 typename BlockHelperDetails::ImplType<MatrixType>::vector_type_1d_view &work,
4974 const int max_num_sweeps,
4975 const typename BlockHelperDetails::ImplType<MatrixType>::magnitude_type tol,
4976 const int check_tol_every) {
4977 IFPACK2_BLOCKHELPER_TIMER(
"BlockTriDi::ApplyInverseJacobi", ApplyInverseJacobi);
4980 using node_memory_space =
typename impl_type::node_memory_space;
4981 using local_ordinal_type =
typename impl_type::local_ordinal_type;
4982 using size_type =
typename impl_type::size_type;
4983 using impl_scalar_type =
typename impl_type::impl_scalar_type;
4984 using magnitude_type =
typename impl_type::magnitude_type;
4985 using local_ordinal_type_1d_view =
typename impl_type::local_ordinal_type_1d_view;
4986 using vector_type_1d_view =
typename impl_type::vector_type_1d_view;
4987 using vector_type_3d_view =
typename impl_type::vector_type_3d_view;
4988 using tpetra_multivector_type =
typename impl_type::tpetra_multivector_type;
4990 using impl_scalar_type_1d_view =
typename impl_type::impl_scalar_type_1d_view;
4994 "Neither Tpetra importer nor Async importer is null.");
4997 "Maximum number of sweeps must be >= 1.");
5000 const bool is_seq_method_requested = !tpetra_importer.is_null();
5001 const bool is_async_importer_active = !async_importer.is_null();
5002 const bool is_norm_manager_active = tol > Kokkos::ArithTraits<magnitude_type>::zero();
5003 const magnitude_type tolerance = tol * tol;
5004 const local_ordinal_type blocksize = btdm.values.extent(1);
5005 const local_ordinal_type num_vectors = Y.getNumVectors();
5006 const local_ordinal_type num_blockrows = interf.part2packrowidx0_back;
5008 const impl_scalar_type zero(0.0);
5010 TEUCHOS_TEST_FOR_EXCEPT_MSG(is_norm_manager_active && is_seq_method_requested,
5011 "The seq method for applyInverseJacobi, "
5012 <<
"which in any case is for developer use only, "
5013 <<
"does not support norm-based termination.");
5014 const bool device_accessible_from_host = Kokkos::SpaceAccessibility<
5015 Kokkos::DefaultHostExecutionSpace, node_memory_space>::accessible;
5017 std::invalid_argument,
5018 "The seq method for applyInverseJacobi, "
5019 <<
"which in any case is for developer use only, "
5020 <<
"only supports memory spaces accessible from host.");
5023 const size_type work_span_required = num_blockrows * num_vectors * blocksize;
5024 if (work.span() < work_span_required)
5025 work = vector_type_1d_view(
"vector workspace 1d view", work_span_required);
5028 const local_ordinal_type W_size = interf.packptr.extent(0) - 1;
5029 if (local_ordinal_type(W.extent(0)) < W_size)
5030 W = impl_scalar_type_1d_view(
"W", W_size);
5032 typename impl_type::impl_scalar_type_2d_view_tpetra remote_multivector;
5034 if (is_seq_method_requested) {
5035 if (Z.getNumVectors() != Y.getNumVectors())
5036 Z = tpetra_multivector_type(tpetra_importer->getTargetMap(), num_vectors,
false);
5038 if (is_async_importer_active) {
5040 async_importer->createDataBuffer(num_vectors);
5041 remote_multivector = async_importer->getRemoteMultiVectorLocalView();
5047 vector_type_3d_view pmv(work.data(), num_blockrows, blocksize, num_vectors);
5048 const auto XX = X.getLocalViewDevice(Tpetra::Access::ReadOnly);
5049 const auto YY = Y.getLocalViewDevice(Tpetra::Access::ReadWrite);
5050 const auto ZZ = Z.getLocalViewDevice(Tpetra::Access::ReadWrite);
5051 if (is_y_zero) Kokkos::deep_copy(YY, zero);
5054 SolveTridiags<MatrixType> solve_tridiags(interf, btdm, pmv,
5055 damping_factor, is_norm_manager_active);
5057 const local_ordinal_type_1d_view dummy_local_ordinal_type_1d_view;
5059 auto A_crs = Teuchos::rcp_dynamic_cast<
const typename impl_type::tpetra_crs_matrix_type>(A);
5060 auto A_bcrs = Teuchos::rcp_dynamic_cast<
const typename impl_type::tpetra_block_crs_matrix_type>(A);
5062 bool hasBlockCrsMatrix = !A_bcrs.is_null();
5065 const auto g = hasBlockCrsMatrix ? A_bcrs->getCrsGraph() : *(A_crs->getCrsGraph());
5067 BlockHelperDetails::ComputeResidualVector<MatrixType>
5068 compute_residual_vector(amd, G->getLocalGraphDevice(), g.getLocalGraphDevice(), blocksize, interf,
5069 is_async_importer_active ? async_importer->dm2cm : dummy_local_ordinal_type_1d_view,
5073 if (is_norm_manager_active)
5074 norm_manager.setCheckFrequency(check_tol_every);
5078 for (; sweep < max_num_sweeps; ++sweep) {
5082 multivector_converter.run(XX);
5084 if (is_seq_method_requested) {
5088 Z.doImport(Y, *tpetra_importer, Tpetra::REPLACE);
5089 compute_residual_vector.run(YY, XX, ZZ);
5092 multivector_converter.run(YY);
5096 if (overlap_communication_and_computation || !is_async_importer_active) {
5097 if (is_async_importer_active) async_importer->asyncSendRecv(YY);
5099 compute_residual_vector.run(pmv, XX, YY, remote_multivector,
true);
5100 if (is_norm_manager_active && norm_manager.checkDone(sweep, tolerance)) {
5101 if (is_async_importer_active) async_importer->cancel();
5104 if (is_async_importer_active) {
5105 async_importer->syncRecv();
5107 compute_residual_vector.run(pmv, XX, YY, remote_multivector,
false);
5110 if (is_async_importer_active)
5111 async_importer->syncExchange(YY);
5112 if (is_norm_manager_active && norm_manager.checkDone(sweep, tolerance))
break;
5114 compute_residual_vector.run(pmv, XX, YY, remote_multivector);
5122 solve_tridiags.run(YY, W);
5125 if (is_norm_manager_active) {
5127 BlockHelperDetails::reduceVector<MatrixType>(W, norm_manager.getBuffer());
5128 if (sweep + 1 == max_num_sweeps) {
5129 norm_manager.ireduce(sweep,
true);
5130 norm_manager.checkDone(sweep + 1, tolerance,
true);
5132 norm_manager.ireduce(sweep);
5140 if (is_norm_manager_active) norm_manager.finalize();
5147 template <
typename MatrixType,
int B>
5148 int applyFusedBlockJacobi_Impl(
5149 const Teuchos::RCP<
const typename BlockHelperDetails::ImplType<MatrixType>::tpetra_import_type> &tpetra_importer,
5150 const Teuchos::RCP<AsyncableImport<MatrixType>> &async_importer,
5151 const bool overlap_communication_and_computation,
5153 const typename BlockHelperDetails::ImplType<MatrixType>::tpetra_multivector_type &X,
5154 typename BlockHelperDetails::ImplType<MatrixType>::tpetra_multivector_type &Y,
5155 typename BlockHelperDetails::ImplType<MatrixType>::impl_scalar_type_1d_view &W,
5157 const BlockHelperDetails::PartInterface<MatrixType> &interf,
5158 const BlockTridiags<MatrixType> &btdm,
5160 typename BlockHelperDetails::ImplType<MatrixType>::impl_scalar_type_1d_view &work,
5165 const int max_num_sweeps,
5166 const typename BlockHelperDetails::ImplType<MatrixType>::magnitude_type tol,
5167 const int check_tol_every) {
5169 using local_ordinal_type =
typename impl_type::local_ordinal_type;
5170 using size_type =
typename impl_type::size_type;
5171 using magnitude_type =
typename impl_type::magnitude_type;
5172 using impl_scalar_type_1d_view =
typename impl_type::impl_scalar_type_1d_view;
5173 using impl_scalar_type_2d_view_tpetra =
typename impl_type::impl_scalar_type_2d_view_tpetra;
5177 "Neither Tpetra importer nor Async importer is null.");
5180 "Maximum number of sweeps must be >= 1.");
5183 const bool is_async_importer_active = !async_importer.is_null();
5184 const bool is_norm_manager_active = tol > Kokkos::ArithTraits<magnitude_type>::zero();
5185 const magnitude_type tolerance = tol * tol;
5186 const local_ordinal_type blocksize = btdm.d_inv.extent(1);
5187 const local_ordinal_type num_vectors = Y.getNumVectors();
5188 const local_ordinal_type num_blockrows = interf.nparts;
5190 typename impl_type::impl_scalar_type_2d_view_tpetra remote_multivector;
5192 if (is_async_importer_active) {
5194 async_importer->createDataBuffer(num_vectors);
5195 remote_multivector = async_importer->getRemoteMultiVectorLocalView();
5199 const auto XX = X.getLocalViewDevice(Tpetra::Access::ReadOnly);
5200 const auto YY = Y.getLocalViewDevice(Tpetra::Access::ReadWrite);
5202 const bool two_pass_residual =
5203 overlap_communication_and_computation && is_async_importer_active;
5208 size_t(num_blockrows) * blocksize * num_vectors != YY.extent(0) * YY.extent(1),
5209 "Local LHS vector (YY) has total size " << YY.extent(0) <<
"x" << YY.extent(1) <<
" = " << YY.extent(0) * YY.extent(1) <<
",\n"
5210 <<
"but expected " << num_blockrows <<
"x" << blocksize <<
"x" << num_vectors <<
" = " << size_t(num_blockrows) * blocksize * num_vectors <<
'\n');
5211 size_type work_required = size_type(num_blockrows) * blocksize * num_vectors;
5212 if (work.extent(0) < work_required) {
5216 Unmanaged<impl_scalar_type_2d_view_tpetra> y_doublebuf(work.data(), num_blockrows * blocksize, num_vectors);
5219 if (W.extent(0) != size_t(num_blockrows))
5223 BlockHelperDetails::ComputeResidualAndSolve_SolveOnly<MatrixType, B>
5224 functor_solve_only(amd, btdm.d_inv, W, blocksize, damping_factor);
5225 BlockHelperDetails::ComputeResidualAndSolve_1Pass<MatrixType, B>
5226 functor_1pass(amd, btdm.d_inv, W, blocksize, damping_factor);
5227 BlockHelperDetails::ComputeResidualAndSolve_2Pass<MatrixType, B>
5228 functor_2pass(amd, btdm.d_inv, W, blocksize, damping_factor);
5231 if (is_norm_manager_active)
5232 norm_manager.setCheckFrequency(check_tol_every);
5237 Unmanaged<impl_scalar_type_2d_view_tpetra> y_buffers[2] = {YY, y_doublebuf};
5242 for (; sweep < max_num_sweeps; ++sweep) {
5245 functor_solve_only.run(XX, y_buffers[1 - current_y]);
5248 if (overlap_communication_and_computation || !is_async_importer_active) {
5249 if (is_async_importer_active) async_importer->asyncSendRecv(y_buffers[current_y]);
5250 if (two_pass_residual) {
5253 functor_2pass.run_pass1(XX, y_buffers[current_y], y_buffers[1 - current_y]);
5257 functor_1pass.run(XX, y_buffers[current_y], remote_multivector, y_buffers[1 - current_y]);
5259 if (is_norm_manager_active && norm_manager.checkDone(sweep, tolerance)) {
5260 if (is_async_importer_active) async_importer->cancel();
5263 if (is_async_importer_active) {
5264 async_importer->syncRecv();
5266 functor_2pass.run_pass2(y_buffers[current_y], remote_multivector, y_buffers[1 - current_y]);
5269 if (is_async_importer_active)
5270 async_importer->syncExchange(y_buffers[current_y]);
5271 if (is_norm_manager_active && norm_manager.checkDone(sweep, tolerance))
break;
5273 functor_1pass.run(XX, y_buffers[current_y], remote_multivector, y_buffers[1 - current_y]);
5278 if (is_norm_manager_active) {
5279 BlockHelperDetails::reduceVector<MatrixType>(W, norm_manager.getBuffer());
5280 if (sweep + 1 == max_num_sweeps) {
5281 norm_manager.ireduce(sweep,
true);
5282 norm_manager.checkDone(sweep + 1, tolerance,
true);
5284 norm_manager.ireduce(sweep);
5289 current_y = 1 - current_y;
5291 if (current_y == 1) {
5293 Kokkos::deep_copy(YY, y_doublebuf);
5297 if (is_norm_manager_active) norm_manager.finalize();
5304 template <
typename MatrixType>
5306 const Teuchos::RCP<
const typename BlockHelperDetails::ImplType<MatrixType>::tpetra_import_type> &tpetra_importer,
5307 const Teuchos::RCP<AsyncableImport<MatrixType>> &async_importer,
5308 const bool overlap_communication_and_computation,
5310 const typename BlockHelperDetails::ImplType<MatrixType>::tpetra_multivector_type &X,
5311 typename BlockHelperDetails::ImplType<MatrixType>::tpetra_multivector_type &Y,
5312 typename BlockHelperDetails::ImplType<MatrixType>::impl_scalar_type_1d_view &W,
5314 const BlockHelperDetails::PartInterface<MatrixType> &interf,
5317 typename BlockHelperDetails::ImplType<MatrixType>::impl_scalar_type_1d_view &work,
5322 const int max_num_sweeps,
5323 const typename BlockHelperDetails::ImplType<MatrixType>::magnitude_type tol,
5324 const int check_tol_every) {
5325 IFPACK2_BLOCKHELPER_TIMER(
"BlockTriDi::ApplyFusedBlockJacobi", ApplyFusedBlockJacobi);
5326 int blocksize = btdm.d_inv.extent(1);
5328 #define BLOCKTRIDICONTAINER_APPLY_FUSED_JACOBI(B) \
5330 sweep = applyFusedBlockJacobi_Impl<MatrixType, B>( \
5331 tpetra_importer, async_importer, overlap_communication_and_computation, \
5332 X, Y, W, interf, btdm, amd, work, \
5333 norm_manager, damping_factor, is_y_zero, \
5334 max_num_sweeps, tol, check_tol_every); \
5337 switch (blocksize) {
5338 case 3: BLOCKTRIDICONTAINER_APPLY_FUSED_JACOBI(3);
5339 case 5: BLOCKTRIDICONTAINER_APPLY_FUSED_JACOBI(5);
5340 case 7: BLOCKTRIDICONTAINER_APPLY_FUSED_JACOBI(7);
5341 case 9: BLOCKTRIDICONTAINER_APPLY_FUSED_JACOBI(9);
5342 case 10: BLOCKTRIDICONTAINER_APPLY_FUSED_JACOBI(10);
5343 case 11: BLOCKTRIDICONTAINER_APPLY_FUSED_JACOBI(11);
5344 case 16: BLOCKTRIDICONTAINER_APPLY_FUSED_JACOBI(16);
5345 case 17: BLOCKTRIDICONTAINER_APPLY_FUSED_JACOBI(17);
5346 case 18: BLOCKTRIDICONTAINER_APPLY_FUSED_JACOBI(18);
5347 default: BLOCKTRIDICONTAINER_APPLY_FUSED_JACOBI(0);
5349 #undef BLOCKTRIDICONTAINER_APPLY_FUSED_JACOBI
5354 template <
typename MatrixType>
5357 using part_interface_type = BlockHelperDetails::PartInterface<MatrixType>;
5358 using block_tridiags_type = BlockTridiags<MatrixType>;
5361 using async_import_type = AsyncableImport<MatrixType>;
5368 bool overlap_communication_and_computation;
5371 mutable typename impl_type::tpetra_multivector_type Z;
5372 mutable typename impl_type::impl_scalar_type_1d_view W;
5375 part_interface_type part_interface;
5376 block_tridiags_type block_tridiags;
5380 bool use_fused_jacobi;
5383 mutable typename impl_type::vector_type_1d_view work;
5385 mutable typename impl_type::impl_scalar_type_1d_view work_flat;
5386 mutable norm_manager_type norm_manager;
Definition: Ifpack2_BlockTriDiContainer_impl.hpp:141
int applyFusedBlockJacobi(const Teuchos::RCP< const typename BlockHelperDetails::ImplType< MatrixType >::tpetra_import_type > &tpetra_importer, const Teuchos::RCP< AsyncableImport< MatrixType >> &async_importer, const bool overlap_communication_and_computation, const typename BlockHelperDetails::ImplType< MatrixType >::tpetra_multivector_type &X, typename BlockHelperDetails::ImplType< MatrixType >::tpetra_multivector_type &Y, typename BlockHelperDetails::ImplType< MatrixType >::impl_scalar_type_1d_view &W, const BlockHelperDetails::PartInterface< MatrixType > &interf, const BlockTridiags< MatrixType > &btdm, const BlockHelperDetails::AmD< MatrixType > &amd, typename BlockHelperDetails::ImplType< MatrixType >::impl_scalar_type_1d_view &work, BlockHelperDetails::NormManager< MatrixType > &norm_manager, const typename BlockHelperDetails::ImplType< MatrixType >::impl_scalar_type &damping_factor, bool is_y_zero, const int max_num_sweeps, const typename BlockHelperDetails::ImplType< MatrixType >::magnitude_type tol, const int check_tol_every)
Definition: Ifpack2_BlockTriDiContainer_impl.hpp:5305
void performNumericPhase(const Teuchos::RCP< const typename BlockHelperDetails::ImplType< MatrixType >::tpetra_row_matrix_type > &A, const Teuchos::RCP< const typename BlockHelperDetails::ImplType< MatrixType >::tpetra_crs_graph_type > &G, const BlockHelperDetails::PartInterface< MatrixType > &interf, BlockTridiags< MatrixType > &btdm, const typename BlockHelperDetails::ImplType< MatrixType >::magnitude_type tiny, bool use_fused_jacobi)
Definition: Ifpack2_BlockTriDiContainer_impl.hpp:3714
void performSymbolicPhase(const Teuchos::RCP< const typename BlockHelperDetails::ImplType< MatrixType >::tpetra_row_matrix_type > &A, const Teuchos::RCP< const typename BlockHelperDetails::ImplType< MatrixType >::tpetra_crs_graph_type > &g, const BlockHelperDetails::PartInterface< MatrixType > &interf, BlockTridiags< MatrixType > &btdm, BlockHelperDetails::AmD< MatrixType > &amd, const bool overlap_communication_and_computation, const Teuchos::RCP< AsyncableImport< MatrixType >> &async_importer, bool useSeqMethod, bool use_fused_jacobi)
Definition: Ifpack2_BlockTriDiContainer_impl.hpp:1923
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
size_t size_type
Definition: Ifpack2_BlockHelper.hpp:274
BlockHelperDetails::PartInterface< MatrixType > createPartInterface(const Teuchos::RCP< const typename BlockHelperDetails::ImplType< MatrixType >::tpetra_row_matrix_type > &A, const Teuchos::RCP< const typename BlockHelperDetails::ImplType< MatrixType >::tpetra_crs_graph_type > &G, const Teuchos::Array< Teuchos::Array< typename BlockHelperDetails::ImplType< MatrixType >::local_ordinal_type >> &partitions, const typename BlockHelperDetails::ImplType< MatrixType >::local_ordinal_type n_subparts_per_part_in)
Definition: Ifpack2_BlockTriDiContainer_impl.hpp:1096
Teuchos::RCP< AsyncableImport< MatrixType > > createBlockCrsAsyncImporter(const Teuchos::RCP< const typename BlockHelperDetails::ImplType< MatrixType >::tpetra_row_matrix_type > &A)
Definition: Ifpack2_BlockTriDiContainer_impl.hpp:889
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
#define TEUCHOS_TEST_FOR_EXCEPT_MSG(throw_exception_test, msg)
BlockTridiags< MatrixType > createBlockTridiags(const BlockHelperDetails::PartInterface< MatrixType > &interf)
Definition: Ifpack2_BlockTriDiContainer_impl.hpp:1671
Kokkos::View< size_type *, device_type > size_type_1d_view
Definition: Ifpack2_BlockHelper.hpp:346
Definition: Ifpack2_BlockHelper.hpp:377
Kokkos::ViewAllocateWithoutInitializing do_not_initialize_tag
Definition: Ifpack2_BlockTriDiContainer_impl.hpp:97
void send(const Packet sendBuffer[], const Ordinal count, const int destRank, const int tag, const Comm< Ordinal > &comm)
Kokkos::Details::ArithTraits< scalar_type >::val_type impl_scalar_type
Definition: Ifpack2_BlockHelper.hpp:283
Definition: Ifpack2_BlockHelper.hpp:211
Definition: Ifpack2_BlockTriDiContainer_impl.hpp:2284
Teuchos::RCP< const typename BlockHelperDetails::ImplType< MatrixType >::tpetra_import_type > createBlockCrsTpetraImporter(const Teuchos::RCP< const typename BlockHelperDetails::ImplType< MatrixType >::tpetra_row_matrix_type > &A)
Definition: Ifpack2_BlockTriDiContainer_impl.hpp:165
RCP< CommRequest< Ordinal > > isend(const ArrayRCP< const Packet > &sendBuffer, const int destRank, const int tag, const Comm< Ordinal > &comm)
#define TEUCHOS_ASSERT(assertion_test)
Definition: Ifpack2_BlockHelper.hpp:270
int applyInverseJacobi(const Teuchos::RCP< const typename BlockHelperDetails::ImplType< MatrixType >::tpetra_row_matrix_type > &A, const Teuchos::RCP< const typename BlockHelperDetails::ImplType< MatrixType >::tpetra_crs_graph_type > &G, const Teuchos::RCP< const typename BlockHelperDetails::ImplType< MatrixType >::tpetra_import_type > &tpetra_importer, const Teuchos::RCP< AsyncableImport< MatrixType >> &async_importer, const bool overlap_communication_and_computation, const typename BlockHelperDetails::ImplType< MatrixType >::tpetra_multivector_type &X, typename BlockHelperDetails::ImplType< MatrixType >::tpetra_multivector_type &Y, typename BlockHelperDetails::ImplType< MatrixType >::tpetra_multivector_type &Z, typename BlockHelperDetails::ImplType< MatrixType >::impl_scalar_type_1d_view &W, const BlockHelperDetails::PartInterface< MatrixType > &interf, const BlockTridiags< MatrixType > &btdm, const BlockHelperDetails::AmD< MatrixType > &amd, typename BlockHelperDetails::ImplType< MatrixType >::vector_type_1d_view &work, BlockHelperDetails::NormManager< MatrixType > &norm_manager, const typename BlockHelperDetails::ImplType< MatrixType >::impl_scalar_type &damping_factor, bool is_y_zero, const int max_num_sweeps, const typename BlockHelperDetails::ImplType< MatrixType >::magnitude_type tol, const int check_tol_every)
Definition: Ifpack2_BlockTriDiContainer_impl.hpp:4954
Definition: Ifpack2_BlockTriDiContainer_impl.hpp:1603
Definition: Ifpack2_BlockComputeResidualVector.hpp:23
Definition: Ifpack2_BlockTriDiContainer_impl.hpp:3764