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