Merge Pull Request #6448 from trilinos/Trilinos/master_merge_20191213…

…_000614

Automatically Merged using Trilinos Master Merge AutoTester
PR Title: Trilinos Master Merge PR Generator: Auto PR created to promote from master_merge_20191213_000614 branch to master
PR Author: trilinos-autotester

packages/ifpack2/doc/UsersGuide/ifpack2guide.bib

@@ -15,6 +15,12 @@ @book{Saad2003
   publisher = {SIAM},
   edition   = {2nd},
 }
+@inproceedings{Saad1999,
+  title     = {Enhanced Parallel Multicolor Preconditioning Techniques for Linear Systems},
+  author    = {Yousef Saad and Masha Sosonkina},
+  booktitle = {PPSC},
+  year      = {1999}
+}
 @article{baker2011multigrid,
   title     = {Multigrid smoothers for ultraparallel computing},
   author    = {Baker, Allison H and Falgout, Robert D and Kolev, Tzanio V and Yang, Ulrike Meier},

packages/ifpack2/doc/UsersGuide/options.tex

@@ -97,9 +97,22 @@ \section{Point relaxation}\label{s:relaxation}
      and communication) to count diagonal entries that are zero, have negative
      real part, or are small in magnitude. This information can be later shown
      in the description.}
+\ccc{relaxation: mtgs cluster size}
+    {int}
+    {1}
+    {Only has an effect if {\tt "relaxation: type"} is {\tt "MT Gauss-Seidel"}
+     or {\tt "MT Symmetric Gauss-Seidel"}. If equal to 1 (default), point
+     coloring parallel Gauss-Seidel is used. This has a faster \texttt{compute()}
+     but may cause the preconditioned solver
+     to converge more slowly. If set to $k > 1$, then multicolor block Gauss-Seidel
+     is used with blocks of size $k$ (see \cite{Saad1999}).
+     In the \texttt{apply()} there is significantly less
+     error due to parallel updates of the LHS vector.}
 \ccc{relaxation: local smoothing indices}
     {Teuchos::ArrayRCP<local\_ordinal>}
     {empty}
