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 
46 #include <Teuchos_Details_MpiTypeTraits.hpp>
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,
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)
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);
95  TEUCHOS_TEST_FOR_EXCEPT_MSG
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 (useSeqMethod)
102  {
103  if (importer.is_null()) // there is no given importer, then create one
104  impl_->tpetra_importer = BlockTriDiContainerDetails::createBlockCrsTpetraImporter<MatrixType>(impl_->A);
105  else
106  impl_->tpetra_importer = importer; // if there is a given importer, use it
107  }
108  else
109  {
110  //Leave tpetra_importer null even if user provided an importer.
111  //It is not used in the performant codepath (!useSeqMethod)
112  impl_->async_importer = BlockTriDiContainerDetails::createBlockCrsAsyncImporter<MatrixType>(impl_->A);
113  }
114 
115  // as a result, there are
116  // 1) tpetra_importer is null , async_importer is null (no need for importer)
117  // 2) tpetra_importer is NOT null , async_importer is null (sequential method is used)
118  // 3) tpetra_importer is null , async_importer is NOT null (async method is used)
119 
120  // temporary disabling
121  impl_->overlap_communication_and_computation = overlapCommAndComp;
122 
123  impl_->Z = typename impl_type::tpetra_multivector_type();
124  impl_->W = typename impl_type::impl_scalar_type_1d_view();
125 
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());
129  }
130 
131  template <typename MatrixType>
132  void
133  BlockTriDiContainer<MatrixType, BlockTriDiContainerDetails::ImplSimdTag>
134  ::clearInternal ()
135  {
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>;
141 
142  impl_->A = Teuchos::null;
143  impl_->tpetra_importer = Teuchos::null;
144  impl_->async_importer = Teuchos::null;
145 
146  impl_->Z = typename impl_type::tpetra_multivector_type();
147  impl_->W = typename impl_type::impl_scalar_type_1d_view();
148 
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();
154 
155  impl_ = Teuchos::null;
156  }
157 
158  template <typename MatrixType>
160  ::BlockTriDiContainer (const Teuchos::RCP<const row_matrix_type>& matrix,
161  const Teuchos::Array<Teuchos::Array<local_ordinal_type> >& partitions,
162  const Teuchos::RCP<const import_type>& importer,
163  bool pointIndexed)
164  : Container<MatrixType>(matrix, partitions, pointIndexed)
165  {
166  const bool useSeqMethod = false;
167  const bool overlapCommAndComp = false;
168  initInternal(matrix, partitions, importer, overlapCommAndComp, useSeqMethod);
169  }
170 
171  template <typename MatrixType>
173  ::BlockTriDiContainer (const Teuchos::RCP<const row_matrix_type>& matrix,
174  const Teuchos::Array<Teuchos::Array<local_ordinal_type> >& partitions,
175  const bool overlapCommAndComp, const bool useSeqMethod)
176  : Container<MatrixType>(matrix, partitions, false)
177  {
178  initInternal(matrix, partitions, Teuchos::null, overlapCommAndComp, useSeqMethod);
179  }
180 
181  template <typename MatrixType>
184  {
185  }
186 
187  template <typename MatrixType>
188  void
190  ::setParameters (const Teuchos::ParameterList& /* List */)
191  {
192  // the solver doesn't currently take any parameters
193  }
194 
195  template <typename MatrixType>
196  void
199  {
200  this->IsInitialized_ = true;
201  // We assume that if you called this method, you intend to recompute
202  // everything.
203  this->IsComputed_ = false;
204  TEUCHOS_ASSERT(!impl_->A.is_null()); // when initInternal is called, A_ must be set
205  {
206  BlockTriDiContainerDetails::performSymbolicPhase<MatrixType>
207  (impl_->A,
208  impl_->part_interface, impl_->block_tridiags,
209  impl_->a_minus_d,
210  impl_->overlap_communication_and_computation);
211  }
212  }
213 
214  template <typename MatrixType>
215  void
218  {
219  this->IsComputed_ = false;
220  if (!this->isInitialized())
221  this->initialize();
222  {
223  BlockTriDiContainerDetails::performNumericPhase<MatrixType>
224  (impl_->A,
225  impl_->part_interface, impl_->block_tridiags,
226  Kokkos::ArithTraits<magnitude_type>::zero());
227  }
228  this->IsComputed_ = true;
229  }
230 
231  template <typename MatrixType>
232  void
235  {
236  clearInternal();
237  this->IsInitialized_ = false;
238  this->IsComputed_ = false;
240  }
241 
242  template <typename MatrixType>
243  void
244  BlockTriDiContainer<MatrixType, BlockTriDiContainerDetails::ImplSimdTag>
245  ::applyInverseJacobi (const mv_type& X, mv_type& Y, scalar_type dampingFactor,
246  bool zeroStartingSolution, int numSweeps) const
247  {
248  const magnitude_type tol = Kokkos::ArithTraits<magnitude_type>::zero();
249  const int check_tol_every = 1;
250 
251  BlockTriDiContainerDetails::applyInverseJacobi<MatrixType>
252  (impl_->A,
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,
258  impl_->work,
259  impl_->norm_manager,
260  dampingFactor,
261  zeroStartingSolution,
262  numSweeps,
263  tol,
264  check_tol_every);
265  }
266 
267  template <typename MatrixType>
268  typename BlockTriDiContainer<MatrixType, BlockTriDiContainerDetails::ImplSimdTag>::ComputeParameters
271  {
272  return ComputeParameters();
273  }
274 
275  template <typename MatrixType>
276  void
278  ::compute (const ComputeParameters& in)
279  {
280  this->IsComputed_ = false;
281  if (!this->isInitialized())
282  this->initialize();
283  {
284  BlockTriDiContainerDetails::performNumericPhase<MatrixType>
285  (impl_->A,
286  impl_->part_interface, impl_->block_tridiags,
287  in.addRadiallyToDiagonal);
288  }
289  this->IsComputed_ = true;
290  }
291 
292  template <typename MatrixType>
296  {
297  ApplyParameters in;
298  in.dampingFactor = Teuchos::ScalarTraits<scalar_type>::one();
299  return in;
300  }
301 
302  template <typename MatrixType>
303  int
305  ::applyInverseJacobi (const mv_type& X, mv_type& Y,
306  const ApplyParameters& in) const
307  {
308  int r_val = 0;
309  {
310  r_val = BlockTriDiContainerDetails::applyInverseJacobi<MatrixType>
311  (impl_->A,
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,
317  impl_->work,
318  impl_->norm_manager,
319  in.dampingFactor,
320  in.zeroStartingSolution,
321  in.maxNumSweeps,
322  in.tolerance,
323  in.checkToleranceEvery);
324  }
325  return r_val;
326  }
327 
328  template <typename MatrixType>
331  ::getNorms0 () const {
332  return impl_->norm_manager.getNorms0();
333  }
334 
335  template <typename MatrixType>
339  return impl_->norm_manager.getNormsFinal();
340  }
341 
342  template <typename MatrixType>
343  void
345  ::apply (ConstHostView /* X */, HostView /* Y */, int /* blockIndex */, Teuchos::ETransp /* mode */,
346  scalar_type /* alpha */, scalar_type /* beta */) const
347  {
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\".");
351  }
352 
353  template <typename MatrixType>
354  void
356  ::weightedApply (ConstHostView /* X */, HostView /* Y */, ConstHostView /* D */, int /* blockIndex */,
357  Teuchos::ETransp /* mode */, scalar_type /* alpha */, scalar_type /* beta */) const
358  {
359  TEUCHOS_TEST_FOR_EXCEPT_MSG(true, "BlockTriDiContainer::weightedApply is not implemented.");
360  }
361 
362  template <typename MatrixType>
363  std::ostream&
365  ::print (std::ostream& os) const
366  {
367  Teuchos::FancyOStream fos(Teuchos::rcp(&os,false));
368  fos.setOutputToRootOnly(0);
369  describe(fos);
370  return os;
371  }
372 
373  template <typename MatrixType>
374  std::string
377  {
378  std::ostringstream oss;
379  oss << Teuchos::Describable::description();
380  if (this->isInitialized()) {
381  if (this->isComputed()) {
382  oss << "{status = initialized, computed";
383  }
384  else {
385  oss << "{status = initialized, not computed";
386  }
387  }
388  else {
389  oss << "{status = not initialized, not computed";
390  }
391 
392  oss << "}";
393  return oss.str();
394  }
395 
396  template <typename MatrixType>
397  void
399  describe (Teuchos::FancyOStream& os,
400  const Teuchos::EVerbosityLevel verbLevel) const
401  {
402  using std::endl;
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
410  << endl;
411  }
412 
413  template <typename MatrixType>
414  std::string
416  ::getName() { return "Ifpack2::BlockTriDiContainer::ImplSimdTag"; }
417 
418 #define IFPACK2_BLOCKTRIDICONTAINER_INSTANT(S,LO,GO,N) \
419  template class Ifpack2::BlockTriDiContainer< Tpetra::RowMatrix<S, LO, GO, N> >;
420 }
421 #endif
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