43 #ifndef IFPACK2_BLOCKTRIDICONTAINER_DEF_HPP
44 #define IFPACK2_BLOCKTRIDICONTAINER_DEF_HPP
48 #include <Tpetra_Distributor.hpp>
49 #include <Tpetra_BlockMultiVector.hpp>
50 #include <Tpetra_BlockCrsMatrix_Helpers.hpp>
52 #include <Kokkos_ArithTraits.hpp>
53 #include <KokkosBatched_Util.hpp>
54 #include <KokkosBatched_Vector.hpp>
55 #include <KokkosBatched_AddRadial_Decl.hpp>
56 #include <KokkosBatched_AddRadial_Impl.hpp>
57 #include <KokkosBatched_Gemm_Decl.hpp>
58 #include <KokkosBatched_Gemm_Serial_Impl.hpp>
59 #include <KokkosBatched_Gemv_Decl.hpp>
60 #include <KokkosBatched_Trsm_Decl.hpp>
61 #include <KokkosBatched_Trsm_Serial_Impl.hpp>
62 #include <KokkosBatched_Trsv_Decl.hpp>
63 #include <KokkosBatched_Trsv_Serial_Impl.hpp>
64 #include <KokkosBatched_LU_Decl.hpp>
65 #include <KokkosBatched_LU_Serial_Impl.hpp>
68 #include "Ifpack2_BlockTriDiContainer_impl.hpp"
79 template <
typename MatrixType>
81 BlockTriDiContainer<MatrixType, BlockTriDiContainerDetails::ImplSimdTag>
84 const bool overlapCommAndComp,
85 const bool useSeqMethod,
88 IFPACK2_BLOCKHELPER_TIMER(
"BlockTriDiContainer::initInternal");
92 IFPACK2_BLOCKHELPER_TIMER(
"BlockTriDiContainer::createImpl");
93 impl_ =
Teuchos::rcp(
new BlockTriDiContainerDetails::ImplObject<MatrixType>());
94 IFPACK2_BLOCKHELPER_TIMER_FENCE(
typename BlockHelperDetails::ImplType<MatrixType>::execution_space)
97 using impl_type = BlockHelperDetails::ImplType<MatrixType>;
101 IFPACK2_BLOCKHELPER_TIMER(
"BlockTriDiContainer::setA");
102 impl_->A = Teuchos::rcp_dynamic_cast<
const block_crs_matrix_type>(matrix);
103 if (impl_->A.is_null()) {
105 (block_size == -1,
"A pointwise matrix and block_size = -1 were given as inputs.");
107 IFPACK2_BLOCKHELPER_TIMER(
"BlockTriDiContainer::setA::convertToBlockCrsMatrix");
108 impl_->A = Tpetra::convertToBlockCrsMatrix(*Teuchos::rcp_dynamic_cast<const crs_matrix_type>(matrix), block_size,
false);
109 IFPACK2_BLOCKHELPER_TIMER_FENCE(
typename BlockHelperDetails::ImplType<MatrixType>::execution_space)
112 IFPACK2_BLOCKHELPER_TIMER_FENCE(
typename BlockHelperDetails::ImplType<MatrixType>::execution_space)
115 impl_->tpetra_importer = Teuchos::null;
116 impl_->async_importer = Teuchos::null;
120 IFPACK2_BLOCKHELPER_TIMER(
"BlockTriDiContainer::createBlockCrsTpetraImporter useSeqMethod");
122 impl_->tpetra_importer = BlockTriDiContainerDetails::createBlockCrsTpetraImporter<MatrixType>(impl_->A);
124 impl_->tpetra_importer = importer;
125 IFPACK2_BLOCKHELPER_TIMER_FENCE(
typename BlockHelperDetails::ImplType<MatrixType>::execution_space)
129 IFPACK2_BLOCKHELPER_TIMER(
"BlockTriDiContainer::createBlockCrsTpetraImporter");
132 impl_->async_importer = BlockTriDiContainerDetails::createBlockCrsAsyncImporter<MatrixType>(impl_->A);
133 IFPACK2_BLOCKHELPER_TIMER_FENCE(
typename BlockHelperDetails::ImplType<MatrixType>::execution_space)
142 impl_->overlap_communication_and_computation = overlapCommAndComp;
145 IFPACK2_BLOCKHELPER_TIMER(
"BlockTriDiContainer::createZ");
146 impl_->Z =
typename impl_type::tpetra_multivector_type();
147 IFPACK2_BLOCKHELPER_TIMER_FENCE(
typename BlockHelperDetails::ImplType<MatrixType>::execution_space)
150 IFPACK2_BLOCKHELPER_TIMER(
"BlockTriDiContainer::createW");
151 impl_->W =
typename impl_type::impl_scalar_type_1d_view();
152 IFPACK2_BLOCKHELPER_TIMER_FENCE(
typename BlockHelperDetails::ImplType<MatrixType>::execution_space)
155 IFPACK2_BLOCKHELPER_TIMER_FENCE(
typename BlockHelperDetails::ImplType<MatrixType>::execution_space)
158 template <
typename MatrixType>
160 BlockTriDiContainer<MatrixType, BlockTriDiContainerDetails::ImplSimdTag>
163 IFPACK2_BLOCKHELPER_TIMER(
"BlockTriDiContainer::clearInternal");
164 using impl_type = BlockHelperDetails::ImplType<MatrixType>;
165 using part_interface_type = BlockHelperDetails::PartInterface<MatrixType>;
166 using block_tridiags_type = BlockTriDiContainerDetails::BlockTridiags<MatrixType>;
167 using amd_type = BlockHelperDetails::AmD<MatrixType>;
168 using norm_manager_type = BlockHelperDetails::NormManager<MatrixType>;
170 impl_->A = Teuchos::null;
171 impl_->tpetra_importer = Teuchos::null;
172 impl_->async_importer = Teuchos::null;
174 impl_->Z =
typename impl_type::tpetra_multivector_type();
175 impl_->W =
typename impl_type::impl_scalar_type_1d_view();
177 impl_->part_interface = part_interface_type();
178 impl_->block_tridiags = block_tridiags_type();
179 impl_->a_minus_d = amd_type();
180 impl_->work =
typename impl_type::vector_type_1d_view();
181 impl_->norm_manager = norm_manager_type();
183 impl_ = Teuchos::null;
184 IFPACK2_BLOCKHELPER_TIMER_FENCE(
typename BlockHelperDetails::ImplType<MatrixType>::execution_space)
187 template <
typename MatrixType>
193 :
Container<MatrixType>(matrix, partitions, pointIndexed), partitions_(partitions)
195 IFPACK2_BLOCKHELPER_TIMER(
"BlockTriDiContainer::BlockTriDiContainer");
196 const bool useSeqMethod =
false;
197 const bool overlapCommAndComp =
false;
198 initInternal(matrix, importer, overlapCommAndComp, useSeqMethod);
199 n_subparts_per_part_ = -1;
200 IFPACK2_BLOCKHELPER_TIMER_FENCE(
typename BlockHelperDetails::ImplType<MatrixType>::execution_space)
203 template <
typename MatrixType>
207 const int n_subparts_per_part,
208 const bool overlapCommAndComp,
209 const bool useSeqMethod,
210 const int block_size)
211 :
Container<MatrixType>(matrix, partitions, false), partitions_(partitions)
213 IFPACK2_BLOCKHELPER_TIMER(
"BlockTriDiContainer::BlockTriDiContainer");
214 initInternal(matrix, Teuchos::null, overlapCommAndComp, useSeqMethod, block_size);
215 n_subparts_per_part_ = n_subparts_per_part;
216 IFPACK2_BLOCKHELPER_TIMER_FENCE(
typename BlockHelperDetails::ImplType<MatrixType>::execution_space)
219 template <
typename MatrixType>
225 template <
typename MatrixType>
230 if (List.
isType<
int>(
"partitioner: subparts per part"))
231 n_subparts_per_part_ = List.
get<
int>(
"partitioner: subparts per part");
234 template <
typename MatrixType>
239 IFPACK2_BLOCKHELPER_TIMER(
"BlockTriDiContainer::initialize");
240 this->IsInitialized_ =
true;
242 IFPACK2_BLOCKHELPER_TIMER(
"BlockTriDiContainer::createPartInterfaceBlockTridiagsNormManager");
243 impl_->part_interface = BlockTriDiContainerDetails::createPartInterface<MatrixType>(impl_->A, partitions_, n_subparts_per_part_);
244 impl_->block_tridiags = BlockTriDiContainerDetails::createBlockTridiags<MatrixType>(impl_->part_interface);
246 IFPACK2_BLOCKHELPER_TIMER_FENCE(
typename BlockHelperDetails::ImplType<MatrixType>::execution_space)
250 this->IsComputed_ =
false;
253 BlockTriDiContainerDetails::performSymbolicPhase<MatrixType>
255 impl_->part_interface, impl_->block_tridiags,
257 impl_->overlap_communication_and_computation);
259 IFPACK2_BLOCKHELPER_TIMER_FENCE(
typename BlockHelperDetails::ImplType<MatrixType>::execution_space)
262 template <
typename MatrixType>
267 IFPACK2_BLOCKHELPER_TIMER(
"BlockTriDiContainer::compute");
268 this->IsComputed_ =
false;
269 if (!this->isInitialized())
272 BlockTriDiContainerDetails::performNumericPhase<MatrixType>
274 impl_->part_interface, impl_->block_tridiags,
275 Kokkos::ArithTraits<magnitude_type>::zero());
277 this->IsComputed_ =
true;
278 IFPACK2_BLOCKHELPER_TIMER_FENCE(
typename BlockHelperDetails::ImplType<MatrixType>::execution_space)
281 template <
typename MatrixType>
286 IFPACK2_BLOCKHELPER_TIMER(
"BlockTriDiContainer::clearBlocks");
288 this->IsInitialized_ =
false;
289 this->IsComputed_ =
false;
291 IFPACK2_BLOCKHELPER_TIMER_FENCE(
typename BlockHelperDetails::ImplType<MatrixType>::execution_space)
294 template <
typename MatrixType>
296 BlockTriDiContainer<MatrixType, BlockTriDiContainerDetails::ImplSimdTag>
297 ::applyInverseJacobi (
const mv_type& X, mv_type& Y, scalar_type dampingFactor,
298 bool zeroStartingSolution,
int numSweeps)
const
300 IFPACK2_BLOCKHELPER_TIMER(
"BlockTriDiContainer::applyInverseJacobi");
301 const magnitude_type tol = Kokkos::ArithTraits<magnitude_type>::zero();
302 const int check_tol_every = 1;
304 BlockTriDiContainerDetails::applyInverseJacobi<MatrixType>
306 impl_->tpetra_importer,
307 impl_->async_importer,
308 impl_->overlap_communication_and_computation,
309 X, Y, impl_->Z, impl_->W,
310 impl_->part_interface, impl_->block_tridiags, impl_->a_minus_d,
314 zeroStartingSolution,
318 IFPACK2_BLOCKHELPER_TIMER_FENCE(
typename BlockHelperDetails::ImplType<MatrixType>::execution_space)
321 template <
typename MatrixType>
322 typename BlockTriDiContainer<MatrixType, BlockTriDiContainerDetails::ImplSimdTag>::ComputeParameters
326 return ComputeParameters();
329 template <
typename MatrixType>
334 IFPACK2_BLOCKHELPER_TIMER(
"BlockTriDiContainer::compute");
335 this->IsComputed_ =
false;
336 if (!this->isInitialized())
339 BlockTriDiContainerDetails::performNumericPhase<MatrixType>
341 impl_->part_interface, impl_->block_tridiags,
342 in.addRadiallyToDiagonal);
344 this->IsComputed_ =
true;
345 IFPACK2_BLOCKHELPER_TIMER_FENCE(
typename BlockHelperDetails::ImplType<MatrixType>::execution_space)
348 template <
typename MatrixType>
358 template <
typename MatrixType>
362 const ApplyParameters& in)
const
364 IFPACK2_BLOCKHELPER_TIMER(
"BlockTriDiContainer::applyInverseJacobi");
367 r_val = BlockTriDiContainerDetails::applyInverseJacobi<MatrixType>
369 impl_->tpetra_importer,
370 impl_->async_importer,
371 impl_->overlap_communication_and_computation,
372 X, Y, impl_->Z, impl_->W,
373 impl_->part_interface, impl_->block_tridiags, impl_->a_minus_d,
377 in.zeroStartingSolution,
380 in.checkToleranceEvery);
382 IFPACK2_BLOCKHELPER_TIMER_FENCE(
typename BlockHelperDetails::ImplType<MatrixType>::execution_space)
386 template <
typename MatrixType>
390 return impl_->norm_manager.getNorms0();
393 template <
typename MatrixType>
397 return impl_->norm_manager.getNormsFinal();
400 template <
typename MatrixType>
404 scalar_type , scalar_type )
const
407 <<
"because you want to use this container's performance-portable Jacobi iteration. In "
408 <<
"that case, set \"relaxation: type\" to \"MT Split Jacobi\" rather than \"Jacobi\".");
411 template <
typename MatrixType>
420 template <
typename MatrixType>
431 template <
typename MatrixType>
436 std::ostringstream oss;
438 if (this->isInitialized()) {
439 if (this->isComputed()) {
440 oss <<
"{status = initialized, computed";
443 oss <<
"{status = initialized, not computed";
447 oss <<
"{status = not initialized, not computed";
454 template <
typename MatrixType>
462 os <<
"================================================================================" << endl
463 <<
"Ifpack2::BlockTriDiContainer" << endl
464 <<
"Number of blocks = " << this->numBlocks_ << endl
465 <<
"isInitialized() = " << this->IsInitialized_ << endl
466 <<
"isComputed() = " << this->IsComputed_ << endl
467 <<
"================================================================================" << endl
471 template <
typename MatrixType>
474 ::getName() {
return "Ifpack2::BlockTriDiContainer::ImplSimdTag"; }
476 #define IFPACK2_BLOCKTRIDICONTAINER_INSTANT(S,LO,GO,N) \
477 template class Ifpack2::BlockTriDiContainer< Tpetra::RowMatrix<S, LO, GO, N> >;
Ifpack2::BlockTriDiContainer class declaration.
T & get(const std::string &name, T def_value)
Store and solve local block tridiagonal linear problems.
Definition: Ifpack2_BlockTriDiContainer_decl.hpp:134
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
#define TEUCHOS_TEST_FOR_EXCEPT_MSG(throw_exception_test, msg)
Definition: Ifpack2_BlockHelper.hpp:388
BlockTriDiContainer(const Teuchos::RCP< const row_matrix_type > &matrix, const Teuchos::Array< Teuchos::Array< local_ordinal_type > > &partitions, const Teuchos::RCP< const import_type > &importer, bool pointIndexed)
Constructor.
Definition: Ifpack2_BlockTriDiContainer_def.hpp:189
virtual std::string description() const
basic_FancyOStream & setOutputToRootOnly(const int rootRank)
Interface for creating and solving a set of local linear problems.
Definition: Ifpack2_Container_decl.hpp:112
bool isType(const std::string &name) const
ComputeParameters createDefaultComputeParameters() const
Create a ComputeParameters struct with default values.
Definition: Ifpack2_BlockTriDiContainer_def.hpp:324
#define TEUCHOS_ASSERT(assertion_test)