43 #ifndef IFPACK2_BLOCKTRIDICONTAINER_DEF_HPP
44 #define IFPACK2_BLOCKTRIDICONTAINER_DEF_HPP
48 #include <Tpetra_Distributor.hpp>
49 #include <Tpetra_BlockMultiVector.hpp>
51 #include <Kokkos_ArithTraits.hpp>
52 #include <KokkosBatched_Util.hpp>
53 #include <KokkosBatched_Vector.hpp>
54 #include <KokkosBatched_AddRadial_Decl.hpp>
55 #include <KokkosBatched_AddRadial_Impl.hpp>
56 #include <KokkosBatched_Gemm_Decl.hpp>
57 #include <KokkosBatched_Gemm_Serial_Impl.hpp>
58 #include <KokkosBatched_Gemv_Decl.hpp>
59 #include <KokkosBatched_Gemv_Serial_Impl.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>
85 const bool overlapCommAndComp,
86 const bool useSeqMethod)
89 impl_ =
Teuchos::rcp(
new BlockTriDiContainerDetails::ImplObject<MatrixType>());
91 using impl_type = BlockTriDiContainerDetails::ImplType<MatrixType>;
92 using block_crs_matrix_type =
typename impl_type::tpetra_block_crs_matrix_type;
94 impl_->A = Teuchos::rcp_dynamic_cast<
const block_crs_matrix_type>(matrix);
96 (impl_->A.is_null(),
"BlockTriDiContainer currently supports Tpetra::BlockCrsMatrix only.");
98 impl_->tpetra_importer = Teuchos::null;
99 impl_->async_importer = Teuchos::null;
103 impl_->tpetra_importer = BlockTriDiContainerDetails::createBlockCrsTpetraImporter<MatrixType>(impl_->A);
105 impl_->async_importer = BlockTriDiContainerDetails::createBlockCrsAsyncImporter<MatrixType>(impl_->A);
107 impl_->tpetra_importer = importer;
115 impl_->overlap_communication_and_computation = overlapCommAndComp;
117 impl_->Z =
typename impl_type::tpetra_multivector_type();
118 impl_->W =
typename impl_type::impl_scalar_type_1d_view();
120 impl_->part_interface = BlockTriDiContainerDetails::createPartInterface<MatrixType>(impl_->A, partitions);
121 impl_->block_tridiags = BlockTriDiContainerDetails::createBlockTridiags<MatrixType>(impl_->part_interface);
122 impl_->norm_manager = BlockTriDiContainerDetails::NormManager<MatrixType>(impl_->A->getComm());
125 template <
typename MatrixType>
127 BlockTriDiContainer<MatrixType, BlockTriDiContainerDetails::ImplSimdTag>
130 using impl_type = BlockTriDiContainerDetails::ImplType<MatrixType>;
131 using part_interface_type = BlockTriDiContainerDetails::PartInterface<MatrixType>;
132 using block_tridiags_type = BlockTriDiContainerDetails::BlockTridiags<MatrixType>;
133 using amd_type = BlockTriDiContainerDetails::AmD<MatrixType>;
134 using norm_manager_type = BlockTriDiContainerDetails::NormManager<MatrixType>;
136 impl_->A = Teuchos::null;
137 impl_->tpetra_importer = Teuchos::null;
138 impl_->async_importer = Teuchos::null;
140 impl_->Z =
typename impl_type::tpetra_multivector_type();
141 impl_->W =
typename impl_type::impl_scalar_type_1d_view();
143 impl_->part_interface = part_interface_type();
144 impl_->block_tridiags = block_tridiags_type();
145 impl_->a_minus_d = amd_type();
146 impl_->work =
typename impl_type::vector_type_1d_view();
147 impl_->norm_manager = norm_manager_type();
149 impl_ = Teuchos::null;
152 template <
typename MatrixType>
157 int OverlapLevel, scalar_type DampingFactor)
158 :
Container<MatrixType>(matrix, partitions, importer, OverlapLevel, DampingFactor)
160 const bool useSeqMethod =
false;
161 const bool overlapCommAndComp =
false;
162 initInternal(matrix, partitions, importer, overlapCommAndComp, useSeqMethod);
165 template <
typename MatrixType>
169 const bool overlapCommAndComp,
const bool useSeqMethod)
170 :
Container<MatrixType>(matrix, partitions, Teuchos::null, 0,
171 Kokkos::ArithTraits<magnitude_type>::one())
173 initInternal(matrix, partitions, Teuchos::null, overlapCommAndComp, useSeqMethod);
176 template <
typename MatrixType>
182 template <
typename MatrixType>
187 return IsInitialized_;
190 template <
typename MatrixType>
198 template <
typename MatrixType>
206 template <
typename MatrixType>
211 IsInitialized_ =
true;
217 BlockTriDiContainerDetails::performSymbolicPhase<MatrixType>
219 impl_->part_interface, impl_->block_tridiags,
221 impl_->overlap_communication_and_computation);
225 template <
typename MatrixType>
231 if (!this->isInitialized())
234 BlockTriDiContainerDetails::performNumericPhase<MatrixType>
236 impl_->part_interface, impl_->block_tridiags,
237 Kokkos::ArithTraits<magnitude_type>::zero());
242 template <
typename MatrixType>
248 IsInitialized_ =
false;
252 template <
typename MatrixType>
258 const magnitude_type tol = Kokkos::ArithTraits<magnitude_type>::zero();
259 const int check_tol_every = 1;
261 BlockTriDiContainerDetails::applyInverseJacobi<MatrixType>
263 impl_->tpetra_importer,
264 impl_->async_importer,
265 impl_->overlap_communication_and_computation,
266 X, Y, impl_->Z, impl_->W,
267 impl_->part_interface, impl_->block_tridiags, impl_->a_minus_d,
270 this->DampingFactor_,
271 zeroStartingSolution,
277 template <
typename MatrixType>
282 return ComputeParameters();
285 template <
typename MatrixType>
291 if (!this->isInitialized())
294 BlockTriDiContainerDetails::performNumericPhase<MatrixType>
296 impl_->part_interface, impl_->block_tridiags,
297 in.addRadiallyToDiagonal);
302 template <
typename MatrixType>
308 in.dampingFactor = this->DampingFactor_;
312 template <
typename MatrixType>
316 const ApplyParameters& in)
const
320 r_val = BlockTriDiContainerDetails::applyInverseJacobi<MatrixType>
322 impl_->tpetra_importer,
323 impl_->async_importer,
324 impl_->overlap_communication_and_computation,
325 X, Y, impl_->Z, impl_->W,
326 impl_->part_interface, impl_->block_tridiags, impl_->a_minus_d,
330 in.zeroStartingSolution,
333 in.checkToleranceEvery);
338 template <
typename MatrixType>
342 const auto p = impl_->norm_manager.getNorms0();
343 return Teuchos::arcp(p, 0, p ? 1 : 0,
false);
346 template <
typename MatrixType>
350 const auto p = impl_->norm_manager.getNormsFinal();
351 return Teuchos::arcp(p, 0, p ? 1 : 0,
false);
354 template <
typename MatrixType>
358 scalar_type , scalar_type )
const
361 <<
"because you want to use this container's performance-portable Jacobi iteration. In "
362 <<
"that case, set \"relaxation: type\" to \"MT Split Jacobi\" rather than \"Jacobi\".");
365 template <
typename MatrixType>
374 template <
typename MatrixType>
385 template <
typename MatrixType>
390 std::ostringstream oss;
392 if (isInitialized()) {
394 oss <<
"{status = initialized, computed";
397 oss <<
"{status = initialized, not computed";
401 oss <<
"{status = not initialized, not computed";
408 template <
typename MatrixType>
416 os <<
"================================================================================" << endl
417 <<
"Ifpack2::BlockTriDiContainer" << endl
418 <<
"Number of blocks = " << this->numBlocks_ << endl
419 <<
"isInitialized() = " << IsInitialized_ << endl
420 <<
"isComputed() = " << IsComputed_ << endl
421 <<
"================================================================================" << endl
425 template <
typename MatrixType>
428 ::getName() {
return "Ifpack2::BlockTriDiContainer::ImplSimdTag"; }
430 #define IFPACK2_BLOCKTRIDICONTAINER_INSTANT(S,LO,GO,N) \
431 template class Ifpack2::BlockTriDiContainer< Tpetra::RowMatrix<S, LO, GO, N> >;
Ifpack2::BlockTriDiContainer class declaration.
Store and solve local block tridiagonal linear problems.
Definition: Ifpack2_BlockTriDiContainer_decl.hpp:134
void applyInverseJacobi(const mv_type &X, mv_type &Y, bool zeroStartingSolution=false, int numSweeps=1) const override
Compute Y := (1 - a) Y + a D^{-1} (X - R*Y).
Definition: Ifpack2_BlockTriDiContainer_def.hpp:255
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
#define TEUCHOS_TEST_FOR_EXCEPT_MSG(throw_exception_test, msg)
virtual std::string description() const
basic_FancyOStream & setOutputToRootOnly(const int rootRank)
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, int OverlapLevel, scalar_type DampingFactor)
Constructor.
Definition: Ifpack2_BlockTriDiContainer_def.hpp:154
Interface for creating and solving a local linear problem.
Definition: Ifpack2_Container.hpp:114
#define TEUCHOS_ASSERT(assertion_test)