Ifpack2 Templated Preconditioning Package  Version 1.0
Ifpack2_BlockRelaxation_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_BLOCKRELAXATION_DEF_HPP
44 #define IFPACK2_BLOCKRELAXATION_DEF_HPP
45 
47 #include "Ifpack2_LinearPartitioner.hpp"
48 #include "Ifpack2_LinePartitioner.hpp"
50 #include "Ifpack2_Details_UserPartitioner_def.hpp"
51 #include <Ifpack2_Parameters.hpp>
52 
53 namespace Ifpack2 {
54 
55 template<class MatrixType,class ContainerType>
58 {
59  if (A.getRawPtr () != A_.getRawPtr ()) { // it's a different matrix
60  IsInitialized_ = false;
61  IsComputed_ = false;
62  Partitioner_ = Teuchos::null;
63  W_ = Teuchos::null;
64 
65  if (! A.is_null ()) {
66  IsParallel_ = (A->getRowMap ()->getComm ()->getSize () != 1);
68  Teuchos::rcp_dynamic_cast<const block_crs_matrix_type> (A);
69  hasBlockCrsMatrix_ = !A_bcrs.is_null();
70  }
71  if (! Container_.is_null ()) {
72  Container_->clearBlocks();
73  }
74  NumLocalBlocks_ = 0;
75 
76  A_ = A;
77  computeImporter();
78  }
79 }
80 
81 template<class MatrixType,class ContainerType>
84 :
85  Time_ (Teuchos::rcp (new Teuchos::Time ("Ifpack2::BlockRelaxation"))),
86  Container_ (Teuchos::null),
87  Partitioner_ (Teuchos::null),
88  PartitionerType_ ("linear"),
89  NumSweeps_ (1),
90  NumLocalBlocks_(0),
91  containerType_ ("TriDi"),
92  PrecType_ (Ifpack2::Details::JACOBI),
93  ZeroStartingSolution_ (true),
94  hasBlockCrsMatrix_ (false),
95  DoBackwardGS_ (false),
96  OverlapLevel_ (0),
97  DampingFactor_ (STS::one ()),
98  IsInitialized_ (false),
99  IsComputed_ (false),
100  NumInitialize_ (0),
101  NumCompute_ (0),
102  NumApply_ (0),
103  InitializeTime_ (0.0),
104  ComputeTime_ (0.0),
105  NumLocalRows_ (0),
106  NumGlobalRows_ (0),
107  NumGlobalNonzeros_ (0),
108  W_ (Teuchos::null),
109  Importer_ (Teuchos::null)
110 {
111  setMatrix(A);
112 }
113 
114 template<class MatrixType,class ContainerType>
117 {}
118 
119 template<class MatrixType,class ContainerType>
123 {
124  Teuchos::RCP<Teuchos::ParameterList> validParams = Teuchos::parameterList ("Ifpack2::BlockRelaxation");
125 
126  validParams->set("relaxation: container", "TriDi");
127  validParams->set("relaxation: backward mode",false);
128  validParams->set("relaxation: type", "Jacobi");
129  validParams->set("relaxation: sweeps", (int)1);
130  validParams->set("relaxation: damping factor", STS::one());
131  validParams->set("relaxation: zero starting solution", true);
132  validParams->set("schwarz: compute condest", false); // mfh 24 Mar 2015: for backwards compatibility ONLY
133  validParams->set("schwarz: combine mode", "ZERO"); // use string mode for this
134  validParams->set("schwarz: use reordering", true);
135  validParams->set("schwarz: filter singletons", false);
136  validParams->set("schwarz: overlap level", (int)0);
137  validParams->set("partitioner: type", "greedy");
138  validParams->set("partitioner: local parts", (int)1);
139  validParams->set("partitioner: overlap", (int)0);
141  validParams->set("partitioner: parts", tmp0);
142  validParams->set("partitioner: maintain sparsity", false);
143 
144  Teuchos::ParameterList dummyList;
145  validParams->set("Amesos2",dummyList);
146  validParams->sublist("Amesos2").disableRecursiveValidation();
147  validParams->set("Amesos2 solver name", "KLU2");
148 
150  validParams->set("partitioner: map", tmp);
151 
152  validParams->set("partitioner: line detection threshold",(double)0.0);
153  validParams->set("partitioner: PDE equations",(int)1);
154  Teuchos::RCP<Tpetra::MultiVector<double,
155  typename MatrixType::local_ordinal_type,
156  typename MatrixType::global_ordinal_type,
157  typename MatrixType::node_type> > dummy;
158  validParams->set("partitioner: coordinates",dummy);
159 
160  return validParams;
161 }
162 
163 template<class MatrixType,class ContainerType>
164 void
167 {
168  // Note that the validation process does not change List.
170  validparams = this->getValidParameters();
171  List.validateParameters (*validparams);
172 
173  if (List.isParameter ("relaxation: container")) {
174  // If the container type isn't a string, this will throw, but it
175  // rightfully should.
176 
177  // If its value does not match the currently registered Container types,
178  // the ContainerFactory will throw with an informative message.
179  containerType_ = List.get<std::string> ("relaxation: container");
180  // Intercept more human-readable aliases for the *TriDi containers
181  if(containerType_ == "Tridiagonal") {
182  containerType_ = "TriDi";
183  }
184  if(containerType_ == "Block Tridiagonal") {
185  containerType_ = "BlockTriDi";
186  }
187  }
188 
189  if (List.isParameter ("relaxation: type")) {
190  const std::string relaxType = List.get<std::string> ("relaxation: type");
191 
192  if (relaxType == "Jacobi") {
193  PrecType_ = Ifpack2::Details::JACOBI;
194  }
195  else if (relaxType == "MT Split Jacobi") {
196  PrecType_ = Ifpack2::Details::MTSPLITJACOBI;
197  }
198  else if (relaxType == "Gauss-Seidel") {
199  PrecType_ = Ifpack2::Details::GS;
200  }
201  else if (relaxType == "Symmetric Gauss-Seidel") {
202  PrecType_ = Ifpack2::Details::SGS;
203  }
204  else {
206  (true, std::invalid_argument, "Ifpack2::BlockRelaxation::"
207  "setParameters: Invalid parameter value \"" << relaxType
208  << "\" for parameter \"relaxation: type\".");
209  }
210  }
211 
212  if (List.isParameter ("relaxation: sweeps")) {
213  NumSweeps_ = List.get<int> ("relaxation: sweeps");
214  }
215 
216  // Users may specify this as a floating-point literal. In that
217  // case, it may have type double, even if scalar_type is something
218  // else. We take care to try the various types that make sense.
219  if (List.isParameter ("relaxation: damping factor")) {
220  if (List.isType<double> ("relaxation: damping factor")) {
221  const double dampingFactor =
222  List.get<double> ("relaxation: damping factor");
223  DampingFactor_ = static_cast<scalar_type> (dampingFactor);
224  }
225  else if (List.isType<scalar_type> ("relaxation: damping factor")) {
226  DampingFactor_ = List.get<scalar_type> ("relaxation: damping factor");
227  }
228  else if (List.isType<magnitude_type> ("relaxation: damping factor")) {
229  const magnitude_type dampingFactor =
230  List.get<magnitude_type> ("relaxation: damping factor");
231  DampingFactor_ = static_cast<scalar_type> (dampingFactor);
232  }
233  else {
235  (true, std::invalid_argument, "Ifpack2::BlockRelaxation::"
236  "setParameters: Parameter \"relaxation: damping factor\" "
237  "has the wrong type.");
238  }
239  }
240 
241  if (List.isParameter ("relaxation: zero starting solution")) {
242  ZeroStartingSolution_ = List.get<bool> ("relaxation: zero starting solution");
243  }
244 
245  if (List.isParameter ("relaxation: backward mode")) {
246  DoBackwardGS_ = List.get<bool> ("relaxation: backward mode");
247  }
248 
249  if (List.isParameter ("partitioner: type")) {
250  PartitionerType_ = List.get<std::string> ("partitioner: type");
251  }
252 
253  // Users may specify this as an int literal, so we need to allow
254  // both int and local_ordinal_type (not necessarily same as int).
255  if (List.isParameter ("partitioner: local parts")) {
256  if (List.isType<local_ordinal_type> ("partitioner: local parts")) {
257  NumLocalBlocks_ = List.get<local_ordinal_type> ("partitioner: local parts");
258  }
259  else if (! std::is_same<int, local_ordinal_type>::value &&
260  List.isType<int> ("partitioner: local parts")) {
261  NumLocalBlocks_ = List.get<int> ("partitioner: local parts");
262  }
263  else {
265  (true, std::invalid_argument, "Ifpack2::BlockRelaxation::"
266  "setParameters: Parameter \"partitioner: local parts\" "
267  "has the wrong type.");
268  }
269  }
270 
271  if (List.isParameter ("partitioner: overlap level")) {
272  if (List.isType<int> ("partitioner: overlap level")) {
273  OverlapLevel_ = List.get<int> ("partitioner: overlap level");
274  }
275  else if (! std::is_same<int, local_ordinal_type>::value &&
276  List.isType<local_ordinal_type> ("partitioner: overlap level")) {
277  OverlapLevel_ = List.get<local_ordinal_type> ("partitioner: overlap level");
278  }
279  else {
281  (true, std::invalid_argument, "Ifpack2::BlockRelaxation::"
282  "setParameters: Parameter \"partitioner: overlap level\" "
283  "has the wrong type.");
284  }
285  }
286 
287  std::string defaultContainer = "TriDi";
288  if(std::is_same<ContainerType, Container<MatrixType> >::value)
289  {
290  //Generic container template parameter, container type specified in List
291  Ifpack2::getParameter(List, "relaxation: container", defaultContainer);
292  }
293  // check parameters
294  if (PrecType_ != Ifpack2::Details::JACOBI) {
295  OverlapLevel_ = 0;
296  }
297  if (NumLocalBlocks_ < static_cast<local_ordinal_type> (0)) {
298  NumLocalBlocks_ = A_->getNodeNumRows() / (-NumLocalBlocks_);
299  }
300  // other checks are performed in Partitioner_
301 
302  // NTS: Sanity check to be removed at a later date when Backward mode is enabled
304  DoBackwardGS_, std::runtime_error,
305  "Ifpack2::BlockRelaxation:setParameters: Setting the \"relaxation: "
306  "backward mode\" parameter to true is not yet supported.");
307 
308  // copy the list as each subblock's constructor will
309  // require it later
310  List_ = List;
311 }
312 
313 template<class MatrixType,class ContainerType>
316 {
318  (A_.is_null (), std::runtime_error, "Ifpack2::BlockRelaxation::getComm: "
319  "The matrix is null. You must call setMatrix() with a nonnull matrix "
320  "before you may call this method.");
321  return A_->getComm ();
322 }
323 
324 template<class MatrixType,class ContainerType>
325 Teuchos::RCP<const Tpetra::RowMatrix<typename MatrixType::scalar_type,
326  typename MatrixType::local_ordinal_type,
327  typename MatrixType::global_ordinal_type,
328  typename MatrixType::node_type> >
330  return A_;
331 }
332 
333 template<class MatrixType,class ContainerType>
334 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,
335  typename MatrixType::global_ordinal_type,
336  typename MatrixType::node_type> >
338 getDomainMap () const
339 {
341  (A_.is_null (), std::runtime_error, "Ifpack2::BlockRelaxation::"
342  "getDomainMap: The matrix is null. You must call setMatrix() with a "
343  "nonnull matrix before you may call this method.");
344  return A_->getDomainMap ();
345 }
346 
347 template<class MatrixType,class ContainerType>
348 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,
349  typename MatrixType::global_ordinal_type,
350  typename MatrixType::node_type> >
352 getRangeMap () const
353 {
355  (A_.is_null (), std::runtime_error, "Ifpack2::BlockRelaxation::"
356  "getRangeMap: The matrix is null. You must call setMatrix() with a "
357  "nonnull matrix before you may call this method.");
358  return A_->getRangeMap ();
359 }
360 
361 template<class MatrixType,class ContainerType>
362 bool
364 hasTransposeApply () const {
365  return true;
366 }
367 
368 template<class MatrixType,class ContainerType>
369 int
372  return NumInitialize_;
373 }
374 
375 template<class MatrixType,class ContainerType>
376 int
379 {
380  return NumCompute_;
381 }
382 
383 template<class MatrixType,class ContainerType>
384 int
386 getNumApply () const
387 {
388  return NumApply_;
389 }
390 
391 template<class MatrixType,class ContainerType>
392 double
395 {
396  return InitializeTime_;
397 }
398 
399 template<class MatrixType,class ContainerType>
400 double
403 {
404  return ComputeTime_;
405 }
406 
407 template<class MatrixType,class ContainerType>
408 double
410 getApplyTime () const
411 {
412  return ApplyTime_;
413 }
414 
415 
416 template<class MatrixType,class ContainerType>
419  A_.is_null (), std::runtime_error, "Ifpack2::BlockRelaxation::getNodeSmootherComplexity: "
420  "The input matrix A is null. Please call setMatrix() with a nonnull "
421  "input matrix, then call compute(), before calling this method.");
422  // Relaxation methods cost roughly one apply + one block-diagonal inverse per iteration
423  // NOTE: This approximates all blocks as dense, which may overstate the cost if you have a sparse (or banded) container.
424  size_t block_nnz = 0;
425  for (local_ordinal_type i = 0; i < NumLocalBlocks_; ++i)
426  block_nnz += Partitioner_->numRowsInPart(i) *Partitioner_->numRowsInPart(i);
427 
428  return block_nnz + A_->getNodeNumEntries();
429 }
430 
431 template<class MatrixType,class ContainerType>
432 void
434 apply (const Tpetra::MultiVector<typename MatrixType::scalar_type,
435  typename MatrixType::local_ordinal_type,
436  typename MatrixType::global_ordinal_type,
437  typename MatrixType::node_type>& X,
438  Tpetra::MultiVector<typename MatrixType::scalar_type,
439  typename MatrixType::local_ordinal_type,
440  typename MatrixType::global_ordinal_type,
441  typename MatrixType::node_type>& Y,
442  Teuchos::ETransp mode,
443  scalar_type alpha,
444  scalar_type beta) const
445 {
447  (A_.is_null (), std::runtime_error, "Ifpack2::BlockRelaxation::apply: "
448  "The matrix is null. You must call setMatrix() with a nonnull matrix, "
449  "then call initialize() and compute() (in that order), before you may "
450  "call this method.");
452  ! isComputed (), std::runtime_error, "Ifpack2::BlockRelaxation::apply: "
453  "isComputed() must be true prior to calling apply.");
455  X.getNumVectors () != Y.getNumVectors (), std::invalid_argument,
456  "Ifpack2::BlockRelaxation::apply: X.getNumVectors() = "
457  << X.getNumVectors () << " != Y.getNumVectors() = "
458  << Y.getNumVectors () << ".");
460  mode != Teuchos::NO_TRANS, std::logic_error, "Ifpack2::BlockRelaxation::"
461  "apply: This method currently only implements the case mode == Teuchos::"
462  "NO_TRANS. Transposed modes are not currently supported.");
464  alpha != Teuchos::ScalarTraits<scalar_type>::one (), std::logic_error,
465  "Ifpack2::BlockRelaxation::apply: This method currently only implements "
466  "the case alpha == 1. You specified alpha = " << alpha << ".");
468  beta != Teuchos::ScalarTraits<scalar_type>::zero (), std::logic_error,
469  "Ifpack2::BlockRelaxation::apply: This method currently only implements "
470  "the case beta == 0. You specified beta = " << beta << ".");
471 
472  Time_->start(true);
473 
474  // If X and Y are pointing to the same memory location,
475  // we need to create an auxiliary vector, Xcopy
476  Teuchos::RCP<const MV> X_copy;
477  {
478  auto X_lcl_host = X.template getLocalView<Kokkos::HostSpace> ();
479  auto Y_lcl_host = Y.template getLocalView<Kokkos::HostSpace> ();
480  if (X_lcl_host.data () == Y_lcl_host.data ()) {
481  X_copy = rcp (new MV (X, Teuchos::Copy));
482  } else {
483  X_copy = rcpFromRef (X);
484  }
485  }
486 
487  if (ZeroStartingSolution_) {
488  Y.putScalar (STS::zero ());
489  }
490 
491  // Flops are updated in each of the following.
492  switch (PrecType_) {
493  case Ifpack2::Details::JACOBI:
494  ApplyInverseJacobi(*X_copy,Y);
495  break;
496  case Ifpack2::Details::GS:
497  ApplyInverseGS(*X_copy,Y);
498  break;
499  case Ifpack2::Details::SGS:
500  ApplyInverseSGS(*X_copy,Y);
501  break;
502  case Ifpack2::Details::MTSPLITJACOBI:
503  Container_->applyInverseJacobi(*X_copy, Y, ZeroStartingSolution_, NumSweeps_);
504  break;
505  default:
507  (true, std::logic_error, "Ifpack2::BlockRelaxation::apply: Invalid "
508  "PrecType_ enum value " << PrecType_ << ". Valid values are Ifpack2::"
509  "Details::JACOBI = " << Ifpack2::Details::JACOBI << ", Ifpack2::Details"
510  "::GS = " << Ifpack2::Details::GS << ", and Ifpack2::Details::SGS = "
511  << Ifpack2::Details::SGS << ". Please report this bug to the Ifpack2 "
512  "developers.");
513  }
514 
515  ++NumApply_;
516  Time_->stop();
517  ApplyTime_ += Time_->totalElapsedTime();
518 }
519 
520 template<class MatrixType,class ContainerType>
521 void
523 applyMat (const Tpetra::MultiVector<typename MatrixType::scalar_type,
524  typename MatrixType::local_ordinal_type,
525  typename MatrixType::global_ordinal_type,
526  typename MatrixType::node_type>& X,
527  Tpetra::MultiVector<typename MatrixType::scalar_type,
528  typename MatrixType::local_ordinal_type,
529  typename MatrixType::global_ordinal_type,
530  typename MatrixType::node_type>& Y,
531  Teuchos::ETransp mode) const
532 {
533  A_->apply (X, Y, mode);
534 }
535 
536 template<class MatrixType,class ContainerType>
537 void
540 {
541  using Teuchos::rcp;
542  typedef Tpetra::RowGraph<local_ordinal_type, global_ordinal_type, node_type>
543  row_graph_type;
544 
546  (A_.is_null (), std::runtime_error, "Ifpack2::BlockRelaxation::initialize: "
547  "The matrix is null. You must call setMatrix() with a nonnull matrix "
548  "before you may call this method.");
549 
550  // Check whether we have a BlockCrsMatrix
552  Teuchos::rcp_dynamic_cast<const block_crs_matrix_type> (A_);
553  hasBlockCrsMatrix_ = !A_bcrs.is_null();
554  if (A_bcrs.is_null ()) {
555  hasBlockCrsMatrix_ = false;
556  }
557  else {
558  hasBlockCrsMatrix_ = true;
559  }
560 
561  IsInitialized_ = false;
562  Time_->start (true);
563 
564  NumLocalRows_ = A_->getNodeNumRows ();
565  NumGlobalRows_ = A_->getGlobalNumRows ();
566  NumGlobalNonzeros_ = A_->getGlobalNumEntries ();
567 
568  // NTS: Will need to add support for Zoltan2 partitions later Also,
569  // will need a better way of handling the Graph typing issue.
570  // Especially with ordinal types w.r.t the container.
571  Partitioner_ = Teuchos::null;
572 
573  if (PartitionerType_ == "linear") {
574  Partitioner_ =
575  rcp (new Ifpack2::LinearPartitioner<row_graph_type> (A_->getGraph ()));
576  } else if (PartitionerType_ == "line") {
577  Partitioner_ =
579  } else if (PartitionerType_ == "user") {
580  Partitioner_ =
581  rcp (new Ifpack2::Details::UserPartitioner<row_graph_type> (A_->getGraph () ) );
582  } else {
583  // We should have checked for this in setParameters(), so it's a
584  // logic_error, not an invalid_argument or runtime_error.
586  (true, std::logic_error, "Ifpack2::BlockRelaxation::initialize: Unknown "
587  "partitioner type " << PartitionerType_ << ". Valid values are "
588  "\"linear\", \"line\", and \"user\".");
589  }
590 
591  // need to partition the graph of A
592  Partitioner_->setParameters (List_);
593  Partitioner_->compute ();
594 
595  // get actual number of partitions
596  NumLocalBlocks_ = Partitioner_->numLocalParts ();
597 
598  // Note: Unlike Ifpack, we'll punt on computing W_ until compute(), which is where
599  // we assume that the type of relaxation has been chosen.
600 
601  if (A_->getComm()->getSize() != 1) {
602  IsParallel_ = true;
603  } else {
604  IsParallel_ = false;
605  }
606 
607  // We should have checked for this in setParameters(), so it's a
608  // logic_error, not an invalid_argument or runtime_error.
610  (NumSweeps_ < 0, std::logic_error, "Ifpack2::BlockRelaxation::initialize: "
611  "NumSweeps_ = " << NumSweeps_ << " < 0.");
612 
613  // Extract the submatrices
614  ExtractSubmatricesStructure ();
615 
616  // Compute the weight vector if we're doing overlapped Jacobi (and
617  // only if we're doing overlapped Jacobi).
618  if (PrecType_ == Ifpack2::Details::JACOBI && OverlapLevel_ > 0) {
620  (hasBlockCrsMatrix_, std::runtime_error,
621  "Ifpack2::BlockRelaxation::initialize: "
622  "We do not support overlapped Jacobi yet for Tpetra::BlockCrsMatrix. Sorry!");
623 
624  // weight of each vertex
625  W_ = rcp (new vector_type (A_->getRowMap ()));
626  W_->putScalar (STS::zero ());
627  Teuchos::ArrayRCP<scalar_type > w_ptr = W_->getDataNonConst(0);
628 
629  for (local_ordinal_type i = 0 ; i < NumLocalBlocks_ ; ++i) {
630  for (size_t j = 0 ; j < Partitioner_->numRowsInPart(i) ; ++j) {
631  // FIXME (mfh 12 Sep 2014) Should this really be int?
632  // Perhaps it should be local_ordinal_type instead.
633  int LID = (*Partitioner_)(i,j);
634  w_ptr[LID]+= STS::one();
635  }
636  }
637  W_->reciprocal (*W_);
638  }
639 
640  ++NumInitialize_;
641  Time_->stop ();
642  InitializeTime_ += Time_->totalElapsedTime ();
643  IsInitialized_ = true;
644 }
645 
646 
647 template<class MatrixType,class ContainerType>
648 void
651 {
652  using Teuchos::rcp;
653 
655  (A_.is_null (), std::runtime_error, "Ifpack2::BlockRelaxation::compute: "
656  "The matrix is null. You must call setMatrix() with a nonnull matrix, "
657  "then call initialize(), before you may call this method.");
658 
659  if (! isInitialized ()) {
660  initialize ();
661  }
662 
663  Time_->start (true);
664 
665  // reset values
666  IsComputed_ = false;
667 
668  Container_->compute(); // compute each block matrix
669 
670  ++NumCompute_;
671  Time_->stop ();
672  ComputeTime_ += Time_->totalElapsedTime();
673  IsComputed_ = true;
674 }
675 
676 template<class MatrixType, class ContainerType>
677 void
680 {
682  (Partitioner_.is_null (), std::runtime_error, "Ifpack2::BlockRelaxation::"
683  "ExtractSubmatricesStructure: Partitioner object is null.");
684 
685  NumLocalBlocks_ = Partitioner_->numLocalParts ();
686  std::string containerType = ContainerType::getName ();
687  if (containerType == "Generic") {
688  // ContainerType is Container<row_matrix_type>. Thus, we need to
689  // get the container name from the parameter list. We use "TriDi"
690  // as the default container type.
691  containerType = containerType_;
692  }
693 
694  Teuchos::Array<Teuchos::Array<local_ordinal_type> > localRows(NumLocalBlocks_);
695  for (local_ordinal_type i = 0; i < NumLocalBlocks_; ++i) {
696  const size_t numRows = Partitioner_->numRowsInPart (i);
697 
698  localRows[i].resize(numRows);
699  // Extract a list of the indices of each partitioner row.
700  for (size_t j = 0; j < numRows; ++j) {
701  localRows[i][j] = (*Partitioner_) (i,j);
702  }
703  }
704  Container_ = ContainerFactory<MatrixType>::build(containerType, A_, localRows, Importer_, OverlapLevel_, DampingFactor_);
705  Container_->setParameters(List_);
706  Container_->initialize();
707 }
708 
709 template<class MatrixType,class ContainerType>
710 void
711 BlockRelaxation<MatrixType,ContainerType>::
712 ApplyInverseJacobi (const MV& X, MV& Y) const
713 {
714  const size_t NumVectors = X.getNumVectors ();
715  typename ContainerType::HostView XView = X.template getLocalView<Kokkos::HostSpace>();
716  typename ContainerType::HostView YView = Y.template getLocalView<Kokkos::HostSpace>();
717  MV AY (Y.getMap (), NumVectors);
718 
719  typename ContainerType::HostView AYView = AY.template getLocalView<Kokkos::HostSpace>();
720  // Initial matvec not needed
721  int starting_iteration = 0;
722  if (OverlapLevel_ > 0)
723  {
724  //Overlapping jacobi, with view of W_
725  typename ContainerType::HostView WView = W_->template getLocalView<Kokkos::HostSpace>();
726  if(ZeroStartingSolution_) {
727  Container_->DoOverlappingJacobi(XView, YView, WView, X.getStride());
728  starting_iteration = 1;
729  }
730  const scalar_type ONE = STS::one();
731  for(int j = starting_iteration; j < NumSweeps_; j++)
732  {
733  applyMat (Y, AY);
734  AY.update (ONE, X, -ONE);
735  Container_->DoOverlappingJacobi (AYView, YView, WView, X.getStride());
736  }
737  }
738  else
739  {
740  //Non-overlapping
741  if(ZeroStartingSolution_)
742  {
743  Container_->DoJacobi (XView, YView, X.getStride());
744  starting_iteration = 1;
745  }
746  const scalar_type ONE = STS::one();
747  for(int j = starting_iteration; j < NumSweeps_; j++)
748  {
749  applyMat (Y, AY);
750  AY.update (ONE, X, -ONE);
751  Container_->DoJacobi (AYView, YView, X.getStride());
752  }
753  }
754 }
755 
756 template<class MatrixType,class ContainerType>
757 void
758 BlockRelaxation<MatrixType,ContainerType>::
759 ApplyInverseGS (const MV& X, MV& Y) const
760 {
761  using Teuchos::Ptr;
762  using Teuchos::ptr;
763  size_t numVecs = X.getNumVectors();
764  //Get view of X (is never modified in this function)
765  typename ContainerType::HostView XView = X.template getLocalView<Kokkos::HostSpace>();
766  typename ContainerType::HostView YView = Y.template getLocalView<Kokkos::HostSpace>();
767  //Pre-import Y, if parallel
768  Ptr<MV> Y2;
769  bool deleteY2 = false;
770  if(IsParallel_)
771  {
772  Y2 = ptr(new MV(Importer_->getTargetMap(), numVecs));
773  deleteY2 = true;
774  }
775  else
776  Y2 = ptr(&Y);
777  if(IsParallel_)
778  {
779  for(int j = 0; j < NumSweeps_; ++j)
780  {
781  //do import once per sweep
782  Y2->doImport(Y, *Importer_, Tpetra::INSERT);
783  typename ContainerType::HostView Y2View = Y2->template getLocalView<Kokkos::HostSpace>();
784  Container_->DoGaussSeidel(XView, YView, Y2View, X.getStride());
785  }
786  }
787  else
788  {
789  typename ContainerType::HostView Y2View = Y2->template getLocalView<Kokkos::HostSpace>();
790  for(int j = 0; j < NumSweeps_; ++j)
791  {
792  Container_->DoGaussSeidel(XView, YView, Y2View, X.getStride());
793  }
794  }
795  if(deleteY2)
796  delete Y2.get();
797 }
798 
799 template<class MatrixType,class ContainerType>
800 void
801 BlockRelaxation<MatrixType,ContainerType>::
802 ApplyInverseSGS (const MV& X, MV& Y) const
803 {
804  using Teuchos::Ptr;
805  using Teuchos::ptr;
806  //Get view of X (is never modified in this function)
807  typename ContainerType::HostView XView = X.template getLocalView<Kokkos::HostSpace>();
808  typename ContainerType::HostView YView = Y.template getLocalView<Kokkos::HostSpace>();
809  //Pre-import Y, if parallel
810  Ptr<MV> Y2;
811  bool deleteY2 = false;
812  if(IsParallel_)
813  {
814  Y2 = ptr(new MV(Importer_->getTargetMap(), X.getNumVectors()));
815  deleteY2 = true;
816  }
817  else
818  Y2 = ptr(&Y);
819  if(IsParallel_)
820  {
821  for(int j = 0; j < NumSweeps_; ++j)
822  {
823  //do import once per sweep
824  Y2->doImport(Y, *Importer_, Tpetra::INSERT);
825  typename ContainerType::HostView Y2View = Y2->template getLocalView<Kokkos::HostSpace>();
826  Container_->DoSGS(XView, YView, Y2View, X.getStride());
827  }
828  }
829  else
830  {
831  typename ContainerType::HostView Y2View = Y2->template getLocalView<Kokkos::HostSpace>();
832  for(int j = 0; j < NumSweeps_; ++j)
833  {
834  Container_->DoSGS(XView, YView, Y2View, X.getStride());
835  }
836  }
837  if(deleteY2)
838  delete Y2.get();
839 }
840 
841 template<class MatrixType,class ContainerType>
842 void BlockRelaxation<MatrixType,ContainerType>::computeImporter () const
843 {
844  using Teuchos::RCP;
845  using Teuchos::rcp;
846  using Teuchos::Ptr;
847  using Teuchos::ArrayView;
848  using Teuchos::Array;
849  using Teuchos::Comm;
850  using Teuchos::rcp_dynamic_cast;
851  if(IsParallel_)
852  {
853  if(hasBlockCrsMatrix_)
854  {
855  const RCP<const block_crs_matrix_type> bcrs =
856  rcp_dynamic_cast<const block_crs_matrix_type>(A_);
857  int bs_ = bcrs->getBlockSize();
858  RCP<const map_type> oldDomainMap = A_->getDomainMap();
859  RCP<const map_type> oldColMap = A_->getColMap();
860  // Because A is a block CRS matrix, import will not do what you think it does
861  // We have to construct the correct maps for it
862  global_size_t numGlobalElements = oldColMap->getGlobalNumElements() * bs_;
863  global_ordinal_type indexBase = oldColMap->getIndexBase();
864  RCP<const Comm<int>> comm = oldColMap->getComm();
865  ArrayView<const global_ordinal_type> oldColElements = oldColMap->getNodeElementList();
866  Array<global_ordinal_type> newColElements(bs_ * oldColElements.size());
867  for(int i = 0; i < oldColElements.size(); i++)
868  {
869  for(int j = 0; j < bs_; j++)
870  newColElements[i * bs_ + j] = oldColElements[i] * bs_ + j;
871  }
872  RCP<map_type> colMap(new map_type(numGlobalElements, newColElements, indexBase, comm));
873  // Create the importer
874  Importer_ = rcp(new import_type(oldDomainMap, colMap));
875  }
876  else if(!A_.is_null())
877  {
878  Importer_ = A_->getGraph()->getImporter();
879  if(Importer_.is_null())
880  Importer_ = rcp(new import_type(A_->getDomainMap(), A_->getColMap()));
881  }
882  }
883  //otherwise, Importer_ is not needed and is left NULL
884 }
885 
886 template<class MatrixType, class ContainerType>
887 std::string
889 description () const
890 {
891  std::ostringstream out;
892 
893  // Output is a valid YAML dictionary in flow style. If you don't
894  // like everything on a single line, you should call describe()
895  // instead.
896  out << "\"Ifpack2::BlockRelaxation\": {";
897  if (this->getObjectLabel () != "") {
898  out << "Label: \"" << this->getObjectLabel () << "\", ";
899  }
900  out << "Initialized: " << (isInitialized () ? "true" : "false") << ", ";
901  out << "Computed: " << (isComputed () ? "true" : "false") << ", ";
902  if (A_.is_null ()) {
903  out << "Matrix: null, ";
904  }
905  else {
906  out << "Matrix: not null"
907  << ", Global matrix dimensions: ["
908  << A_->getGlobalNumRows () << ", " << A_->getGlobalNumCols () << "], ";
909  }
910 
911  // It's useful to print this instance's relaxation method. If you
912  // want more info than that, call describe() instead.
913  out << "\"relaxation: type\": ";
914  if (PrecType_ == Ifpack2::Details::JACOBI) {
915  out << "Block Jacobi";
916  } else if (PrecType_ == Ifpack2::Details::GS) {
917  out << "Block Gauss-Seidel";
918  } else if (PrecType_ == Ifpack2::Details::SGS) {
919  out << "Block Symmetric Gauss-Seidel";
920  } else {
921  out << "INVALID";
922  }
923 
924  out << ", overlap: " << OverlapLevel_;
925 
926  out << ", " << "sweeps: " << NumSweeps_ << ", "
927  << "damping factor: " << DampingFactor_ << ", ";
928 
929  std::string containerType = ContainerType::getName();
930  out << "block container: " << ((containerType == "Generic") ? containerType_ : containerType);
931 
932  out << "}";
933  return out.str();
934 }
935 
936 template<class MatrixType,class ContainerType>
937 void
940  const Teuchos::EVerbosityLevel verbLevel) const
941 {
942  using std::endl;
943  using std::setw;
945  using Teuchos::VERB_DEFAULT;
946  using Teuchos::VERB_NONE;
947  using Teuchos::VERB_LOW;
948  using Teuchos::VERB_MEDIUM;
949  using Teuchos::VERB_HIGH;
950  using Teuchos::VERB_EXTREME;
951 
952  Teuchos::EVerbosityLevel vl = verbLevel;
953  if (vl == VERB_DEFAULT) vl = VERB_LOW;
954  const int myImageID = A_->getComm()->getRank();
955 
956  // Convention requires that describe() begin with a tab.
957  Teuchos::OSTab tab (out);
958 
959  // none: print nothing
960  // low: print O(1) info from node 0
961  // medium:
962  // high:
963  // extreme:
964  if (vl != VERB_NONE && myImageID == 0) {
965  out << "Ifpack2::BlockRelaxation<"
966  << TypeNameTraits<MatrixType>::name () << ", "
967  << TypeNameTraits<ContainerType>::name () << " >:";
968  Teuchos::OSTab tab1 (out);
969 
970  if (this->getObjectLabel () != "") {
971  out << "label: \"" << this->getObjectLabel () << "\"" << endl;
972  }
973  out << "initialized: " << (isInitialized () ? "true" : "false") << endl
974  << "computed: " << (isComputed () ? "true" : "false") << endl;
975 
976  std::string precType;
977  if (PrecType_ == Ifpack2::Details::JACOBI) {
978  precType = "Block Jacobi";
979  } else if (PrecType_ == Ifpack2::Details::GS) {
980  precType = "Block Gauss-Seidel";
981  } else if (PrecType_ == Ifpack2::Details::SGS) {
982  precType = "Block Symmetric Gauss-Seidel";
983  }
984  out << "type: " << precType << endl
985  << "global number of rows: " << A_->getGlobalNumRows () << endl
986  << "global number of columns: " << A_->getGlobalNumCols () << endl
987  << "number of sweeps: " << NumSweeps_ << endl
988  << "damping factor: " << DampingFactor_ << endl
989  << "backwards mode: "
990  << ((PrecType_ == Ifpack2::Details::GS && DoBackwardGS_) ? "true" : "false")
991  << endl
992  << "zero starting solution: "
993  << (ZeroStartingSolution_ ? "true" : "false") << endl;
994  }
995 }
996 
997 }//namespace Ifpack2
998 
999 
1000 // Macro that does explicit template instantiation (ETI) for
1001 // Ifpack2::BlockRelaxation. S, LO, GO, N correspond to the four
1002 // template parameters of Ifpack2::Preconditioner and
1003 // Tpetra::RowMatrix.
1004 //
1005 // We only instantiate for MatrixType = Tpetra::RowMatrix. There's no
1006 // need to instantiate for Tpetra::CrsMatrix too. All Ifpack2
1007 // preconditioners can and should do dynamic casts if they need a type
1008 // more specific than RowMatrix. This keeps build time short and
1009 // library and executable sizes small.
1010 
1011 #ifdef HAVE_IFPACK2_EXPLICIT_INSTANTIATION
1012 
1013 #define IFPACK2_BLOCKRELAXATION_INSTANT(S,LO,GO,N) \
1014  template \
1015  class Ifpack2::BlockRelaxation< \
1016  Tpetra::RowMatrix<S, LO, GO, N>, \
1017  Ifpack2::Container<Tpetra::RowMatrix<S, LO, GO, N> > >;
1018 
1019 #endif // HAVE_IFPACK2_EXPLICIT_INSTANTIATION
1020 
1021 #endif // IFPACK2_BLOCKRELAXATION_DEF_HPP
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
The communicator over which the input matrix is distributed.
Definition: Ifpack2_BlockRelaxation_def.hpp:315
Ifpack2::BlockRelaxation class declaration.
double getApplyTime() const
Returns the time spent in apply().
Definition: Ifpack2_BlockRelaxation_def.hpp:410
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
Definition: Ifpack2_BlockRelaxation_def.hpp:417
VERB_DEFAULT
ParameterList & disableRecursiveValidation()
double getComputeTime() const
Returns the time spent in compute().
Definition: Ifpack2_BlockRelaxation_def.hpp:402
int getNumInitialize() const
Returns the number of calls to initialize().
Definition: Ifpack2_BlockRelaxation_def.hpp:371
VERB_LOW
VERB_NONE
T & get(const std::string &name, T def_value)
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const
Print the object with some verbosity level to an FancyOStream object.
Definition: Ifpack2_BlockRelaxation_def.hpp:939
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
double getInitializeTime() const
Returns the time spent in initialize().
Definition: Ifpack2_BlockRelaxation_def.hpp:394
static Teuchos::RCP< BaseContainer > build(std::string containerType, const Teuchos::RCP< const MatrixType > &A, const Teuchos::Array< Teuchos::Array< local_ordinal_type >> &localRows, const Teuchos::RCP< const import_type > importer, int OverlapLevel, scalar_type DampingFactor)
Build a specialization of Ifpack2::Container given a key that has been registered.
Definition: Ifpack2_ContainerFactory_def.hpp:89
Teuchos::RCP< const row_matrix_type > getMatrix() const
The input matrix of this preconditioner&#39;s constructor.
Definition: Ifpack2_BlockRelaxation_def.hpp:329
Teuchos::RCP< const map_type > getDomainMap() const
Returns the Tpetra::Map object associated with the domain of this operator.
Definition: Ifpack2_BlockRelaxation_def.hpp:338
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
VERB_EXTREME
void getParameter(const Teuchos::ParameterList &params, const std::string &name, T &value)
Set a value from a ParameterList if a parameter with the specified name exists.
Definition: Ifpack2_Parameters.hpp:59
void apply(const MV &X, MV &Y, Teuchos::ETransp mode=Teuchos::NO_TRANS, scalar_type alpha=Teuchos::ScalarTraits< scalar_type >::one(), scalar_type beta=Teuchos::ScalarTraits< scalar_type >::zero()) const
Applies the preconditioner to X, returns the result in Y.
Definition: Ifpack2_BlockRelaxation_def.hpp:434
std::string description() const
A one-line description of this object.
Definition: Ifpack2_BlockRelaxation_def.hpp:889
Teuchos::RCP< const map_type > getRangeMap() const
Returns the Tpetra::Map object associated with the range of this operator.
Definition: Ifpack2_BlockRelaxation_def.hpp:352
void setParameters(const Teuchos::ParameterList &params)
Sets all the parameters for the preconditioner.
Definition: Ifpack2_BlockRelaxation_def.hpp:166
bool is_null() const
void compute()
compute the preconditioner for the specified matrix, diagonal perturbation thresholds and relaxation ...
Definition: Ifpack2_BlockRelaxation_def.hpp:650
VERB_MEDIUM
Block relaxation preconditioners (or smoothers) for Tpetra::RowMatrix and Tpetra::CrsMatrix sparse ma...
Definition: Ifpack2_BlockRelaxation_decl.hpp:83
Ifpack2 implementation details.
Declaration of a user-defined partitioner in which the user can define a partition of the graph...
void validateParameters(ParameterList const &validParamList, int const depth=1000, EValidateUsed const validateUsed=VALIDATE_USED_ENABLED, EValidateDefaults const validateDefaults=VALIDATE_DEFAULTS_ENABLED) const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
VERB_HIGH
void applyMat(const MV &X, MV &Y, Teuchos::ETransp mode=Teuchos::NO_TRANS) const
Applies the matrix to a Tpetra::MultiVector.
Definition: Ifpack2_BlockRelaxation_def.hpp:523
Teuchos::ScalarTraits< scalar_type >::magnitudeType magnitude_type
The type of the magnitude (absolute value) of a matrix entry.
Definition: Ifpack2_BlockRelaxation_decl.hpp:109
virtual void setMatrix(const Teuchos::RCP< const row_matrix_type > &A)
Change the matrix to be preconditioned.
Definition: Ifpack2_BlockRelaxation_def.hpp:57
virtual ~BlockRelaxation()
Destructor.
Definition: Ifpack2_BlockRelaxation_def.hpp:116
bool isType(const std::string &name) const
Partition in which the user can define a nonoverlapping partition of the graph in any way they choose...
Definition: Ifpack2_Details_UserPartitioner_decl.hpp:69
void initialize()
Initialize.
Definition: Ifpack2_BlockRelaxation_def.hpp:539
MatrixType::local_ordinal_type local_ordinal_type
The type of local indices in the input MatrixType.
Definition: Ifpack2_BlockRelaxation_decl.hpp:101
int getNumCompute() const
Returns the number of calls to compute().
Definition: Ifpack2_BlockRelaxation_def.hpp:378
int getNumApply() const
Returns the number of calls to apply().
Definition: Ifpack2_BlockRelaxation_def.hpp:386
MatrixType::scalar_type scalar_type
The type of the entries of the input MatrixType.
Definition: Ifpack2_BlockRelaxation_decl.hpp:98
Interface for creating and solving a local linear problem.
Definition: Ifpack2_Container.hpp:114
T * getRawPtr() const
bool isParameter(const std::string &name) const
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:72
void getValidParameters(Teuchos::ParameterList &params)
Fills a list which contains all the parameters possibly used by Ifpack2.
Definition: Ifpack2_Parameters.cpp:50
Ifpack2::LinePartitioner: A class to define partitions into a set of lines.
Definition: Ifpack2_LinePartitioner_decl.hpp:77
A class to define linear partitions.
Definition: Ifpack2_LinearPartitioner_decl.hpp:60
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Return a list of all the parameters that this class accepts.
Definition: Ifpack2_BlockRelaxation_def.hpp:122
BlockRelaxation(const Teuchos::RCP< const row_matrix_type > &Matrix)
Constructor.
Definition: Ifpack2_BlockRelaxation_def.hpp:83