Ifpack2 Templated Preconditioning Package  Version 1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
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.");
454  TEUCHOS_TEST_FOR_EXCEPTION(
455  X.getNumVectors () != Y.getNumVectors (), std::invalid_argument,
456  "Ifpack2::BlockRelaxation::apply: X.getNumVectors() = "
457  << X.getNumVectors () << " != Y.getNumVectors() = "
458  << Y.getNumVectors () << ".");
459  TEUCHOS_TEST_FOR_EXCEPTION(
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.");
463  TEUCHOS_TEST_FOR_EXCEPTION(
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 << ".");
467  TEUCHOS_TEST_FOR_EXCEPTION(
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 #ifdef KOKKOS_ENABLE_DEPRECATED_CODE
479  auto X_lcl_host = X.template getLocalView<Kokkos::HostSpace> ();
480  auto Y_lcl_host = Y.template getLocalView<Kokkos::HostSpace> ();
481 #else
482  auto X_lcl_host = X.getLocalViewHost ();
483  auto Y_lcl_host = Y.getLocalViewHost ();
484 #endif
485  if (X_lcl_host.data () == Y_lcl_host.data ()) {
486  X_copy = rcp (new MV (X, Teuchos::Copy));
487  } else {
488  X_copy = rcpFromRef (X);
489  }
490  }
491 
492  if (ZeroStartingSolution_) {
493  Y.putScalar (STS::zero ());
494  }
495 
496  // Flops are updated in each of the following.
497  switch (PrecType_) {
498  case Ifpack2::Details::JACOBI:
499  ApplyInverseJacobi(*X_copy,Y);
500  break;
501  case Ifpack2::Details::GS:
502  ApplyInverseGS(*X_copy,Y);
503  break;
504  case Ifpack2::Details::SGS:
505  ApplyInverseSGS(*X_copy,Y);
506  break;
507  case Ifpack2::Details::MTSPLITJACOBI:
508  Container_->applyInverseJacobi(*X_copy, Y, ZeroStartingSolution_, NumSweeps_);
509  break;
510  default:
511  TEUCHOS_TEST_FOR_EXCEPTION
512  (true, std::logic_error, "Ifpack2::BlockRelaxation::apply: Invalid "
513  "PrecType_ enum value " << PrecType_ << ". Valid values are Ifpack2::"
514  "Details::JACOBI = " << Ifpack2::Details::JACOBI << ", Ifpack2::Details"
515  "::GS = " << Ifpack2::Details::GS << ", and Ifpack2::Details::SGS = "
516  << Ifpack2::Details::SGS << ". Please report this bug to the Ifpack2 "
517  "developers.");
518  }
519 
520  ++NumApply_;
521  Time_->stop();
522  ApplyTime_ += Time_->totalElapsedTime();
523 }
524 
525 template<class MatrixType,class ContainerType>
526 void
528 applyMat (const Tpetra::MultiVector<typename MatrixType::scalar_type,
529  typename MatrixType::local_ordinal_type,
530  typename MatrixType::global_ordinal_type,
531  typename MatrixType::node_type>& X,
532  Tpetra::MultiVector<typename MatrixType::scalar_type,
533  typename MatrixType::local_ordinal_type,
534  typename MatrixType::global_ordinal_type,
535  typename MatrixType::node_type>& Y,
536  Teuchos::ETransp mode) const
537 {
538  A_->apply (X, Y, mode);
539 }
540 
541 template<class MatrixType,class ContainerType>
542 void
545 {
546  using Teuchos::rcp;
547  typedef Tpetra::RowGraph<local_ordinal_type, global_ordinal_type, node_type>
548  row_graph_type;
549 
551  (A_.is_null (), std::runtime_error, "Ifpack2::BlockRelaxation::initialize: "
552  "The matrix is null. You must call setMatrix() with a nonnull matrix "
553  "before you may call this method.");
554 
555  // Check whether we have a BlockCrsMatrix
557  Teuchos::rcp_dynamic_cast<const block_crs_matrix_type> (A_);
558  hasBlockCrsMatrix_ = !A_bcrs.is_null();
559  if (A_bcrs.is_null ()) {
560  hasBlockCrsMatrix_ = false;
561  }
562  else {
563  hasBlockCrsMatrix_ = true;
564  }
565 
566  IsInitialized_ = false;
567  Time_->start (true);
568 
569  NumLocalRows_ = A_->getNodeNumRows ();
570  NumGlobalRows_ = A_->getGlobalNumRows ();
571  NumGlobalNonzeros_ = A_->getGlobalNumEntries ();
572 
573  // NTS: Will need to add support for Zoltan2 partitions later Also,
574  // will need a better way of handling the Graph typing issue.
575  // Especially with ordinal types w.r.t the container.
576  Partitioner_ = Teuchos::null;
577 
578  if (PartitionerType_ == "linear") {
579  Partitioner_ =
580  rcp (new Ifpack2::LinearPartitioner<row_graph_type> (A_->getGraph ()));
581  } else if (PartitionerType_ == "line") {
582  Partitioner_ =
584  } else if (PartitionerType_ == "user") {
585  Partitioner_ =
586  rcp (new Ifpack2::Details::UserPartitioner<row_graph_type> (A_->getGraph () ) );
587  } else {
588  // We should have checked for this in setParameters(), so it's a
589  // logic_error, not an invalid_argument or runtime_error.
591  (true, std::logic_error, "Ifpack2::BlockRelaxation::initialize: Unknown "
592  "partitioner type " << PartitionerType_ << ". Valid values are "
593  "\"linear\", \"line\", and \"user\".");
594  }
595 
596  // need to partition the graph of A
597  Partitioner_->setParameters (List_);
598  Partitioner_->compute ();
599 
600  // get actual number of partitions
601  NumLocalBlocks_ = Partitioner_->numLocalParts ();
602 
603  // Note: Unlike Ifpack, we'll punt on computing W_ until compute(), which is where
604  // we assume that the type of relaxation has been chosen.
605 
606  if (A_->getComm()->getSize() != 1) {
607  IsParallel_ = true;
608  } else {
609  IsParallel_ = false;
610  }
611 
612  // We should have checked for this in setParameters(), so it's a
613  // logic_error, not an invalid_argument or runtime_error.
615  (NumSweeps_ < 0, std::logic_error, "Ifpack2::BlockRelaxation::initialize: "
616  "NumSweeps_ = " << NumSweeps_ << " < 0.");
617 
618  // Extract the submatrices
619  ExtractSubmatricesStructure ();
620 
621  // Compute the weight vector if we're doing overlapped Jacobi (and
622  // only if we're doing overlapped Jacobi).
623  if (PrecType_ == Ifpack2::Details::JACOBI && OverlapLevel_ > 0) {
625  (hasBlockCrsMatrix_, std::runtime_error,
626  "Ifpack2::BlockRelaxation::initialize: "
627  "We do not support overlapped Jacobi yet for Tpetra::BlockCrsMatrix. Sorry!");
628 
629  // weight of each vertex
630  W_ = rcp (new vector_type (A_->getRowMap ()));
631  W_->putScalar (STS::zero ());
632  Teuchos::ArrayRCP<scalar_type > w_ptr = W_->getDataNonConst(0);
633 
634  for (local_ordinal_type i = 0 ; i < NumLocalBlocks_ ; ++i) {
635  for (size_t j = 0 ; j < Partitioner_->numRowsInPart(i) ; ++j) {
636  // FIXME (mfh 12 Sep 2014) Should this really be int?
637  // Perhaps it should be local_ordinal_type instead.
638  int LID = (*Partitioner_)(i,j);
639  w_ptr[LID]+= STS::one();
640  }
641  }
642  W_->reciprocal (*W_);
643  }
644 
645  ++NumInitialize_;
646  Time_->stop ();
647  InitializeTime_ += Time_->totalElapsedTime ();
648  IsInitialized_ = true;
649 }
650 
651 
652 template<class MatrixType,class ContainerType>
653 void
656 {
657  using Teuchos::rcp;
658 
660  (A_.is_null (), std::runtime_error, "Ifpack2::BlockRelaxation::compute: "
661  "The matrix is null. You must call setMatrix() with a nonnull matrix, "
662  "then call initialize(), before you may call this method.");
663 
664  if (! isInitialized ()) {
665  initialize ();
666  }
667 
668  Time_->start (true);
669 
670  // reset values
671  IsComputed_ = false;
672 
673  Container_->compute(); // compute each block matrix
674 
675  ++NumCompute_;
676  Time_->stop ();
677  ComputeTime_ += Time_->totalElapsedTime();
678  IsComputed_ = true;
679 }
680 
681 template<class MatrixType, class ContainerType>
682 void
685 {
687  (Partitioner_.is_null (), std::runtime_error, "Ifpack2::BlockRelaxation::"
688  "ExtractSubmatricesStructure: Partitioner object is null.");
689 
690  NumLocalBlocks_ = Partitioner_->numLocalParts ();
691  std::string containerType = ContainerType::getName ();
692  if (containerType == "Generic") {
693  // ContainerType is Container<row_matrix_type>. Thus, we need to
694  // get the container name from the parameter list. We use "TriDi"
695  // as the default container type.
696  containerType = containerType_;
697  }
698 
699  Teuchos::Array<Teuchos::Array<local_ordinal_type> > localRows(NumLocalBlocks_);
700  for (local_ordinal_type i = 0; i < NumLocalBlocks_; ++i) {
701  const size_t numRows = Partitioner_->numRowsInPart (i);
702 
703  localRows[i].resize(numRows);
704  // Extract a list of the indices of each partitioner row.
705  for (size_t j = 0; j < numRows; ++j) {
706  localRows[i][j] = (*Partitioner_) (i,j);
707  }
708  }
709  Container_ = ContainerFactory<MatrixType>::build(containerType, A_, localRows, Importer_, OverlapLevel_, DampingFactor_);
710  Container_->setParameters(List_);
711  Container_->initialize();
712 }
713 
714 template<class MatrixType,class ContainerType>
715 void
716 BlockRelaxation<MatrixType,ContainerType>::
717 ApplyInverseJacobi (const MV& X, MV& Y) const
718 {
719  const size_t NumVectors = X.getNumVectors ();
720 #ifdef KOKKOS_ENABLE_DEPRECATED_CODE
721  typename ContainerType::HostView XView = X.template getLocalView<Kokkos::HostSpace> ();
722  typename ContainerType::HostView YView = Y.template getLocalView<Kokkos::HostSpace> ();
723 #else
724  auto XView = X.getLocalViewHost ();
725  auto YView = Y.getLocalViewHost ();
726 #endif
727  MV AY (Y.getMap (), NumVectors);
728 
729 #ifdef KOKKOS_ENABLE_DEPRECATED_CODE
730  typename ContainerType::HostView AYView = AY.template getLocalView<Kokkos::HostSpace> ();
731 #else
732  auto AYView = AY.getLocalViewHost ();
733 #endif
734  // Initial matvec not needed
735  int starting_iteration = 0;
736  if (OverlapLevel_ > 0)
737  {
738  //Overlapping jacobi, with view of W_
739 #ifdef KOKKOS_ENABLE_DEPRECATED_CODE
740  typename ContainerType::HostView WView = W_->template getLocalView<Kokkos::HostSpace> ();
741 #else
742  auto WView = W_->getLocalViewHost ();
743 #endif
744  if(ZeroStartingSolution_) {
745  Container_->DoOverlappingJacobi(XView, YView, WView, X.getStride());
746  starting_iteration = 1;
747  }
748  const scalar_type ONE = STS::one();
749  for(int j = starting_iteration; j < NumSweeps_; j++)
750  {
751  applyMat (Y, AY);
752  AY.update (ONE, X, -ONE);
753  Container_->DoOverlappingJacobi (AYView, YView, WView, X.getStride());
754  }
755  }
756  else
757  {
758  //Non-overlapping
759  if(ZeroStartingSolution_)
760  {
761  Container_->DoJacobi (XView, YView, X.getStride());
762  starting_iteration = 1;
763  }
764  const scalar_type ONE = STS::one();
765  for(int j = starting_iteration; j < NumSweeps_; j++)
766  {
767  applyMat (Y, AY);
768  AY.update (ONE, X, -ONE);
769  Container_->DoJacobi (AYView, YView, X.getStride());
770  }
771  }
772 }
773 
774 template<class MatrixType,class ContainerType>
775 void
776 BlockRelaxation<MatrixType,ContainerType>::
777 ApplyInverseGS (const MV& X, MV& Y) const
778 {
779  using Teuchos::Ptr;
780  using Teuchos::ptr;
781  size_t numVecs = X.getNumVectors();
782  //Get view of X (is never modified in this function)
783 #ifdef KOKKOS_ENABLE_DEPRECATED_CODE
784  typename ContainerType::HostView XView = X.template getLocalView<Kokkos::HostSpace> ();
785  typename ContainerType::HostView YView = Y.template getLocalView<Kokkos::HostSpace> ();
786 #else
787  auto XView = X.getLocalViewHost ();
788  auto YView = Y.getLocalViewHost ();
789 #endif
790  //Pre-import Y, if parallel
791  Ptr<MV> Y2;
792  bool deleteY2 = false;
793  if(IsParallel_)
794  {
795  Y2 = ptr(new MV(Importer_->getTargetMap(), numVecs));
796  deleteY2 = true;
797  }
798  else
799  Y2 = ptr(&Y);
800  if(IsParallel_)
801  {
802  for(int j = 0; j < NumSweeps_; ++j)
803  {
804  //do import once per sweep
805  Y2->doImport(Y, *Importer_, Tpetra::INSERT);
806 #ifdef KOKKOS_ENABLE_DEPRECATED_CODE
807  typename ContainerType::HostView Y2View = Y2->template getLocalView<Kokkos::HostSpace> ();
808 #else
809  auto Y2View = Y2->getLocalViewHost ();
810 #endif
811  Container_->DoGaussSeidel(XView, YView, Y2View, X.getStride());
812  }
813  }
814  else
815  {
816 #ifdef KOKKOS_ENABLE_DEPRECATED_CODE
817  typename ContainerType::HostView Y2View = Y2->template getLocalView<Kokkos::HostSpace> ();
818 #else
819  auto Y2View = Y2->getLocalViewHost ();
820 #endif
821  for(int j = 0; j < NumSweeps_; ++j)
822  {
823  Container_->DoGaussSeidel(XView, YView, Y2View, X.getStride());
824  }
825  }
826  if(deleteY2)
827  delete Y2.get();
828 }
829 
830 template<class MatrixType,class ContainerType>
831 void
832 BlockRelaxation<MatrixType,ContainerType>::
833 ApplyInverseSGS (const MV& X, MV& Y) const
834 {
835  using Teuchos::Ptr;
836  using Teuchos::ptr;
837  //Get view of X (is never modified in this function)
838 #ifdef KOKKOS_ENABLE_DEPRECATED_CODE
839  typename ContainerType::HostView XView = X.template getLocalView<Kokkos::HostSpace> ();
840  typename ContainerType::HostView YView = Y.template getLocalView<Kokkos::HostSpace> ();
841 #else
842  auto XView = X.getLocalViewHost ();
843  auto YView = Y.getLocalViewHost ();
844 #endif
845  //Pre-import Y, if parallel
846  Ptr<MV> Y2;
847  bool deleteY2 = false;
848  if(IsParallel_)
849  {
850  Y2 = ptr(new MV(Importer_->getTargetMap(), X.getNumVectors()));
851  deleteY2 = true;
852  }
853  else
854  Y2 = ptr(&Y);
855  if(IsParallel_)
856  {
857  for(int j = 0; j < NumSweeps_; ++j)
858  {
859  //do import once per sweep
860  Y2->doImport(Y, *Importer_, Tpetra::INSERT);
861 #ifdef KOKKOS_ENABLE_DEPRECATED_CODE
862  typename ContainerType::HostView Y2View = Y2->template getLocalView<Kokkos::HostSpace> ();
863 #else
864  auto Y2View = Y2->getLocalViewHost ();
865 #endif
866  Container_->DoSGS(XView, YView, Y2View, X.getStride());
867  }
868  }
869  else
870  {
871 #ifdef KOKKOS_ENABLE_DEPRECATED_CODE
872  typename ContainerType::HostView Y2View = Y2->template getLocalView<Kokkos::HostSpace> ();
873 #else
874  auto Y2View = Y2->getLocalViewHost ();
875 #endif
876  for(int j = 0; j < NumSweeps_; ++j)
877  {
878  Container_->DoSGS(XView, YView, Y2View, X.getStride());
879  }
880  }
881  if(deleteY2)
882  delete Y2.get();
883 }
884 
885 template<class MatrixType,class ContainerType>
886 void BlockRelaxation<MatrixType,ContainerType>::computeImporter () const
887 {
888  using Teuchos::RCP;
889  using Teuchos::rcp;
890  using Teuchos::Ptr;
891  using Teuchos::ArrayView;
892  using Teuchos::Array;
893  using Teuchos::Comm;
894  using Teuchos::rcp_dynamic_cast;
895  if(IsParallel_)
896  {
897  if(hasBlockCrsMatrix_)
898  {
899  const RCP<const block_crs_matrix_type> bcrs =
900  rcp_dynamic_cast<const block_crs_matrix_type>(A_);
901  int bs_ = bcrs->getBlockSize();
902  RCP<const map_type> oldDomainMap = A_->getDomainMap();
903  RCP<const map_type> oldColMap = A_->getColMap();
904  // Because A is a block CRS matrix, import will not do what you think it does
905  // We have to construct the correct maps for it
906  global_size_t numGlobalElements = oldColMap->getGlobalNumElements() * bs_;
907  global_ordinal_type indexBase = oldColMap->getIndexBase();
908  RCP<const Comm<int>> comm = oldColMap->getComm();
909  ArrayView<const global_ordinal_type> oldColElements = oldColMap->getNodeElementList();
910  Array<global_ordinal_type> newColElements(bs_ * oldColElements.size());
911  for(int i = 0; i < oldColElements.size(); i++)
912  {
913  for(int j = 0; j < bs_; j++)
914  newColElements[i * bs_ + j] = oldColElements[i] * bs_ + j;
915  }
916  RCP<map_type> colMap(new map_type(numGlobalElements, newColElements, indexBase, comm));
917  // Create the importer
918  Importer_ = rcp(new import_type(oldDomainMap, colMap));
919  }
920  else if(!A_.is_null())
921  {
922  Importer_ = A_->getGraph()->getImporter();
923  if(Importer_.is_null())
924  Importer_ = rcp(new import_type(A_->getDomainMap(), A_->getColMap()));
925  }
926  }
927  //otherwise, Importer_ is not needed and is left NULL
928 }
929 
930 template<class MatrixType, class ContainerType>
931 std::string
933 description () const
934 {
935  std::ostringstream out;
936 
937  // Output is a valid YAML dictionary in flow style. If you don't
938  // like everything on a single line, you should call describe()
939  // instead.
940  out << "\"Ifpack2::BlockRelaxation\": {";
941  if (this->getObjectLabel () != "") {
942  out << "Label: \"" << this->getObjectLabel () << "\", ";
943  }
944  // out << "Initialized: " << (isInitialized () ? "true" : "false") << ", ";
945  // out << "Computed: " << (isComputed () ? "true" : "false") << ", ";
946  if (A_.is_null ()) {
947  out << "Matrix: null, ";
948  }
949  else {
950  // out << "Matrix: not null"
951  // << ", Global matrix dimensions: ["
952  out << "Global matrix dimensions: ["
953  << A_->getGlobalNumRows () << ", " << A_->getGlobalNumCols () << "], ";
954  }
955 
956  // It's useful to print this instance's relaxation method. If you
957  // want more info than that, call describe() instead.
958  out << "\"relaxation: type\": ";
959  if (PrecType_ == Ifpack2::Details::JACOBI) {
960  out << "Block Jacobi";
961  } else if (PrecType_ == Ifpack2::Details::GS) {
962  out << "Block Gauss-Seidel";
963  } else if (PrecType_ == Ifpack2::Details::SGS) {
964  out << "Block Symmetric Gauss-Seidel";
965  } else {
966  out << "INVALID";
967  }
968 
969  out << ", overlap: " << OverlapLevel_;
970 
971  out << ", " << "sweeps: " << NumSweeps_ << ", "
972  << "damping factor: " << DampingFactor_ << ", ";
973 
974  std::string containerType = ContainerType::getName();
975  out << "block container: " << ((containerType == "Generic") ? containerType_ : containerType);
976  if(List_.isParameter("partitioner: PDE equations"))
977  out << ", dofs/node: "<<List_.get<int>("partitioner: PDE equations");
978 
979 
980  out << "}";
981  return out.str();
982 }
983 
984 template<class MatrixType,class ContainerType>
985 void
988  const Teuchos::EVerbosityLevel verbLevel) const
989 {
990  using std::endl;
991  using std::setw;
993  using Teuchos::VERB_DEFAULT;
994  using Teuchos::VERB_NONE;
995  using Teuchos::VERB_LOW;
996  using Teuchos::VERB_MEDIUM;
997  using Teuchos::VERB_HIGH;
998  using Teuchos::VERB_EXTREME;
999 
1000  Teuchos::EVerbosityLevel vl = verbLevel;
1001  if (vl == VERB_DEFAULT) vl = VERB_LOW;
1002  const int myImageID = A_->getComm()->getRank();
1003 
1004  // Convention requires that describe() begin with a tab.
1005  Teuchos::OSTab tab (out);
1006 
1007  // none: print nothing
1008  // low: print O(1) info from node 0
1009  // medium:
1010  // high:
1011  // extreme:
1012  if (vl != VERB_NONE && myImageID == 0) {
1013  out << "Ifpack2::BlockRelaxation<"
1014  << TypeNameTraits<MatrixType>::name () << ", "
1015  << TypeNameTraits<ContainerType>::name () << " >:";
1016  Teuchos::OSTab tab1 (out);
1017 
1018  if (this->getObjectLabel () != "") {
1019  out << "label: \"" << this->getObjectLabel () << "\"" << endl;
1020  }
1021  out << "initialized: " << (isInitialized () ? "true" : "false") << endl
1022  << "computed: " << (isComputed () ? "true" : "false") << endl;
1023 
1024  std::string precType;
1025  if (PrecType_ == Ifpack2::Details::JACOBI) {
1026  precType = "Block Jacobi";
1027  } else if (PrecType_ == Ifpack2::Details::GS) {
1028  precType = "Block Gauss-Seidel";
1029  } else if (PrecType_ == Ifpack2::Details::SGS) {
1030  precType = "Block Symmetric Gauss-Seidel";
1031  }
1032  out << "type: " << precType << endl
1033  << "global number of rows: " << A_->getGlobalNumRows () << endl
1034  << "global number of columns: " << A_->getGlobalNumCols () << endl
1035  << "number of sweeps: " << NumSweeps_ << endl
1036  << "damping factor: " << DampingFactor_ << endl
1037  << "backwards mode: "
1038  << ((PrecType_ == Ifpack2::Details::GS && DoBackwardGS_) ? "true" : "false")
1039  << endl
1040  << "zero starting solution: "
1041  << (ZeroStartingSolution_ ? "true" : "false") << endl;
1042  }
1043 }
1044 
1045 }//namespace Ifpack2
1046 
1047 
1048 // Macro that does explicit template instantiation (ETI) for
1049 // Ifpack2::BlockRelaxation. S, LO, GO, N correspond to the four
1050 // template parameters of Ifpack2::Preconditioner and
1051 // Tpetra::RowMatrix.
1052 //
1053 // We only instantiate for MatrixType = Tpetra::RowMatrix. There's no
1054 // need to instantiate for Tpetra::CrsMatrix too. All Ifpack2
1055 // preconditioners can and should do dynamic casts if they need a type
1056 // more specific than RowMatrix. This keeps build time short and
1057 // library and executable sizes small.
1058 
1059 #ifdef HAVE_IFPACK2_EXPLICIT_INSTANTIATION
1060 
1061 #define IFPACK2_BLOCKRELAXATION_INSTANT(S,LO,GO,N) \
1062  template \
1063  class Ifpack2::BlockRelaxation< \
1064  Tpetra::RowMatrix<S, LO, GO, N>, \
1065  Ifpack2::Container<Tpetra::RowMatrix<S, LO, GO, N> > >;
1066 
1067 #endif // HAVE_IFPACK2_EXPLICIT_INSTANTIATION
1068 
1069 #endif // IFPACK2_BLOCKRELAXATION_DEF_HPP
Ifpack2::BlockRelaxation class declaration.
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
The communicator over which the input matrix is distributed.
Definition: Ifpack2_BlockRelaxation_def.hpp:315
ParameterList & disableRecursiveValidation()
std::string description() const
A one-line description of this object.
Definition: Ifpack2_BlockRelaxation_def.hpp:933
double getApplyTime() const
Returns the time spent in apply().
Definition: Ifpack2_BlockRelaxation_def.hpp:410
T & get(const std::string &name, T def_value)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
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
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
int getNumApply() const
Returns the number of calls to apply().
Definition: Ifpack2_BlockRelaxation_def.hpp:386
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
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
int getNumCompute() const
Returns the number of calls to compute().
Definition: Ifpack2_BlockRelaxation_def.hpp:378
void setParameters(const Teuchos::ParameterList &params)
Sets all the parameters for the preconditioner.
Definition: Ifpack2_BlockRelaxation_def.hpp:166
void compute()
compute the preconditioner for the specified matrix, diagonal perturbation thresholds and relaxation ...
Definition: Ifpack2_BlockRelaxation_def.hpp:655
Block relaxation preconditioners (or smoothers) for Tpetra::RowMatrix and Tpetra::CrsMatrix sparse ma...
Definition: Ifpack2_BlockRelaxation_decl.hpp:83
double getInitializeTime() const
Returns the time spent in initialize().
Definition: Ifpack2_BlockRelaxation_def.hpp:394
Declaration of a user-defined partitioner in which the user can define a partition of the graph...
bool isParameter(const std::string &name) const
T * getRawPtr() const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
double getComputeTime() const
Returns the time spent in compute().
Definition: Ifpack2_BlockRelaxation_def.hpp:402
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
Definition: Ifpack2_BlockRelaxation_def.hpp:417
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
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
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
int getNumInitialize() const
Returns the number of calls to initialize().
Definition: Ifpack2_BlockRelaxation_def.hpp:371
void initialize()
Initialize.
Definition: Ifpack2_BlockRelaxation_def.hpp:544
void validateParameters(ParameterList const &validParamList, int const depth=1000, EValidateUsed const validateUsed=VALIDATE_USED_ENABLED, EValidateDefaults const validateDefaults=VALIDATE_DEFAULTS_ENABLED) const
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:987
MatrixType::local_ordinal_type local_ordinal_type
The type of local indices in the input MatrixType.
Definition: Ifpack2_BlockRelaxation_decl.hpp:101
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:528
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
bool isType(const std::string &name) const
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
Teuchos::RCP< const row_matrix_type > getMatrix() const
The input matrix of this preconditioner&#39;s constructor.
Definition: Ifpack2_BlockRelaxation_def.hpp:329
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
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
BlockRelaxation(const Teuchos::RCP< const row_matrix_type > &Matrix)
Constructor.
Definition: Ifpack2_BlockRelaxation_def.hpp:83
bool is_null() const