Ifpack2 Templated Preconditioning Package  Version 1.0
Ifpack2_BlockRelaxation_def.hpp
1 /*@HEADER
2 // ***********************************************************************
3 //
4 // Ifpack2: Tempated 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_ ("Dense"),
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>
120 void
123 {
124  //Teuchos::ParameterList validparams;
125  //Ifpack2::getValidParameters (validparams);
126  //List.validateParameters (validparams);
127 
128  // Don't change ANY parameter unless the user specified it
129  // explicitly, thus indicating that the user wants to change it.
130  // This avoids the (not really sensible) defaults in
131  // Ifpack2::getValidParameters.
132 
133  if (List.isParameter ("relaxation: container")) {
134  // If the container type isn't a string, this will throw, but it
135  // rightfully should.
136  containerType_ = List.get<std::string> ("relaxation: container");
137  }
138 
139  if (List.isParameter ("relaxation: type")) {
140  const std::string relaxType = List.get<std::string> ("relaxation: type");
141 
142  if (relaxType == "Jacobi") {
143  PrecType_ = Ifpack2::Details::JACOBI;
144  }
145  else if (relaxType == "Gauss-Seidel") {
146  PrecType_ = Ifpack2::Details::GS;
147  }
148  else if (relaxType == "Symmetric Gauss-Seidel") {
149  PrecType_ = Ifpack2::Details::SGS;
150  }
151  else {
153  (true, std::invalid_argument, "Ifpack2::BlockRelaxation::"
154  "setParameters: Invalid parameter value \"" << relaxType
155  << "\" for parameter \"relaxation: type\".");
156  }
157  }
158 
159  if (List.isParameter ("relaxation: sweeps")) {
160  NumSweeps_ = List.get<int> ("relaxation: sweeps");
161  }
162 
163  // Users may specify this as a floating-point literal. In that
164  // case, it may have type double, even if scalar_type is something
165  // else. We take care to try the various types that make sense.
166  if (List.isParameter ("relaxation: damping factor")) {
167  if (List.isType<double> ("relaxation: damping factor")) {
168  const double dampingFactor =
169  List.get<double> ("relaxation: damping factor");
170  DampingFactor_ = static_cast<scalar_type> (dampingFactor);
171  }
172  else if (List.isType<scalar_type> ("relaxation: damping factor")) {
173  DampingFactor_ = List.get<scalar_type> ("relaxation: damping factor");
174  }
175  else if (List.isType<magnitude_type> ("relaxation: damping factor")) {
176  const magnitude_type dampingFactor =
177  List.get<magnitude_type> ("relaxation: damping factor");
178  DampingFactor_ = static_cast<scalar_type> (dampingFactor);
179  }
180  else {
182  (true, std::invalid_argument, "Ifpack2::BlockRelaxation::"
183  "setParameters: Parameter \"relaxation: damping factor\" "
184  "has the wrong type.");
185  }
186  }
187 
188  if (List.isParameter ("relaxation: zero starting solution")) {
189  ZeroStartingSolution_ = List.get<bool> ("relaxation: zero starting solution");
190  }
191 
192  if (List.isParameter ("relaxation: backward mode")) {
193  DoBackwardGS_ = List.get<bool> ("relaxation: backward mode");
194  }
195 
196  if (List.isParameter ("partitioner: type")) {
197  PartitionerType_ = List.get<std::string> ("partitioner: type");
198  }
199 
200  // Users may specify this as an int literal, so we need to allow
201  // both int and local_ordinal_type (not necessarily same as int).
202  if (List.isParameter ("partitioner: local parts")) {
203  if (List.isType<local_ordinal_type> ("partitioner: local parts")) {
204  NumLocalBlocks_ = List.get<local_ordinal_type> ("partitioner: local parts");
205  }
206  else if (! std::is_same<int, local_ordinal_type>::value &&
207  List.isType<int> ("partitioner: local parts")) {
208  NumLocalBlocks_ = List.get<int> ("partitioner: local parts");
209  }
210  else {
212  (true, std::invalid_argument, "Ifpack2::BlockRelaxation::"
213  "setParameters: Parameter \"partitioner: local parts\" "
214  "has the wrong type.");
215  }
216  }
217 
218  if (List.isParameter ("partitioner: overlap level")) {
219  if (List.isType<int> ("partitioner: overlap level")) {
220  OverlapLevel_ = List.get<int> ("partitioner: overlap level");
221  }
222  else if (! std::is_same<int, local_ordinal_type>::value &&
223  List.isType<local_ordinal_type> ("partitioner: overlap level")) {
224  OverlapLevel_ = List.get<local_ordinal_type> ("partitioner: overlap level");
225  }
226  else {
228  (true, std::invalid_argument, "Ifpack2::BlockRelaxation::"
229  "setParameters: Parameter \"partitioner: overlap level\" "
230  "has the wrong type.");
231  }
232  }
233 
234  std::string defaultContainer = "Dense";
235  if(std::is_same<ContainerType, Container<MatrixType> >::value)
236  {
237  //Generic container template parameter, container type specified in List
238  Ifpack2::getParameter(List, "relaxation: container", defaultContainer);
239  }
240  // check parameters
241  if (PrecType_ != Ifpack2::Details::JACOBI) {
242  OverlapLevel_ = 0;
243  }
244  if (NumLocalBlocks_ < static_cast<local_ordinal_type> (0)) {
245  NumLocalBlocks_ = A_->getNodeNumRows() / (-NumLocalBlocks_);
246  }
247  // other checks are performed in Partitioner_
248 
249  // NTS: Sanity check to be removed at a later date when Backward mode is enabled
251  DoBackwardGS_, std::runtime_error,
252  "Ifpack2::BlockRelaxation:setParameters: Setting the \"relaxation: "
253  "backward mode\" parameter to true is not yet supported.");
254 
255  // copy the list as each subblock's constructor will
256  // require it later
257  List_ = List;
258 }
259 
260 template<class MatrixType,class ContainerType>
263 {
265  (A_.is_null (), std::runtime_error, "Ifpack2::BlockRelaxation::getComm: "
266  "The matrix is null. You must call setMatrix() with a nonnull matrix "
267  "before you may call this method.");
268  return A_->getComm ();
269 }
270 
271 template<class MatrixType,class ContainerType>
272 Teuchos::RCP<const Tpetra::RowMatrix<typename MatrixType::scalar_type,
273  typename MatrixType::local_ordinal_type,
274  typename MatrixType::global_ordinal_type,
275  typename MatrixType::node_type> >
277  return A_;
278 }
279 
280 template<class MatrixType,class ContainerType>
281 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,
282  typename MatrixType::global_ordinal_type,
283  typename MatrixType::node_type> >
285 getDomainMap () const
286 {
288  (A_.is_null (), std::runtime_error, "Ifpack2::BlockRelaxation::"
289  "getDomainMap: The matrix is null. You must call setMatrix() with a "
290  "nonnull matrix before you may call this method.");
291  return A_->getDomainMap ();
292 }
293 
294 template<class MatrixType,class ContainerType>
295 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,
296  typename MatrixType::global_ordinal_type,
297  typename MatrixType::node_type> >
299 getRangeMap () const
300 {
302  (A_.is_null (), std::runtime_error, "Ifpack2::BlockRelaxation::"
303  "getRangeMap: The matrix is null. You must call setMatrix() with a "
304  "nonnull matrix before you may call this method.");
305  return A_->getRangeMap ();
306 }
307 
308 template<class MatrixType,class ContainerType>
309 bool
311 hasTransposeApply () const {
312  return true;
313 }
314 
315 template<class MatrixType,class ContainerType>
316 int
319  return NumInitialize_;
320 }
321 
322 template<class MatrixType,class ContainerType>
323 int
326 {
327  return NumCompute_;
328 }
329 
330 template<class MatrixType,class ContainerType>
331 int
333 getNumApply () const
334 {
335  return NumApply_;
336 }
337 
338 template<class MatrixType,class ContainerType>
339 double
342 {
343  return InitializeTime_;
344 }
345 
346 template<class MatrixType,class ContainerType>
347 double
350 {
351  return ComputeTime_;
352 }
353 
354 template<class MatrixType,class ContainerType>
355 double
357 getApplyTime () const
358 {
359  return ApplyTime_;
360 }
361 
362 template<class MatrixType,class ContainerType>
363 void
365 apply (const Tpetra::MultiVector<typename MatrixType::scalar_type,
366  typename MatrixType::local_ordinal_type,
367  typename MatrixType::global_ordinal_type,
368  typename MatrixType::node_type>& X,
369  Tpetra::MultiVector<typename MatrixType::scalar_type,
370  typename MatrixType::local_ordinal_type,
371  typename MatrixType::global_ordinal_type,
372  typename MatrixType::node_type>& Y,
373  Teuchos::ETransp mode,
374  scalar_type alpha,
375  scalar_type beta) const
376 {
378  (A_.is_null (), std::runtime_error, "Ifpack2::BlockRelaxation::apply: "
379  "The matrix is null. You must call setMatrix() with a nonnull matrix, "
380  "then call initialize() and compute() (in that order), before you may "
381  "call this method.");
383  ! isComputed (), std::runtime_error, "Ifpack2::BlockRelaxation::apply: "
384  "isComputed() must be true prior to calling apply.");
386  X.getNumVectors () != Y.getNumVectors (), std::invalid_argument,
387  "Ifpack2::BlockRelaxation::apply: X.getNumVectors() = "
388  << X.getNumVectors () << " != Y.getNumVectors() = "
389  << Y.getNumVectors () << ".");
391  mode != Teuchos::NO_TRANS, std::logic_error, "Ifpack2::BlockRelaxation::"
392  "apply: This method currently only implements the case mode == Teuchos::"
393  "NO_TRANS. Transposed modes are not currently supported.");
395  alpha != Teuchos::ScalarTraits<scalar_type>::one (), std::logic_error,
396  "Ifpack2::BlockRelaxation::apply: This method currently only implements "
397  "the case alpha == 1. You specified alpha = " << alpha << ".");
399  beta != Teuchos::ScalarTraits<scalar_type>::zero (), std::logic_error,
400  "Ifpack2::BlockRelaxation::apply: This method currently only implements "
401  "the case beta == 0. You specified beta = " << beta << ".");
402 
403  Time_->start(true);
404 
405  // If X and Y are pointing to the same memory location,
406  // we need to create an auxiliary vector, Xcopy
407  Teuchos::RCP<const MV> X_copy;
408  {
409  auto X_lcl_host = X.template getLocalView<Kokkos::HostSpace> ();
410  auto Y_lcl_host = Y.template getLocalView<Kokkos::HostSpace> ();
411  if (X_lcl_host.ptr_on_device () == Y_lcl_host.ptr_on_device ()) {
412  X_copy = rcp (new MV (X, Teuchos::Copy));
413  } else {
414  X_copy = rcpFromRef (X);
415  }
416  }
417 
418  if (ZeroStartingSolution_) {
419  Y.putScalar (STS::zero ());
420  }
421 
422  // Flops are updated in each of the following.
423  switch (PrecType_) {
424  case Ifpack2::Details::JACOBI:
425  ApplyInverseJacobi(*X_copy,Y);
426  break;
427  case Ifpack2::Details::GS:
428  ApplyInverseGS(*X_copy,Y);
429  break;
430  case Ifpack2::Details::SGS:
431  ApplyInverseSGS(*X_copy,Y);
432  break;
433  default:
435  (true, std::logic_error, "Ifpack2::BlockRelaxation::apply: Invalid "
436  "PrecType_ enum value " << PrecType_ << ". Valid values are Ifpack2::"
437  "Details::JACOBI = " << Ifpack2::Details::JACOBI << ", Ifpack2::Details"
438  "::GS = " << Ifpack2::Details::GS << ", and Ifpack2::Details::SGS = "
439  << Ifpack2::Details::SGS << ". Please report this bug to the Ifpack2 "
440  "developers.");
441  }
442 
443  ++NumApply_;
444  Time_->stop();
445  ApplyTime_ += Time_->totalElapsedTime();
446 }
447 
448 template<class MatrixType,class ContainerType>
449 void
451 applyMat (const Tpetra::MultiVector<typename MatrixType::scalar_type,
452  typename MatrixType::local_ordinal_type,
453  typename MatrixType::global_ordinal_type,
454  typename MatrixType::node_type>& X,
455  Tpetra::MultiVector<typename MatrixType::scalar_type,
456  typename MatrixType::local_ordinal_type,
457  typename MatrixType::global_ordinal_type,
458  typename MatrixType::node_type>& Y,
459  Teuchos::ETransp mode) const
460 {
461  A_->apply (X, Y, mode);
462 }
463 
464 template<class MatrixType,class ContainerType>
465 void
468 {
469  using Teuchos::rcp;
470  typedef Tpetra::RowGraph<local_ordinal_type, global_ordinal_type, node_type>
471  row_graph_type;
472 
474  (A_.is_null (), std::runtime_error, "Ifpack2::BlockRelaxation::initialize: "
475  "The matrix is null. You must call setMatrix() with a nonnull matrix "
476  "before you may call this method.");
477 
478  // Check whether we have a BlockCrsMatrix
480  Teuchos::rcp_dynamic_cast<const block_crs_matrix_type> (A_);
481  hasBlockCrsMatrix_ = !A_bcrs.is_null();
482  if (A_bcrs.is_null ()) {
483  hasBlockCrsMatrix_ = false;
484  }
485  else {
486  hasBlockCrsMatrix_ = true;
487  }
488 
489  IsInitialized_ = false;
490  Time_->start (true);
491 
492  NumLocalRows_ = A_->getNodeNumRows ();
493  NumGlobalRows_ = A_->getGlobalNumRows ();
494  NumGlobalNonzeros_ = A_->getGlobalNumEntries ();
495 
496  // NTS: Will need to add support for Zoltan2 partitions later Also,
497  // will need a better way of handling the Graph typing issue.
498  // Especially with ordinal types w.r.t the container.
499  Partitioner_ = Teuchos::null;
500 
501  if (PartitionerType_ == "linear") {
502  Partitioner_ =
503  rcp (new Ifpack2::LinearPartitioner<row_graph_type> (A_->getGraph ()));
504  } else if (PartitionerType_ == "line") {
505  Partitioner_ =
507  } else if (PartitionerType_ == "user") {
508  Partitioner_ =
509  rcp (new Ifpack2::Details::UserPartitioner<row_graph_type> (A_->getGraph () ) );
510  } else {
511  // We should have checked for this in setParameters(), so it's a
512  // logic_error, not an invalid_argument or runtime_error.
514  (true, std::logic_error, "Ifpack2::BlockRelaxation::initialize: Unknown "
515  "partitioner type " << PartitionerType_ << ". Valid values are "
516  "\"linear\", \"line\", and \"user\".");
517  }
518 
519  // need to partition the graph of A
520  Partitioner_->setParameters (List_);
521  Partitioner_->compute ();
522 
523  // get actual number of partitions
524  NumLocalBlocks_ = Partitioner_->numLocalParts ();
525 
526  // Note: Unlike Ifpack, we'll punt on computing W_ until compute(), which is where
527  // we assume that the type of relaxation has been chosen.
528 
529  if (A_->getComm()->getSize() != 1) {
530  IsParallel_ = true;
531  } else {
532  IsParallel_ = false;
533  }
534 
535  ++NumInitialize_;
536  Time_->stop ();
537  InitializeTime_ += Time_->totalElapsedTime ();
538  IsInitialized_ = true;
539 }
540 
541 
542 template<class MatrixType,class ContainerType>
543 void
546 {
547  // typedef Tpetra::Map<local_ordinal_type,global_ordinal_type,node_type> map_type; // unused
548  // typedef Tpetra::Import<local_ordinal_type,global_ordinal_type, node_type> import_type; // unused
549  using Teuchos::Ptr;
550  using Teuchos::RCP;
551  using Teuchos::rcp;
552  using Teuchos::Array;
553  using Teuchos::ArrayView;
554 
555  Time_->start (true);
556 
557  // reset values
558  IsComputed_ = false;
559 
560  // Extract the submatrices
561  ExtractSubmatrices ();
562 
563  // Compute the weight vector if we're doing overlapped Jacobi (and
564  // only if we're doing overlapped Jacobi).
566  (PrecType_ == Ifpack2::Details::JACOBI && OverlapLevel_ > 0, std::runtime_error,
567  "Ifpack2::BlockRelaxation::computeBlockCrs: "
568  "We do not support overlapped Jacobi yet for Tpetra::BlockCrsMatrix. Sorry!");
569 
570  ++NumCompute_;
571  Time_->stop ();
572  ComputeTime_ += Time_->totalElapsedTime();
573  IsComputed_ = true;
574 }
575 
576 
577 template<class MatrixType,class ContainerType>
578 void
581 {
582  using Teuchos::rcp;
583  typedef Tpetra::Vector<scalar_type,
585  // typedef Tpetra::Import<local_ordinal_type,
586  // global_ordinal_type, node_type> import_type; // unused
587 
589  (A_.is_null (), std::runtime_error, "Ifpack2::BlockRelaxation::compute: "
590  "The matrix is null. You must call setMatrix() with a nonnull matrix, "
591  "then call initialize(), before you may call this method.");
592 
593  // We should have checked for this in setParameters(), so it's a
594  // logic_error, not an invalid_argument or runtime_error.
596  (NumSweeps_ < 0, std::logic_error, "Ifpack2::BlockRelaxation::compute: "
597  "NumSweeps_ = " << NumSweeps_ << " < 0.");
598 
599  if (! isInitialized ()) {
600  initialize ();
601  }
602 
603  if (hasBlockCrsMatrix_) {
604  computeBlockCrs ();
605  return;
606  }
607 
608  Time_->start (true);
609 
610  // reset values
611  IsComputed_ = false;
612 
613  // Extract the submatrices
614  ExtractSubmatrices ();
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) {
619  // weight of each vertex
620  W_ = rcp (new vector_type (A_->getRowMap ()));
621  W_->putScalar (STS::zero ());
622  Teuchos::ArrayRCP<scalar_type > w_ptr = W_->getDataNonConst(0);
623 
624  for (local_ordinal_type i = 0 ; i < NumLocalBlocks_ ; ++i) {
625  for (size_t j = 0 ; j < Partitioner_->numRowsInPart(i) ; ++j) {
626  // FIXME (mfh 12 Sep 2014) Should this really be int?
627  // Perhaps it should be local_ordinal_type instead.
628  int LID = (*Partitioner_)(i,j);
629  w_ptr[LID]+= STS::one();
630  }
631  }
632  W_->reciprocal (*W_);
633  }
634  ++NumCompute_;
635  Time_->stop ();
636  ComputeTime_ += Time_->totalElapsedTime();
637  IsComputed_ = true;
638 }
639 
640 template<class MatrixType, class ContainerType>
641 void
644 {
646  (Partitioner_.is_null (), std::runtime_error, "Ifpack2::BlockRelaxation::"
647  "ExtractSubmatrices: Partitioner object is null.");
648 
649  NumLocalBlocks_ = Partitioner_->numLocalParts ();
650  std::string containerType = ContainerType::getName ();
651  if (containerType == "Generic") {
652  // ContainerType is Container<row_matrix_type>. Thus, we need to
653  // get the container name from the parameter list. We use "Dense"
654  // as the default container type.
655  containerType = containerType_;
656  }
657 
658  Teuchos::Array<Teuchos::Array<local_ordinal_type> > localRows(NumLocalBlocks_);
659  for (local_ordinal_type i = 0; i < NumLocalBlocks_; ++i) {
660  const size_t numRows = Partitioner_->numRowsInPart (i);
661 
662  localRows[i].resize(numRows);
663  // Extract a list of the indices of each partitioner row.
664  for (size_t j = 0; j < numRows; ++j) {
665  localRows[i][j] = (*Partitioner_) (i,j);
666  }
667  }
668  Container_ = Teuchos::rcp_static_cast<Container<MatrixType>>
669  (Details::createContainer<row_matrix_type> (containerType, A_, localRows, Importer_,
670  OverlapLevel_, DampingFactor_));
671  Container_->setParameters(List_);
672  Container_->initialize();
673  Container_->compute(); //initialize + compute each block matrix
674 }
675 
676 template<class MatrixType,class ContainerType>
677 void
678 BlockRelaxation<MatrixType,ContainerType>::
679 ApplyInverseJacobi (const MV& X, MV& Y) const
680 {
681  const size_t NumVectors = X.getNumVectors ();
682  typename ContainerType::HostView XView = X.template getLocalView<Kokkos::HostSpace>();
683  typename ContainerType::HostView YView = Y.template getLocalView<Kokkos::HostSpace>();
684  MV AY (Y.getMap (), NumVectors);
685 
686  typename ContainerType::HostView AYView = AY.template getLocalView<Kokkos::HostSpace>();
687  // Initial matvec not needed
688  int starting_iteration = 0;
689  if (OverlapLevel_ > 0)
690  {
691  //Overlapping jacobi, with view of W_
692  typename ContainerType::HostView WView = W_->template getLocalView<Kokkos::HostSpace>();
693  if(ZeroStartingSolution_) {
694  Container_->DoOverlappingJacobi(XView, YView, WView, X.getStride());
695  starting_iteration = 1;
696  }
697  const scalar_type ONE = STS::one();
698  for(int j = starting_iteration; j < NumSweeps_; j++)
699  {
700  applyMat (Y, AY);
701  AY.update (ONE, X, -ONE);
702  Container_->DoOverlappingJacobi (AYView, YView, WView, X.getStride());
703  }
704  }
705  else
706  {
707  //Non-overlapping
708  if(ZeroStartingSolution_)
709  {
710  Container_->DoJacobi (XView, YView, X.getStride());
711  starting_iteration = 1;
712  }
713  const scalar_type ONE = STS::one();
714  for(int j = starting_iteration; j < NumSweeps_; j++)
715  {
716  applyMat (Y, AY);
717  AY.update (ONE, X, -ONE);
718  Container_->DoJacobi (AYView, YView, X.getStride());
719  }
720  }
721 }
722 
723 template<class MatrixType,class ContainerType>
724 void
725 BlockRelaxation<MatrixType,ContainerType>::
726 ApplyInverseGS (const MV& X, MV& Y) const
727 {
728  using Teuchos::Ptr;
729  using Teuchos::ptr;
730  size_t numVecs = X.getNumVectors();
731  //Get view of X (is never modified in this function)
732  typename ContainerType::HostView XView = X.template getLocalView<Kokkos::HostSpace>();
733  typename ContainerType::HostView YView = Y.template getLocalView<Kokkos::HostSpace>();
734  //Pre-import Y, if parallel
735  Ptr<MV> Y2;
736  bool deleteY2 = false;
737  if(IsParallel_)
738  {
739  Y2 = ptr(new MV(Importer_->getTargetMap(), numVecs));
740  deleteY2 = true;
741  }
742  else
743  Y2 = ptr(&Y);
744  if(IsParallel_)
745  {
746  for(int j = 0; j < NumSweeps_; ++j)
747  {
748  //do import once per sweep
749  Y2->doImport(Y, *Importer_, Tpetra::INSERT);
750  typename ContainerType::HostView Y2View = Y2->template getLocalView<Kokkos::HostSpace>();
751  Container_->DoGaussSeidel(XView, YView, Y2View, X.getStride());
752  }
753  }
754  else
755  {
756  typename ContainerType::HostView Y2View = Y2->template getLocalView<Kokkos::HostSpace>();
757  for(int j = 0; j < NumSweeps_; ++j)
758  {
759  Container_->DoGaussSeidel(XView, YView, Y2View, X.getStride());
760  }
761  }
762  if(deleteY2)
763  delete Y2.get();
764 }
765 
766 template<class MatrixType,class ContainerType>
767 void
768 BlockRelaxation<MatrixType,ContainerType>::
769 ApplyInverseSGS (const MV& X, MV& Y) const
770 {
771  using Teuchos::Ptr;
772  using Teuchos::ptr;
773  //Get view of X (is never modified in this function)
774  typename ContainerType::HostView XView = X.template getLocalView<Kokkos::HostSpace>();
775  typename ContainerType::HostView YView = Y.template getLocalView<Kokkos::HostSpace>();
776  //Pre-import Y, if parallel
777  Ptr<MV> Y2;
778  bool deleteY2 = false;
779  if(IsParallel_)
780  {
781  Y2 = ptr(new MV(Importer_->getTargetMap(), X.getNumVectors()));
782  deleteY2 = true;
783  }
784  else
785  Y2 = ptr(&Y);
786  if(IsParallel_)
787  {
788  for(int j = 0; j < NumSweeps_; ++j)
789  {
790  //do import once per sweep
791  Y2->doImport(Y, *Importer_, Tpetra::INSERT);
792  typename ContainerType::HostView Y2View = Y2->template getLocalView<Kokkos::HostSpace>();
793  Container_->DoSGS(XView, YView, Y2View, X.getStride());
794  }
795  }
796  else
797  {
798  typename ContainerType::HostView Y2View = Y2->template getLocalView<Kokkos::HostSpace>();
799  for(int j = 0; j < NumSweeps_; ++j)
800  {
801  Container_->DoSGS(XView, YView, Y2View, X.getStride());
802  }
803  }
804  if(deleteY2)
805  delete Y2.get();
806 }
807 
808 template<class MatrixType,class ContainerType>
809 void BlockRelaxation<MatrixType,ContainerType>::computeImporter () const
810 {
811  using Teuchos::RCP;
812  using Teuchos::rcp;
813  using Teuchos::Ptr;
814  using Teuchos::ArrayView;
815  using Teuchos::Array;
816  using Teuchos::Comm;
817  using Teuchos::rcp_dynamic_cast;
818  if(IsParallel_)
819  {
820  if(hasBlockCrsMatrix_)
821  {
822  const RCP<const block_crs_matrix_type> bcrs =
823  rcp_dynamic_cast<const block_crs_matrix_type>(A_);
824  int bs_ = bcrs->getBlockSize();
825  RCP<const map_type> oldDomainMap = A_->getDomainMap();
826  RCP<const map_type> oldColMap = A_->getColMap();
827  // Because A is a block CRS matrix, import will not do what you think it does
828  // We have to construct the correct maps for it
829  global_size_t numGlobalElements = oldColMap->getGlobalNumElements() * bs_;
830  global_ordinal_type indexBase = oldColMap->getIndexBase();
831  RCP<const Comm<int>> comm = oldColMap->getComm();
832  ArrayView<const global_ordinal_type> oldColElements = oldColMap->getNodeElementList();
833  Array<global_ordinal_type> newColElements(bs_ * oldColElements.size());
834  for(int i = 0; i < oldColElements.size(); i++)
835  {
836  for(int j = 0; j < bs_; j++)
837  newColElements[i * bs_ + j] = oldColElements[i] * bs_ + j;
838  }
839  RCP<map_type> colMap(new map_type(numGlobalElements, newColElements, indexBase, comm));
840  // Create the importer
841  Importer_ = rcp(new import_type(oldDomainMap, colMap));
842  }
843  else if(!A_.is_null())
844  {
845  Importer_ = A_->getGraph()->getImporter();
846  if(Importer_.is_null())
847  Importer_ = rcp(new import_type(A_->getDomainMap(), A_->getColMap()));
848  }
849  }
850  //otherwise, Importer_ is not needed and is left NULL
851 }
852 
853 template<class MatrixType, class ContainerType>
854 std::string
856 description () const
857 {
858  std::ostringstream out;
859 
860  // Output is a valid YAML dictionary in flow style. If you don't
861  // like everything on a single line, you should call describe()
862  // instead.
863  out << "\"Ifpack2::BlockRelaxation\": {";
864  if (this->getObjectLabel () != "") {
865  out << "Label: \"" << this->getObjectLabel () << "\", ";
866  }
867  out << "Initialized: " << (isInitialized () ? "true" : "false") << ", ";
868  out << "Computed: " << (isComputed () ? "true" : "false") << ", ";
869  if (A_.is_null ()) {
870  out << "Matrix: null, ";
871  }
872  else {
873  out << "Matrix: not null"
874  << ", Global matrix dimensions: ["
875  << A_->getGlobalNumRows () << ", " << A_->getGlobalNumCols () << "], ";
876  }
877 
878  // It's useful to print this instance's relaxation method. If you
879  // want more info than that, call describe() instead.
880  out << "\"relaxation: type\": ";
881  if (PrecType_ == Ifpack2::Details::JACOBI) {
882  out << "Block Jacobi";
883  } else if (PrecType_ == Ifpack2::Details::GS) {
884  out << "Block Gauss-Seidel";
885  } else if (PrecType_ == Ifpack2::Details::SGS) {
886  out << "Block Symmetric Gauss-Seidel";
887  } else {
888  out << "INVALID";
889  }
890 
891  out << ", " << "sweeps: " << NumSweeps_ << ", "
892  << "damping factor: " << DampingFactor_ << ", ";
893 
894  out << "}";
895  return out.str();
896 }
897 
898 template<class MatrixType,class ContainerType>
899 void
902  const Teuchos::EVerbosityLevel verbLevel) const
903 {
904  using std::endl;
905  using std::setw;
907  using Teuchos::VERB_DEFAULT;
908  using Teuchos::VERB_NONE;
909  using Teuchos::VERB_LOW;
910  using Teuchos::VERB_MEDIUM;
911  using Teuchos::VERB_HIGH;
912  using Teuchos::VERB_EXTREME;
913 
914  Teuchos::EVerbosityLevel vl = verbLevel;
915  if (vl == VERB_DEFAULT) vl = VERB_LOW;
916  const int myImageID = A_->getComm()->getRank();
917 
918  // Convention requires that describe() begin with a tab.
919  Teuchos::OSTab tab (out);
920 
921  // none: print nothing
922  // low: print O(1) info from node 0
923  // medium:
924  // high:
925  // extreme:
926  if (vl != VERB_NONE && myImageID == 0) {
927  out << "Ifpack2::BlockRelaxation<"
928  << TypeNameTraits<MatrixType>::name () << ", "
929  << TypeNameTraits<ContainerType>::name () << " >:";
930  Teuchos::OSTab tab1 (out);
931 
932  if (this->getObjectLabel () != "") {
933  out << "label: \"" << this->getObjectLabel () << "\"" << endl;
934  }
935  out << "initialized: " << (isInitialized () ? "true" : "false") << endl
936  << "computed: " << (isComputed () ? "true" : "false") << endl;
937 
938  std::string precType;
939  if (PrecType_ == Ifpack2::Details::JACOBI) {
940  precType = "Block Jacobi";
941  } else if (PrecType_ == Ifpack2::Details::GS) {
942  precType = "Block Gauss-Seidel";
943  } else if (PrecType_ == Ifpack2::Details::SGS) {
944  precType = "Block Symmetric Gauss-Seidel";
945  }
946  out << "type: " << precType << endl
947  << "global number of rows: " << A_->getGlobalNumRows () << endl
948  << "global number of columns: " << A_->getGlobalNumCols () << endl
949  << "number of sweeps: " << NumSweeps_ << endl
950  << "damping factor: " << DampingFactor_ << endl
951  << "backwards mode: "
952  << ((PrecType_ == Ifpack2::Details::GS && DoBackwardGS_) ? "true" : "false")
953  << endl
954  << "zero starting solution: "
955  << (ZeroStartingSolution_ ? "true" : "false") << endl;
956  }
957 }
958 
959 }//namespace Ifpack2
960 
961 
962 // Macro that does explicit template instantiation (ETI) for
963 // Ifpack2::BlockRelaxation. S, LO, GO, N correspond to the four
964 // template parameters of Ifpack2::Preconditioner and
965 // Tpetra::RowMatrix.
966 //
967 // We only instantiate for MatrixType = Tpetra::RowMatrix. There's no
968 // need to instantiate for Tpetra::CrsMatrix too. All Ifpack2
969 // preconditioners can and should do dynamic casts if they need a type
970 // more specific than RowMatrix. This keeps build time short and
971 // library and executable sizes small.
972 
973 #ifdef HAVE_IFPACK2_EXPLICIT_INSTANTIATION
974 
975 #define IFPACK2_BLOCKRELAXATION_INSTANT(S,LO,GO,N) \
976  template \
977  class Ifpack2::BlockRelaxation< \
978  Tpetra::RowMatrix<S, LO, GO, N>, \
979  Ifpack2::Container<Tpetra::RowMatrix<S, LO, GO, N> > >;
980 
981 #endif // HAVE_IFPACK2_EXPLICIT_INSTANTIATION
982 
983 #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:262
Ifpack2::BlockRelaxation class declaration.
double getApplyTime() const
Returns the time spent in apply().
Definition: Ifpack2_BlockRelaxation_def.hpp:357
MatrixType::global_ordinal_type global_ordinal_type
The type of global indices in the input MatrixType.
Definition: Ifpack2_BlockRelaxation_decl.hpp:104
VERB_DEFAULT
double getComputeTime() const
Returns the time spent in compute().
Definition: Ifpack2_BlockRelaxation_def.hpp:349
int getNumInitialize() const
Returns the number of calls to initialize().
Definition: Ifpack2_BlockRelaxation_def.hpp:318
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:901
double getInitializeTime() const
Returns the time spent in initialize().
Definition: Ifpack2_BlockRelaxation_def.hpp:341
Teuchos::RCP< const row_matrix_type > getMatrix() const
The input matrix of this preconditioner&#39;s constructor.
Definition: Ifpack2_BlockRelaxation_def.hpp:276
Teuchos::RCP< const map_type > getDomainMap() const
Returns the Tpetra::Map object associated with the domain of this operator.
Definition: Ifpack2_BlockRelaxation_def.hpp:285
#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:365
std::string description() const
A one-line description of this object.
Definition: Ifpack2_BlockRelaxation_def.hpp:856
Teuchos::RCP< const map_type > getRangeMap() const
Returns the Tpetra::Map object associated with the range of this operator.
Definition: Ifpack2_BlockRelaxation_def.hpp:299
void setParameters(const Teuchos::ParameterList &params)
Sets all the parameters for the preconditioner.
Definition: Ifpack2_BlockRelaxation_def.hpp:122
bool is_null() const
void compute()
compute the preconditioner for the specified matrix, diagonal perturbation thresholds and relaxation ...
Definition: Ifpack2_BlockRelaxation_def.hpp:580
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...
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:451
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:467
MatrixType::node_type node_type
Node type of the input MatrixType.
Definition: Ifpack2_BlockRelaxation_decl.hpp:106
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:325
int getNumApply() const
Returns the number of calls to apply().
Definition: Ifpack2_BlockRelaxation_def.hpp:333
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
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:72
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
BlockRelaxation(const Teuchos::RCP< const row_matrix_type > &Matrix)
Constructor.
Definition: Ifpack2_BlockRelaxation_def.hpp:83