+    {}
+
 %Teuchos::ArrayRCP MatrixType::local_ordinal_type}{\texttt{Teuchos::null}}
     {A given method will only relax on the local indices listed in the
      \texttt{ArrayRCP}, in the order that they are listed. This can be used to

packages/ifpack2/src/Ifpack2_Parameters.cpp

@@ -64,97 +64,98 @@ void getValidParameters(Teuchos::ParameterList& params)
   // ============================================================ //
 
   // Ifpack2_IlukGraph.hpp
-  params.set("fact: iluk level-of-fill", (int)1);
-  params.set("fact: iluk level-of-overlap", (int)0);
+  params.set("fact: iluk level-of-fill", 1);
+  params.set("fact: iluk level-of-overlap", 0);
 
   // Ifpack2_Chebyshev
-  params.set("chebyshev: max eigenvalue",STS::nan());
-  params.set("chebyshev: ratio eigenvalue",STS::nan());
-  params.set("chebyshev: min eigenvalue",(double)30.0);
-  params.set("chebyshev: degree",(int)1);
-  params.set("chebyshev: eigenvalue max iterations",(int)10);
-  params.set("chebyshev: assume matrix does not change",false);
+  params.set("chebyshev: max eigenvalue", STS::nan());
+  params.set("chebyshev: ratio eigenvalue", STS::nan());
+  params.set("chebyshev: min eigenvalue", 30.0);
+  params.set("chebyshev: degree", 1);
+  params.set("chebyshev: eigenvalue max iterations", 10);
+  params.set("chebyshev: assume matrix does not change", false);
   // params.set("chebyshev: operator inv diagonal",Teuchos::null);
-  params.set("chebyshev: min diagonal value",STS::eps());
-  params.set("chebyshev: zero starting solution",true);
+  params.set("chebyshev: min diagonal value", STS::eps());
+  params.set("chebyshev: zero starting solution", true);
 
   // Ifpack2_Amesos.cpp
   params.set("amesos: solver type", "Amesos_Klu");
 
   // Ifpack2_IC.cpp
-  params.set("fact: level-of-fill", (int)1);
-  params.set("fact: absolute threshold", (double)0.0);
-  params.set("fact: relative threshold", (double)0.0);
-  params.set("fact: drop tolerance", (double)0.0);
+  params.set("fact: level-of-fill", 1);
+  params.set("fact: absolute threshold", 0.0);
+  params.set("fact: relative threshold", 0.0);
+  params.set("fact: drop tolerance", 0.0);
 
   // Ifpack2_ICT.cpp
-  params.set("fact: ict level-of-fill", (double)1.0);
-  params.set("fact: absolute threshold", (double)0.0);
-  params.set("fact: relative threshold", (double)1.0);
-  params.set("fact: relax value", (double)0.0);
-  params.set("fact: drop tolerance", (double)0.0);
+  params.set("fact: ict level-of-fill", 1.0);
+  params.set("fact: absolute threshold", 0.0);
+  params.set("fact: relative threshold", 1.0);
+  params.set("fact: relax value", 0.0);
+  params.set("fact: drop tolerance", 0.0);
 
   // Ifpack2_ILU.cpp
-  params.set("fact: level-of-fill", (int)0);
-  params.set("fact: absolute threshold", (double)0.0);
-  params.set("fact: relative threshold", (double)1.0);
-  params.set("fact: relax value", (double)0.0);
+  params.set("fact: level-of-fill", 0);
+  params.set("fact: absolute threshold", 0.0);
+  params.set("fact: relative threshold", 1.0);
+  params.set("fact: relax value", 0.0);
 
   // Ifpack2_ILUT.cpp
-  params.set("fact: ilut level-of-fill", (double)1.0);
-  params.set("fact: absolute threshold", (double)0.0);
-  params.set("fact: relative threshold", (double)1.0);
-  params.set("fact: relax value", (double)0.0);
+  params.set("fact: ilut level-of-fill", 1.0);
+  params.set("fact: absolute threshold", 0.0);
+  params.set("fact: relative threshold", 1.0);
+  params.set("fact: relax value", 0.0);
 
   // Ifpack2_LocalSparseTriangularSolver.cpp
   params.set("trisolver: type", "Internal");
-  params.set("trisolver: block size", (int)1);
+  params.set("trisolver: block size", 1);
   params.set("trisolver: reverse U", false);
 
   // Overlapping partitioner
-  params.set("partitioner: local parts", (int)1);
-  params.set("partitioner: overlap", (int)0);
-  params.set("partitioner: print level", (int)0);
+  params.set("partitioner: local parts", 1);
+  params.set("partitioner: overlap", 0);
+  params.set("partitioner: print level", 0);
 
   // Ifpack2_Relaxation.cpp
   params.set("relaxation: container", "TriDi");
   params.set("relaxation: type", "Jacobi");
-  params.set("relaxation: sweeps", (int)1);
+  params.set("relaxation: sweeps", 1);
   params.set("relaxation: direction", "forward");
-  params.set("relaxation: damping factor", (double)1.0);
-  params.set("relaxation: min diagonal value", (double)1.0);
+  params.set("relaxation: damping factor", 1.0);
+  params.set("relaxation: min diagonal value", 1.0);
   params.set("relaxation: zero starting solution", true);
-  params.set("relaxation: backward mode",false);
-  params.set("relaxation: use l1",false);
-  params.set("relaxation: l1 eta",(double)1.5);
+  params.set("relaxation: backward mode", false);
+  params.set("relaxation: use l1", false);
+  params.set("relaxation: l1 eta", 1.5);
   params.set("relaxation: banded container superdiagonals", -1);
   params.set("relaxation: banded container subdiagonals", -1);
+  params.set("relaxation: mtgs cluster size", 1);
 
   // Ifpack2_SPARSKIT.cpp
   // ap 25 May 2016: all SPARSKIT for backwards compatibility ONLY
-  params.set("fact: sparskit: lfil", (int)0);
-  params.set("fact: sparskit: tol", (double)0.0);
-  params.set("fact: sparskit: droptol", (double)0.0);
-  params.set("fact: sparskit: permtol", (double)0.1);
-  params.set("fact: sparskit: alph", (double)0.0);
-  params.set("fact: sparskit: mbloc", (int)(-1));
-  params.set("fact: sparskit: type", ("ILUT"));
+  params.set("fact: sparskit: lfil",    0);
+  params.set("fact: sparskit: tol",     0.0);
+  params.set("fact: sparskit: droptol", 0.0);
+  params.set("fact: sparskit: permtol", 0.1);
+  params.set("fact: sparskit: alph",    0.0);
+  params.set("fact: sparskit: mbloc",   -1);
+  params.set("fact: sparskit: type",    "ILUT");
 
   // Additive Schwarz preconditioner
   params.set("schwarz: compute condest", false); // mfh 24 Mar 2015: for backwards compatibility ONLY
   params.set("schwarz: combine mode", "ZERO"); // use string mode for this
   params.set("schwarz: use reordering", true);
   params.set("schwarz: filter singletons", false);
-  params.set("schwarz: overlap level", (int)0);
+  params.set("schwarz: overlap level", 0);
 
   // Ifpack2_BlockRelaxation.hpp
   // params.set("relaxation: type", "Jacobi"); // already set
   // params.set("relaxation: sweeps", 1); // already set
   // params.get("relaxation: damping factor", 1.0); // already set
   // params.get("relaxation: zero starting solution", true); // already set
   params.set("partitioner: type", "greedy");
-  params.set("partitioner: local parts", (int)1);
-  params.set("partitioner: overlap", (int)0);
+  params.set("partitioner: local parts", 1);
+  params.set("partitioner: overlap", 0);
   Teuchos::Array<Teuchos::ArrayRCP<int>> tmp0;
   params.set("partitioner: parts", tmp0);
   params.set("partitioner: maintain sparsity", false);
@@ -175,8 +176,8 @@ void getValidParameters(Teuchos::ParameterList& params)
   params.set("partitioner: map", tmp);
 
   // Ifpack2_LinePartitioner.hpp (FIXME)
-  params.set("partitioner: line detection threshold",(double)0.0);
-  params.set("partitioner: PDE equations",(int)1);
+  params.set("partitioner: line detection threshold", 0.0);
+  params.set("partitioner: PDE equations", 1);
   Teuchos::RCP<Tpetra::MultiVector<> > dummy;
   params.set("partitioner: coordinates",dummy);
 }

