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_Trsm_Decl.hpp>
60#include <KokkosBatched_Trsm_Serial_Impl.hpp>
61#include <KokkosBatched_Trsv_Decl.hpp>
62#include <KokkosBatched_Trsv_Serial_Impl.hpp>
63#include <KokkosBatched_LU_Decl.hpp>
64#include <KokkosBatched_LU_Serial_Impl.hpp>
67#include "Ifpack2_BlockTriDiContainer_impl.hpp"
78 template <
typename MatrixType>
80 BlockTriDiContainer<MatrixType, BlockTriDiContainerDetails::ImplSimdTag>
81 ::initInternal (
const Teuchos::RCP<const row_matrix_type>& matrix,
82 const Teuchos::Array<Teuchos::Array<local_ordinal_type> >& partitions,
83 const Teuchos::RCP<const import_type>& importer,
84 const bool overlapCommAndComp,
85 const bool useSeqMethod)
88 impl_ = Teuchos::rcp(
new BlockTriDiContainerDetails::ImplObject<MatrixType>());
90 using impl_type = BlockTriDiContainerDetails::ImplType<MatrixType>;
93 impl_->A = Teuchos::rcp_dynamic_cast<const block_crs_matrix_type>(matrix);
94 TEUCHOS_TEST_FOR_EXCEPT_MSG
95 (impl_->A.is_null(),
"BlockTriDiContainer currently supports Tpetra::BlockCrsMatrix only.");
97 impl_->tpetra_importer = Teuchos::null;
98 impl_->async_importer = Teuchos::null;
102 if (importer.is_null())
103 impl_->tpetra_importer = BlockTriDiContainerDetails::createBlockCrsTpetraImporter<MatrixType>(impl_->A);
105 impl_->tpetra_importer = importer;
111 impl_->async_importer = BlockTriDiContainerDetails::createBlockCrsAsyncImporter<MatrixType>(impl_->A);
120 impl_->overlap_communication_and_computation = overlapCommAndComp;
122 impl_->Z =
typename impl_type::tpetra_multivector_type();
123 impl_->W =
typename impl_type::impl_scalar_type_1d_view();
125 impl_->part_interface = BlockTriDiContainerDetails::createPartInterface<MatrixType>(impl_->A, partitions);
126 impl_->block_tridiags = BlockTriDiContainerDetails::createBlockTridiags<MatrixType>(impl_->part_interface);
127 impl_->norm_manager = BlockTriDiContainerDetails::NormManager<MatrixType>(impl_->A->getComm());
130 template <
typename MatrixType>
132 BlockTriDiContainer<MatrixType, BlockTriDiContainerDetails::ImplSimdTag>
135 using impl_type = BlockTriDiContainerDetails::ImplType<MatrixType>;
136 using part_interface_type = BlockTriDiContainerDetails::PartInterface<MatrixType>;
137 using block_tridiags_type = BlockTriDiContainerDetails::BlockTridiags<MatrixType>;
138 using amd_type = BlockTriDiContainerDetails::AmD<MatrixType>;
139 using norm_manager_type = BlockTriDiContainerDetails::NormManager<MatrixType>;
141 impl_->A = Teuchos::null;
142 impl_->tpetra_importer = Teuchos::null;
143 impl_->async_importer = Teuchos::null;
145 impl_->Z =
typename impl_type::tpetra_multivector_type();
146 impl_->W =
typename impl_type::impl_scalar_type_1d_view();
148 impl_->part_interface = part_interface_type();
149 impl_->block_tridiags = block_tridiags_type();
150 impl_->a_minus_d = amd_type();
151 impl_->work =
typename impl_type::vector_type_1d_view();
152 impl_->norm_manager = norm_manager_type();
154 impl_ = Teuchos::null;
157 template <
typename MatrixType>
158 BlockTriDiContainer<MatrixType, BlockTriDiContainerDetails::ImplSimdTag>
159 ::BlockTriDiContainer (
const Teuchos::RCP<const row_matrix_type>& matrix,
160 const Teuchos::Array<Teuchos::Array<local_ordinal_type> >& partitions,
161 const Teuchos::RCP<const import_type>& importer,
163 :
Container<MatrixType>(matrix, partitions, pointIndexed)
165 const bool useSeqMethod =
false;
166 const bool overlapCommAndComp =
false;
167 initInternal(matrix, partitions, importer, overlapCommAndComp, useSeqMethod);
170 template <
typename MatrixType>
173 const Teuchos::Array<Teuchos::Array<local_ordinal_type> >& partitions,
174 const bool overlapCommAndComp,
const bool useSeqMethod)
175 :
Container<MatrixType>(matrix, partitions, false)
177 initInternal(matrix, partitions, Teuchos::null, overlapCommAndComp, useSeqMethod);
180 template <
typename MatrixType>
186 template <
typename MatrixType>
194 template <
typename MatrixType>
199 this->IsInitialized_ =
true;
202 this->IsComputed_ =
false;
203 TEUCHOS_ASSERT(!impl_->A.is_null());
205 BlockTriDiContainerDetails::performSymbolicPhase<MatrixType>
207 impl_->part_interface, impl_->block_tridiags,
209 impl_->overlap_communication_and_computation);
213 template <
typename MatrixType>
218 this->IsComputed_ =
false;
219 if (!this->isInitialized())
222 BlockTriDiContainerDetails::performNumericPhase<MatrixType>
224 impl_->part_interface, impl_->block_tridiags,
225 Kokkos::ArithTraits<magnitude_type>::zero());
227 this->IsComputed_ =
true;
230 template <
typename MatrixType>
236 this->IsInitialized_ =
false;
237 this->IsComputed_ =
false;
241 template <
typename MatrixType>
243 BlockTriDiContainer<MatrixType, BlockTriDiContainerDetails::ImplSimdTag>
244 ::applyInverseJacobi (
const mv_type& X, mv_type& Y, scalar_type dampingFactor,
245 bool zeroStartingSolution,
int numSweeps)
const
247 const magnitude_type tol = Kokkos::ArithTraits<magnitude_type>::zero();
248 const int check_tol_every = 1;
250 BlockTriDiContainerDetails::applyInverseJacobi<MatrixType>
252 impl_->tpetra_importer,
253 impl_->async_importer,
254 impl_->overlap_communication_and_computation,
255 X, Y, impl_->Z, impl_->W,
256 impl_->part_interface, impl_->block_tridiags, impl_->a_minus_d,
260 zeroStartingSolution,
266 template <
typename MatrixType>
271 return ComputeParameters();
274 template <
typename MatrixType>
279 this->IsComputed_ =
false;
280 if (!this->isInitialized())
283 BlockTriDiContainerDetails::performNumericPhase<MatrixType>
285 impl_->part_interface, impl_->block_tridiags,
286 in.addRadiallyToDiagonal);
288 this->IsComputed_ =
true;
291 template <
typename MatrixType>
297 in.dampingFactor = Teuchos::ScalarTraits<scalar_type>::one();
301 template <
typename MatrixType>
305 const ApplyParameters& in)
const
309 r_val = BlockTriDiContainerDetails::applyInverseJacobi<MatrixType>
311 impl_->tpetra_importer,
312 impl_->async_importer,
313 impl_->overlap_communication_and_computation,
314 X, Y, impl_->Z, impl_->W,
315 impl_->part_interface, impl_->block_tridiags, impl_->a_minus_d,
319 in.zeroStartingSolution,
322 in.checkToleranceEvery);
327 template <
typename MatrixType>
331 return impl_->norm_manager.getNorms0();
334 template <
typename MatrixType>
338 return impl_->norm_manager.getNormsFinal();
341 template <
typename MatrixType>
344 ::apply (ConstHostView , HostView ,
int , Teuchos::ETransp ,
345 scalar_type , scalar_type )
const
347 TEUCHOS_TEST_FOR_EXCEPT_MSG(
true,
"BlockTriDiContainer::apply is not implemented. You may have reached this message "
348 <<
"because you want to use this container's performance-portable Jacobi iteration. In "
349 <<
"that case, set \"relaxation: type\" to \"MT Split Jacobi\" rather than \"Jacobi\".");
352 template <
typename MatrixType>
356 Teuchos::ETransp , scalar_type , scalar_type )
const
358 TEUCHOS_TEST_FOR_EXCEPT_MSG(
true,
"BlockTriDiContainer::weightedApply is not implemented.");
361 template <
typename MatrixType>
366 Teuchos::FancyOStream fos(Teuchos::rcp(&os,
false));
367 fos.setOutputToRootOnly(0);
372 template <
typename MatrixType>
377 std::ostringstream oss;
378 oss << Teuchos::Describable::description();
379 if (this->isInitialized()) {
380 if (this->isComputed()) {
381 oss <<
"{status = initialized, computed";
384 oss <<
"{status = initialized, not computed";
388 oss <<
"{status = not initialized, not computed";
395 template <
typename MatrixType>
398 describe (Teuchos::FancyOStream& os,
399 const Teuchos::EVerbosityLevel verbLevel)
const
402 if(verbLevel==Teuchos::VERB_NONE)
return;
403 os <<
"================================================================================" << endl
404 <<
"Ifpack2::BlockTriDiContainer" << endl
405 <<
"Number of blocks = " << this->numBlocks_ << endl
406 <<
"isInitialized() = " << this->IsInitialized_ << endl
407 <<
"isComputed() = " << this->IsComputed_ << endl
408 <<
"================================================================================" << endl
412 template <
typename MatrixType>
415 ::getName() {
return "Ifpack2::BlockTriDiContainer::ImplSimdTag"; }
417#define IFPACK2_BLOCKTRIDICONTAINER_INSTANT(S,LO,GO,N) \
418 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
Interface for creating and solving a set of local linear problems.
Definition: Ifpack2_Container_decl.hpp:113
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:74