Ifpack2 Templated Preconditioning Package Version 1.0
Loading...
Searching...
No Matches
Ifpack2_BlockTriDiContainer_def.hpp
1// @HEADER
2// *****************************************************************************
3// Ifpack2: Templated Object-Oriented Algebraic Preconditioner Package
4//
5// Copyright 2009 NTESS and the Ifpack2 contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9
10#ifndef IFPACK2_BLOCKTRIDICONTAINER_DEF_HPP
11#define IFPACK2_BLOCKTRIDICONTAINER_DEF_HPP
12
13#include <Teuchos_Details_MpiTypeTraits.hpp>
14
15#include <Tpetra_Distributor.hpp>
16#include <Tpetra_BlockMultiVector.hpp>
17#include <Tpetra_BlockCrsMatrix_Helpers.hpp>
18
19#include <Kokkos_ArithTraits.hpp>
20#include <KokkosBatched_Util.hpp>
21#include <KokkosBatched_Vector.hpp>
22#include <KokkosBatched_AddRadial_Decl.hpp>
23#include <KokkosBatched_AddRadial_Impl.hpp>
24#include <KokkosBatched_Gemm_Decl.hpp>
25#include <KokkosBatched_Gemm_Serial_Impl.hpp>
26#include <KokkosBatched_Gemv_Decl.hpp>
27#include <KokkosBatched_Trsm_Decl.hpp>
28#include <KokkosBatched_Trsm_Serial_Impl.hpp>
29#include <KokkosBatched_Trsv_Decl.hpp>
30#include <KokkosBatched_Trsv_Serial_Impl.hpp>
31#include <KokkosBatched_LU_Decl.hpp>
32#include <KokkosBatched_LU_Serial_Impl.hpp>
33
34#include "Ifpack2_BlockTriDiContainer_decl.hpp"
35#include "Ifpack2_BlockTriDiContainer_impl.hpp"
36
37#include <memory>
38
39
40namespace Ifpack2 {
41
45
46 template <typename MatrixType>
47 void
48 BlockTriDiContainer<MatrixType, BlockTriDiContainerDetails::ImplSimdTag>
49 ::initInternal (const Teuchos::RCP<const row_matrix_type>& matrix,
50 const Teuchos::RCP<const import_type>& importer,
51 const bool overlapCommAndComp,
52 const bool useSeqMethod,
53 const int block_size,
54 const bool explicitConversion)
55 {
56 IFPACK2_BLOCKHELPER_TIMER_WITH_FENCE("BlockTriDiContainer::initInternal", initInternal, typename BlockHelperDetails::ImplType<MatrixType>::execution_space);
57
58 // create pointer of impl
59 {
60 IFPACK2_BLOCKHELPER_TIMER("BlockTriDiContainer::createImpl", createImpl);
61 impl_ = Teuchos::rcp(new BlockTriDiContainerDetails::ImplObject<MatrixType>());
62 IFPACK2_BLOCKHELPER_TIMER_FENCE(typename BlockHelperDetails::ImplType<MatrixType>::execution_space)
63 }
64
65 using impl_type = BlockHelperDetails::ImplType<MatrixType>;
66 // using block_crs_matrix_type = typename impl_type::tpetra_block_crs_matrix_type;
67
68 {
69 IFPACK2_BLOCKHELPER_TIMER("BlockTriDiContainer::setA", setA);
70 if (explicitConversion) {
71 impl_->A = Teuchos::rcp_dynamic_cast<const block_crs_matrix_type>(matrix);
72 if (impl_->A.is_null()) {
73 TEUCHOS_TEST_FOR_EXCEPT_MSG
74 (block_size == -1, "A pointwise matrix and block_size = -1 were given as inputs.");
75 {
76 IFPACK2_BLOCKHELPER_TIMER("BlockTriDiContainer::setA::convertToBlockCrsMatrix", convertToBlockCrsMatrix);
77 impl_->A = Tpetra::convertToBlockCrsMatrix(*Teuchos::rcp_dynamic_cast<const crs_matrix_type>(matrix), block_size, true);
78 IFPACK2_BLOCKHELPER_TIMER_FENCE(typename BlockHelperDetails::ImplType<MatrixType>::execution_space)
79 }
80 }
81 }
82 else {
83 impl_->A = matrix;
84 }
85 IFPACK2_BLOCKHELPER_TIMER_FENCE(typename BlockHelperDetails::ImplType<MatrixType>::execution_space)
86 }
87
88 impl_->tpetra_importer = Teuchos::null;
89 impl_->async_importer = Teuchos::null;
90
91 if (useSeqMethod)
92 {
93 IFPACK2_BLOCKHELPER_TIMER("BlockTriDiContainer::createBlockCrsTpetraImporter useSeqMethod", useSeqMethod);
94 if (importer.is_null()) // there is no given importer, then create one
95 impl_->tpetra_importer = BlockTriDiContainerDetails::createBlockCrsTpetraImporter<MatrixType>(impl_->A);
96 else
97 impl_->tpetra_importer = importer; // if there is a given importer, use it
98 IFPACK2_BLOCKHELPER_TIMER_FENCE(typename BlockHelperDetails::ImplType<MatrixType>::execution_space)
99 }
100 else
101 {
102 IFPACK2_BLOCKHELPER_TIMER("BlockTriDiContainer::createBlockCrsTpetraImporter", createBlockCrsTpetraImporter);
103 //Leave tpetra_importer null even if user provided an importer.
104 //It is not used in the performant codepath (!useSeqMethod)
105 impl_->async_importer = BlockTriDiContainerDetails::createBlockCrsAsyncImporter<MatrixType>(impl_->A);
106 IFPACK2_BLOCKHELPER_TIMER_FENCE(typename BlockHelperDetails::ImplType<MatrixType>::execution_space)
107 }
108
109 // as a result, there are
110 // 1) tpetra_importer is null , async_importer is null (no need for importer)
111 // 2) tpetra_importer is NOT null , async_importer is null (sequential method is used)
112 // 3) tpetra_importer is null , async_importer is NOT null (async method is used)
113
114 // temporary disabling
115 impl_->overlap_communication_and_computation = overlapCommAndComp;
116
117 {
118 IFPACK2_BLOCKHELPER_TIMER("BlockTriDiContainer::createZ", createZ);
119 impl_->Z = typename impl_type::tpetra_multivector_type();
120 IFPACK2_BLOCKHELPER_TIMER_FENCE(typename BlockHelperDetails::ImplType<MatrixType>::execution_space)
121 }
122 {
123 IFPACK2_BLOCKHELPER_TIMER("BlockTriDiContainer::createW", createW);
124 impl_->W = typename impl_type::impl_scalar_type_1d_view();
125 IFPACK2_BLOCKHELPER_TIMER_FENCE(typename BlockHelperDetails::ImplType<MatrixType>::execution_space)
126 }
127
128 IFPACK2_BLOCKHELPER_TIMER_FENCE(typename BlockHelperDetails::ImplType<MatrixType>::execution_space)
129 }
130
131 template <typename MatrixType>
132 void
133 BlockTriDiContainer<MatrixType, BlockTriDiContainerDetails::ImplSimdTag>
134 ::clearInternal ()
135 {
136 IFPACK2_BLOCKHELPER_TIMER("BlockTriDiContainer::clearInternal", clearInternal);
137 using impl_type = BlockHelperDetails::ImplType<MatrixType>;
138 using part_interface_type = BlockHelperDetails::PartInterface<MatrixType>;
139 using block_tridiags_type = BlockTriDiContainerDetails::BlockTridiags<MatrixType>;
140 using amd_type = BlockHelperDetails::AmD<MatrixType>;
141 using norm_manager_type = BlockHelperDetails::NormManager<MatrixType>;
142
143 impl_->A = Teuchos::null;
144 impl_->tpetra_importer = Teuchos::null;
145 impl_->async_importer = Teuchos::null;
146
147 impl_->Z = typename impl_type::tpetra_multivector_type();
148 impl_->W = typename impl_type::impl_scalar_type_1d_view();
149
150 impl_->part_interface = part_interface_type();
151 impl_->block_tridiags = block_tridiags_type();
152 impl_->a_minus_d = amd_type();
153 impl_->work = typename impl_type::vector_type_1d_view();
154 impl_->norm_manager = norm_manager_type();
155
156 impl_ = Teuchos::null;
157 IFPACK2_BLOCKHELPER_TIMER_FENCE(typename BlockHelperDetails::ImplType<MatrixType>::execution_space)
158 }
159
160 template <typename MatrixType>
161 BlockTriDiContainer<MatrixType, BlockTriDiContainerDetails::ImplSimdTag>
162 ::BlockTriDiContainer (const Teuchos::RCP<const row_matrix_type>& matrix,
163 const Teuchos::Array<Teuchos::Array<local_ordinal_type> >& partitions,
164 const Teuchos::RCP<const import_type>& importer,
165 bool pointIndexed)
166 : Container<MatrixType>(matrix, partitions, pointIndexed), partitions_(partitions)
167 {
168 IFPACK2_BLOCKHELPER_TIMER("BlockTriDiContainer::BlockTriDiContainer", BlockTriDiContainer);
169 const bool useSeqMethod = false;
170 const bool overlapCommAndComp = false;
171 initInternal(matrix, importer, overlapCommAndComp, useSeqMethod);
172 n_subparts_per_part_ = -1;
173 block_size_ = -1;
174 IFPACK2_BLOCKHELPER_TIMER_FENCE(typename BlockHelperDetails::ImplType<MatrixType>::execution_space)
175 }
176
177 template <typename MatrixType>
178 BlockTriDiContainer<MatrixType, BlockTriDiContainerDetails::ImplSimdTag>
179 ::BlockTriDiContainer (const Teuchos::RCP<const row_matrix_type>& matrix,
180 const Teuchos::Array<Teuchos::Array<local_ordinal_type> >& partitions,
181 const int n_subparts_per_part,
182 const bool overlapCommAndComp,
183 const bool useSeqMethod,
184 const int block_size,
185 const bool explicitConversion)
186 : Container<MatrixType>(matrix, partitions, false), partitions_(partitions)
187 {
188 IFPACK2_BLOCKHELPER_TIMER("BlockTriDiContainer::BlockTriDiContainer", BlockTriDiContainer);
189 initInternal(matrix, Teuchos::null, overlapCommAndComp, useSeqMethod, block_size, explicitConversion);
190 n_subparts_per_part_ = n_subparts_per_part;
191 block_size_ = block_size;
192 IFPACK2_BLOCKHELPER_TIMER_FENCE(typename BlockHelperDetails::ImplType<MatrixType>::execution_space)
193 }
194
195 template <typename MatrixType>
196 BlockTriDiContainer<MatrixType, BlockTriDiContainerDetails::ImplSimdTag>
200
201 template <typename MatrixType>
202 void
203 BlockTriDiContainer<MatrixType, BlockTriDiContainerDetails::ImplSimdTag>
204 ::setParameters (const Teuchos::ParameterList& List)
205 {
206 if (List.isType<int>("partitioner: subparts per part"))
207 n_subparts_per_part_ = List.get<int>("partitioner: subparts per part");
208 if (List.isType<int>("partitioner: block size"))
209 block_size_ = List.get<int>("partitioner: block size");
210 }
211
212 template <typename MatrixType>
213 void
214 BlockTriDiContainer<MatrixType, BlockTriDiContainerDetails::ImplSimdTag>
215 ::initialize ()
216 {
217 IFPACK2_BLOCKHELPER_TIMER("BlockTriDiContainer::initialize", initialize);
218 this->IsInitialized_ = true;
219 {
220 auto bA = Teuchos::rcp_dynamic_cast<const block_crs_matrix_type>(impl_->A);
221 if (bA.is_null()) {
222 TEUCHOS_TEST_FOR_EXCEPT_MSG
223 (block_size_ == -1, "A pointwise matrix and block_size = -1 were given as inputs.");
224 {
225 IFPACK2_BLOCKHELPER_TIMER("BlockTriDiContainer::initialize::getBlockCrsGraph", getBlockCrsGraph);
226 auto A = Teuchos::rcp_dynamic_cast<const crs_matrix_type>(impl_->A);
227 impl_->blockGraph = Tpetra::getBlockCrsGraph(*A, block_size_, true);
228 IFPACK2_BLOCKHELPER_TIMER_FENCE(typename BlockHelperDetails::ImplType<MatrixType>::execution_space)
229 }
230 }
231 else {
232 impl_->blockGraph = Teuchos::rcpFromRef(bA->getCrsGraph());
233 }
234 }
235
236 {
237 IFPACK2_BLOCKHELPER_TIMER("BlockTriDiContainer::createPartInterfaceBlockTridiagsNormManager", createPartInterfaceBlockTridiagsNormManager);
238 impl_->part_interface = BlockTriDiContainerDetails::createPartInterface<MatrixType>(impl_->A, impl_->blockGraph, partitions_, n_subparts_per_part_);
239 impl_->block_tridiags = BlockTriDiContainerDetails::createBlockTridiags<MatrixType>(impl_->part_interface);
240 impl_->norm_manager = BlockHelperDetails::NormManager<MatrixType>(impl_->A->getComm());
241 IFPACK2_BLOCKHELPER_TIMER_FENCE(typename BlockHelperDetails::ImplType<MatrixType>::execution_space)
242 }
243
244 // We assume that if you called this method, you intend to recompute
245 // everything.
246 this->IsComputed_ = false;
247 TEUCHOS_ASSERT(!impl_->A.is_null()); // when initInternal is called, A_ must be set
248 {
249 bool useSeqMethod = !impl_->tpetra_importer.is_null();
250 BlockTriDiContainerDetails::performSymbolicPhase<MatrixType>
251 (impl_->A,
252 impl_->blockGraph,
253 impl_->part_interface,
254 impl_->block_tridiags,
255 impl_->a_minus_d,
256 impl_->overlap_communication_and_computation,
257 impl_->async_importer, useSeqMethod);
258 }
259 IFPACK2_BLOCKHELPER_TIMER_FENCE(typename BlockHelperDetails::ImplType<MatrixType>::execution_space)
260 }
261
262 template <typename MatrixType>
263 void
264 BlockTriDiContainer<MatrixType, BlockTriDiContainerDetails::ImplSimdTag>
265 ::compute ()
266 {
267 IFPACK2_BLOCKHELPER_TIMER("BlockTriDiContainer::compute", compute);
268 this->IsComputed_ = false;
269 if (!this->isInitialized())
270 this->initialize();
271 {
272 BlockTriDiContainerDetails::performNumericPhase<MatrixType>
273 (impl_->A,
274 impl_->blockGraph,
275 impl_->part_interface, impl_->block_tridiags,
276 Kokkos::ArithTraits<magnitude_type>::zero());
277 }
278 this->IsComputed_ = true;
279 IFPACK2_BLOCKHELPER_TIMER_FENCE(typename BlockHelperDetails::ImplType<MatrixType>::execution_space)
280 }
281
282 template <typename MatrixType>
283 void
284 BlockTriDiContainer<MatrixType, BlockTriDiContainerDetails::ImplSimdTag>
285 ::clearBlocks ()
286 {
287 IFPACK2_BLOCKHELPER_TIMER("BlockTriDiContainer::clearBlocks", clearBlocks);
288 clearInternal();
289 this->IsInitialized_ = false;
290 this->IsComputed_ = false;
291 Container<MatrixType>::clearBlocks();
292 IFPACK2_BLOCKHELPER_TIMER_FENCE(typename BlockHelperDetails::ImplType<MatrixType>::execution_space)
293 }
294
295 template <typename MatrixType>
296 void
297 BlockTriDiContainer<MatrixType, BlockTriDiContainerDetails::ImplSimdTag>
298 ::applyInverseJacobi (const mv_type& X, mv_type& Y, scalar_type dampingFactor,
299 bool zeroStartingSolution, int numSweeps) const
300 {
301 IFPACK2_BLOCKHELPER_TIMER("BlockTriDiContainer::applyInverseJacobi", applyInverseJacobi);
302 const magnitude_type tol = Kokkos::ArithTraits<magnitude_type>::zero();
303 const int check_tol_every = 1;
304
305 BlockTriDiContainerDetails::applyInverseJacobi<MatrixType>
306 (impl_->A,
307 impl_->blockGraph,
308 impl_->tpetra_importer,
309 impl_->async_importer,
310 impl_->overlap_communication_and_computation,
311 X, Y, impl_->Z, impl_->W,
312 impl_->part_interface, impl_->block_tridiags, impl_->a_minus_d,
313 impl_->work,
314 impl_->norm_manager,
315 dampingFactor,
316 zeroStartingSolution,
317 numSweeps,
318 tol,
319 check_tol_every);
320 IFPACK2_BLOCKHELPER_TIMER_FENCE(typename BlockHelperDetails::ImplType<MatrixType>::execution_space)
321 }
322
323 template <typename MatrixType>
324 typename BlockTriDiContainer<MatrixType, BlockTriDiContainerDetails::ImplSimdTag>::ComputeParameters
325 BlockTriDiContainer<MatrixType, BlockTriDiContainerDetails::ImplSimdTag>
326 ::createDefaultComputeParameters () const
327 {
328 return ComputeParameters();
329 }
330
331 template <typename MatrixType>
332 void
333 BlockTriDiContainer<MatrixType, BlockTriDiContainerDetails::ImplSimdTag>
334 ::compute (const ComputeParameters& in)
335 {
336 IFPACK2_BLOCKHELPER_TIMER("BlockTriDiContainer::compute", compute);
337 this->IsComputed_ = false;
338 if (!this->isInitialized())
339 this->initialize();
340 {
341 BlockTriDiContainerDetails::performNumericPhase<MatrixType>
342 (impl_->A,
343 impl_->blockGraph,
344 impl_->part_interface, impl_->block_tridiags,
345 in.addRadiallyToDiagonal);
346 }
347 this->IsComputed_ = true;
348 IFPACK2_BLOCKHELPER_TIMER_FENCE(typename BlockHelperDetails::ImplType<MatrixType>::execution_space)
349 }
350
351 template <typename MatrixType>
352 typename BlockTriDiContainer<MatrixType, BlockTriDiContainerDetails::ImplSimdTag>::ApplyParameters
353 BlockTriDiContainer<MatrixType, BlockTriDiContainerDetails::ImplSimdTag>
354 ::createDefaultApplyParameters () const
355 {
356 ApplyParameters in;
357 in.dampingFactor = Teuchos::ScalarTraits<scalar_type>::one();
358 return in;
359 }
360
361 template <typename MatrixType>
362 int
363 BlockTriDiContainer<MatrixType, BlockTriDiContainerDetails::ImplSimdTag>
364 ::applyInverseJacobi (const mv_type& X, mv_type& Y,
365 const ApplyParameters& in) const
366 {
367 IFPACK2_BLOCKHELPER_TIMER("BlockTriDiContainer::applyInverseJacobi", applyInverseJacobi);
368 int r_val = 0;
369 {
370 r_val = BlockTriDiContainerDetails::applyInverseJacobi<MatrixType>
371 (impl_->A,
372 impl_->blockGraph,
373 impl_->tpetra_importer,
374 impl_->async_importer,
375 impl_->overlap_communication_and_computation,
376 X, Y, impl_->Z, impl_->W,
377 impl_->part_interface, impl_->block_tridiags, impl_->a_minus_d,
378 impl_->work,
379 impl_->norm_manager,
380 in.dampingFactor,
381 in.zeroStartingSolution,
382 in.maxNumSweeps,
383 in.tolerance,
384 in.checkToleranceEvery);
385 }
386 IFPACK2_BLOCKHELPER_TIMER_FENCE(typename BlockHelperDetails::ImplType<MatrixType>::execution_space)
387 return r_val;
388 }
389
390 template <typename MatrixType>
391 const typename BlockTriDiContainer<MatrixType, BlockTriDiContainerDetails::ImplSimdTag>::magnitude_type
392 BlockTriDiContainer<MatrixType, BlockTriDiContainerDetails::ImplSimdTag>
393 ::getNorms0 () const {
394 return impl_->norm_manager.getNorms0();
395 }
396
397 template <typename MatrixType>
398 const typename BlockTriDiContainer<MatrixType, BlockTriDiContainerDetails::ImplSimdTag>::magnitude_type
399 BlockTriDiContainer<MatrixType, BlockTriDiContainerDetails::ImplSimdTag>
400 ::getNormsFinal () const {
401 return impl_->norm_manager.getNormsFinal();
402 }
403
404 template <typename MatrixType>
405 void
406 BlockTriDiContainer<MatrixType, BlockTriDiContainerDetails::ImplSimdTag>
407 ::apply (ConstHostView /* X */, HostView /* Y */, int /* blockIndex */, Teuchos::ETransp /* mode */,
408 scalar_type /* alpha */, scalar_type /* beta */) const
409 {
410 TEUCHOS_TEST_FOR_EXCEPT_MSG(true, "BlockTriDiContainer::apply is not implemented. You may have reached this message "
411 << "because you want to use this container's performance-portable Jacobi iteration. In "
412 << "that case, set \"relaxation: type\" to \"MT Split Jacobi\" rather than \"Jacobi\".");
413 }
414
415 template <typename MatrixType>
416 void
417 BlockTriDiContainer<MatrixType, BlockTriDiContainerDetails::ImplSimdTag>
418 ::weightedApply (ConstHostView /* X */, HostView /* Y */, ConstHostView /* D */, int /* blockIndex */,
419 Teuchos::ETransp /* mode */, scalar_type /* alpha */, scalar_type /* beta */) const
420 {
421 TEUCHOS_TEST_FOR_EXCEPT_MSG(true, "BlockTriDiContainer::weightedApply is not implemented.");
422 }
423
424 template <typename MatrixType>
425 std::ostream&
426 BlockTriDiContainer<MatrixType, BlockTriDiContainerDetails::ImplSimdTag>
427 ::print (std::ostream& os) const
428 {
429 Teuchos::FancyOStream fos(Teuchos::rcp(&os,false));
430 fos.setOutputToRootOnly(0);
431 describe(fos);
432 return os;
433 }
434
435 template <typename MatrixType>
436 std::string
437 BlockTriDiContainer<MatrixType, BlockTriDiContainerDetails::ImplSimdTag>
438 ::description () const
439 {
440 std::ostringstream oss;
441 oss << Teuchos::Describable::description();
442 if (this->isInitialized()) {
443 if (this->isComputed()) {
444 oss << "{status = initialized, computed";
445 }
446 else {
447 oss << "{status = initialized, not computed";
448 }
449 }
450 else {
451 oss << "{status = not initialized, not computed";
452 }
453
454 oss << "}";
455 return oss.str();
456 }
457
458 template <typename MatrixType>
459 void
460 BlockTriDiContainer<MatrixType, BlockTriDiContainerDetails::ImplSimdTag>::
461 describe (Teuchos::FancyOStream& os,
462 const Teuchos::EVerbosityLevel verbLevel) const
463 {
464 using std::endl;
465 if(verbLevel==Teuchos::VERB_NONE) return;
466 os << "================================================================================" << endl
467 << "Ifpack2::BlockTriDiContainer" << endl
468 << "Number of blocks = " << this->numBlocks_ << endl
469 << "isInitialized() = " << this->IsInitialized_ << endl
470 << "isComputed() = " << this->IsComputed_ << endl
471 << "================================================================================" << endl
472 << endl;
473 }
474
475 template <typename MatrixType>
476 std::string
477 BlockTriDiContainer<MatrixType, BlockTriDiContainerDetails::ImplSimdTag>
478 ::getName() { return "Ifpack2::BlockTriDiContainer::ImplSimdTag"; }
479
480#define IFPACK2_BLOCKTRIDICONTAINER_INSTANT(S,LO,GO,N) \
481 template class Ifpack2::BlockTriDiContainer< Tpetra::RowMatrix<S, LO, GO, N> >;
482}
483#endif
Ifpack2_BlockTriDiContainer_decl.hpp
Ifpack2::BlockTriDiContainer class declaration.
Ifpack2::BlockTriDiContainer
Store and solve local block tridiagonal linear problems.
Definition Ifpack2_BlockTriDiContainer_decl.hpp:101
Ifpack2::Container
Interface for creating and solving a set of local linear problems.
Definition Ifpack2_Container_decl.hpp:80
Ifpack2
Preconditioners and smoothers for Tpetra sparse matrices.
Definition Ifpack2_AdditiveSchwarz_decl.hpp:41
Ifpack2::BlockHelperDetails::NormManager
Definition Ifpack2_BlockHelper.hpp:351
Generated by doxygen 1.9.8