packages/ifpack2/src/Ifpack2_Relaxation_decl.hpp

@@ -761,71 +761,73 @@ class Relaxation :
   Teuchos::RCP<block_multivector_type> yBlockColumnPointMap_;
 
   //! How many times to apply the relaxation per apply() call.
-  int NumSweeps_;
+  int NumSweeps_ = 1;
   //! Which relaxation method to use.
-  Details::RelaxationType PrecType_;
+  Details::RelaxationType PrecType_ = Ifpack2::Details::JACOBI;
   //! Damping factor
-  scalar_type DampingFactor_;
+  scalar_type DampingFactor_ = STS::one();
   //! If \c true, more than 1 processor is currently used.
   bool IsParallel_;
   //! If \c true, the starting solution is always the zero vector.
-  bool ZeroStartingSolution_;
+  bool ZeroStartingSolution_ = true;
   //! If true, do backward-mode Gauss-Seidel.
-  bool DoBackwardGS_;
+  bool DoBackwardGS_ = false;
   //! If true, do the L1 version of Jacobi, Gauss-Seidel, or symmetric Gauss-Seidel.
-  bool DoL1Method_;
+  bool DoL1Method_ = false;
   //! Eta parameter for modified L1 method
-  magnitude_type L1Eta_;
+  magnitude_type L1Eta_ = Teuchos::as<magnitude_type>(1.5);
   //! Minimum diagonal value
-  scalar_type MinDiagonalValue_;
+  scalar_type MinDiagonalValue_ = STS::zero();
   //! Whether to fix up zero or tiny diagonal entries.
