43 #ifndef IFPACK2_BLOCKTRIDICONTAINER_DEF_HPP 44 #define IFPACK2_BLOCKTRIDICONTAINER_DEF_HPP 46 #include <Teuchos_Details_MpiTypeTraits.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>
82 ::initInternal (
const Teuchos::RCP<const row_matrix_type>& matrix,
83 const Teuchos::Array<Teuchos::Array<local_ordinal_type> >& partitions,
84 const Teuchos::RCP<const import_type>& importer,
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);
95 TEUCHOS_TEST_FOR_EXCEPT_MSG
96 (impl_->A.is_null(),
"BlockTriDiContainer currently supports Tpetra::BlockCrsMatrix only.");
98 impl_->tpetra_importer = Teuchos::null;
99 impl_->async_importer = Teuchos::null;
103 if (importer.is_null())
104 impl_->tpetra_importer = BlockTriDiContainerDetails::createBlockCrsTpetraImporter<MatrixType>(impl_->A);
106 impl_->tpetra_importer = importer;
112 impl_->async_importer = BlockTriDiContainerDetails::createBlockCrsAsyncImporter<MatrixType>(impl_->A);
121 impl_->overlap_communication_and_computation = overlapCommAndComp;
123 impl_->Z =
typename impl_type::tpetra_multivector_type();
124 impl_->W =
typename impl_type::impl_scalar_type_1d_view();
126 impl_->part_interface = BlockTriDiContainerDetails::createPartInterface<MatrixType>(impl_->A, partitions);
127 impl_->block_tridiags = BlockTriDiContainerDetails::createBlockTridiags<MatrixType>(impl_->part_interface);
128 impl_->norm_manager = BlockTriDiContainerDetails::NormManager<MatrixType>(impl_->A->getComm());
131 template <
typename MatrixType>
133 BlockTriDiContainer<MatrixType, BlockTriDiContainerDetails::ImplSimdTag>
136 using impl_type = BlockTriDiContainerDetails::ImplType<MatrixType>;
137 using part_interface_type = BlockTriDiContainerDetails::PartInterface<MatrixType>;
138 using block_tridiags_type = BlockTriDiContainerDetails::BlockTridiags<MatrixType>;
139 using amd_type = BlockTriDiContainerDetails::AmD<MatrixType>;
140 using norm_manager_type = BlockTriDiContainerDetails::NormManager<MatrixType>;
142 impl_->A = Teuchos::null;
143 impl_->tpetra_importer = Teuchos::null;
144 impl_->async_importer = Teuchos::null;
146 impl_->Z =
typename impl_type::tpetra_multivector_type();
147 impl_->W =
typename impl_type::impl_scalar_type_1d_view();
149 impl_->part_interface = part_interface_type();
150 impl_->block_tridiags = block_tridiags_type();
151 impl_->a_minus_d = amd_type();
152 impl_->work =
typename impl_type::vector_type_1d_view();
153 impl_->norm_manager = norm_manager_type();
155 impl_ = Teuchos::null;
158 template <
typename MatrixType>
161 const Teuchos::Array<Teuchos::Array<local_ordinal_type> >& partitions,
162 const Teuchos::RCP<const import_type>& importer,
164 :
Container<MatrixType>(matrix, partitions, pointIndexed)
166 const bool useSeqMethod =
false;
167 const bool overlapCommAndComp =
false;
168 initInternal(matrix, partitions, importer, overlapCommAndComp, useSeqMethod);
171 template <
typename MatrixType>
174 const Teuchos::Array<Teuchos::Array<local_ordinal_type> >& partitions,
175 const bool overlapCommAndComp,
const bool useSeqMethod)
176 :
Container<MatrixType>(matrix, partitions, false)
178 initInternal(matrix, partitions, Teuchos::null, overlapCommAndComp, useSeqMethod);
181 template <
typename MatrixType>
187 template <
typename MatrixType>
195 template <
typename MatrixType>
200 this->IsInitialized_ =
true;
203 this->IsComputed_ =
false;
204 TEUCHOS_ASSERT(!impl_->A.is_null());
206 BlockTriDiContainerDetails::performSymbolicPhase<MatrixType>
208 impl_->part_interface, impl_->block_tridiags,
210 impl_->overlap_communication_and_computation);
214 template <
typename MatrixType>
219 this->IsComputed_ =
false;
220 if (!this->isInitialized())
223 BlockTriDiContainerDetails::performNumericPhase<MatrixType>
225 impl_->part_interface, impl_->block_tridiags,
226 Kokkos::ArithTraits<magnitude_type>::zero());
228 this->IsComputed_ =
true;
231 template <
typename MatrixType>
237 this->IsInitialized_ =
false;
238 this->IsComputed_ =
false;
242 template <
typename MatrixType>
244 BlockTriDiContainer<MatrixType, BlockTriDiContainerDetails::ImplSimdTag>
245 ::applyInverseJacobi (
const mv_type& X, mv_type& Y, scalar_type dampingFactor,
246 bool zeroStartingSolution,
int numSweeps)
const 248 const magnitude_type tol = Kokkos::ArithTraits<magnitude_type>::zero();
249 const int check_tol_every = 1;
251 BlockTriDiContainerDetails::applyInverseJacobi<MatrixType>
253 impl_->tpetra_importer,
254 impl_->async_importer,
255 impl_->overlap_communication_and_computation,
256 X, Y, impl_->Z, impl_->W,
257 impl_->part_interface, impl_->block_tridiags, impl_->a_minus_d,
261 zeroStartingSolution,
267 template <
typename MatrixType>
268 typename BlockTriDiContainer<MatrixType, BlockTriDiContainerDetails::ImplSimdTag>::ComputeParameters
272 return ComputeParameters();
275 template <
typename MatrixType>
280 this->IsComputed_ =
false;
281 if (!this->isInitialized())
284 BlockTriDiContainerDetails::performNumericPhase<MatrixType>
286 impl_->part_interface, impl_->block_tridiags,
287 in.addRadiallyToDiagonal);
289 this->IsComputed_ =
true;
292 template <
typename MatrixType>
298 in.dampingFactor = Teuchos::ScalarTraits<scalar_type>::one();
302 template <
typename MatrixType>
306 const ApplyParameters& in)
const 310 r_val = BlockTriDiContainerDetails::applyInverseJacobi<MatrixType>
312 impl_->tpetra_importer,
313 impl_->async_importer,
314 impl_->overlap_communication_and_computation,
315 X, Y, impl_->Z, impl_->W,
316 impl_->part_interface, impl_->block_tridiags, impl_->a_minus_d,
320 in.zeroStartingSolution,
323 in.checkToleranceEvery);
328 template <
typename MatrixType>
332 return impl_->norm_manager.getNorms0();
335 template <
typename MatrixType>
339 return impl_->norm_manager.getNormsFinal();
342 template <
typename MatrixType>
345 ::apply (ConstHostView , HostView ,
int , Teuchos::ETransp ,
346 scalar_type , scalar_type )
const 348 TEUCHOS_TEST_FOR_EXCEPT_MSG(
true,
"BlockTriDiContainer::apply is not implemented. You may have reached this message " 349 <<
"because you want to use this container's performance-portable Jacobi iteration. In " 350 <<
"that case, set \"relaxation: type\" to \"MT Split Jacobi\" rather than \"Jacobi\".");
353 template <
typename MatrixType>
357 Teuchos::ETransp , scalar_type , scalar_type )
const 359 TEUCHOS_TEST_FOR_EXCEPT_MSG(
true,
"BlockTriDiContainer::weightedApply is not implemented.");
362 template <
typename MatrixType>
367 Teuchos::FancyOStream fos(Teuchos::rcp(&os,
false));
368 fos.setOutputToRootOnly(0);
373 template <
typename MatrixType>
378 std::ostringstream oss;
379 oss << Teuchos::Describable::description();
380 if (this->isInitialized()) {
381 if (this->isComputed()) {
382 oss <<
"{status = initialized, computed";
385 oss <<
"{status = initialized, not computed";
389 oss <<
"{status = not initialized, not computed";
396 template <
typename MatrixType>
400 const Teuchos::EVerbosityLevel verbLevel)
const 403 if(verbLevel==Teuchos::VERB_NONE)
return;
404 os <<
"================================================================================" << endl
405 <<
"Ifpack2::BlockTriDiContainer" << endl
406 <<
"Number of blocks = " << this->numBlocks_ << endl
407 <<
"isInitialized() = " << this->IsInitialized_ << endl
408 <<
"isComputed() = " << this->IsComputed_ << endl
409 <<
"================================================================================" << endl
413 template <
typename MatrixType>
416 ::getName() {
return "Ifpack2::BlockTriDiContainer::ImplSimdTag"; }
418 #define IFPACK2_BLOCKTRIDICONTAINER_INSTANT(S,LO,GO,N) \ 419 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
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:160
Interface for creating and solving a set of local linear problems.
Definition: Ifpack2_Container_decl.hpp:112
ComputeParameters createDefaultComputeParameters() const
Create a ComputeParameters struct with default values.
Definition: Ifpack2_BlockTriDiContainer_def.hpp:270
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:73