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>
57 setMatrix (const Teuchos::RCP<const row_matrix_type>& A)
58 {
59  if (A.getRawPtr () != A_.getRawPtr ()) { // it's a different matrix
60  IsInitialized_ = false;
61  IsComputed_ = false;
62  Partitioner_ = Teuchos::null;
63  Importer_ = Teuchos::null;
64  W_ = Teuchos::null;
65  hasBlockCrsMatrix_ = false;
66 
67  if (! A.is_null ()) {
68  IsParallel_ = (A->getRowMap ()->getComm ()->getSize () > 1);
69  }
70 
71  std::vector<Teuchos::RCP<ContainerType> > emptyVec;
72  std::swap (Containers_, emptyVec);
73  NumLocalBlocks_ = 0;
74 
75  A_ = A;
76  }
77 }
78 
79 template<class MatrixType,class ContainerType>
81 BlockRelaxation (const Teuchos::RCP<const row_matrix_type>& A)
82 : A_ (A),
83  Time_ (Teuchos::rcp (new Teuchos::Time ("Ifpack2::BlockRelaxation"))),
84  OverlapLevel_ (0),
85  PartitionerType_ ("linear"),
86  NumSweeps_ (1),
87  NumLocalBlocks_(0),
88  PrecType_ (Ifpack2::Details::JACOBI),
89  DampingFactor_ (STS::one ()),
90  IsParallel_ (false),
91  ZeroStartingSolution_ (true),
92  DoBackwardGS_ (false),
93  IsInitialized_ (false),
94  IsComputed_ (false),
95  NumInitialize_ (0),
96  NumCompute_ (0),
97  NumApply_ (0),
98  InitializeTime_ (0.0),
99  ComputeTime_ (0.0),
100  ApplyTime_ (0.0),
101  ComputeFlops_ (0.0),
102  ApplyFlops_ (0.0),
103  NumMyRows_ (0),
104  NumGlobalRows_ (0),
105  NumGlobalNonzeros_ (0),
106  hasBlockCrsMatrix_ (false)
107 {}
108 
109 template<class MatrixType,class ContainerType>
112 {}
113 
114 template<class MatrixType,class ContainerType>
115 void
117 setParameters (const Teuchos::ParameterList& List)
118 {
119  Teuchos::ParameterList validparams;
120  Ifpack2::getValidParameters (validparams);
121  List.validateParameters (validparams);
122 
123  std::string PT;
124  if (PrecType_ == Ifpack2::Details::JACOBI) {
125  PT = "Jacobi";
126  } else if (PrecType_ == Ifpack2::Details::GS) {
127  PT = "Gauss-Seidel";
128  } else if (PrecType_ == Ifpack2::Details::SGS) {
129  PT = "Symmetric Gauss-Seidel";
130  }
131 
132  Ifpack2::getParameter (List, "relaxation: type", PT);
133 
134  if (PT == "Jacobi") {
135  PrecType_ = Ifpack2::Details::JACOBI;
136  }
137  else if (PT == "Gauss-Seidel") {
138  PrecType_ = Ifpack2::Details::GS;
139  }
140  else if (PT == "Symmetric Gauss-Seidel") {
141  PrecType_ = Ifpack2::Details::SGS;
142  }
143  else {
144  TEUCHOS_TEST_FOR_EXCEPTION(
145  true, std::invalid_argument, "Ifpack2::BlockRelaxation::setParameters: "
146  "Invalid parameter value \"" << PT << "\" for parameter \"relaxation: "
147  "type\".");
148  }
149 
150  Ifpack2::getParameter (List, "relaxation: sweeps",NumSweeps_);
151  Ifpack2::getParameter (List, "relaxation: damping factor", DampingFactor_);
152  Ifpack2::getParameter (List, "relaxation: zero starting solution", ZeroStartingSolution_);
153  Ifpack2::getParameter (List, "relaxation: backward mode",DoBackwardGS_);
154  Ifpack2::getParameter (List, "partitioner: type",PartitionerType_);
155  Ifpack2::getParameter (List, "partitioner: local parts",NumLocalBlocks_);
156  Ifpack2::getParameter (List, "partitioner: overlap",OverlapLevel_);
157 
158  // check parameters
159  if (PrecType_ != Ifpack2::Details::JACOBI) {
160  OverlapLevel_ = 0;
161  }
162  if (NumLocalBlocks_ < 0) {
163  NumLocalBlocks_ = A_->getNodeNumRows() / (-NumLocalBlocks_);
164  }
165  // other checks are performed in Partitioner_
166 
167  // NTS: Sanity check to be removed at a later date when Backward mode is enabled
168  TEUCHOS_TEST_FOR_EXCEPTION(
169  DoBackwardGS_, std::runtime_error,
170  "Ifpack2::BlockRelaxation:setParameters: Setting the \"relaxation: "
171  "backward mode\" parameter to true is not yet supported.");
172 
173  // copy the list as each subblock's constructor will
174  // require it later
175  List_ = List;
176 }
177 
178 template<class MatrixType,class ContainerType>
179 Teuchos::RCP<const Teuchos::Comm<int> >
181 {
182  TEUCHOS_TEST_FOR_EXCEPTION
183  (A_.is_null (), std::runtime_error, "Ifpack2::BlockRelaxation::getComm: "
184  "The matrix is null. You must call setMatrix() with a nonnull matrix "
185  "before you may call this method.");
186  return A_->getComm ();
187 }
188 
189 template<class MatrixType,class ContainerType>
190 Teuchos::RCP<const Tpetra::RowMatrix<typename MatrixType::scalar_type,
191  typename MatrixType::local_ordinal_type,
192  typename MatrixType::global_ordinal_type,
193  typename MatrixType::node_type> >
195  return A_;
196 }
197 
198 template<class MatrixType,class ContainerType>
199 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,
200  typename MatrixType::global_ordinal_type,
201  typename MatrixType::node_type> >
203 getDomainMap () const
204 {
205  TEUCHOS_TEST_FOR_EXCEPTION
206  (A_.is_null (), std::runtime_error, "Ifpack2::BlockRelaxation::"
207  "getDomainMap: The matrix is null. You must call setMatrix() with a "
208  "nonnull matrix before you may call this method.");
209  return A_->getDomainMap ();
210 }
211 
212 template<class MatrixType,class ContainerType>
213 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,
214  typename MatrixType::global_ordinal_type,
215  typename MatrixType::node_type> >
217 getRangeMap () const
218 {
219  TEUCHOS_TEST_FOR_EXCEPTION
220  (A_.is_null (), std::runtime_error, "Ifpack2::BlockRelaxation::"
221  "getRangeMap: The matrix is null. You must call setMatrix() with a "
222  "nonnull matrix before you may call this method.");
223  return A_->getRangeMap ();
224 }
225 
226 template<class MatrixType,class ContainerType>
227 bool
229 hasTransposeApply () const {
230  return true;
231 }
232 
233 template<class MatrixType,class ContainerType>
234 int
237  return NumInitialize_;
238 }
239 
240 template<class MatrixType,class ContainerType>
241 int
244 {
245  return NumCompute_;
246 }
247 
248 template<class MatrixType,class ContainerType>
249 int
251 getNumApply () const
252 {
253  return NumApply_;
254 }
255 
256 template<class MatrixType,class ContainerType>
257 double
260 {
261  return InitializeTime_;
262 }
263 
264 template<class MatrixType,class ContainerType>
265 double
268 {
269  return ComputeTime_;
270 }
271 
272 template<class MatrixType,class ContainerType>
273 double
275 getApplyTime () const
276 {
277  return ApplyTime_;
278 }
279 
280 template<class MatrixType,class ContainerType>
281 double
284 {
285  return ComputeFlops_;
286 }
287 
288 template<class MatrixType,class ContainerType>
289 double
292 {
293  return ApplyFlops_;
294 }
295 
296 template<class MatrixType,class ContainerType>
297 void
299 apply (const Tpetra::MultiVector<typename MatrixType::scalar_type,
300  typename MatrixType::local_ordinal_type,
301  typename MatrixType::global_ordinal_type,
302  typename MatrixType::node_type>& X,
303  Tpetra::MultiVector<typename MatrixType::scalar_type,
304  typename MatrixType::local_ordinal_type,
305  typename MatrixType::global_ordinal_type,
306  typename MatrixType::node_type>& Y,
307  Teuchos::ETransp mode,
308  scalar_type alpha,
309  scalar_type beta) const
310 {
311  TEUCHOS_TEST_FOR_EXCEPTION
312  (A_.is_null (), std::runtime_error, "Ifpack2::BlockRelaxation::apply: "
313  "The matrix is null. You must call setMatrix() with a nonnull matrix, "
314  "then call initialize() and compute() (in that order), before you may "
315  "call this method.");
316  TEUCHOS_TEST_FOR_EXCEPTION(
317  ! isComputed (), std::runtime_error, "Ifpack2::BlockRelaxation::apply: "
318  "isComputed() must be true prior to calling apply.");
319  TEUCHOS_TEST_FOR_EXCEPTION(
320  X.getNumVectors () != Y.getNumVectors (), std::invalid_argument,
321  "Ifpack2::BlockRelaxation::apply: X.getNumVectors() = "
322  << X.getNumVectors () << " != Y.getNumVectors() = "
323  << Y.getNumVectors () << ".");
324  TEUCHOS_TEST_FOR_EXCEPTION(
325  mode != Teuchos::NO_TRANS, std::logic_error, "Ifpack2::BlockRelaxation::"
326  "apply: This method currently only implements the case mode == Teuchos::"
327  "NO_TRANS. Transposed modes are not currently supported.");
328  TEUCHOS_TEST_FOR_EXCEPTION(
329  alpha != Teuchos::ScalarTraits<scalar_type>::one (), std::logic_error,
330  "Ifpack2::BlockRelaxation::apply: This method currently only implements "
331  "the case alpha == 1. You specified alpha = " << alpha << ".");
332  TEUCHOS_TEST_FOR_EXCEPTION(
333  beta != Teuchos::ScalarTraits<scalar_type>::zero (), std::logic_error,
334  "Ifpack2::BlockRelaxation::apply: This method currently only implements "
335  "the case beta == 0. You specified beta = " << beta << ".");
336 
337  Time_->start(true);
338 
339  // If X and Y are pointing to the same memory location,
340  // we need to create an auxiliary vector, Xcopy
341  Teuchos::RCP<const MV> X_copy;
342  {
343  auto X_lcl_host = X.template getLocalView<Kokkos::HostSpace> ();
344  auto Y_lcl_host = Y.template getLocalView<Kokkos::HostSpace> ();
345  if (X_lcl_host.ptr_on_device () == Y_lcl_host.ptr_on_device ()) {
346  X_copy = rcp (new MV (X, Teuchos::Copy));
347  } else {
348  X_copy = rcpFromRef (X);
349  }
350  }
351 
352  if (ZeroStartingSolution_) {
353  Y.putScalar (STS::zero ());
354  }
355 
356  // Flops are updated in each of the following.
357  switch (PrecType_) {
358  case Ifpack2::Details::JACOBI:
359  ApplyInverseJacobi(*X_copy,Y);
360  break;
361  case Ifpack2::Details::GS:
362  ApplyInverseGS(*X_copy,Y);
363  break;
364  case Ifpack2::Details::SGS:
365  ApplyInverseSGS(*X_copy,Y);
366  break;
367  default:
368  TEUCHOS_TEST_FOR_EXCEPTION
369  (true, std::logic_error, "Ifpack2::BlockRelaxation::apply: Invalid "
370  "PrecType_ enum value " << PrecType_ << ". Valid values are Ifpack2::"
371  "Details::JACOBI = " << Ifpack2::Details::JACOBI << ", Ifpack2::Details"
372  "::GS = " << Ifpack2::Details::GS << ", and Ifpack2::Details::SGS = "
373  << Ifpack2::Details::SGS << ". Please report this bug to the Ifpack2 "
374  "developers.");
375  }
376 
377  ++NumApply_;
378  Time_->stop();
379  ApplyTime_ += Time_->totalElapsedTime();
380 }
381 
382 template<class MatrixType,class ContainerType>
383 void
385 applyMat (const Tpetra::MultiVector<typename MatrixType::scalar_type,
386  typename MatrixType::local_ordinal_type,
387  typename MatrixType::global_ordinal_type,
388  typename MatrixType::node_type>& X,
389  Tpetra::MultiVector<typename MatrixType::scalar_type,
390  typename MatrixType::local_ordinal_type,
391  typename MatrixType::global_ordinal_type,
392  typename MatrixType::node_type>& Y,
393  Teuchos::ETransp mode) const
394 {
395  A_->apply (X, Y, mode);
396 }
397 
398 template<class MatrixType,class ContainerType>
399 void
402 {
403  using Teuchos::rcp;
404  typedef Tpetra::RowGraph<local_ordinal_type, global_ordinal_type, node_type>
405  row_graph_type;
406 
407  TEUCHOS_TEST_FOR_EXCEPTION
408  (A_.is_null (), std::runtime_error, "Ifpack2::BlockRelaxation::initialize: "
409  "The matrix is null. You must call setMatrix() with a nonnull matrix "
410  "before you may call this method.");
411 
412  // Check whether we have a BlockCrsMatrix
413  Teuchos::RCP<const block_crs_matrix_type> A_bcrs =
414  Teuchos::rcp_dynamic_cast<const block_crs_matrix_type> (A_);
415  if (A_bcrs.is_null ()) {
416  hasBlockCrsMatrix_ = false;
417  }
418  else {
419  hasBlockCrsMatrix_ = true;
420  }
421 
422  IsInitialized_ = false;
423  Time_->start (true);
424 
425  NumMyRows_ = A_->getNodeNumRows ();
426  NumGlobalRows_ = A_->getGlobalNumRows ();
427  NumGlobalNonzeros_ = A_->getGlobalNumEntries ();
428 
429  // NTS: Will need to add support for Zoltan2 partitions later Also,
430  // will need a better way of handling the Graph typing issue.
431  // Especially with ordinal types w.r.t the container.
432 
433  if (PartitionerType_ == "linear") {
434  Partitioner_ =
435  rcp (new Ifpack2::LinearPartitioner<row_graph_type> (A_->getGraph ()));
436  } else if (PartitionerType_ == "line") {
437  Partitioner_ =
439  } else if (PartitionerType_ == "user") {
440  Partitioner_ =
441  rcp (new Ifpack2::Details::UserPartitioner<row_graph_type> (A_->getGraph () ) );
442  } else {
443  // We should have checked for this in setParameters(), so it's a
444  // logic_error, not an invalid_argument or runtime_error.
445  TEUCHOS_TEST_FOR_EXCEPTION
446  (true, std::logic_error, "Ifpack2::BlockRelaxation::initialize: Unknown "
447  "partitioner type " << PartitionerType_ << ". Valid values are "
448  "\"linear\", \"line\", and \"user\".");
449  }
450 
451  // need to partition the graph of A
452  Partitioner_->setParameters (List_);
453  Partitioner_->compute ();
454 
455  // get actual number of partitions
456  NumLocalBlocks_ = Partitioner_->numLocalParts ();
457 
458  // Note: Unlike Ifpack, we'll punt on computing W_ until compute(), which is where
459  // we assume that the type of relaxation has been chosen.
460 
461  if (A_->getComm()->getSize() != 1) {
462  IsParallel_ = true;
463  } else {
464  IsParallel_ = false;
465  }
466 
467  ++NumInitialize_;
468  Time_->stop ();
469  InitializeTime_ += Time_->totalElapsedTime ();
470  IsInitialized_ = true;
471 }
472 
473 
474 template<class MatrixType,class ContainerType>
475 void
478 {
479  typedef Tpetra::Map<local_ordinal_type,global_ordinal_type,node_type> map_type;
480  typedef Tpetra::Import<local_ordinal_type,global_ordinal_type, node_type> import_type;
481  using Teuchos::RCP;
482  using Teuchos::rcp;
483  using Teuchos::Array;
484  using Teuchos::ArrayView;
485 
486  Time_->start (true);
487 
488  // reset values
489  IsComputed_ = false;
490 
491  // Extract the submatrices
492  ExtractSubmatrices ();
493 
494  // Compute the weight vector if we're doing overlapped Jacobi (and
495  // only if we're doing overlapped Jacobi).
496  TEUCHOS_TEST_FOR_EXCEPTION
497  (PrecType_ == Ifpack2::Details::JACOBI && OverlapLevel_ > 0, std::runtime_error,
498  "Ifpack2::BlockRelaxation::computeBlockCrs: "
499  "We do not support overlapped Jacobi yet for Tpetra::BlockCrsMatrix. Sorry!");
500 
501  // We need to import data from external processors. Here I create a
502  // Tpetra::Import object if needed (stealing from A_ if possible)
503  // Marzio's comment:
504  // Note that I am doing some strange stuff to set the components of Y
505  // from Y2 (to save some time).
506  //
507  if (IsParallel_ && (PrecType_ == Ifpack2::Details::GS ||
508  PrecType_ == Ifpack2::Details::SGS)) {
509  // Get the block size
510  RCP<const block_crs_matrix_type> A_bcrs =
511  Teuchos::rcp_dynamic_cast<const block_crs_matrix_type> (A_);
512  int bs = A_bcrs->getBlockSize();
513 
514  // Get the maps describing the block data
515  // Each "element" is actually a block of numbers
516  RCP<const map_type> oldDomainMap = A_->getDomainMap();
517  RCP<const map_type> oldColMap = A_->getColMap();
518 
519  // Because A is a block CRS matrix, import will not do what you think it does
520  // We have to construct the correct maps for it
521  global_size_t numGlobalElements = oldColMap->getGlobalNumElements()*bs;
522  global_ordinal_type indexBase = oldColMap->getIndexBase();
523  RCP<const Teuchos::Comm<int> >comm = oldColMap->getComm();
524  ArrayView<const global_ordinal_type> oldColElements = oldColMap->getNodeElementList();
525  Array<global_ordinal_type> newColElements(bs*oldColElements.size());
526 
527  for(int i=0; i<oldColElements.size(); i++) {
528  for(int j=0; j<bs; j++) {
529  newColElements[i*bs+j] = oldColElements[i]*bs+j;
530  }
531  }
532  RCP<map_type> colMap = rcp(new map_type(numGlobalElements,newColElements,indexBase,comm));
533 
534  // Create the importer
535  Importer_ = rcp (new import_type (oldDomainMap, colMap));
536  }
537 
538  ++NumCompute_;
539  Time_->stop ();
540  ComputeTime_ += Time_->totalElapsedTime();
541  IsComputed_ = true;
542 }
543 
544 
545 template<class MatrixType,class ContainerType>
546 void
549 {
550  using Teuchos::rcp;
551  typedef Tpetra::Vector<scalar_type,
553  typedef Tpetra::Import<local_ordinal_type,
554  global_ordinal_type, node_type> import_type;
555 
556  TEUCHOS_TEST_FOR_EXCEPTION
557  (A_.is_null (), std::runtime_error, "Ifpack2::BlockRelaxation::compute: "
558  "The matrix is null. You must call setMatrix() with a nonnull matrix, "
559  "then call initialize(), before you may call this method.");
560 
561  // We should have checked for this in setParameters(), so it's a
562  // logic_error, not an invalid_argument or runtime_error.
563  TEUCHOS_TEST_FOR_EXCEPTION
564  (NumSweeps_ < 0, std::logic_error, "Ifpack2::BlockRelaxation::compute: "
565  "NumSweeps_ = " << NumSweeps_ << " < 0.");
566 
567  if (! isInitialized ()) {
568  initialize ();
569  }
570 
571  if (hasBlockCrsMatrix_) {
572  computeBlockCrs ();
573  return;
574  }
575 
576  Time_->start (true);
577 
578  // reset values
579  IsComputed_ = false;
580 
581  // Extract the submatrices
582  ExtractSubmatrices ();
583 
584  // Compute the weight vector if we're doing overlapped Jacobi (and
585  // only if we're doing overlapped Jacobi).
586  if (PrecType_ == Ifpack2::Details::JACOBI && OverlapLevel_ > 0) {
587  // weight of each vertex
588  W_ = rcp (new vector_type (A_->getRowMap ()));
589  W_->putScalar (STS::zero ());
590  Teuchos::ArrayRCP<scalar_type > w_ptr = W_->getDataNonConst(0);
591 
592  for (local_ordinal_type i = 0 ; i < NumLocalBlocks_ ; ++i) {
593  for (size_t j = 0 ; j < Partitioner_->numRowsInPart(i) ; ++j) {
594  // FIXME (mfh 12 Sep 2014) Should this really be int?
595  // Perhaps it should be local_ordinal_type instead.
596  int LID = (*Partitioner_)(i,j);
597  w_ptr[LID]+= STS::one();
598  }
599  }
600  W_->reciprocal (*W_);
601  }
602 
603  // We need to import data from external processors. Here I create a
604  // Tpetra::Import object if needed (stealing from A_ if possible)
605  // Marzio's comment:
606  // Note that I am doing some strange stuff to set the components of Y
607  // from Y2 (to save some time).
608  //
609  if (IsParallel_ && (PrecType_ == Ifpack2::Details::GS ||
610  PrecType_ == Ifpack2::Details::SGS)) {
611  Importer_ = A_->getGraph ()->getImporter ();
612  if (Importer_.is_null ()) {
613  Importer_ = rcp (new import_type (A_->getDomainMap (), A_->getColMap ()));
614  }
615  }
616 
617  ++NumCompute_;
618  Time_->stop ();
619  ComputeTime_ += Time_->totalElapsedTime();
620  IsComputed_ = true;
621 }
622 
623 template<class MatrixType,class ContainerType>
624 void
627 {
628  TEUCHOS_TEST_FOR_EXCEPTION
629  (Partitioner_.is_null (), std::runtime_error, "Ifpack2::BlockRelaxation::"
630  "ExtractSubmatrices: Partitioner object is null.");
631 
632  NumLocalBlocks_ = Partitioner_->numLocalParts ();
633  Containers_.resize (NumLocalBlocks_);
634 
635  for (local_ordinal_type i = 0; i < NumLocalBlocks_; ++i) {
636  const size_t numRows = Partitioner_->numRowsInPart (i);
637 
638  // Extract a list of the indices of each partitioner row.
639  Teuchos::Array<local_ordinal_type> localRows (numRows);
640  for (size_t j = 0; j < numRows; ++j) {
641  localRows[j] = (*Partitioner_) (i,j);
642  }
643  if(numRows>1 || hasBlockCrsMatrix_) { // only do for non-singletons
644  Containers_[i] = Teuchos::rcp (new ContainerType (A_, localRows ()));
645  Containers_[i]->setParameters (List_);
646  Containers_[i]->initialize ();
647  Containers_[i]->compute ();
648  }
649  }
650 }
651 
652 template<class MatrixType,class ContainerType>
653 void
655 ApplyInverseJacobi (const MV& X, MV& Y) const
656 {
657  const size_t NumVectors = X.getNumVectors ();
658  MV AY (Y.getMap (), NumVectors);
659 
660  // Initial matvec not needed
661  int starting_iteration = 0;
662  if (ZeroStartingSolution_) {
663  DoJacobi (X, Y);
664  starting_iteration = 1;
665  }
666 
667  const scalar_type ONE = STS::one ();
668  for (int j = starting_iteration; j < NumSweeps_; ++j) {
669  applyMat (Y, AY);
670  AY.update (ONE, X, -ONE);
671  DoJacobi (AY, Y);
672 
673  // Flops for matrix apply & update
674  ApplyFlops_ += NumVectors * (2 * NumGlobalNonzeros_ + 2 * NumGlobalRows_);
675  }
676 
677 }
678 
679 template<class MatrixType,class ContainerType>
680 void
682 DoJacobi (const MV& X, MV& Y) const
683 {
684  const size_t NumVectors = X.getNumVectors ();
685  const scalar_type one = STS::one ();
686  // Note: Flop counts copied naively from Ifpack.
687 
688  if (OverlapLevel_ == 0) {
689  // Non-overlapping Jacobi
690  for (local_ordinal_type i = 0; i < NumLocalBlocks_; ++i) {
691  // may happen that a partition is empty
692  if( Partitioner_->numRowsInPart (i) != 1 || hasBlockCrsMatrix_) {
693  if(Containers_[i]->getNumRows () == 0 ) continue;
694  Containers_[i]->apply (X, Y, Teuchos::NO_TRANS, DampingFactor_, one);
695  ApplyFlops_ += NumVectors * 2 * NumGlobalRows_;
696  }
697  else { // singleton, can't access Containers_[i] as it was never filled and may be null.
698  local_ordinal_type LRID = (*Partitioner_)(i,0); // by definition, a singleton 1 row in block.
699  getMatDiag();
700  Teuchos::ArrayRCP< const scalar_type > Diag = DiagRCP_->getData();
701  scalar_type d = Diag[LRID];
702  for(unsigned int nv = 0;nv < NumVectors ; ++nv ) {
703  Teuchos::ArrayRCP< const scalar_type > xRCP = X.getData(nv);
704  scalar_type x = xRCP[LRID];
705  Teuchos::ArrayRCP< scalar_type > yRCP = Y.getDataNonConst(nv);
706 
707  scalar_type newy= x/d;
708  yRCP[LRID]= newy;
709  }
710  }
711  }
712  }
713  else {
714  // Overlapping Jacobi
715  for (local_ordinal_type i = 0 ; i < NumLocalBlocks_ ; i++) {
716  // may happen that a partition is empty
717  if(Containers_[i]->getNumRows() == 0) continue;
718  if ( Partitioner_->numRowsInPart (i) != 1 ) {
719  try {
720  Containers_[i]->weightedApply(X,Y,*W_,Teuchos::NO_TRANS,DampingFactor_,one);
721  } catch (std::exception& e) {
722  std::cerr << "BlockRelaxation::DoJacobi: Containers_[" << i
723  << "]->weightedApply() threw an exception: " << e.what ()
724  << std::endl;
725  throw;
726  }
727  } // end Partitioner_->numRowsInPart (i) != 1
728 
729  // NOTE: do not count (for simplicity) the flops due to overlapping rows
730  ApplyFlops_ += NumVectors * 4 * NumGlobalRows_;
731  }
732  }
733 }
734 
735 template<class MatrixType,class ContainerType>
736 void
738 ApplyInverseGS (const MV& X, MV& Y) const
739 {
740  MV Xcopy (X, Teuchos::Copy);
741  for (int j = 0; j < NumSweeps_; ++j) {
742  DoGaussSeidel (Xcopy, Y);
743  if (j != NumSweeps_ - 1) {
744  Tpetra::deep_copy (Xcopy, X);
745  }
746  }
747 }
748 
749 template<class MatrixType,class ContainerType>
750 void
752 DoGaussSeidel (MV& X, MV& Y) const
753 {
754  using Teuchos::Array;
755  using Teuchos::ArrayRCP;
756  using Teuchos::ArrayView;
757  using Teuchos::RCP;
758  using Teuchos::rcp;
759  using Teuchos::rcpFromRef;
760 
761  // Note: Flop counts copied naively from Ifpack.
762 
763  const scalar_type one = STS::one ();
764  const size_t Length = A_->getNodeMaxNumRowEntries();
765  const size_t NumVectors = X.getNumVectors();
766  Array<scalar_type> Values;
767  Array<local_ordinal_type> Indices;
768  Indices.resize (Length);
769 
770  if(hasBlockCrsMatrix_)
771  {
772  Teuchos::RCP<const block_crs_matrix_type> A_bcrs =
773  Teuchos::rcp_dynamic_cast<const block_crs_matrix_type> (A_);
774  int bs = A_bcrs->getBlockSize();
775  Values.resize (bs*bs*Length);
776  }
777  else
778  Values.resize (Length);
779 
780  // an additonal vector is needed by parallel computations
781  // (note that applications through Ifpack2_AdditiveSchwarz
782  // are always seen are serial)
783  RCP<MV> Y2;
784  if (IsParallel_) {
785  Y2 = rcp (new MV (Importer_->getTargetMap (), NumVectors));
786  } else {
787  Y2 = rcpFromRef (Y);
788  }
789 
790  // I think I decided I need two extra vectors:
791  // One to store the sum of the corrections (initialized to zero)
792  // One to store the temporary residual (doesn't matter if it is zeroed or not)
793  // My apologies for making the names clear and meaningful. (X=RHS, Y=guess?! Nice.)
794  MV Residual (X.getMap (), NumVectors, false);
795 
796  ArrayRCP<ArrayRCP<scalar_type> > x_ptr = X.get2dViewNonConst();
797  ArrayRCP<ArrayRCP<scalar_type> > y_ptr = Y.get2dViewNonConst();
798  ArrayRCP<ArrayRCP<scalar_type> > y2_ptr = Y2->get2dViewNonConst();
799  ArrayRCP<ArrayRCP<scalar_type> > residual_ptr = Residual.get2dViewNonConst();
800 
801  // data exchange is here, once per sweep
802  if (IsParallel_) Y2->doImport(Y,*Importer_,Tpetra::INSERT);
803 
804  for (local_ordinal_type i = 0; i < NumLocalBlocks_; ++i) {
805  if( Partitioner_->numRowsInPart (i) != 1 || hasBlockCrsMatrix_) {
806  if (Containers_[i]->getNumRows () == 0) continue;
807  // update from previous block
808  ArrayView<const local_ordinal_type> localRows =
809  Containers_[i]->getLocalRows ();
810  const size_t localNumRows = Containers_[i]->getNumRows ();
811  for (size_t j = 0; j < localNumRows; ++j) {
812  const local_ordinal_type LID = localRows[j]; // Containers_[i]->ID (j);
813  size_t NumEntries;
814  A_->getLocalRowCopy (LID, Indices (), Values (), NumEntries);
815 
816  for (size_t m = 0; m < NumVectors; ++m) {
817  ArrayView<const scalar_type> x_local = (x_ptr())[m]();
818  ArrayView<scalar_type> y2_local = (y2_ptr())[m]();
819  ArrayView<scalar_type> r_local = (residual_ptr())[m]();
820 
821  if(hasBlockCrsMatrix_) {
822  Teuchos::RCP<const block_crs_matrix_type> A_bcrs =
823  Teuchos::rcp_dynamic_cast<const block_crs_matrix_type> (A_);
824  int bs = A_bcrs->getBlockSize();
825  for (int localR = 0; localR < bs; localR++)
826  r_local[LID*bs+localR] = x_local[LID*bs+localR];
827  for (size_t k = 0; k < NumEntries; ++k) {
828  const local_ordinal_type col = Indices[k];
829  for (int localR = 0; localR < bs; localR++) {
830  for(int localC = 0; localC < bs; localC++) {
831 // std::cerr << "r[" << LID*bs+localR << "] -= Values["
832 // << k*bs*bs+localR+localC*bs << "] * y2["
833 // << col*bs+localC << "]\n";
834  r_local[LID*bs+localR] -= Values[k*bs*bs+localR+localC*bs] * y2_local[col*bs+localC];
835  }
836  }
837  }
838  }
839  else {
840  r_local[LID] = x_local[LID];
841  for (size_t k = 0; k < NumEntries; ++k) {
842  const local_ordinal_type col = Indices[k];
843  r_local[LID] -= Values[k] * y2_local[col];
844  }
845  }
846  }
847  }
848 
849 // Teuchos::RCP<Teuchos::FancyOStream> wrappedStream = Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cout));
850 // Residual.describe (*wrappedStream, Teuchos::VERB_EXTREME);
851 
852  // solve with this block
853  //
854  // Note: I'm abusing the ordering information, knowing that X/Y
855  // and Y2 have the same ordering for on-proc unknowns.
856  //
857  // Note: Add flop counts for inverse apply
858  Containers_[i]->apply (Residual, *Y2, Teuchos::NO_TRANS,
859  DampingFactor_,one);
860 
861 // Y2->describe (*wrappedStream, Teuchos::VERB_EXTREME);
862 
863  // operations for all getrow's
864  ApplyFlops_ += NumVectors * (2 * NumGlobalNonzeros_ + 2 * NumGlobalRows_);
865  }
866  else { // singleton, can't access Containers_[i] as it was never filled and may be null.
867  // a singleton calculation is exact, all residuals should be zero.
868  local_ordinal_type LRID = (*Partitioner_)(i,0); // by definition, a singleton 1 row in block.
869  getMatDiag();
870  Teuchos::ArrayRCP< const scalar_type > Diag = DiagRCP_->getData();
871  scalar_type d = Diag[LRID];
872  ArrayRCP<ArrayRCP<scalar_type> > y2_ptr2 = Y2->get2dViewNonConst();
873  for(unsigned int nv = 0;nv < NumVectors ; ++nv ) {
874  Teuchos::ArrayRCP< const scalar_type > xRCP = X.getData(nv);
875  scalar_type x = xRCP[LRID];
876  ArrayView<scalar_type> y2_local = (y2_ptr2())[nv]();
877  scalar_type newy= x/d;
878  y2_local[LRID]= newy;
879  }
880  } // end else
881  } // end for NumLocalBlocks_
882  // Attention: this is delicate... Not all combinations
883  // of Y2 and Y will always work (tough for ML it should be ok)
884  if (IsParallel_) {
885  if(hasBlockCrsMatrix_) {
886  Teuchos::RCP<const block_crs_matrix_type> A_bcrs =
887  Teuchos::rcp_dynamic_cast<const block_crs_matrix_type> (A_);
888  int bs = A_bcrs->getBlockSize();
889  for (size_t m = 0; m < NumVectors; ++m) {
890  ArrayView<scalar_type> y2_local = (y2_ptr())[m]();
891  ArrayView<scalar_type> y_local = (y_ptr())[m]();
892  for (size_t i = 0; i < NumMyRows_*bs; ++i) {
893  y_local[i] = y2_local[i];
894  }
895  }
896  }
897  else {
898  for (size_t m = 0; m < NumVectors; ++m) {
899  ArrayView<scalar_type> y2_local = (y2_ptr())[m]();
900  ArrayView<scalar_type> y_local = (y_ptr())[m]();
901  for (size_t i = 0; i < NumMyRows_; ++i) {
902  y_local[i] = y2_local[i];
903  }
904  }
905  }
906  }
907 }
908 
909 template<class MatrixType,class ContainerType>
910 void
912 ApplyInverseSGS (const MV& X, MV& Y) const
913 {
914  MV Xcopy (X, Teuchos::Copy);
915  for (int j = 0; j < NumSweeps_; ++j) {
916  DoSGS (Xcopy, Y);
917  if (j != NumSweeps_ - 1) {
918  Tpetra::deep_copy (Xcopy, X);
919  }
920  }
921 }
922 
923 template<class MatrixType,class ContainerType>
924 void
926 DoSGS (MV& X, MV& Y) const
927 {
928  using Teuchos::Array;
929  using Teuchos::ArrayRCP;
930  using Teuchos::ArrayView;
931  using Teuchos::RCP;
932  using Teuchos::rcp;
933  using Teuchos::rcpFromRef;
934 
935  const scalar_type one = STS::one ();
936  const size_t Length = A_->getNodeMaxNumRowEntries();
937  const size_t NumVectors = X.getNumVectors();
938  Array<scalar_type> Values;
939  Array<local_ordinal_type> Indices;
940  Indices.resize(Length);
941 
942  if(hasBlockCrsMatrix_)
943  {
944  Teuchos::RCP<const block_crs_matrix_type> A_bcrs =
945  Teuchos::rcp_dynamic_cast<const block_crs_matrix_type> (A_);
946  int bs = A_bcrs->getBlockSize();
947  Values.resize (bs*bs*Length);
948  }
949  else
950  Values.resize (Length);
951 
952  // an additonal vector is needed by parallel computations
953  // (note that applications through Ifpack2_AdditiveSchwarz
954  // are always seen are serial)
955  RCP<MV> Y2;
956  if (IsParallel_) {
957  Y2 = rcp (new MV (Importer_->getTargetMap (), NumVectors));
958  } else {
959  Y2 = rcpFromRef (Y);
960  }
961 
962  // I think I decided I need two extra vectors:
963  // One to store the sum of the corrections (initialized to zero)
964  // One to store the temporary residual (doesn't matter if it is zeroed or not)
965  // My apologies for making the names clear and meaningful. (X=RHS, Y=guess?! Nice.)
966  MV Residual (X.getMap (), NumVectors, false);
967 
968  ArrayRCP<ArrayRCP<scalar_type> > x_ptr = X.get2dViewNonConst();
969  ArrayRCP<ArrayRCP<scalar_type> > y_ptr = Y.get2dViewNonConst();
970  ArrayRCP<ArrayRCP<scalar_type> > y2_ptr = Y2->get2dViewNonConst();
971  ArrayRCP<ArrayRCP<scalar_type> > residual_ptr = Residual.get2dViewNonConst();
972 
973  // data exchange is here, once per sweep
974  if (IsParallel_) {
975  Y2->doImport (Y, *Importer_, Tpetra::INSERT);
976  }
977 
978  // Forward Sweep
979  for (local_ordinal_type i = 0; i < NumLocalBlocks_; ++i) {
980  if( Partitioner_->numRowsInPart (i) != 1 || hasBlockCrsMatrix_) {
981  if (Containers_[i]->getNumRows () == 0) {
982  continue; // Skip empty partitions
983  }
984  // update from previous block
985  ArrayView<const local_ordinal_type> localRows =
986  Containers_[i]->getLocalRows ();
987  for (size_t j = 0; j < Containers_[i]->getNumRows (); ++j) {
988  const local_ordinal_type LID = localRows[j]; // Containers_[i]->ID (j);
989  size_t NumEntries;
990  A_->getLocalRowCopy (LID, Indices (), Values (), NumEntries);
991 
992  //set tmpresid = initresid - A*correction
993  for (size_t m = 0; m < NumVectors; ++m) {
994  ArrayView<const scalar_type> x_local = (x_ptr())[m]();
995  ArrayView<scalar_type> y2_local = (y2_ptr())[m]();
996  ArrayView<scalar_type> r_local = (residual_ptr())[m]();
997 
998  if(hasBlockCrsMatrix_) {
999  Teuchos::RCP<const block_crs_matrix_type> A_bcrs =
1000  Teuchos::rcp_dynamic_cast<const block_crs_matrix_type> (A_);
1001  int bs = A_bcrs->getBlockSize();
1002  for (int localR = 0; localR < bs; localR++)
1003  r_local[LID*bs+localR] = x_local[LID*bs+localR];
1004  for (size_t k = 0; k < NumEntries; ++k) {
1005  const local_ordinal_type col = Indices[k];
1006  for (int localR = 0; localR < bs; localR++) {
1007  for(int localC = 0; localC < bs; localC++) {
1008  r_local[LID*bs+localR] -= Values[k*bs*bs+localR+localC*bs] * y2_local[col*bs+localC];
1009  }
1010  }
1011  }
1012  }
1013  else {
1014  r_local[LID] = x_local[LID];
1015  for (size_t k = 0 ; k < NumEntries ; k++) {
1016  local_ordinal_type col = Indices[k];
1017  r_local[LID] -= Values[k] * y2_local[col];
1018  }
1019  }
1020  }
1021  }
1022  // solve with this block
1023  //
1024  // Note: I'm abusing the ordering information, knowing that X/Y
1025  // and Y2 have the same ordering for on-proc unknowns.
1026  //
1027  // Note: Add flop counts for inverse apply
1028  Containers_[i]->apply (Residual, *Y2, Teuchos::NO_TRANS,
1029  DampingFactor_, one);
1030 
1031  // operations for all getrow's
1032  ApplyFlops_ += NumVectors * (2 * NumGlobalNonzeros_ + 2 * NumGlobalRows_);
1033 
1034  }
1035  else { // singleton, can't access Containers_[i] as it was never filled and may be null.
1036  local_ordinal_type LRID = (*Partitioner_)(i,0); // by definition, a singleton 1 row in block.
1037  getMatDiag();
1038  Teuchos::ArrayRCP< const scalar_type > Diag = DiagRCP_->getData();
1039  scalar_type d = Diag[LRID];
1040  for(unsigned int nv = 0;nv < NumVectors ; ++nv ) {
1041  Teuchos::ArrayRCP< const scalar_type > xRCP = X.getData(nv);
1042  scalar_type x = xRCP[LRID];
1043  Teuchos::ArrayRCP< scalar_type > yRCP = Y.getDataNonConst(nv);
1044 
1045  scalar_type newy= x/d;
1046  yRCP[LRID]= newy;
1047  }
1048  } // end else
1049  } // end forward sweep over NumLocalBlocks
1050 
1051  // Reverse Sweep
1052  //
1053  // mfh 12 July 2013: The unusual iteration bounds, and the use of
1054  // i-1 rather than i in the loop body, ensure correctness even if
1055  // local_ordinal_type is unsigned. "i = NumLocalBlocks_-1; i >= 0;
1056  // i--" will loop forever if local_ordinal_type is unsigned, because
1057  // unsigned integers are (trivially) always nonnegative.
1058  for (local_ordinal_type i = NumLocalBlocks_; i > 0; --i) {
1059  if( hasBlockCrsMatrix_ || Partitioner_->numRowsInPart (i) != 1 ) {
1060  if (Containers_[i-1]->getNumRows () == 0) continue;
1061 
1062  // update from previous block
1063  ArrayView<const local_ordinal_type> localRows =
1064  Containers_[i-1]->getLocalRows ();
1065  for (size_t j = 0; j < Containers_[i-1]->getNumRows (); ++j) {
1066  const local_ordinal_type LID = localRows[j]; // Containers_[i-1]->ID (j);
1067  size_t NumEntries;
1068  A_->getLocalRowCopy (LID, Indices (), Values (), NumEntries);
1069 
1070  //set tmpresid = initresid - A*correction
1071  for (size_t m = 0; m < NumVectors; ++m) {
1072  ArrayView<const scalar_type> x_local = (x_ptr())[m]();
1073  ArrayView<scalar_type> y2_local = (y2_ptr())[m]();
1074  ArrayView<scalar_type> r_local = (residual_ptr())[m]();
1075 
1076  if(hasBlockCrsMatrix_) {
1077  Teuchos::RCP<const block_crs_matrix_type> A_bcrs =
1078  Teuchos::rcp_dynamic_cast<const block_crs_matrix_type> (A_);
1079  int bs = A_bcrs->getBlockSize();
1080  for (int localR = 0; localR < bs; localR++)
1081  r_local[LID*bs+localR] = x_local[LID*bs+localR];
1082  for (size_t k = 0; k < NumEntries; ++k) {
1083  const local_ordinal_type col = Indices[k];
1084  for (int localR = 0; localR < bs; localR++) {
1085  for(int localC = 0; localC < bs; localC++) {
1086  r_local[LID*bs+localR] -= Values[k*bs*bs+localR+localC*bs] * y2_local[col*bs+localC];
1087  }
1088  }
1089  }
1090  }
1091  else {
1092  r_local [LID] = x_local[LID];
1093  for (size_t k = 0; k < NumEntries; ++k) {
1094  local_ordinal_type col = Indices[k];
1095  r_local[LID] -= Values[k] * y2_local[col];
1096  }
1097  }
1098  }
1099  }
1100 
1101  // solve with this block
1102  //
1103  // Note: I'm abusing the ordering information, knowing that X/Y
1104  // and Y2 have the same ordering for on-proc unknowns.
1105  //
1106  // Note: Add flop counts for inverse apply
1107  Containers_[i-1]->apply (Residual, *Y2, Teuchos::NO_TRANS,
1108  DampingFactor_, one);
1109 
1110  // operations for all getrow's
1111  ApplyFlops_ += NumVectors * (2 * NumGlobalNonzeros_ + 2 * NumGlobalRows_);
1112  } // end Partitioner_->numRowsInPart (i) != 1 ) {
1113  // else do nothing, as by definition with a singleton, the residuals are zero.
1114  } //end reverse sweep
1115 
1116  // Attention: this is delicate... Not all combinations
1117  // of Y2 and Y will always work (though for ML it should be ok)
1118  if (IsParallel_) {
1119  if(hasBlockCrsMatrix_) {
1120  Teuchos::RCP<const block_crs_matrix_type> A_bcrs =
1121  Teuchos::rcp_dynamic_cast<const block_crs_matrix_type> (A_);
1122  int bs = A_bcrs->getBlockSize();
1123  for (size_t m = 0; m < NumVectors; ++m) {
1124  ArrayView<scalar_type> y2_local = (y2_ptr())[m]();
1125  ArrayView<scalar_type> y_local = (y_ptr())[m]();
1126  for (size_t i = 0; i < NumMyRows_*bs; ++i) {
1127  y_local[i] = y2_local[i];
1128  }
1129  }
1130  }
1131  else {
1132  for (size_t m = 0; m < NumVectors; ++m) {
1133  ArrayView<scalar_type> y2_local = (y2_ptr())[m]();
1134  ArrayView<scalar_type> y_local = (y_ptr())[m]();
1135  for (size_t i = 0; i < NumMyRows_; ++i) {
1136  y_local[i] = y2_local[i];
1137  }
1138  }
1139  }
1140  }
1141 }
1142 
1143 template<class MatrixType, class ContainerType>
1144 std::string
1146 description () const
1147 {
1148  std::ostringstream out;
1149 
1150  // Output is a valid YAML dictionary in flow style. If you don't
1151  // like everything on a single line, you should call describe()
1152  // instead.
1153  out << "\"Ifpack2::BlockRelaxation\": {";
1154  if (this->getObjectLabel () != "") {
1155  out << "Label: \"" << this->getObjectLabel () << "\", ";
1156  }
1157  out << "Initialized: " << (isInitialized () ? "true" : "false") << ", ";
1158  out << "Computed: " << (isComputed () ? "true" : "false") << ", ";
1159 
1160  if (A_.is_null ()) {
1161  out << "Matrix: null, ";
1162  }
1163  else {
1164  out << "Matrix: not null"
1165  << ", Global matrix dimensions: ["
1166  << A_->getGlobalNumRows () << ", " << A_->getGlobalNumCols () << "], ";
1167  }
1168 
1169  // It's useful to print this instance's relaxation method. If you
1170  // want more info than that, call describe() instead.
1171  out << "\"relaxation: type\": ";
1172  if (PrecType_ == Ifpack2::Details::JACOBI) {
1173  out << "Block Jacobi";
1174  } else if (PrecType_ == Ifpack2::Details::GS) {
1175  out << "Block Gauss-Seidel";
1176  } else if (PrecType_ == Ifpack2::Details::SGS) {
1177  out << "Block Symmetric Gauss-Seidel";
1178  } else {
1179  out << "INVALID";
1180  }
1181 
1182  out << ", " << "sweeps: " << NumSweeps_ << ", "
1183  << "damping factor: " << DampingFactor_ << ", ";
1184 
1185  out << "}";
1186  return out.str();
1187 }
1188 
1189 template<class MatrixType,class ContainerType>
1190 void
1192 describe (Teuchos::FancyOStream& out,
1193  const Teuchos::EVerbosityLevel verbLevel) const
1194 {
1195  using std::endl;
1196  using std::setw;
1197  using Teuchos::TypeNameTraits;
1198  using Teuchos::VERB_DEFAULT;
1199  using Teuchos::VERB_NONE;
1200  using Teuchos::VERB_LOW;
1201  using Teuchos::VERB_MEDIUM;
1202  using Teuchos::VERB_HIGH;
1203  using Teuchos::VERB_EXTREME;
1204 
1205  Teuchos::EVerbosityLevel vl = verbLevel;
1206  if (vl == VERB_DEFAULT) vl = VERB_LOW;
1207  const int myImageID = A_->getComm()->getRank();
1208 
1209  // Convention requires that describe() begin with a tab.
1210  Teuchos::OSTab tab (out);
1211 
1212  // none: print nothing
1213  // low: print O(1) info from node 0
1214  // medium:
1215  // high:
1216  // extreme:
1217  if (vl != VERB_NONE && myImageID == 0) {
1218  out << "Ifpack2::BlockRelaxation<"
1219  << TypeNameTraits<MatrixType>::name () << ", "
1220  << TypeNameTraits<ContainerType>::name () << " >:";
1221  Teuchos::OSTab tab1 (out);
1222 
1223  if (this->getObjectLabel () != "") {
1224  out << "label: \"" << this->getObjectLabel () << "\"" << endl;
1225  }
1226  out << "initialized: " << (isInitialized () ? "true" : "false") << endl
1227  << "computed: " << (isComputed () ? "true" : "false") << endl;
1228 
1229  std::string precType;
1230  if (PrecType_ == Ifpack2::Details::JACOBI) {
1231  precType = "Block Jacobi";
1232  } else if (PrecType_ == Ifpack2::Details::GS) {
1233  precType = "Block Gauss-Seidel";
1234  } else if (PrecType_ == Ifpack2::Details::SGS) {
1235  precType = "Block Symmetric Gauss-Seidel";
1236  }
1237  out << "type: " << precType << endl
1238  << "global number of rows: " << A_->getGlobalNumRows () << endl
1239  << "global number of columns: " << A_->getGlobalNumCols () << endl
1240  << "number of sweeps: " << NumSweeps_ << endl
1241  << "damping factor: " << DampingFactor_ << endl
1242  << "backwards mode: "
1243  << ((PrecType_ == Ifpack2::Details::GS && DoBackwardGS_) ? "true" : "false")
1244  << endl
1245  << "zero starting solution: "
1246  << (ZeroStartingSolution_ ? "true" : "false") << endl;
1247 
1248  out << "===============================================================================" << endl;
1249  out << "Phase # calls Total Time (s) Total MFlops MFlops/s " << endl;
1250  out << "------------ ------- --------------- --------------- ---------------" << endl;
1251  out << setw(12) << "initialize()" << setw(5) << getNumInitialize() << " " << setw(15) << getInitializeTime() << endl;
1252  out << setw(12) << "compute()" << setw(5) << getNumCompute() << " " << setw(15) << getComputeTime() << " "
1253  << setw(15) << getComputeFlops() << " "
1254  << setw(15) << (getComputeTime() != 0.0 ? getComputeFlops() / getComputeTime() * 1.0e-6 : 0.0) << endl;
1255  out << setw(12) << "apply()" << setw(5) << getNumApply() << " " << setw(15) << getApplyTime() << " "
1256  << setw(15) << getApplyFlops() << " "
1257  << setw(15) << (getApplyTime() != 0.0 ? getApplyFlops() / getApplyTime() * 1.0e-6 : 0.0) << endl;
1258  out << "===============================================================================" << endl;
1259  out << endl;
1260  }
1261 }
1262 
1263 
1264 template<class MatrixType,class ContainerType>
1265 void
1267 {
1268  // TODO amk: Is this map correct for BlockCRSMatrix?
1269  if(DiagRCP_ == Teuchos::null) {
1270  DiagRCP_ = Teuchos::rcp(new vector_type(A_->getDomainMap ()));
1271  A_->getLocalDiagCopy (*DiagRCP_);
1272  }
1273 }
1274 
1275 }//namespace Ifpack2
1276 
1277 
1278 #ifdef HAVE_IFPACK2_EXPLICIT_INSTANTIATION
1279 
1280 // For ETI
1285 #include "Ifpack2_ILUT_decl.hpp"
1286 
1287 // There's no need to instantiate for CrsMatrix too. All Ifpack2
1288 // preconditioners can and should do dynamic casts if they need a type
1289 // more specific than RowMatrix.
1290 
1291 #define IFPACK2_BLOCKRELAXATION_INSTANT(S,LO,GO,N) \
1292  template \
1293  class Ifpack2::BlockRelaxation< \
1294  Tpetra::RowMatrix<S, LO, GO, N>, \
1295  Ifpack2::SparseContainer< \
1296  Tpetra::RowMatrix<S, LO, GO, N>, \
1297  Ifpack2::ILUT< ::Tpetra::RowMatrix<S,LO,GO,N> > > >; \
1298  template \
1299  class Ifpack2::BlockRelaxation< \
1300  Tpetra::RowMatrix<S, LO, GO, N>, \
1301  Ifpack2::DenseContainer< \
1302  Tpetra::RowMatrix<S, LO, GO, N>, \
1303  S > >; \
1304  template \
1305  class Ifpack2::BlockRelaxation< \
1306  Tpetra::RowMatrix<S, LO, GO, N>, \
1307  Ifpack2::TriDiContainer< \
1308  Tpetra::RowMatrix<S, LO, GO, N>, \
1309  S > >; \
1310  template \
1311  class Ifpack2::BlockRelaxation< \
1312  Tpetra::RowMatrix<S, LO, GO, N>, \
1313  Ifpack2::BandedContainer< \
1314  Tpetra::RowMatrix<S, LO, GO, N>, \
1315  S > >;
1316 
1317 #endif // HAVE_IFPACK2_EXPLICIT_INSTANTIATION
1318 
1319 #endif // IFPACK2_BLOCKRELAXATION_DEF_HPP
double getApplyFlops() const
Returns the number of flops for the application of the preconditioner.
Definition: Ifpack2_BlockRelaxation_def.hpp:291
void applyMat(const Tpetra::MultiVector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > &X, Tpetra::MultiVector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > &Y, Teuchos::ETransp mode=Teuchos::NO_TRANS) const
Applies the matrix to a Tpetra::MultiVector.
Definition: Ifpack2_BlockRelaxation_def.hpp:385
Ifpack2::BlockRelaxation class declaration.
Ifpack2::TriDiContainer class declaration.
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
The communicator over which the input matrix is distributed.
Definition: Ifpack2_BlockRelaxation_def.hpp:180
MatrixType::global_ordinal_type global_ordinal_type
The type of global indices in the input MatrixType.
Definition: Ifpack2_BlockRelaxation_decl.hpp:102
std::string description() const
A one-line description of this object.
Definition: Ifpack2_BlockRelaxation_def.hpp:1146
double getApplyTime() const
Returns the time spent in apply().
Definition: Ifpack2_BlockRelaxation_def.hpp:275
int getNumApply() const
Returns the number of calls to apply().
Definition: Ifpack2_BlockRelaxation_def.hpp:251
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 Tpetra::Map< local_ordinal_type, global_ordinal_type, node_type > > getRangeMap() const
Returns the Tpetra::Map object associated with the range of this operator.
Definition: Ifpack2_BlockRelaxation_def.hpp:217
int getNumCompute() const
Returns the number of calls to compute().
Definition: Ifpack2_BlockRelaxation_def.hpp:243
void setParameters(const Teuchos::ParameterList &params)
Sets all the parameters for the preconditioner.
Definition: Ifpack2_BlockRelaxation_def.hpp:117
void compute()
compute the preconditioner for the specified matrix, diagonal perturbation thresholds and relaxation ...
Definition: Ifpack2_BlockRelaxation_def.hpp:548
Block relaxation preconditioners (or smoothers) for Tpetra::RowMatrix and Tpetra::CrsMatrix sparse ma...
Definition: Ifpack2_BlockRelaxation_decl.hpp:81
Ifpack2 implementation details.
double getInitializeTime() const
Returns the time spent in initialize().
Definition: Ifpack2_BlockRelaxation_def.hpp:259
Declaration of a user-defined partitioner in which the user can define a nonoverlapping partition of ...
double getComputeFlops() const
Returns the number of flops in the computation phase.
Definition: Ifpack2_BlockRelaxation_def.hpp:283
double getComputeTime() const
Returns the time spent in compute().
Definition: Ifpack2_BlockRelaxation_def.hpp:267
bool isComputed() const
Return true if compute() has been called.
Definition: Ifpack2_BlockRelaxation_decl.hpp:202
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:111
Teuchos::RCP< const Tpetra::Map< local_ordinal_type, global_ordinal_type, node_type > > getDomainMap() const
Returns the Tpetra::Map object associated with the domain of this operator.
Definition: Ifpack2_BlockRelaxation_def.hpp:203
Partition in which the user can define a nonoverlapping partition of the graph in any way they choose...
Definition: Ifpack2_Details_UserPartitioner_decl.hpp:66
bool isInitialized() const
Returns true if the preconditioner has been successfully initialized.
Definition: Ifpack2_BlockRelaxation_decl.hpp:192
int getNumInitialize() const
Returns the number of calls to initialize().
Definition: Ifpack2_BlockRelaxation_def.hpp:236
void initialize()
Initialize.
Definition: Ifpack2_BlockRelaxation_def.hpp:401
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:1192
Ifpack2::BandedContainer class declaration.
Ifpack2::DenseContainer class declaration.
MatrixType::node_type node_type
Node type of the input MatrixType.
Definition: Ifpack2_BlockRelaxation_decl.hpp:105
MatrixType::local_ordinal_type local_ordinal_type
The type of local indices in the input MatrixType.
Definition: Ifpack2_BlockRelaxation_decl.hpp:99
Definition: Ifpack2_Details_Amesos2Wrapper_decl.hpp:63
MatrixType::scalar_type scalar_type
The type of the entries of the input MatrixType.
Definition: Ifpack2_BlockRelaxation_decl.hpp:96
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:72
Declaration of ILUT preconditioner.
Teuchos::RCP< const Tpetra::RowMatrix< scalar_type, local_ordinal_type, global_ordinal_type, node_type > > getMatrix() const
The input matrix of this preconditioner&#39;s constructor.
Definition: Ifpack2_BlockRelaxation_def.hpp:194
void getValidParameters(Teuchos::ParameterList &params)
Fills a list which contains all the parameters possibly used by Ifpack2.
Definition: Ifpack2_Parameters.cpp:50
Ifpack2::SparseContainer class declaration.
Definition: Ifpack2_LinePartitioner_decl.hpp:75
A class to define linear partitions.
Definition: Ifpack2_LinearPartitioner_decl.hpp:54
void apply(const Tpetra::MultiVector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > &X, Tpetra::MultiVector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > &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:299
BlockRelaxation(const Teuchos::RCP< const row_matrix_type > &Matrix)
Constructor.
Definition: Ifpack2_BlockRelaxation_def.hpp:81