-  bool fixTinyDiagEntries_;
+  bool fixTinyDiagEntries_ = false;
   //! Whether to spend extra effort and all-reduces checking diagonal entries.
-  bool checkDiagEntries_;
+  bool checkDiagEntries_ = false;
+  //! For MTSGS, the cluster size (use point coloring if equal to 1)
+  int clusterSize_ = 1;
 
   //!Wheter the provided matrix is structurally symmetric or not.
-  bool is_matrix_structurally_symmetric_;
+  bool is_matrix_structurally_symmetric_ = false;
 
   //!Whether to write the given input file
-  bool ifpack2_dump_matrix_;
+  bool ifpack2_dump_matrix_ = false;
 
 
   //! If \c true, the preconditioner has been initialized successfully.
-  bool isInitialized_;
+  bool isInitialized_ = false;
   //! If \c true, the preconditioner has been computed successfully.
-  bool IsComputed_;
+  bool IsComputed_ = false;
   //! The number of successful calls to initialize().
-  int NumInitialize_;
+  int NumInitialize_ = 0;
   //! the number of successful calls to compute().
-  int NumCompute_;
+  int NumCompute_ = 0;
   //! The number of successful calls to apply().
-  mutable int NumApply_;
+  mutable int NumApply_ = 0;
   //! Total time in seconds for all successful calls to initialize().
-  double InitializeTime_;
+  double InitializeTime_ = 0.0;
   //! Total time in seconds for all successful calls to compute().
-  double ComputeTime_;
+  double ComputeTime_ = 0.0;
   //! Total time in seconds for all successful calls to apply().
-  mutable double ApplyTime_;
+  mutable double ApplyTime_ = 0.0;
   //! The total number of floating-point operations for all successful calls to compute().
-  double ComputeFlops_;
+  double ComputeFlops_ = 0.0;
   //! The total number of floating-point operations for all successful calls to apply().
-  mutable double ApplyFlops_;
+  mutable double ApplyFlops_ = 0.0;
 
   //! Global magnitude of the diagonal entry with the minimum magnitude.
-  magnitude_type globalMinMagDiagEntryMag_;
+  magnitude_type globalMinMagDiagEntryMag_ = STM::zero();
   //! Global magnitude of the diagonal entry with the maximum magnitude.
-  magnitude_type globalMaxMagDiagEntryMag_;
+  magnitude_type globalMaxMagDiagEntryMag_ = STM::zero();
   //! Global number of small (in magnitude) diagonal entries detected by compute().
-  size_t globalNumSmallDiagEntries_;
+  size_t globalNumSmallDiagEntries_ = 0;
   //! Global number of zero diagonal entries detected by compute().
-  size_t globalNumZeroDiagEntries_;
+  size_t globalNumZeroDiagEntries_ = 0;
   //! Global number of negative (real part) diagonal entries detected by compute().
-  size_t globalNumNegDiagEntries_;
+  size_t globalNumNegDiagEntries_ = 0;
   /// \brief Absolute two-norm difference between computed and actual inverse diagonal.
   ///
   /// "Actual inverse diagonal" means the result of 1/diagonal,
   /// without any protection against zero or small diagonal entries.
-  magnitude_type globalDiagNormDiff_;
+  magnitude_type globalDiagNormDiff_ = STM::zero();
 
   /// \brief Precomputed offsets of local diagonal entries of the matrix.
   ///
@@ -839,9 +841,9 @@ class Relaxation :
   /// We need this flag because it is not enough just to test if
   /// diagOffsets_ has size zero.  It is perfectly legitimate for the
   /// matrix to have zero rows on the calling process.
-  bool savedDiagOffsets_;
+  bool savedDiagOffsets_ = false;
 
-  bool hasBlockCrsMatrix_;
+  bool hasBlockCrsMatrix_ = false;
 
   /// \brief In case of local/reordered smoothing, the unknowns to use
   Teuchos::ArrayRCP<local_ordinal_type> localSmoothingIndices_;