Ifpack2 Templated Preconditioning Package  Version 1.0
Ifpack2_BlockTriDiContainer_def.hpp
1 /*@HEADER
2 // ***********************************************************************
3 //
4 // Ifpack2: Templated Object-Oriented Algebraic Preconditioner Package
5 // Copyright (2009) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ***********************************************************************
40 //@HEADER
41 */
42 
43 #ifndef IFPACK2_BLOCKTRIDICONTAINER_DEF_HPP
44 #define IFPACK2_BLOCKTRIDICONTAINER_DEF_HPP
45 
47 
48 #include <Tpetra_Distributor.hpp>
49 #include <Tpetra_BlockMultiVector.hpp>
50 
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>
66 
68 #include "Ifpack2_BlockTriDiContainer_impl.hpp"
69 
70 #include <memory>
71 
72 
73 namespace Ifpack2 {
74 
78 
79  template <typename MatrixType>
80  void
81  BlockTriDiContainer<MatrixType, BlockTriDiContainerDetails::ImplSimdTag>
82  ::initInternal (const Teuchos::RCP<const row_matrix_type>& matrix,
84  const Teuchos::RCP<const import_type>& importer,
85  const bool overlapCommAndComp,
86  const bool useSeqMethod)
87  {
88  // create pointer of impl
89  impl_ = Teuchos::rcp(new BlockTriDiContainerDetails::ImplObject<MatrixType>());
90 
91  using impl_type = BlockTriDiContainerDetails::ImplType<MatrixType>;
92  using block_crs_matrix_type = typename impl_type::tpetra_block_crs_matrix_type;
93 
94  impl_->A = Teuchos::rcp_dynamic_cast<const block_crs_matrix_type>(matrix);
96  (impl_->A.is_null(), "BlockTriDiContainer currently supports Tpetra::BlockCrsMatrix only.");
97 
98  impl_->tpetra_importer = Teuchos::null;
99  impl_->async_importer = Teuchos::null;
100 
101  if (importer.is_null()) // there is no given importer, then create one
102  if (useSeqMethod)
103  impl_->tpetra_importer = BlockTriDiContainerDetails::createBlockCrsTpetraImporter<MatrixType>(impl_->A);
104  else
105  impl_->async_importer = BlockTriDiContainerDetails::createBlockCrsAsyncImporter<MatrixType>(impl_->A);
106  else
107  impl_->tpetra_importer = importer; // if there is a given importer, use it
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  impl_->Z = typename impl_type::tpetra_multivector_type();
118 
119  impl_->part_interface = BlockTriDiContainerDetails::createPartInterface<MatrixType>(impl_->A, partitions);
120  impl_->block_tridiags = BlockTriDiContainerDetails::createBlockTridiags<MatrixType>(impl_->part_interface);
121  impl_->norm_manager = BlockTriDiContainerDetails::NormManager<MatrixType>(impl_->A->getComm());
122  }
123 
124  template <typename MatrixType>
125  void
126  BlockTriDiContainer<MatrixType, BlockTriDiContainerDetails::ImplSimdTag>
127  ::clearInternal ()
128  {
129  using impl_type = BlockTriDiContainerDetails::ImplType<MatrixType>;
130  using part_interface_type = BlockTriDiContainerDetails::PartInterface<MatrixType>;
131  using block_tridiags_type = BlockTriDiContainerDetails::BlockTridiags<MatrixType>;
132  using amd_type = BlockTriDiContainerDetails::AmD<MatrixType>;
133  using norm_manager_type = BlockTriDiContainerDetails::NormManager<MatrixType>;
134 
135  impl_->A = Teuchos::null;
136  impl_->tpetra_importer = Teuchos::null;
137  impl_->async_importer = Teuchos::null;
138 
139  impl_->Z = typename impl_type::tpetra_multivector_type();
140 
141  impl_->part_interface = part_interface_type();
142  impl_->block_tridiags = block_tridiags_type();
143  impl_->a_minus_d = amd_type();
144  impl_->work = typename impl_type::vector_type_1d_view();
145  impl_->norm_manager = norm_manager_type();
146 
147  impl_ = Teuchos::null;
148  }
149 
150  template <typename MatrixType>
154  const Teuchos::RCP<const import_type>& importer,
155  int OverlapLevel, scalar_type DampingFactor)
156  : Container<MatrixType>(matrix, partitions, importer, OverlapLevel, DampingFactor)
157  {
158  const bool useSeqMethod = false;
159  const bool overlapCommAndComp = false;
160  initInternal(matrix, partitions, importer, overlapCommAndComp, useSeqMethod);
161  }
162 
163  template <typename MatrixType>
167  const bool overlapCommAndComp, const bool useSeqMethod)
168  : Container<MatrixType>(matrix, partitions, Teuchos::null, 0,
169  Kokkos::ArithTraits<magnitude_type>::one())
170  {
171  initInternal(matrix, partitions, Teuchos::null, overlapCommAndComp, useSeqMethod);
172  }
173 
174  template <typename MatrixType>
177  {
178  }
179 
180  template <typename MatrixType>
181  bool
184  {
185  return IsInitialized_;
186  }
187 
188  template <typename MatrixType>
189  bool
192  {
193  return IsComputed_;
194  }
195 
196  template <typename MatrixType>
197  void
200  {
201  // the solver doesn't currently take any parameters
202  }
203 
204  template <typename MatrixType>
205  void
208  {
209  IsInitialized_ = true;
210  // We assume that if you called this method, you intend to recompute
211  // everything.
212  IsComputed_ = false;
213  TEUCHOS_ASSERT(!impl_->A.is_null()); // when initInternal is called, A_ must be set
214  {
215  BlockTriDiContainerDetails::performSymbolicPhase<MatrixType>
216  (impl_->A,
217  impl_->part_interface, impl_->block_tridiags,
218  impl_->a_minus_d,
219  impl_->overlap_communication_and_computation);
220  }
221  }
222 
223  template <typename MatrixType>
224  void
227  {
228  IsComputed_ = false;
229  if (!this->isInitialized())
230  this->initialize();
231  {
232  BlockTriDiContainerDetails::performNumericPhase<MatrixType>
233  (impl_->A,
234  impl_->part_interface, impl_->block_tridiags,
235  Kokkos::ArithTraits<magnitude_type>::zero());
236  }
237  IsComputed_ = true;
238  }
239 
240  template <typename MatrixType>
241  void
244  {
245  clearInternal();
246  IsInitialized_ = false;
247  IsComputed_ = false;
248  }
249 
250  template <typename MatrixType>
251  void
253  ::applyInverseJacobi (const mv_type& X, mv_type& Y, bool zeroStartingSolution,
254  int numSweeps) const
255  {
256  const magnitude_type tol = Kokkos::ArithTraits<magnitude_type>::zero();
257  const int check_tol_every = 1;
258 
259  BlockTriDiContainerDetails::applyInverseJacobi<MatrixType>
260  (impl_->A,
261  impl_->tpetra_importer,
262  impl_->async_importer,
263  impl_->overlap_communication_and_computation,
264  X, Y, impl_->Z,
265  impl_->part_interface, impl_->block_tridiags, impl_->a_minus_d,
266  impl_->work,
267  impl_->norm_manager,
268  this->DampingFactor_,
269  zeroStartingSolution,
270  numSweeps,
271  tol,
272  check_tol_every);
273  }
274 
275  template <typename MatrixType>
279  {
280  return ComputeParameters();
281  }
282 
283  template <typename MatrixType>
284  void
286  ::compute (const ComputeParameters& in)
287  {
288  IsComputed_ = false;
289  if (!this->isInitialized())
290  this->initialize();
291  {
292  BlockTriDiContainerDetails::performNumericPhase<MatrixType>
293  (impl_->A,
294  impl_->part_interface, impl_->block_tridiags,
295  in.addRadiallyToDiagonal);
296  }
297  IsComputed_ = true;
298  }
299 
300  template <typename MatrixType>
304  {
305  ApplyParameters in;
306  in.dampingFactor = this->DampingFactor_;
307  return in;
308  }
309 
310  template <typename MatrixType>
311  int
313  ::applyInverseJacobi (const mv_type& X, mv_type& Y,
314  const ApplyParameters& in) const
315  {
316  int r_val = 0;
317  {
318  r_val = BlockTriDiContainerDetails::applyInverseJacobi<MatrixType>
319  (impl_->A,
320  impl_->tpetra_importer,
321  impl_->async_importer,
322  impl_->overlap_communication_and_computation,
323  X, Y, impl_->Z,
324  impl_->part_interface, impl_->block_tridiags, impl_->a_minus_d,
325  impl_->work,
326  impl_->norm_manager,
327  in.dampingFactor,
328  in.zeroStartingSolution,
329  in.maxNumSweeps,
330  in.tolerance,
331  in.checkToleranceEvery);
332  }
333  return r_val;
334  }
335 
336  template <typename MatrixType>
339  ::getNorms0 () const {
340  const auto p = impl_->norm_manager.getNorms0();
341  return Teuchos::arcp(p, 0, p ? impl_->norm_manager.getNumVectors() : 0, false);
342  }
343 
344  template <typename MatrixType>
348  const auto p = impl_->norm_manager.getNormsFinal();
349  return Teuchos::arcp(p, 0, p ? impl_->norm_manager.getNumVectors() : 0, false);
350  }
351 
352  template <typename MatrixType>
353  void
355  ::apply (HostView& /* X */, HostView& /* Y */, int /* blockIndex */, int /* stride */, Teuchos::ETransp /* mode */,
356  scalar_type /* alpha */, scalar_type /* beta */) const
357  {
358  TEUCHOS_TEST_FOR_EXCEPT_MSG(true, "BlockTriDiContainer::apply is not implemented. You may have reached this message "
359  << "because you want to use this container's performance-portable Jacobi iteration. In "
360  << "that case, set \"relaxation: type\" to \"MT Split Jacobi\" rather than \"Jacobi\".");
361  }
362 
363  template <typename MatrixType>
364  void
366  ::weightedApply (HostView& /* X */, HostView& /* Y */, HostView& /* D */, int /* blockIndex */, int /* stride */,
367  Teuchos::ETransp /* mode */, scalar_type /* alpha */, scalar_type /* beta */) const
368  {
369  TEUCHOS_TEST_FOR_EXCEPT_MSG(true, "BlockTriDiContainer::weightedApply is not implemented.");
370  }
371 
372  template <typename MatrixType>
373  std::ostream&
375  ::print (std::ostream& os) const
376  {
377  Teuchos::FancyOStream fos(Teuchos::rcp(&os,false));
378  fos.setOutputToRootOnly(0);
379  describe(fos);
380  return os;
381  }
382 
383  template <typename MatrixType>
384  std::string
387  {
388  std::ostringstream oss;
390  if (isInitialized()) {
391  if (isComputed()) {
392  oss << "{status = initialized, computed";
393  }
394  else {
395  oss << "{status = initialized, not computed";
396  }
397  }
398  else {
399  oss << "{status = not initialized, not computed";
400  }
401 
402  oss << "}";
403  return oss.str();
404  }
405 
406  template <typename MatrixType>
407  void
410  const Teuchos::EVerbosityLevel verbLevel) const
411  {
412  using std::endl;
413  if(verbLevel==Teuchos::VERB_NONE) return;
414  os << "================================================================================" << endl
415  << "Ifpack2::BlockTriDiContainer" << endl
416  << "Number of blocks = " << this->numBlocks_ << endl
417  << "isInitialized() = " << IsInitialized_ << endl
418  << "isComputed() = " << IsComputed_ << endl
419  << "================================================================================" << endl
420  << endl;
421  }
422 
423  template <typename MatrixType>
424  std::string
426  ::getName() { return "Ifpack2::BlockTriDiContainer::ImplSimdTag"; }
427 
428 #define IFPACK2_BLOCKTRIDICONTAINER_INSTANT(S,LO,GO,N) \
429  template class Ifpack2::BlockTriDiContainer< Tpetra::RowMatrix<S, LO, GO, N> >;
430 }
431 #endif
Ifpack2::BlockTriDiContainer class declaration.
VERB_NONE
Store and solve local block tridiagonal linear problems.
Definition: Ifpack2_BlockTriDiContainer_decl.hpp:134
virtual std::string description() const
bool is_null() const
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:253
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
#define TEUCHOS_TEST_FOR_EXCEPT_MSG(throw_exception_test, msg)
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:152
Interface for creating and solving a local linear problem.
Definition: Ifpack2_Container.hpp:114
#define TEUCHOS_ASSERT(assertion_test)
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:72