41#ifndef IFPACK2_RELAXATION_DEF_HPP
42#define IFPACK2_RELAXATION_DEF_HPP
44#include "Teuchos_StandardParameterEntryValidators.hpp"
45#include "Teuchos_TimeMonitor.hpp"
46#include "Tpetra_CrsMatrix.hpp"
47#include "Tpetra_BlockCrsMatrix.hpp"
48#include "Tpetra_BlockView.hpp"
50#include "Ifpack2_Details_getCrsMatrix.hpp"
51#include "MatrixMarket_Tpetra.hpp"
52#include "Tpetra_Details_residual.hpp"
56#include "KokkosSparse_gauss_seidel.hpp"
60 class NonnegativeIntValidator :
public Teuchos::ParameterEntryValidator {
63 NonnegativeIntValidator () {}
66 Teuchos::ParameterEntryValidator::ValidStringsList validStringValues ()
const {
72 validate (
const Teuchos::ParameterEntry& entry,
73 const std::string& paramName,
74 const std::string& sublistName)
const
77 Teuchos::any anyVal = entry.getAny (
true);
78 const std::string entryName = entry.getAny (
false).typeName ();
80 TEUCHOS_TEST_FOR_EXCEPTION(
81 anyVal.type () != typeid (
int),
82 Teuchos::Exceptions::InvalidParameterType,
83 "Parameter \"" << paramName <<
"\" in sublist \"" << sublistName
84 <<
"\" has the wrong type." << endl <<
"Parameter: " << paramName
85 << endl <<
"Type specified: " << entryName << endl
86 <<
"Type required: int" << endl);
88 const int val = Teuchos::any_cast<int> (anyVal);
89 TEUCHOS_TEST_FOR_EXCEPTION(
90 val < 0, Teuchos::Exceptions::InvalidParameterValue,
91 "Parameter \"" << paramName <<
"\" in sublist \"" << sublistName
92 <<
"\" is negative." << endl <<
"Parameter: " << paramName
93 << endl <<
"Value specified: " << val << endl
94 <<
"Required range: [0, INT_MAX]" << endl);
98 const std::string getXMLTypeName ()
const {
99 return "NonnegativeIntValidator";
104 printDoc (
const std::string& docString,
105 std::ostream &out)
const
107 Teuchos::StrUtils::printLines (out,
"# ", docString);
108 out <<
"#\tValidator Used: " << std::endl;
109 out <<
"#\t\tNonnegativeIntValidator" << std::endl;
116 template<class Scalar, const bool isOrdinal=Teuchos::ScalarTraits<Scalar>::isOrdinal>
120 static const Scalar eps ();
124 template<
class Scalar>
125 class SmallTraits<Scalar, true> {
127 static const Scalar eps () {
128 return Teuchos::ScalarTraits<Scalar>::one ();
133 template<
class Scalar>
134 class SmallTraits<Scalar, false> {
136 static const Scalar eps () {
137 return Teuchos::ScalarTraits<Scalar>::eps ();
142 template<
class Scalar,
143 const bool isComplex = Teuchos::ScalarTraits<Scalar>::isComplex>
144 struct RealTraits {};
146 template<
class Scalar>
147 struct RealTraits<Scalar, false> {
148 using val_type = Scalar;
149 using mag_type = Scalar;
150 static KOKKOS_INLINE_FUNCTION mag_type real (
const val_type& z) {
155 template<
class Scalar>
156 struct RealTraits<Scalar, true> {
157 using val_type = Scalar;
158 using mag_type =
typename Teuchos::ScalarTraits<Scalar>::magnitudeType;
159 static KOKKOS_INLINE_FUNCTION mag_type real (
const val_type& z) {
160 return Kokkos::ArithTraits<val_type>::real (z);
164 template<
class Scalar>
165 KOKKOS_INLINE_FUNCTION
typename RealTraits<Scalar>::mag_type
166 getRealValue (
const Scalar& z) {
167 return RealTraits<Scalar>::real (z);
174template<
class MatrixType>
176Relaxation<MatrixType>::
177updateCachedMultiVector (
const Teuchos::RCP<
const Tpetra::Map<local_ordinal_type, global_ordinal_type, node_type>>& map,
178 size_t numVecs)
const
182 if (cachedMV_.is_null () ||
183 map.get () != cachedMV_->getMap ().get () ||
184 cachedMV_->getNumVectors () != numVecs) {
185 using MV = Tpetra::MultiVector<scalar_type, local_ordinal_type,
186 global_ordinal_type, node_type>;
187 cachedMV_ = Teuchos::rcp (
new MV (map, numVecs,
false));
191template<
class MatrixType>
193setMatrix(
const Teuchos::RCP<const row_matrix_type>& A)
195 if (A.getRawPtr() != A_.getRawPtr()) {
196 Importer_ = Teuchos::null;
197 pointImporter_ = Teuchos::null;
198 Diagonal_ = Teuchos::null;
199 isInitialized_ =
false;
201 diagOffsets_ = Kokkos::View<size_t*, device_type>();
202 savedDiagOffsets_ =
false;
203 hasBlockCrsMatrix_ =
false;
204 serialGaussSeidel_ = Teuchos::null;
205 if (! A.is_null ()) {
206 IsParallel_ = (A->getRowMap ()->getComm ()->getSize () > 1);
212template<
class MatrixType>
214Relaxation (
const Teuchos::RCP<const row_matrix_type>& A)
216 IsParallel_ ((A.is_null () || A->getRowMap ().is_null () || A->getRowMap ()->getComm ().is_null ()) ?
218 A->getRowMap ()->getComm ()->getSize () > 1)
220 this->setObjectLabel (
"Ifpack2::Relaxation");
224template<
class MatrixType>
225Teuchos::RCP<const Teuchos::ParameterList>
228 using Teuchos::Array;
229 using Teuchos::ParameterList;
230 using Teuchos::parameterList;
233 using Teuchos::rcp_const_cast;
234 using Teuchos::rcp_implicit_cast;
235 using Teuchos::setStringToIntegralParameter;
236 typedef Teuchos::ParameterEntryValidator PEV;
238 if (validParams_.is_null ()) {
239 RCP<ParameterList> pl = parameterList (
"Ifpack2::Relaxation");
243 Array<std::string> precTypes (8);
244 precTypes[0] =
"Jacobi";
245 precTypes[1] =
"Gauss-Seidel";
246 precTypes[2] =
"Symmetric Gauss-Seidel";
247 precTypes[3] =
"MT Gauss-Seidel";
248 precTypes[4] =
"MT Symmetric Gauss-Seidel";
249 precTypes[5] =
"Richardson";
250 precTypes[6] =
"Two-stage Gauss-Seidel";
251 precTypes[7] =
"Two-stage Symmetric Gauss-Seidel";
252 Array<Details::RelaxationType> precTypeEnums (8);
253 precTypeEnums[0] = Details::JACOBI;
254 precTypeEnums[1] = Details::GS;
255 precTypeEnums[2] = Details::SGS;
256 precTypeEnums[3] = Details::MTGS;
257 precTypeEnums[4] = Details::MTSGS;
258 precTypeEnums[5] = Details::RICHARDSON;
259 precTypeEnums[6] = Details::GS2;
260 precTypeEnums[7] = Details::SGS2;
261 const std::string defaultPrecType (
"Jacobi");
262 setStringToIntegralParameter<Details::RelaxationType> (
"relaxation: type",
263 defaultPrecType,
"Relaxation method", precTypes (), precTypeEnums (),
266 const int numSweeps = 1;
267 RCP<PEV> numSweepsValidator =
268 rcp_implicit_cast<PEV> (rcp (
new NonnegativeIntValidator));
269 pl->set (
"relaxation: sweeps", numSweeps,
"Number of relaxation sweeps",
270 rcp_const_cast<const PEV> (numSweepsValidator));
273 const int numOuterSweeps = 1;
274 RCP<PEV> numOuterSweepsValidator =
275 rcp_implicit_cast<PEV> (rcp (
new NonnegativeIntValidator));
276 pl->set (
"relaxation: outer sweeps", numOuterSweeps,
"Number of outer local relaxation sweeps for two-stage GS",
277 rcp_const_cast<const PEV> (numOuterSweepsValidator));
279 const int numInnerSweeps = 1;
280 RCP<PEV> numInnerSweepsValidator =
281 rcp_implicit_cast<PEV> (rcp (
new NonnegativeIntValidator));
282 pl->set (
"relaxation: inner sweeps", numInnerSweeps,
"Number of inner local relaxation sweeps for two-stage GS",
283 rcp_const_cast<const PEV> (numInnerSweepsValidator));
285 const scalar_type innerDampingFactor = STS::one ();
286 pl->set (
"relaxation: inner damping factor", innerDampingFactor,
"Damping factor for the inner sweep of two-stage GS");
288 const bool innerSpTrsv =
false;
289 pl->set (
"relaxation: inner sparse-triangular solve", innerSpTrsv,
"Specify whether to use sptrsv instead of JR iterations for two-stage GS");
291 const bool compactForm =
false;
292 pl->set (
"relaxation: compact form", compactForm,
"Specify whether to use compact form of recurrence for two-stage GS");
295 pl->set (
"relaxation: damping factor", dampingFactor);
297 const bool zeroStartingSolution =
true;
298 pl->set (
"relaxation: zero starting solution", zeroStartingSolution);
300 const bool doBackwardGS =
false;
301 pl->set (
"relaxation: backward mode", doBackwardGS);
303 const bool doL1Method =
false;
304 pl->set (
"relaxation: use l1", doL1Method);
306 const magnitude_type l1eta = (STM::one() + STM::one() + STM::one()) /
307 (STM::one() + STM::one());
308 pl->set (
"relaxation: l1 eta", l1eta);
311 pl->set (
"relaxation: min diagonal value", minDiagonalValue);
313 const bool fixTinyDiagEntries =
false;
314 pl->set (
"relaxation: fix tiny diagonal entries", fixTinyDiagEntries);
316 const bool checkDiagEntries =
false;
317 pl->set (
"relaxation: check diagonal entries", checkDiagEntries);
319 Teuchos::ArrayRCP<local_ordinal_type> localSmoothingIndices = Teuchos::null;
320 pl->set(
"relaxation: local smoothing indices", localSmoothingIndices);
322 const bool is_matrix_structurally_symmetric =
false;
323 pl->set(
"relaxation: symmetric matrix structure", is_matrix_structurally_symmetric);
325 const bool ifpack2_dump_matrix =
false;
326 pl->set(
"relaxation: ifpack2 dump matrix", ifpack2_dump_matrix);
328 const int cluster_size = 1;
329 pl->set(
"relaxation: mtgs cluster size", cluster_size);
331 pl->set(
"relaxation: mtgs coloring algorithm",
"Default");
333 const int long_row_threshold = 0;
334 pl->set(
"relaxation: long row threshold", long_row_threshold);
336 validParams_ = rcp_const_cast<const ParameterList> (pl);
342template<
class MatrixType>
345 using Teuchos::getIntegralValue;
346 using Teuchos::ParameterList;
348 typedef scalar_type ST;
350 if (pl.isType<
double>(
"relaxation: damping factor")) {
353 ST df = pl.get<
double>(
"relaxation: damping factor");
354 pl.remove(
"relaxation: damping factor");
355 pl.set(
"relaxation: damping factor",df);
360 const Details::RelaxationType precType =
361 getIntegralValue<Details::RelaxationType> (pl,
"relaxation: type");
362 const int numSweeps = pl.get<
int> (
"relaxation: sweeps");
363 const ST dampingFactor = pl.get<ST> (
"relaxation: damping factor");
364 const bool zeroStartSol = pl.get<
bool> (
"relaxation: zero starting solution");
365 const bool doBackwardGS = pl.get<
bool> (
"relaxation: backward mode");
366 const bool doL1Method = pl.get<
bool> (
"relaxation: use l1");
367 const magnitude_type l1Eta = pl.get<magnitude_type> (
"relaxation: l1 eta");
368 const ST minDiagonalValue = pl.get<ST> (
"relaxation: min diagonal value");
369 const bool fixTinyDiagEntries = pl.get<
bool> (
"relaxation: fix tiny diagonal entries");
370 const bool checkDiagEntries = pl.get<
bool> (
"relaxation: check diagonal entries");
371 const bool is_matrix_structurally_symmetric = pl.get<
bool> (
"relaxation: symmetric matrix structure");
372 const bool ifpack2_dump_matrix = pl.get<
bool> (
"relaxation: ifpack2 dump matrix");
373 int cluster_size = 1;
374 if(pl.isParameter (
"relaxation: mtgs cluster size"))
375 cluster_size = pl.get<
int> (
"relaxation: mtgs cluster size");
376 int long_row_threshold = 0;
377 if(pl.isParameter (
"relaxation: long row threshold"))
378 long_row_threshold = pl.get<
int> (
"relaxation: long row threshold");
379 std::string color_algo_name = pl.get<std::string>(
"relaxation: mtgs coloring algorithm");
381 for(
char& c : color_algo_name)
383 if(color_algo_name ==
"default")
384 this->mtColoringAlgorithm_ = KokkosGraph::COLORING_DEFAULT;
385 else if(color_algo_name ==
"serial")
386 this->mtColoringAlgorithm_ = KokkosGraph::COLORING_SERIAL;
387 else if(color_algo_name ==
"vb")
388 this->mtColoringAlgorithm_ = KokkosGraph::COLORING_VB;
389 else if(color_algo_name ==
"vbbit")
390 this->mtColoringAlgorithm_ = KokkosGraph::COLORING_VBBIT;
391 else if(color_algo_name ==
"vbcs")
392 this->mtColoringAlgorithm_ = KokkosGraph::COLORING_VBCS;
393 else if(color_algo_name ==
"vbd")
394 this->mtColoringAlgorithm_ = KokkosGraph::COLORING_VBD;
395 else if(color_algo_name ==
"vbdbit")
396 this->mtColoringAlgorithm_ = KokkosGraph::COLORING_VBDBIT;
397 else if(color_algo_name ==
"eb")
398 this->mtColoringAlgorithm_ = KokkosGraph::COLORING_EB;
401 std::ostringstream msg;
402 msg <<
"Ifpack2::Relaxation: 'relaxation: mtgs coloring algorithm' = '" << color_algo_name <<
"' is not valid.\n";
403 msg <<
"Choices (not case sensitive) are: Default, Serial, VB, VBBIT, VBCS, VBD, VBDBIT, EB.";
404 TEUCHOS_TEST_FOR_EXCEPTION(
405 true, std::invalid_argument, msg.str());
408 Teuchos::ArrayRCP<local_ordinal_type> localSmoothingIndices = pl.get<Teuchos::ArrayRCP<local_ordinal_type> >(
"relaxation: local smoothing indices");
411 if (!std::is_same<double, ST>::value && pl.isType<
double>(
"relaxation: inner damping factor")) {
414 ST df = pl.get<
double>(
"relaxation: inner damping factor");
415 pl.remove(
"relaxation: inner damping factor");
416 pl.set(
"relaxation: inner damping factor",df);
419 if (long_row_threshold > 0) {
420 TEUCHOS_TEST_FOR_EXCEPTION(
421 cluster_size != 1, std::invalid_argument,
"Ifpack2::Relaxation: "
422 "Requested long row MTGS/MTSGS algorithm and cluster GS/SGS, but those are not compatible.");
423 TEUCHOS_TEST_FOR_EXCEPTION(
424 precType != Details::RelaxationType::MTGS && precType != Details::RelaxationType::MTSGS,
425 std::invalid_argument,
"Ifpack2::Relaxation: "
426 "Requested long row MTGS/MTSGS algorithm, but this is only compatible with preconditioner types "
427 "'MT Gauss-Seidel' and 'MT Symmetric Gauss-Seidel'.");
430 const ST innerDampingFactor = pl.get<ST> (
"relaxation: inner damping factor");
431 const int numInnerSweeps = pl.get<
int> (
"relaxation: inner sweeps");
432 const int numOuterSweeps = pl.get<
int> (
"relaxation: outer sweeps");
433 const bool innerSpTrsv = pl.get<
bool> (
"relaxation: inner sparse-triangular solve");
434 const bool compactForm = pl.get<
bool> (
"relaxation: compact form");
437 PrecType_ = precType;
438 NumSweeps_ = numSweeps;
439 DampingFactor_ = dampingFactor;
440 ZeroStartingSolution_ = zeroStartSol;
441 DoBackwardGS_ = doBackwardGS;
442 DoL1Method_ = doL1Method;
444 MinDiagonalValue_ = minDiagonalValue;
445 fixTinyDiagEntries_ = fixTinyDiagEntries;
446 checkDiagEntries_ = checkDiagEntries;
447 clusterSize_ = cluster_size;
448 longRowThreshold_ = long_row_threshold;
449 is_matrix_structurally_symmetric_ = is_matrix_structurally_symmetric;
450 ifpack2_dump_matrix_ = ifpack2_dump_matrix;
451 localSmoothingIndices_ = localSmoothingIndices;
453 NumInnerSweeps_ = numInnerSweeps;
454 NumOuterSweeps_ = numOuterSweeps;
455 InnerSpTrsv_ = innerSpTrsv;
456 InnerDampingFactor_ = innerDampingFactor;
457 CompactForm_ = compactForm;
461template<
class MatrixType>
466 this->setParametersImpl(
const_cast<Teuchos::ParameterList&
>(pl));
470template<
class MatrixType>
471Teuchos::RCP<const Teuchos::Comm<int> >
473 TEUCHOS_TEST_FOR_EXCEPTION(
474 A_.is_null (), std::runtime_error,
"Ifpack2::Relaxation::getComm: "
475 "The input matrix A is null. Please call setMatrix() with a nonnull "
476 "input matrix before calling this method.");
477 return A_->getRowMap ()->getComm ();
481template<
class MatrixType>
482Teuchos::RCP<const typename Relaxation<MatrixType>::row_matrix_type>
488template<
class MatrixType>
489Teuchos::RCP<
const Tpetra::Map<
typename MatrixType::local_ordinal_type,
490 typename MatrixType::global_ordinal_type,
491 typename MatrixType::node_type> >
493 TEUCHOS_TEST_FOR_EXCEPTION(
494 A_.is_null (), std::runtime_error,
"Ifpack2::Relaxation::getDomainMap: "
495 "The input matrix A is null. Please call setMatrix() with a nonnull "
496 "input matrix before calling this method.");
497 return A_->getDomainMap ();
501template<
class MatrixType>
502Teuchos::RCP<
const Tpetra::Map<
typename MatrixType::local_ordinal_type,
503 typename MatrixType::global_ordinal_type,
504 typename MatrixType::node_type> >
506 TEUCHOS_TEST_FOR_EXCEPTION(
507 A_.is_null (), std::runtime_error,
"Ifpack2::Relaxation::getRangeMap: "
508 "The input matrix A is null. Please call setMatrix() with a nonnull "
509 "input matrix before calling this method.");
510 return A_->getRangeMap ();
514template<
class MatrixType>
520template<
class MatrixType>
522 return(NumInitialize_);
526template<
class MatrixType>
532template<
class MatrixType>
538template<
class MatrixType>
540 return(InitializeTime_);
544template<
class MatrixType>
546 return(ComputeTime_);
550template<
class MatrixType>
556template<
class MatrixType>
558 return(ComputeFlops_);
562template<
class MatrixType>
569template<
class MatrixType>
571 TEUCHOS_TEST_FOR_EXCEPTION(
572 A_.is_null (), std::runtime_error,
"Ifpack2::Relaxation::getNodeSmootherComplexity: "
573 "The input matrix A is null. Please call setMatrix() with a nonnull "
574 "input matrix, then call compute(), before calling this method.");
576 return A_->getLocalNumRows() + A_->getLocalNumEntries();
580template<
class MatrixType>
583apply (
const Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& X,
584 Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& Y,
592 using Teuchos::rcpFromRef;
595 TEUCHOS_TEST_FOR_EXCEPTION(
596 A_.is_null (), std::runtime_error,
"Ifpack2::Relaxation::apply: "
597 "The input matrix A is null. Please call setMatrix() with a nonnull "
598 "input matrix, then call compute(), before calling this method.");
599 TEUCHOS_TEST_FOR_EXCEPTION(
602 "Ifpack2::Relaxation::apply: You must call compute() on this Ifpack2 "
603 "preconditioner instance before you may call apply(). You may call "
604 "isComputed() to find out if compute() has been called already.");
605 TEUCHOS_TEST_FOR_EXCEPTION(
606 X.getNumVectors() != Y.getNumVectors(),
608 "Ifpack2::Relaxation::apply: X and Y have different numbers of columns. "
609 "X has " << X.getNumVectors() <<
" columns, but Y has "
610 << Y.getNumVectors() <<
" columns.");
611 TEUCHOS_TEST_FOR_EXCEPTION(
612 beta != STS::zero (), std::logic_error,
613 "Ifpack2::Relaxation::apply: beta = " << beta <<
" != 0 case not "
616 const std::string timerName (
"Ifpack2::Relaxation::apply");
617 Teuchos::RCP<Teuchos::Time> timer = Teuchos::TimeMonitor::lookupCounter (timerName);
618 if (timer.is_null ()) {
619 timer = Teuchos::TimeMonitor::getNewCounter (timerName);
622 double startTime = timer->wallTime();
624 Teuchos::TimeMonitor timeMon (*timer);
626 if (alpha == STS::zero ()) {
628 Y.putScalar (STS::zero ());
636 Xcopy = rcp (
new MV (X, Teuchos::Copy));
638 Xcopy = rcpFromRef (X);
646 case Ifpack2::Details::JACOBI:
647 ApplyInverseJacobi(*Xcopy,Y);
649 case Ifpack2::Details::GS:
650 ApplyInverseSerialGS(*Xcopy, Y, DoBackwardGS_ ? Tpetra::Backward : Tpetra::Forward);
652 case Ifpack2::Details::SGS:
653 ApplyInverseSerialGS(*Xcopy, Y, Tpetra::Symmetric);
655 case Ifpack2::Details::MTGS:
656 case Ifpack2::Details::GS2:
657 ApplyInverseMTGS_CrsMatrix(*Xcopy, Y, DoBackwardGS_ ? Tpetra::Backward : Tpetra::Forward);
659 case Ifpack2::Details::MTSGS:
660 case Ifpack2::Details::SGS2:
661 ApplyInverseMTGS_CrsMatrix(*Xcopy, Y, Tpetra::Symmetric);
663 case Ifpack2::Details::RICHARDSON:
664 ApplyInverseRichardson(*Xcopy,Y);
668 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
669 "Ifpack2::Relaxation::apply: Invalid preconditioner type enum value "
670 << PrecType_ <<
". Please report this bug to the Ifpack2 developers.");
672 if (alpha != STS::one ()) {
674 const double numGlobalRows = as<double> (A_->getGlobalNumRows ());
675 const double numVectors = as<double> (Y.getNumVectors ());
676 ApplyFlops_ += numGlobalRows * numVectors;
680 ApplyTime_ += (timer->wallTime() - startTime);
685template<
class MatrixType>
688applyMat (
const Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& X,
689 Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& Y,
690 Teuchos::ETransp mode)
const
692 TEUCHOS_TEST_FOR_EXCEPTION(
693 A_.is_null (), std::runtime_error,
"Ifpack2::Relaxation::applyMat: "
694 "The input matrix A is null. Please call setMatrix() with a nonnull "
695 "input matrix, then call compute(), before calling this method.");
696 TEUCHOS_TEST_FOR_EXCEPTION(
697 ! isComputed (), std::runtime_error,
"Ifpack2::Relaxation::applyMat: "
698 "isComputed() must be true before you may call applyMat(). "
699 "Please call compute() before calling this method.");
700 TEUCHOS_TEST_FOR_EXCEPTION(
701 X.getNumVectors () != Y.getNumVectors (), std::invalid_argument,
702 "Ifpack2::Relaxation::applyMat: X.getNumVectors() = " << X.getNumVectors ()
703 <<
" != Y.getNumVectors() = " << Y.getNumVectors () <<
".");
704 A_->apply (X, Y, mode);
708template<
class MatrixType>
711 const char methodName[] =
"Ifpack2::Relaxation::initialize";
713 TEUCHOS_TEST_FOR_EXCEPTION
714 (A_.is_null (), std::runtime_error, methodName <<
": The "
715 "input matrix A is null. Please call setMatrix() with "
716 "a nonnull input matrix before calling this method.");
718 Teuchos::RCP<Teuchos::Time> timer =
719 Teuchos::TimeMonitor::getNewCounter (methodName);
721 double startTime = timer->wallTime();
724 Teuchos::TimeMonitor timeMon (*timer);
725 isInitialized_ =
false;
728 auto rowMap = A_->getRowMap ();
729 auto comm = rowMap.is_null () ? Teuchos::null : rowMap->getComm ();
730 IsParallel_ = ! comm.is_null () && comm->getSize () > 1;
740 Importer_ = IsParallel_ ? A_->getGraph ()->getImporter () :
744 Teuchos::RCP<const block_crs_matrix_type> A_bcrs =
745 Teuchos::rcp_dynamic_cast<const block_crs_matrix_type> (A_);
746 hasBlockCrsMatrix_ = ! A_bcrs.is_null ();
749 serialGaussSeidel_ = Teuchos::null;
751 if (PrecType_ == Details::MTGS || PrecType_ == Details::MTSGS ||
752 PrecType_ == Details::GS2 || PrecType_ == Details::SGS2) {
753 auto crsMat = Details::getCrsMatrix(A_);
754 TEUCHOS_TEST_FOR_EXCEPTION
755 (crsMat.is_null(), std::logic_error, methodName <<
": "
756 "Multithreaded Gauss-Seidel methods currently only work "
757 "when the input matrix is a Tpetra::CrsMatrix.");
762 if (ifpack2_dump_matrix_) {
763 static int sequence_number = 0;
764 const std::string file_name =
"Ifpack2_MT_GS_" +
765 std::to_string (sequence_number++) +
".mtx";
766 using writer_type = Tpetra::MatrixMarket::Writer<crs_matrix_type>;
767 writer_type::writeSparseFile (file_name, crsMat);
770 this->mtKernelHandle_ = Teuchos::rcp (
new mt_kernel_handle_type ());
771 if (mtKernelHandle_->get_gs_handle () ==
nullptr) {
772 if (PrecType_ == Details::GS2 || PrecType_ == Details::SGS2)
773 mtKernelHandle_->create_gs_handle (KokkosSparse::GS_TWOSTAGE);
774 else if(this->clusterSize_ == 1) {
775 mtKernelHandle_->create_gs_handle (KokkosSparse::GS_DEFAULT, this->mtColoringAlgorithm_);
776 mtKernelHandle_->get_point_gs_handle()->set_long_row_threshold(longRowThreshold_);
779 mtKernelHandle_->create_gs_handle (KokkosSparse::CLUSTER_DEFAULT, this->clusterSize_, this->mtColoringAlgorithm_);
781 local_matrix_device_type kcsr = crsMat->getLocalMatrixDevice ();
782 if (PrecType_ == Details::GS2 || PrecType_ == Details::SGS2) {
784 mtKernelHandle_->set_gs_set_num_inner_sweeps (NumInnerSweeps_);
785 mtKernelHandle_->set_gs_set_num_outer_sweeps (NumOuterSweeps_);
786 mtKernelHandle_->set_gs_set_inner_damp_factor (InnerDampingFactor_);
787 mtKernelHandle_->set_gs_twostage (!InnerSpTrsv_, A_->getLocalNumRows ());
788 mtKernelHandle_->set_gs_twostage_compact_form (CompactForm_);
791 KokkosSparse::Experimental::gauss_seidel_symbolic(
792 mtKernelHandle_.getRawPtr (),
793 A_->getLocalNumRows (),
794 A_->getLocalNumCols (),
797 is_matrix_structurally_symmetric_);
801 InitializeTime_ += (timer->wallTime() - startTime);
803 isInitialized_ =
true;
807template <
typename BlockDiagView>
808struct InvertDiagBlocks {
809 typedef typename BlockDiagView::size_type Size;
812 typedef Kokkos::MemoryTraits<Kokkos::Unmanaged> Unmanaged;
813 template <
typename View>
814 using UnmanagedView = Kokkos::View<
typename View::data_type,
typename View::array_layout,
815 typename View::device_type, Unmanaged>;
817 typedef typename BlockDiagView::non_const_value_type Scalar;
818 typedef typename BlockDiagView::device_type Device;
819 typedef Kokkos::View<Scalar**, Kokkos::LayoutRight, Device> RWrk;
820 typedef Kokkos::View<int**, Kokkos::LayoutRight, Device> IWrk;
822 UnmanagedView<BlockDiagView> block_diag_;
826 UnmanagedView<RWrk> rwrk_;
828 UnmanagedView<IWrk> iwrk_;
831 InvertDiagBlocks (BlockDiagView& block_diag)
832 : block_diag_(block_diag)
834 const auto blksz = block_diag.extent(1);
835 Kokkos::resize(rwrk_buf_, block_diag_.extent(0), blksz);
837 Kokkos::resize(iwrk_buf_, block_diag_.extent(0), blksz);
841 KOKKOS_INLINE_FUNCTION
842 void operator() (
const Size i,
int& jinfo)
const {
843 auto D_cur = Kokkos::subview(block_diag_, i, Kokkos::ALL(), Kokkos::ALL());
844 auto ipiv = Kokkos::subview(iwrk_, i, Kokkos::ALL());
845 auto work = Kokkos::subview(rwrk_, i, Kokkos::ALL());
847 Tpetra::GETF2(D_cur, ipiv, info);
852 Tpetra::GETRI(D_cur, ipiv, work, info);
858template<
class MatrixType>
859void Relaxation<MatrixType>::computeBlockCrs ()
862 using Teuchos::Array;
863 using Teuchos::ArrayRCP;
864 using Teuchos::ArrayView;
869 using Teuchos::REDUCE_MAX;
870 using Teuchos::REDUCE_MIN;
871 using Teuchos::REDUCE_SUM;
872 using Teuchos::rcp_dynamic_cast;
873 using Teuchos::reduceAll;
874 typedef local_ordinal_type LO;
876 const std::string timerName (
"Ifpack2::Relaxation::computeBlockCrs");
877 Teuchos::RCP<Teuchos::Time> timer = Teuchos::TimeMonitor::lookupCounter (timerName);
878 if (timer.is_null ()) {
879 timer = Teuchos::TimeMonitor::getNewCounter (timerName);
881 double startTime = timer->wallTime();
883 Teuchos::TimeMonitor timeMon (*timer);
885 TEUCHOS_TEST_FOR_EXCEPTION(
886 A_.is_null (), std::runtime_error,
"Ifpack2::Relaxation::"
887 "computeBlockCrs: The input matrix A is null. Please call setMatrix() "
888 "with a nonnull input matrix, then call initialize(), before calling "
890 auto blockCrsA = rcp_dynamic_cast<const block_crs_matrix_type> (A_);
891 TEUCHOS_TEST_FOR_EXCEPTION(
892 blockCrsA.is_null(), std::logic_error,
"Ifpack2::Relaxation::"
893 "computeBlockCrs: A_ is not a BlockCrsMatrix, but it should be if we "
894 "got this far. Please report this bug to the Ifpack2 developers.");
896 const scalar_type one = STS::one ();
901 const LO lclNumMeshRows =
902 blockCrsA->getCrsGraph ().getLocalNumRows ();
903 const LO blockSize = blockCrsA->getBlockSize ();
905 if (! savedDiagOffsets_) {
906 blockDiag_ = block_diag_type ();
907 blockDiag_ = block_diag_type (
"Ifpack2::Relaxation::blockDiag_",
908 lclNumMeshRows, blockSize, blockSize);
909 if (Teuchos::as<LO>(diagOffsets_.extent (0) ) < lclNumMeshRows) {
912 diagOffsets_ = Kokkos::View<size_t*, device_type> ();
913 diagOffsets_ = Kokkos::View<size_t*, device_type> (
"offsets", lclNumMeshRows);
915 blockCrsA->getCrsGraph ().getLocalDiagOffsets (diagOffsets_);
916 TEUCHOS_TEST_FOR_EXCEPTION
917 (
static_cast<size_t> (diagOffsets_.extent (0)) !=
918 static_cast<size_t> (blockDiag_.extent (0)),
919 std::logic_error,
"diagOffsets_.extent(0) = " <<
920 diagOffsets_.extent (0) <<
" != blockDiag_.extent(0) = "
921 << blockDiag_.extent (0) <<
922 ". Please report this bug to the Ifpack2 developers.");
923 savedDiagOffsets_ =
true;
925 blockCrsA->getLocalDiagCopy (blockDiag_, diagOffsets_);
932 unmanaged_block_diag_type blockDiag = blockDiag_;
934 if (DoL1Method_ && IsParallel_) {
935 const scalar_type two = one + one;
936 const size_t maxLength = A_->getLocalMaxNumRowEntries ();
937 nonconst_local_inds_host_view_type indices (
"indices",maxLength);
938 nonconst_values_host_view_type values_ (
"values",maxLength * blockSize * blockSize);
939 size_t numEntries = 0;
941 for (LO i = 0; i < lclNumMeshRows; ++i) {
943 blockCrsA->getLocalRowCopy (i, indices, values_, numEntries);
944 scalar_type * values =
reinterpret_cast<scalar_type*
>(values_.data());
946 auto diagBlock = Kokkos::subview (blockDiag, i, ALL (), ALL ());
947 for (LO subRow = 0; subRow < blockSize; ++subRow) {
948 magnitude_type diagonal_boost = STM::zero ();
949 for (
size_t k = 0 ; k < numEntries ; ++k) {
950 if (indices[k] >= lclNumMeshRows) {
951 const size_t offset = blockSize*blockSize*k + subRow*blockSize;
952 for (LO subCol = 0; subCol < blockSize; ++subCol) {
953 diagonal_boost += STS::magnitude (values[offset+subCol] / two);
957 if (STS::magnitude (diagBlock(subRow, subRow)) < L1Eta_ * diagonal_boost) {
958 diagBlock(subRow, subRow) += diagonal_boost;
966 Impl::InvertDiagBlocks<unmanaged_block_diag_type> idb(blockDiag);
967 typedef typename unmanaged_block_diag_type::execution_space exec_space;
968 typedef Kokkos::RangePolicy<exec_space, LO> range_type;
970 Kokkos::parallel_reduce (range_type (0, lclNumMeshRows), idb, info);
975 TEUCHOS_TEST_FOR_EXCEPTION
976 (info > 0, std::runtime_error,
"GETF2 or GETRI failed on = " << info
977 <<
" diagonal blocks.");
982#ifdef HAVE_IFPACK2_DEBUG
983 const int numResults = 2;
985 int lclResults[2], gblResults[2];
986 lclResults[0] = info;
987 lclResults[1] = -info;
990 reduceAll<int, int> (* (A_->getGraph ()->getComm ()), REDUCE_MIN,
991 numResults, lclResults, gblResults);
992 TEUCHOS_TEST_FOR_EXCEPTION
993 (gblResults[0] != 0 || gblResults[1] != 0, std::runtime_error,
994 "Ifpack2::Relaxation::compute: When processing the input "
995 "Tpetra::BlockCrsMatrix, one or more diagonal block LU factorizations "
996 "failed on one or more (MPI) processes.");
998 serialGaussSeidel_ = rcp(
new SerialGaussSeidel(blockCrsA, blockDiag_, localSmoothingIndices_, DampingFactor_));
1001 ComputeTime_ += (timer->wallTime() - startTime);
1006template<
class MatrixType>
1009 using Teuchos::Array;
1010 using Teuchos::ArrayRCP;
1011 using Teuchos::ArrayView;
1013 using Teuchos::Comm;
1016 using Teuchos::REDUCE_MAX;
1017 using Teuchos::REDUCE_MIN;
1018 using Teuchos::REDUCE_SUM;
1019 using Teuchos::reduceAll;
1023 using IST =
typename vector_type::impl_scalar_type;
1024 using KAT = Kokkos::ArithTraits<IST>;
1026 const char methodName[] =
"Ifpack2::Relaxation::compute";
1030 TEUCHOS_TEST_FOR_EXCEPTION
1031 (A_.is_null (), std::runtime_error, methodName <<
": "
1032 "The input matrix A is null. Please call setMatrix() with a nonnull "
1033 "input matrix, then call initialize(), before calling this method.");
1036 if (! isInitialized ()) {
1040 if (hasBlockCrsMatrix_) {
1045 Teuchos::RCP<Teuchos::Time> timer =
1046 Teuchos::TimeMonitor::getNewCounter (methodName);
1047 double startTime = timer->wallTime();
1050 Teuchos::TimeMonitor timeMon (*timer);
1059 IST oneOverMinDiagVal = KAT::one () /
static_cast<IST
> (SmallTraits<scalar_type>::eps ());
1060 if ( MinDiagonalValue_ != zero)
1061 oneOverMinDiagVal = KAT::one () /
static_cast<IST
> (MinDiagonalValue_);
1064 const magnitude_type minDiagValMag = STS::magnitude (MinDiagonalValue_);
1066 const LO numMyRows =
static_cast<LO
> (A_->getLocalNumRows ());
1068 TEUCHOS_TEST_FOR_EXCEPTION
1069 (NumSweeps_ < 0, std::logic_error, methodName
1070 <<
": NumSweeps_ = " << NumSweeps_ <<
" < 0. "
1071 "Please report this bug to the Ifpack2 developers.");
1072 IsComputed_ =
false;
1074 if (Diagonal_.is_null()) {
1077 Diagonal_ = rcp (
new vector_type (A_->getRowMap (),
false));
1080 if (checkDiagEntries_) {
1086 size_t numSmallDiagEntries = 0;
1087 size_t numZeroDiagEntries = 0;
1088 size_t numNegDiagEntries = 0;
1095 A_->getLocalDiagCopy (*Diagonal_);
1096 std::unique_ptr<vector_type> origDiag;
1097 origDiag = std::unique_ptr<vector_type> (
new vector_type (*Diagonal_, Teuchos::Copy));
1098 auto diag2d = Diagonal_->getLocalViewHost(Tpetra::Access::ReadWrite);
1099 auto diag = Kokkos::subview(diag2d, Kokkos::ALL(), 0);
1106 if (numMyRows != 0) {
1108 minMagDiagEntryMag = d_0_mag;
1109 maxMagDiagEntryMag = d_0_mag;
1118 for (LO i = 0; i < numMyRows; ++i) {
1119 const IST d_i = diag(i);
1123 const auto d_i_real = getRealValue (d_i);
1127 if (d_i_real < STM::zero ()) {
1128 ++numNegDiagEntries;
1130 if (d_i_mag < minMagDiagEntryMag) {
1131 minMagDiagEntryMag = d_i_mag;
1133 if (d_i_mag > maxMagDiagEntryMag) {
1134 maxMagDiagEntryMag = d_i_mag;
1137 if (fixTinyDiagEntries_) {
1139 if (d_i_mag <= minDiagValMag) {
1140 ++numSmallDiagEntries;
1141 if (d_i_mag == STM::zero ()) {
1142 ++numZeroDiagEntries;
1144 diag(i) = oneOverMinDiagVal;
1147 diag(i) = KAT::one () / d_i;
1152 if (d_i_mag <= minDiagValMag) {
1153 ++numSmallDiagEntries;
1154 if (d_i_mag == STM::zero ()) {
1155 ++numZeroDiagEntries;
1158 diag(i) = KAT::one () / d_i;
1176 const Comm<int>& comm = * (A_->getRowMap ()->getComm ());
1179 Array<magnitude_type> localVals (2);
1180 localVals[0] = minMagDiagEntryMag;
1182 localVals[1] = -maxMagDiagEntryMag;
1183 Array<magnitude_type> globalVals (2);
1184 reduceAll<int, magnitude_type> (comm, REDUCE_MIN, 2,
1185 localVals.getRawPtr (),
1186 globalVals.getRawPtr ());
1187 globalMinMagDiagEntryMag_ = globalVals[0];
1188 globalMaxMagDiagEntryMag_ = -globalVals[1];
1190 Array<size_t> localCounts (3);
1191 localCounts[0] = numSmallDiagEntries;
1192 localCounts[1] = numZeroDiagEntries;
1193 localCounts[2] = numNegDiagEntries;
1194 Array<size_t> globalCounts (3);
1195 reduceAll<int, size_t> (comm, REDUCE_SUM, 3,
1196 localCounts.getRawPtr (),
1197 globalCounts.getRawPtr ());
1198 globalNumSmallDiagEntries_ = globalCounts[0];
1199 globalNumZeroDiagEntries_ = globalCounts[1];
1200 globalNumNegDiagEntries_ = globalCounts[2];
1204 vector_type diff (A_->getRowMap ());
1205 diff.reciprocal (*origDiag);
1206 diff.update (-one, *Diagonal_, one);
1207 globalDiagNormDiff_ = diff.norm2 ();
1214 bool debugAgainstSlowPath =
false;
1216 auto crsMat = Details::getCrsMatrix(A_);
1218 if (crsMat.get() && crsMat->isFillComplete ()) {
1223 if (invDiagKernel_.is_null())
1226 invDiagKernel_->setMatrix(crsMat);
1227 invDiagKernel_->compute(*Diagonal_,
1228 DoL1Method_ && IsParallel_, L1Eta_,
1229 fixTinyDiagEntries_, minDiagValMag);
1230 savedDiagOffsets_ =
true;
1234#ifdef HAVE_IFPACK2_DEBUG
1235 debugAgainstSlowPath =
true;
1239 if (crsMat.is_null() || ! crsMat->isFillComplete () || debugAgainstSlowPath) {
1246 RCP<vector_type> Diagonal;
1247 if (debugAgainstSlowPath)
1248 Diagonal = rcp(
new vector_type(A_->getRowMap ()));
1250 Diagonal = Diagonal_;
1252 A_->getLocalDiagCopy (*Diagonal);
1262 if (DoL1Method_ && IsParallel_) {
1264 auto diag = Diagonal->getLocalViewHost(Tpetra::Access::ReadWrite);
1266 const size_t maxLength = A_row.getLocalMaxNumRowEntries ();
1267 nonconst_local_inds_host_view_type indices(
"indices",maxLength);
1268 nonconst_values_host_view_type values(
"values",maxLength);
1271 for (LO i = 0; i < numMyRows; ++i) {
1272 A_row.getLocalRowCopy (i, indices, values, numEntries);
1274 for (
size_t k = 0 ; k < numEntries; ++k) {
1275 if (indices[k] >= numMyRows) {
1276 diagonal_boost += STS::magnitude (values[k] / two);
1279 if (KAT::magnitude (diag(i, 0)) < L1Eta_ * diagonal_boost) {
1280 diag(i, 0) += diagonal_boost;
1288 if (fixTinyDiagEntries_) {
1292 auto localDiag = Diagonal->getLocalViewDevice(Tpetra::Access::ReadWrite);
1293 Kokkos::parallel_for(Kokkos::RangePolicy<MyExecSpace>(0, localDiag.extent(0)),
1295 auto d_i = localDiag(i, 0);
1298 if (d_i_mag <= minDiagValMag) {
1299 d_i = oneOverMinDiagVal;
1304 d_i = IST (KAT::one () / d_i);
1306 localDiag(i, 0) = d_i;
1310 Diagonal->reciprocal (*Diagonal);
1313 if (debugAgainstSlowPath) {
1315 Diagonal->update (STS::one (), *Diagonal_, -STS::one ());
1319 TEUCHOS_TEST_FOR_EXCEPTION
1320 (err > 100*STM::eps(), std::logic_error, methodName <<
": "
1321 <<
"\"fast-path\" diagonal computation failed. "
1322 "\\|D1 - D2\\|_inf = " << err <<
".");
1329 if (STS::isComplex) {
1330 ComputeFlops_ += 4.0 * numMyRows;
1333 ComputeFlops_ += numMyRows;
1336 if (PrecType_ == Ifpack2::Details::MTGS ||
1337 PrecType_ == Ifpack2::Details::MTSGS ||
1338 PrecType_ == Ifpack2::Details::GS2 ||
1339 PrecType_ == Ifpack2::Details::SGS2) {
1342 using scalar_view_t =
typename local_matrix_device_type::values_type;
1344 TEUCHOS_TEST_FOR_EXCEPTION
1345 (crsMat.is_null(), std::logic_error, methodName <<
": "
1346 "Multithreaded Gauss-Seidel methods currently only work "
1347 "when the input matrix is a Tpetra::CrsMatrix.");
1348 local_matrix_device_type kcsr = crsMat->getLocalMatrixDevice ();
1352 auto diagView_2d = Diagonal_->getLocalViewDevice (Tpetra::Access::ReadWrite);
1353 scalar_view_t diagView_1d = Kokkos::subview (diagView_2d, Kokkos::ALL (), 0);
1354 KokkosSparse::Experimental::gauss_seidel_numeric(
1355 mtKernelHandle_.getRawPtr (),
1356 A_->getLocalNumRows (),
1357 A_->getLocalNumCols (),
1362 is_matrix_structurally_symmetric_);
1364 else if(PrecType_ == Ifpack2::Details::GS || PrecType_ == Ifpack2::Details::SGS) {
1366 serialGaussSeidel_ = rcp(
new SerialGaussSeidel(crsMat, Diagonal_, localSmoothingIndices_, DampingFactor_));
1368 serialGaussSeidel_ = rcp(
new SerialGaussSeidel(A_, Diagonal_, localSmoothingIndices_, DampingFactor_));
1372 ComputeTime_ += (timer->wallTime() - startTime);
1378template<
class MatrixType>
1381ApplyInverseRichardson (
const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
1382 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y)
const
1385 const double numGlobalRows = as<double> (A_->getGlobalNumRows ());
1386 const double numVectors = as<double> (X.getNumVectors ());
1387 if (ZeroStartingSolution_) {
1391 Y.scale(DampingFactor_,X);
1397 double flopUpdate = 0.0;
1398 if (DampingFactor_ == STS::one ()) {
1400 flopUpdate = numGlobalRows * numVectors;
1404 flopUpdate = numGlobalRows * numVectors;
1406 ApplyFlops_ += flopUpdate;
1407 if (NumSweeps_ == 1) {
1413 const int startSweep = ZeroStartingSolution_ ? 1 : 0;
1416 updateCachedMultiVector(Y.getMap(),as<size_t>(numVectors));
1418 for (
int j = startSweep; j < NumSweeps_; ++j) {
1420 Tpetra::Details::residual(*A_,Y,X,*cachedMV_);
1421 Y.update(DampingFactor_,*cachedMV_,STS::one());
1435 const double numGlobalNonzeros = as<double> (A_->getGlobalNumEntries ());
1436 const double dampingFlops = (DampingFactor_ == STS::one ()) ? 0.0 : 1.0;
1437 ApplyFlops_ += as<double> (NumSweeps_ - startSweep) * numVectors *
1438 (2.0 * numGlobalNonzeros + dampingFlops);
1442template<
class MatrixType>
1444Relaxation<MatrixType>::
1445ApplyInverseJacobi (
const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
1446 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y)
const
1449 if (hasBlockCrsMatrix_) {
1450 ApplyInverseJacobi_BlockCrsMatrix (X, Y);
1454 const double numGlobalRows = as<double> (A_->getGlobalNumRows ());
1455 const double numVectors = as<double> (X.getNumVectors ());
1456 if (ZeroStartingSolution_) {
1461 Y.elementWiseMultiply (DampingFactor_, *Diagonal_, X, STS::zero ());
1467 double flopUpdate = 0.0;
1468 if (DampingFactor_ == STS::one ()) {
1470 flopUpdate = numGlobalRows * numVectors;
1474 flopUpdate = 2.0 * numGlobalRows * numVectors;
1476 ApplyFlops_ += flopUpdate;
1477 if (NumSweeps_ == 1) {
1483 const int startSweep = ZeroStartingSolution_ ? 1 : 0;
1486 updateCachedMultiVector(Y.getMap(),as<size_t>(numVectors));
1488 for (
int j = startSweep; j < NumSweeps_; ++j) {
1490 Tpetra::Details::residual(*A_,Y,X,*cachedMV_);
1491 Y.elementWiseMultiply (DampingFactor_, *Diagonal_, *cachedMV_, STS::one ());
1505 const double numGlobalNonzeros = as<double> (A_->getGlobalNumEntries ());
1506 const double dampingFlops = (DampingFactor_ == STS::one ()) ? 0.0 : 1.0;
1507 ApplyFlops_ += as<double> (NumSweeps_ - startSweep) * numVectors *
1508 (2.0 * numGlobalRows + 2.0 * numGlobalNonzeros + dampingFlops);
1511template<
class MatrixType>
1513Relaxation<MatrixType>::
1514ApplyInverseJacobi_BlockCrsMatrix (
const Tpetra::MultiVector<scalar_type,
1516 global_ordinal_type,
1518 Tpetra::MultiVector<scalar_type,
1520 global_ordinal_type,
1521 node_type>& Y)
const
1523 using Tpetra::BlockMultiVector;
1524 using BMV = BlockMultiVector<scalar_type, local_ordinal_type,
1525 global_ordinal_type, node_type>;
1527 const block_crs_matrix_type* blockMatConst =
1528 dynamic_cast<const block_crs_matrix_type*
> (A_.getRawPtr ());
1529 TEUCHOS_TEST_FOR_EXCEPTION
1530 (blockMatConst ==
nullptr, std::logic_error,
"This method should "
1531 "never be called if the matrix A_ is not a BlockCrsMatrix. "
1532 "Please report this bug to the Ifpack2 developers.");
1537 block_crs_matrix_type* blockMat =
1538 const_cast<block_crs_matrix_type*
> (blockMatConst);
1540 auto meshRowMap = blockMat->getRowMap ();
1541 auto meshColMap = blockMat->getColMap ();
1542 const local_ordinal_type blockSize = blockMat->getBlockSize ();
1544 BMV X_blk (X, *meshColMap, blockSize);
1545 BMV Y_blk (Y, *meshRowMap, blockSize);
1547 if (ZeroStartingSolution_) {
1551 Y_blk.blockWiseMultiply (DampingFactor_, blockDiag_, X_blk);
1552 if (NumSweeps_ == 1) {
1557 auto pointRowMap = Y.getMap ();
1558 const size_t numVecs = X.getNumVectors ();
1562 BMV A_times_Y (*meshRowMap, *pointRowMap, blockSize, numVecs);
1566 const int startSweep = ZeroStartingSolution_ ? 1 : 0;
1568 for (
int j = startSweep; j < NumSweeps_; ++j) {
1569 blockMat->applyBlock (Y_blk, A_times_Y);
1571 Y_blk.blockJacobiUpdate (DampingFactor_, blockDiag_,
1572 X_blk, A_times_Y, STS::one ());
1576template<
class MatrixType>
1578Relaxation<MatrixType>::
1579ApplyInverseSerialGS (
const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
1580 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
1581 Tpetra::ESweepDirection direction)
const
1583 using this_type = Relaxation<MatrixType>;
1593 auto blockCrsMat = Teuchos::rcp_dynamic_cast<const block_crs_matrix_type> (A_);
1594 auto crsMat = Details::getCrsMatrix(A_);
1595 if (blockCrsMat.get()) {
1596 const_cast<this_type&
> (*this).ApplyInverseSerialGS_BlockCrsMatrix (*blockCrsMat, X, Y, direction);
1598 else if (crsMat.get()) {
1599 ApplyInverseSerialGS_CrsMatrix (*crsMat, X, Y, direction);
1602 ApplyInverseSerialGS_RowMatrix (X, Y, direction);
1607template<
class MatrixType>
1609Relaxation<MatrixType>::
1610ApplyInverseSerialGS_RowMatrix (
const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
1611 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
1612 Tpetra::ESweepDirection direction)
const {
1613 using Teuchos::Array;
1614 using Teuchos::ArrayRCP;
1615 using Teuchos::ArrayView;
1619 using Teuchos::rcpFromRef;
1620 typedef Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> MV;
1625 if (ZeroStartingSolution_) {
1626 Y.putScalar (STS::zero ());
1629 size_t NumVectors = X.getNumVectors();
1633 if (Importer_.is_null ()) {
1634 updateCachedMultiVector (Y.getMap (), NumVectors);
1637 updateCachedMultiVector (Importer_->getTargetMap (), NumVectors);
1642 Y2 = rcpFromRef (Y);
1645 for (
int j = 0; j < NumSweeps_; ++j) {
1648 if (Importer_.is_null ()) {
1654 Y2->doImport (Y, *Importer_, Tpetra::INSERT);
1657 serialGaussSeidel_->apply(*Y2, X, direction);
1661 Tpetra::deep_copy (Y, *Y2);
1666 const double dampingFlops = (DampingFactor_ == STS::one()) ? 0.0 : 1.0;
1667 const double numVectors = as<double> (X.getNumVectors ());
1668 const double numGlobalRows = as<double> (A_->getGlobalNumRows ());
1669 const double numGlobalNonzeros = as<double> (A_->getGlobalNumEntries ());
1670 ApplyFlops_ += 2.0 * NumSweeps_ * numVectors *
1671 (2.0 * numGlobalRows + 2.0 * numGlobalNonzeros + dampingFlops);
1675template<
class MatrixType>
1677Relaxation<MatrixType>::
1678ApplyInverseSerialGS_CrsMatrix(
const crs_matrix_type& A,
1679 const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& B,
1680 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
1681 Tpetra::ESweepDirection direction)
const
1683 using Teuchos::null;
1686 using Teuchos::rcpFromRef;
1687 using Teuchos::rcp_const_cast;
1688 typedef scalar_type Scalar;
1689 const char prefix[] =
"Ifpack2::Relaxation::SerialGS: ";
1690 const scalar_type ZERO = Teuchos::ScalarTraits<Scalar>::zero ();
1692 TEUCHOS_TEST_FOR_EXCEPTION(
1693 ! A.isFillComplete (), std::runtime_error,
1694 prefix <<
"The matrix is not fill complete.");
1695 TEUCHOS_TEST_FOR_EXCEPTION(
1696 NumSweeps_ < 0, std::invalid_argument,
1697 prefix <<
"The number of sweeps must be nonnegative, "
1698 "but you provided numSweeps = " << NumSweeps_ <<
" < 0.");
1700 if (NumSweeps_ == 0) {
1704 RCP<const import_type> importer = A.getGraph ()->getImporter ();
1706 RCP<const map_type> domainMap = A.getDomainMap ();
1707 RCP<const map_type> rangeMap = A.getRangeMap ();
1708 RCP<const map_type> rowMap = A.getGraph ()->getRowMap ();
1709 RCP<const map_type> colMap = A.getGraph ()->getColMap ();
1711#ifdef HAVE_IFPACK2_DEBUG
1716 TEUCHOS_TEST_FOR_EXCEPTION(
1717 ! X.getMap ()->isSameAs (*domainMap), std::runtime_error,
1718 "Tpetra::CrsMatrix::gaussSeidelCopy requires that the input "
1719 "multivector X be in the domain Map of the matrix.");
1720 TEUCHOS_TEST_FOR_EXCEPTION(
1721 ! B.getMap ()->isSameAs (*rangeMap), std::runtime_error,
1722 "Tpetra::CrsMatrix::gaussSeidelCopy requires that the input "
1723 "B be in the range Map of the matrix.");
1724 TEUCHOS_TEST_FOR_EXCEPTION(
1725 ! Diagonal_->getMap ()->isSameAs (*rowMap), std::runtime_error,
1726 "Tpetra::CrsMatrix::gaussSeidelCopy requires that the input "
1727 "D be in the row Map of the matrix.");
1728 TEUCHOS_TEST_FOR_EXCEPTION(
1729 ! rowMap->isSameAs (*rangeMap), std::runtime_error,
1730 "Tpetra::CrsMatrix::gaussSeidelCopy requires that the row Map and the "
1731 "range Map be the same (in the sense of Tpetra::Map::isSameAs).");
1732 TEUCHOS_TEST_FOR_EXCEPTION(
1733 ! domainMap->isSameAs (*rangeMap), std::runtime_error,
1734 "Tpetra::CrsMatrix::gaussSeidelCopy requires that the domain Map and "
1735 "the range Map of the matrix be the same.");
1747 RCP<multivector_type> X_colMap;
1748 RCP<multivector_type> X_domainMap;
1749 bool copyBackOutput =
false;
1750 if (importer.is_null ()) {
1751 X_colMap = Teuchos::rcpFromRef (X);
1752 X_domainMap = Teuchos::rcpFromRef (X);
1758 if (ZeroStartingSolution_) {
1759 X_colMap->putScalar (ZERO);
1764 updateCachedMultiVector(colMap, X.getNumVectors());
1765 X_colMap = cachedMV_;
1766 X_domainMap = X_colMap->offsetViewNonConst (domainMap, 0);
1768 if (ZeroStartingSolution_) {
1770 X_colMap->putScalar (ZERO);
1780 X_colMap->doImport (X, *importer, Tpetra::INSERT);
1782 copyBackOutput =
true;
1785 for (
int sweep = 0; sweep < NumSweeps_; ++sweep) {
1786 if (! importer.is_null () && sweep > 0) {
1789 X_colMap->doImport (*X_domainMap, *importer, Tpetra::INSERT);
1792 serialGaussSeidel_->apply(*X_colMap, B, direction);
1795 if (copyBackOutput) {
1797 deep_copy (X , *X_domainMap);
1798 }
catch (std::exception& e) {
1799 TEUCHOS_TEST_FOR_EXCEPTION(
1800 true, std::runtime_error, prefix <<
"deep_copy(X, *X_domainMap) "
1801 "threw an exception: " << e.what ());
1816 const double dampingFlops = (DampingFactor_ == STS::one()) ? 0.0 : 1.0;
1817 const double numVectors = X.getNumVectors ();
1818 const double numGlobalRows = A_->getGlobalNumRows ();
1819 const double numGlobalNonzeros = A_->getGlobalNumEntries ();
1820 ApplyFlops_ += NumSweeps_ * numVectors *
1821 (2.0 * numGlobalRows + 2.0 * numGlobalNonzeros + dampingFlops);
1824template<
class MatrixType>
1826Relaxation<MatrixType>::
1827ApplyInverseSerialGS_BlockCrsMatrix (
const block_crs_matrix_type& A,
1828 const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
1829 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
1830 Tpetra::ESweepDirection direction)
1832 using Tpetra::INSERT;
1835 using Teuchos::rcpFromRef;
1844 block_multivector_type yBlock(Y, *(A.getGraph ()->getDomainMap()), A.getBlockSize());
1845 const block_multivector_type xBlock(X, *(A.getColMap ()), A.getBlockSize());
1847 bool performImport =
false;
1848 RCP<block_multivector_type> yBlockCol;
1849 if (Importer_.is_null()) {
1850 yBlockCol = rcpFromRef(yBlock);
1853 if (yBlockColumnPointMap_.is_null() ||
1854 yBlockColumnPointMap_->getNumVectors() != yBlock.getNumVectors() ||
1855 yBlockColumnPointMap_->getBlockSize() != yBlock.getBlockSize()) {
1856 yBlockColumnPointMap_ =
1857 rcp(
new block_multivector_type(*(A.getColMap()), A.getBlockSize(),
1858 static_cast<local_ordinal_type
>(yBlock.getNumVectors())));
1860 yBlockCol = yBlockColumnPointMap_;
1861 if (pointImporter_.is_null()) {
1862 auto srcMap = rcp(
new map_type(yBlock.getPointMap()));
1863 auto tgtMap = rcp(
new map_type(yBlockCol->getPointMap()));
1864 pointImporter_ = rcp(
new import_type(srcMap, tgtMap));
1866 performImport =
true;
1869 multivector_type yBlock_mv;
1870 multivector_type yBlockCol_mv;
1871 RCP<const multivector_type> yBlockColPointDomain;
1872 if (performImport) {
1873 yBlock_mv = yBlock.getMultiVectorView();
1874 yBlockCol_mv = yBlockCol->getMultiVectorView();
1875 yBlockColPointDomain =
1876 yBlockCol_mv.offsetView(A.getDomainMap(), 0);
1879 if (ZeroStartingSolution_) {
1880 yBlockCol->putScalar(STS::zero ());
1882 else if (performImport) {
1883 yBlockCol_mv.doImport(yBlock_mv, *pointImporter_, Tpetra::INSERT);
1886 for (
int sweep = 0; sweep < NumSweeps_; ++sweep) {
1887 if (performImport && sweep > 0) {
1888 yBlockCol_mv.doImport(yBlock_mv, *pointImporter_, Tpetra::INSERT);
1890 serialGaussSeidel_->applyBlock(*yBlockCol, xBlock, direction);
1891 if (performImport) {
1892 Tpetra::deep_copy(Y, *yBlockColPointDomain);
1897template<
class MatrixType>
1899Relaxation<MatrixType>::
1900ApplyInverseMTGS_CrsMatrix(
1901 const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& B,
1902 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
1903 Tpetra::ESweepDirection direction)
const
1905 using Teuchos::null;
1908 using Teuchos::rcpFromRef;
1909 using Teuchos::rcp_const_cast;
1912 typedef scalar_type Scalar;
1914 const char prefix[] =
"Ifpack2::Relaxation::(reordered)MTGaussSeidel: ";
1915 const Scalar ZERO = Teuchos::ScalarTraits<Scalar>::zero ();
1917 auto crsMat = Details::getCrsMatrix(A_);
1918 TEUCHOS_TEST_FOR_EXCEPTION
1919 (crsMat.is_null(), std::logic_error,
"Ifpack2::Relaxation::apply: "
1920 "Multithreaded Gauss-Seidel methods currently only work when the "
1921 "input matrix is a Tpetra::CrsMatrix.");
1924 TEUCHOS_TEST_FOR_EXCEPTION
1925 (! localSmoothingIndices_.is_null (), std::logic_error,
1926 "Ifpack2's implementation of Multithreaded Gauss-Seidel does not "
1927 "implement the case where the user supplies an iteration order. "
1928 "This error used to appear as \"MT GaussSeidel ignores the given "
1930 "I tried to add more explanation, but I didn't implement \"MT "
1931 "GaussSeidel\" [sic]. "
1932 "You'll have to ask the person who did.");
1934 TEUCHOS_TEST_FOR_EXCEPTION
1935 (! crsMat->isFillComplete (), std::runtime_error, prefix <<
"The "
1936 "input CrsMatrix is not fill complete. Please call fillComplete "
1937 "on the matrix, then do setup again, before calling apply(). "
1938 "\"Do setup\" means that if only the matrix's values have changed "
1939 "since last setup, you need only call compute(). If the matrix's "
1940 "structure may also have changed, you must first call initialize(), "
1941 "then call compute(). If you have not set up this preconditioner "
1942 "for this matrix before, you must first call initialize(), then "
1944 TEUCHOS_TEST_FOR_EXCEPTION
1945 (NumSweeps_ < 0, std::logic_error, prefix <<
": NumSweeps_ = "
1946 << NumSweeps_ <<
" < 0. Please report this bug to the Ifpack2 "
1948 if (NumSweeps_ == 0) {
1952 RCP<const import_type> importer = crsMat->getGraph ()->getImporter ();
1953 TEUCHOS_TEST_FOR_EXCEPTION(
1954 ! crsMat->getGraph ()->getExporter ().is_null (), std::runtime_error,
1955 "This method's implementation currently requires that the matrix's row, "
1956 "domain, and range Maps be the same. This cannot be the case, because "
1957 "the matrix has a nontrivial Export object.");
1959 RCP<const map_type> domainMap = crsMat->getDomainMap ();
1960 RCP<const map_type> rangeMap = crsMat->getRangeMap ();
1961 RCP<const map_type> rowMap = crsMat->getGraph ()->getRowMap ();
1962 RCP<const map_type> colMap = crsMat->getGraph ()->getColMap ();
1964#ifdef HAVE_IFPACK2_DEBUG
1969 TEUCHOS_TEST_FOR_EXCEPTION(
1970 ! X.getMap ()->isSameAs (*domainMap), std::runtime_error,
1971 "Ifpack2::Relaxation::MTGaussSeidel requires that the input "
1972 "multivector X be in the domain Map of the matrix.");
1973 TEUCHOS_TEST_FOR_EXCEPTION(
1974 ! B.getMap ()->isSameAs (*rangeMap), std::runtime_error,
1975 "Ifpack2::Relaxation::MTGaussSeidel requires that the input "
1976 "B be in the range Map of the matrix.");
1977 TEUCHOS_TEST_FOR_EXCEPTION(
1978 ! rowMap->isSameAs (*rangeMap), std::runtime_error,
1979 "Ifpack2::Relaxation::MTGaussSeidel requires that the row Map and the "
1980 "range Map be the same (in the sense of Tpetra::Map::isSameAs).");
1981 TEUCHOS_TEST_FOR_EXCEPTION(
1982 ! domainMap->isSameAs (*rangeMap), std::runtime_error,
1983 "Ifpack2::Relaxation::MTGaussSeidel requires that the domain Map and "
1984 "the range Map of the matrix be the same.");
2000 RCP<multivector_type> X_colMap;
2001 RCP<multivector_type> X_domainMap;
2002 bool copyBackOutput =
false;
2003 if (importer.is_null ()) {
2004 if (X.isConstantStride ()) {
2005 X_colMap = rcpFromRef (X);
2006 X_domainMap = rcpFromRef (X);
2013 if (ZeroStartingSolution_) {
2014 X_colMap->putScalar (ZERO);
2023 updateCachedMultiVector(colMap,as<size_t>(X.getNumVectors()));
2024 X_colMap = cachedMV_;
2028 X_domainMap = X_colMap;
2029 if (! ZeroStartingSolution_) {
2031 deep_copy (*X_domainMap , X);
2032 }
catch (std::exception& e) {
2033 std::ostringstream os;
2034 os <<
"Ifpack2::Relaxation::MTGaussSeidel: "
2035 "deep_copy(*X_domainMap, X) threw an exception: "
2036 << e.what () <<
".";
2037 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::runtime_error, e.what ());
2040 copyBackOutput =
true;
2054 updateCachedMultiVector(colMap,as<size_t>(X.getNumVectors()));
2055 X_colMap = cachedMV_;
2057 X_domainMap = X_colMap->offsetViewNonConst (domainMap, 0);
2059 if (ZeroStartingSolution_) {
2061 X_colMap->putScalar (ZERO);
2071 X_colMap->doImport (X, *importer, Tpetra::CombineMode::INSERT);
2073 copyBackOutput =
true;
2079 RCP<const multivector_type> B_in;
2080 if (B.isConstantStride ()) {
2081 B_in = rcpFromRef (B);
2087 RCP<multivector_type> B_in_nonconst = rcp (
new multivector_type(rowMap, B.getNumVectors()));
2089 deep_copy (*B_in_nonconst, B);
2090 }
catch (std::exception& e) {
2091 std::ostringstream os;
2092 os <<
"Ifpack2::Relaxation::MTGaussSeidel: "
2093 "deep_copy(*B_in_nonconst, B) threw an exception: "
2094 << e.what () <<
".";
2095 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::runtime_error, e.what ());
2097 B_in = rcp_const_cast<const multivector_type> (B_in_nonconst);
2111 local_matrix_device_type kcsr = crsMat->getLocalMatrixDevice ();
2113 bool update_y_vector =
true;
2115 bool zero_x_vector =
false;
2117 for (
int sweep = 0; sweep < NumSweeps_; ++sweep) {
2118 if (! importer.is_null () && sweep > 0) {
2121 X_colMap->doImport (*X_domainMap, *importer, Tpetra::CombineMode::INSERT);
2124 if (direction == Tpetra::Symmetric) {
2125 KokkosSparse::Experimental::symmetric_gauss_seidel_apply
2126 (mtKernelHandle_.getRawPtr(), A_->getLocalNumRows(), A_->getLocalNumCols(),
2127 kcsr.graph.row_map, kcsr.graph.entries, kcsr.values,
2128 X_colMap->getLocalViewDevice(Tpetra::Access::ReadWrite),
2129 B_in->getLocalViewDevice(Tpetra::Access::ReadOnly),
2130 zero_x_vector, update_y_vector, DampingFactor_, 1);
2132 else if (direction == Tpetra::Forward) {
2133 KokkosSparse::Experimental::forward_sweep_gauss_seidel_apply
2134 (mtKernelHandle_.getRawPtr(), A_->getLocalNumRows(), A_->getLocalNumCols(),
2135 kcsr.graph.row_map,kcsr.graph.entries, kcsr.values,
2136 X_colMap->getLocalViewDevice(Tpetra::Access::ReadWrite),
2137 B_in->getLocalViewDevice(Tpetra::Access::ReadOnly),
2138 zero_x_vector, update_y_vector, DampingFactor_, 1);
2140 else if (direction == Tpetra::Backward) {
2141 KokkosSparse::Experimental::backward_sweep_gauss_seidel_apply
2142 (mtKernelHandle_.getRawPtr(), A_->getLocalNumRows(), A_->getLocalNumCols(),
2143 kcsr.graph.row_map,kcsr.graph.entries, kcsr.values,
2144 X_colMap->getLocalViewDevice(Tpetra::Access::ReadWrite),
2145 B_in->getLocalViewDevice(Tpetra::Access::ReadOnly),
2146 zero_x_vector, update_y_vector, DampingFactor_, 1);
2149 TEUCHOS_TEST_FOR_EXCEPTION(
2150 true, std::invalid_argument,
2151 prefix <<
"The 'direction' enum does not have any of its valid "
2152 "values: Forward, Backward, or Symmetric.");
2154 update_y_vector =
false;
2157 if (copyBackOutput) {
2159 deep_copy (X , *X_domainMap);
2160 }
catch (std::exception& e) {
2161 TEUCHOS_TEST_FOR_EXCEPTION(
2162 true, std::runtime_error, prefix <<
"deep_copy(X, *X_domainMap) "
2163 "threw an exception: " << e.what ());
2167 const double dampingFlops = (DampingFactor_ == STS::one ()) ? 0.0 : 1.0;
2168 const double numVectors = as<double> (X.getNumVectors ());
2169 const double numGlobalRows = as<double> (A_->getGlobalNumRows ());
2170 const double numGlobalNonzeros = as<double> (A_->getGlobalNumEntries ());
2171 double ApplyFlops = NumSweeps_ * numVectors *
2172 (2.0 * numGlobalRows + 2.0 * numGlobalNonzeros + dampingFlops);
2173 if (direction == Tpetra::Symmetric)
2175 ApplyFlops_ += ApplyFlops;
2178template<
class MatrixType>
2181 std::ostringstream os;
2186 os <<
"\"Ifpack2::Relaxation\": {";
2188 os <<
"Initialized: " << (isInitialized () ?
"true" :
"false") <<
", "
2189 <<
"Computed: " << (isComputed () ?
"true" :
"false") <<
", ";
2195 if (PrecType_ == Ifpack2::Details::JACOBI) {
2197 }
else if (PrecType_ == Ifpack2::Details::GS) {
2198 os <<
"Gauss-Seidel";
2199 }
else if (PrecType_ == Ifpack2::Details::SGS) {
2200 os <<
"Symmetric Gauss-Seidel";
2201 }
else if (PrecType_ == Ifpack2::Details::MTGS) {
2202 os <<
"MT Gauss-Seidel";
2203 }
else if (PrecType_ == Ifpack2::Details::MTSGS) {
2204 os <<
"MT Symmetric Gauss-Seidel";
2205 }
else if (PrecType_ == Ifpack2::Details::GS2) {
2206 os <<
"Two-stage Gauss-Seidel";
2207 }
else if (PrecType_ == Ifpack2::Details::SGS2) {
2208 os <<
"Two-stage Symmetric Gauss-Seidel";
2213 if(hasBlockCrsMatrix_)
2216 os <<
", " <<
"sweeps: " << NumSweeps_ <<
", "
2217 <<
"damping factor: " << DampingFactor_ <<
", ";
2219 if (PrecType_ == Ifpack2::Details::MTGS || PrecType_ == Ifpack2::Details::MTSGS) {
2220 os <<
"\"relaxation: mtgs cluster size\": " << clusterSize_ <<
", ";
2221 os <<
"\"relaxation: long row threshold\": " << longRowThreshold_ <<
", ";
2222 os <<
"\"relaxation: symmetric matrix structure\": " << (is_matrix_structurally_symmetric_ ?
"true" :
"false") <<
", ";
2223 os <<
"\"relaxation: relaxation: mtgs coloring algorithm\": ";
2224 switch(mtColoringAlgorithm_)
2226 case KokkosGraph::COLORING_DEFAULT:
2227 os <<
"DEFAULT";
break;
2228 case KokkosGraph::COLORING_SERIAL:
2229 os <<
"SERIAL";
break;
2230 case KokkosGraph::COLORING_VB:
2232 case KokkosGraph::COLORING_VBBIT:
2233 os <<
"VBBIT";
break;
2234 case KokkosGraph::COLORING_VBCS:
2235 os <<
"VBCS";
break;
2236 case KokkosGraph::COLORING_VBD:
2238 case KokkosGraph::COLORING_VBDBIT:
2239 os <<
"VBDBIT";
break;
2240 case KokkosGraph::COLORING_EB:
2248 if (PrecType_ == Ifpack2::Details::GS2 ||
2249 PrecType_ == Ifpack2::Details::SGS2) {
2250 os <<
"outer sweeps: " << NumOuterSweeps_ <<
", "
2251 <<
"inner sweeps: " << NumInnerSweeps_ <<
", "
2252 <<
"inner damping factor: " << InnerDampingFactor_ <<
", ";
2256 os <<
"use l1: " << DoL1Method_ <<
", "
2257 <<
"l1 eta: " << L1Eta_ <<
", ";
2260 if (A_.is_null ()) {
2261 os <<
"Matrix: null";
2264 os <<
"Global matrix dimensions: ["
2265 << A_->getGlobalNumRows () <<
", " << A_->getGlobalNumCols () <<
"]"
2266 <<
", Global nnz: " << A_->getGlobalNumEntries();
2274template<
class MatrixType>
2277describe (Teuchos::FancyOStream &out,
2278 const Teuchos::EVerbosityLevel verbLevel)
const
2280 using Teuchos::OSTab;
2281 using Teuchos::TypeNameTraits;
2282 using Teuchos::VERB_DEFAULT;
2283 using Teuchos::VERB_NONE;
2284 using Teuchos::VERB_LOW;
2285 using Teuchos::VERB_MEDIUM;
2286 using Teuchos::VERB_HIGH;
2287 using Teuchos::VERB_EXTREME;
2290 const Teuchos::EVerbosityLevel vl =
2291 (verbLevel == VERB_DEFAULT) ? VERB_LOW : verbLevel;
2293 const int myRank = this->getComm ()->getRank ();
2301 if (vl != VERB_NONE && myRank == 0) {
2305 out <<
"\"Ifpack2::Relaxation\":" << endl;
2307 out <<
"MatrixType: \"" << TypeNameTraits<MatrixType>::name () <<
"\""
2309 if (this->getObjectLabel () !=
"") {
2310 out <<
"Label: " << this->getObjectLabel () << endl;
2312 out <<
"Initialized: " << (isInitialized () ?
"true" :
"false") << endl
2313 <<
"Computed: " << (isComputed () ?
"true" :
"false") << endl
2314 <<
"Parameters: " << endl;
2317 out <<
"\"relaxation: type\": ";
2318 if (PrecType_ == Ifpack2::Details::JACOBI) {
2320 }
else if (PrecType_ == Ifpack2::Details::GS) {
2321 out <<
"Gauss-Seidel";
2322 }
else if (PrecType_ == Ifpack2::Details::SGS) {
2323 out <<
"Symmetric Gauss-Seidel";
2324 }
else if (PrecType_ == Ifpack2::Details::MTGS) {
2325 out <<
"MT Gauss-Seidel";
2326 }
else if (PrecType_ == Ifpack2::Details::MTSGS) {
2327 out <<
"MT Symmetric Gauss-Seidel";
2328 }
else if (PrecType_ == Ifpack2::Details::GS2) {
2329 out <<
"Two-stage Gauss-Seidel";
2330 }
else if (PrecType_ == Ifpack2::Details::SGS2) {
2331 out <<
"Two-stage Symmetric Gauss-Seidel";
2338 <<
"\"relaxation: sweeps\": " << NumSweeps_ << endl
2339 <<
"\"relaxation: damping factor\": " << DampingFactor_ << endl
2340 <<
"\"relaxation: min diagonal value\": " << MinDiagonalValue_ << endl
2341 <<
"\"relaxation: zero starting solution\": " << ZeroStartingSolution_ << endl
2342 <<
"\"relaxation: backward mode\": " << DoBackwardGS_ << endl
2343 <<
"\"relaxation: use l1\": " << DoL1Method_ << endl
2344 <<
"\"relaxation: l1 eta\": " << L1Eta_ << endl;
2345 if (PrecType_ == Ifpack2::Details::MTGS || PrecType_ == Ifpack2::Details::MTSGS) {
2346 out <<
"\"relaxation: mtgs cluster size\": " << clusterSize_ << endl;
2347 out <<
"\"relaxation: long row threshold\": " << longRowThreshold_ << endl;
2348 out <<
"\"relaxation: symmetric matrix structure\": " << (is_matrix_structurally_symmetric_ ?
"true" :
"false") << endl;
2349 out <<
"\"relaxation: relaxation: mtgs coloring algorithm\": ";
2350 switch(mtColoringAlgorithm_)
2352 case KokkosGraph::COLORING_DEFAULT:
2353 out <<
"DEFAULT";
break;
2354 case KokkosGraph::COLORING_SERIAL:
2355 out <<
"SERIAL";
break;
2356 case KokkosGraph::COLORING_VB:
2358 case KokkosGraph::COLORING_VBBIT:
2359 out <<
"VBBIT";
break;
2360 case KokkosGraph::COLORING_VBCS:
2361 out <<
"VBCS";
break;
2362 case KokkosGraph::COLORING_VBD:
2363 out <<
"VBD";
break;
2364 case KokkosGraph::COLORING_VBDBIT:
2365 out <<
"VBDBIT";
break;
2366 case KokkosGraph::COLORING_EB:
2373 if (PrecType_ == Ifpack2::Details::GS2 || PrecType_ == Ifpack2::Details::SGS2) {
2374 out <<
"\"relaxation: inner damping factor\": " << InnerDampingFactor_ << endl;
2375 out <<
"\"relaxation: outer sweeps\" : " << NumOuterSweeps_ << endl;
2376 out <<
"\"relaxation: inner sweeps\" : " << NumInnerSweeps_ << endl;
2379 out <<
"Computed quantities:" << endl;
2382 out <<
"Global number of rows: " << A_->getGlobalNumRows () << endl
2383 <<
"Global number of columns: " << A_->getGlobalNumCols () << endl;
2385 if (checkDiagEntries_ && isComputed ()) {
2386 out <<
"Properties of input diagonal entries:" << endl;
2389 out <<
"Magnitude of minimum-magnitude diagonal entry: "
2390 << globalMinMagDiagEntryMag_ << endl
2391 <<
"Magnitude of maximum-magnitude diagonal entry: "
2392 << globalMaxMagDiagEntryMag_ << endl
2393 <<
"Number of diagonal entries with small magnitude: "
2394 << globalNumSmallDiagEntries_ << endl
2395 <<
"Number of zero diagonal entries: "
2396 << globalNumZeroDiagEntries_ << endl
2397 <<
"Number of diagonal entries with negative real part: "
2398 << globalNumNegDiagEntries_ << endl
2399 <<
"Abs 2-norm diff between computed and actual inverse "
2400 <<
"diagonal: " << globalDiagNormDiff_ << endl;
2403 if (isComputed ()) {
2404 out <<
"Saved diagonal offsets: "
2405 << (savedDiagOffsets_ ?
"true" :
"false") << endl;
2407 out <<
"Call counts and total times (in seconds): " << endl;
2410 out <<
"initialize: " << endl;
2413 out <<
"Call count: " << NumInitialize_ << endl;
2414 out <<
"Total time: " << InitializeTime_ << endl;
2416 out <<
"compute: " << endl;
2419 out <<
"Call count: " << NumCompute_ << endl;
2420 out <<
"Total time: " << ComputeTime_ << endl;
2422 out <<
"apply: " << endl;
2425 out <<
"Call count: " << NumApply_ << endl;
2426 out <<
"Total time: " << ApplyTime_ << endl;
2435#define IFPACK2_RELAXATION_INSTANT(S,LO,GO,N) \
2436 template class Ifpack2::Relaxation< Tpetra::RowMatrix<S, LO, GO, N> >;
File for utility functions.
Compute scaled damped residual for Chebyshev.
Definition: Ifpack2_Details_InverseDiagonalKernel_decl.hpp:77
Relaxation preconditioners for Tpetra::RowMatrix and Tpetra::CrsMatrix sparse matrices.
Definition: Ifpack2_Relaxation_decl.hpp:248
Relaxation(const Teuchos::RCP< const row_matrix_type > &A)
Constructor.
Definition: Ifpack2_Relaxation_def.hpp:214
MatrixType::global_ordinal_type global_ordinal_type
The type of global indices in the input MatrixType.
Definition: Ifpack2_Relaxation_decl.hpp:260
void setParameters(const Teuchos::ParameterList ¶ms)
Set the relaxation / preconditioner parameters.
Definition: Ifpack2_Relaxation_def.hpp:462
MatrixType::node_type node_type
The Node type used by the input MatrixType.
Definition: Ifpack2_Relaxation_decl.hpp:263
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
The communicator over which the matrix and vectors are distributed.
Definition: Ifpack2_Relaxation_def.hpp:472
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_Relaxation_def.hpp:492
double getApplyTime() const
Total time in seconds spent in all calls to apply().
Definition: Ifpack2_Relaxation_def.hpp:551
int getNumCompute() const
Total number of calls to compute().
Definition: Ifpack2_Relaxation_def.hpp:527
Teuchos::RCP< const row_matrix_type > getMatrix() const
The matrix to be preconditioned.
Definition: Ifpack2_Relaxation_def.hpp:483
void compute()
Compute the preconditioner ("numeric setup");.
Definition: Ifpack2_Relaxation_def.hpp:1007
void initialize()
Initialize the preconditioner ("symbolic setup").
Definition: Ifpack2_Relaxation_def.hpp:709
double getComputeTime() const
Total time in seconds spent in all calls to compute().
Definition: Ifpack2_Relaxation_def.hpp:545
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const
Print the object's attributes to the given output stream.
Definition: Ifpack2_Relaxation_def.hpp:2277
std::string description() const
A simple one-line description of this object.
Definition: Ifpack2_Relaxation_def.hpp:2179
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Return a list of all the parameters that this class accepts.
Definition: Ifpack2_Relaxation_def.hpp:226
double getInitializeTime() const
Total time in seconds spent in all calls to initialize().
Definition: Ifpack2_Relaxation_def.hpp:539
virtual void setMatrix(const Teuchos::RCP< const row_matrix_type > &A)
Change the matrix to be preconditioned.
Definition: Ifpack2_Relaxation_def.hpp:193
MatrixType::scalar_type scalar_type
The type of the entries of the input MatrixType.
Definition: Ifpack2_Relaxation_decl.hpp:254
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
Apply the preconditioner to X, returning the result in Y.
Definition: Ifpack2_Relaxation_def.hpp:583
MatrixType::local_ordinal_type local_ordinal_type
The type of local indices in the input MatrixType.
Definition: Ifpack2_Relaxation_decl.hpp:257
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
Apply the input matrix to X, returning the result in Y.
Definition: Ifpack2_Relaxation_def.hpp:688
Tpetra::RowMatrix< scalar_type, local_ordinal_type, global_ordinal_type, node_type > row_matrix_type
Tpetra::RowMatrix specialization used by this class.
Definition: Ifpack2_Relaxation_decl.hpp:273
Teuchos::ScalarTraits< scalar_type >::magnitudeType magnitude_type
The type of the magnitude (absolute value) of a matrix entry.
Definition: Ifpack2_Relaxation_decl.hpp:269
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
Definition: Ifpack2_Relaxation_def.hpp:570
int getNumApply() const
Total number of calls to apply().
Definition: Ifpack2_Relaxation_def.hpp:533
int getNumInitialize() const
Total number of calls to initialize().
Definition: Ifpack2_Relaxation_def.hpp:521
double getComputeFlops() const
Total number of floating-point operations over all calls to compute().
Definition: Ifpack2_Relaxation_def.hpp:557
bool hasTransposeApply() const
Whether apply() and applyMat() let you apply the transpose or conjugate transpose.
Definition: Ifpack2_Relaxation_def.hpp:515
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_Relaxation_def.hpp:505
double getApplyFlops() const
Total number of floating-point operations over all calls to apply().
Definition: Ifpack2_Relaxation_def.hpp:563
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:74
void getValidParameters(Teuchos::ParameterList ¶ms)
Fills a list which contains all the parameters possibly used by Ifpack2.
Definition: Ifpack2_Parameters.cpp:51