Ifpack2: Improve LocalSparseTriangularSolver unit test

@trilinos/ifpack2 Test PASSES.

Add a new unit test for LocalSparseTriangularSolver.  It exercises
lower and upper triangular solves with a matrix whose LU factorization
is known and exact (in binary floating-point arithmetic), and a
resulting linear system whose solution is exact (again in binary
floating-point arithmetic).

Exactness means we don't need an error tolerance.  (This assumes that
the triangular solve is exact.  It's possible to do approximate
triangular solves.  We ignore this case for now.)

This relates to #558 and other outstanding issues.  Given that
LocalSparseTriangularSolver behaves correctly, we can now focus on
classes that use it.  For example, it would make sense to extend the
"arrow matrix" test case in this new test, to include "remote" entries
that Ifpack2::LocalFilter should filter out.

I also instrumented Ifpack2::LocalSparseTriangularSolver with optional
debug output.  It takes a Teuchos::FancyOStream, which helps us get
debug output in sync with unit test output.

Build/Test Cases Summary
Enabled Packages: Ifpack2
Disabled Packages: FEI,PyTrilinos,Moertel,STK,SEACAS,ThreadPool,OptiPack,Rythmos,Intrepid,ROL
0) MPI_DEBUG => passed: passed=29,notpassed=0 (7.80 min)
1) SERIAL_RELEASE => passed: passed=25,notpassed=0 (6.10 min)

packages/ifpack2/src/Ifpack2_LocalSparseTriangularSolver_decl.hpp

@@ -45,6 +45,7 @@
 
 #include "Ifpack2_Preconditioner.hpp"
 #include "Ifpack2_Details_CanChangeMatrix.hpp"
+#include "Teuchos_FancyOStream.hpp"
 #include <type_traits>
 
 #ifndef DOXYGEN_SHOULD_SKIP_THIS
@@ -143,6 +144,17 @@ class LocalSparseTriangularSolver :
  /// those properties.
  LocalSparseTriangularSolver (const Teuchos::RCP<const row_matrix_type>& A);
 
+ /// \brief Constructor that takes an optional debug output stream.
+ ///
+ /// \param A [in] The input sparse matrix. Though its type is
+ /// Tpetra::RowMatrix for consistency with other Ifpack2 solvers,
+ /// this must be a Tpetra::CrsMatrix specialization.
+ ///
+ /// \param out [in/out] Optional debug output stream. If nonnull,
+ /// this solver will print copious debug output to the stream.
+ LocalSparseTriangularSolver (const Teuchos::RCP<const row_matrix_type>& A,
+ const Teuchos::RCP<Teuchos::FancyOStream>& out);
+
  //! Destructor (virtual for memory safety).
  virtual ~LocalSparseTriangularSolver ();
 
@@ -296,11 +308,14 @@ class LocalSparseTriangularSolver :
 private:
  //! The original input matrix.
  Teuchos::RCP<const row_matrix_type> A_;
+ //! Debug output stream; may be null (not used in that case)
+ Teuchos::RCP<Teuchos::FancyOStream> out_;
  //! The original input matrix, as a Tpetra::CrsMatrix.
  Teuchos::RCP<const Tpetra::CrsMatrix<scalar_type,
  local_ordinal_type,
  global_ordinal_type,
  node_type, false> > A_crs_;
+
  typedef Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type> MV;
  mutable Teuchos::RCP<MV> X_colMap_;
  mutable Teuchos::RCP<MV> Y_rowMap_;

packages/ifpack2/src/Ifpack2_LocalSparseTriangularSolver_def.hpp

@@ -74,6 +74,38 @@ LocalSparseTriangularSolver (const Teuchos::RCP<const row_matrix_type>& A) :
  }
 }
 
+template<class MatrixType>
+LocalSparseTriangularSolver<MatrixType>::
+LocalSparseTriangularSolver (const Teuchos::RCP<const row_matrix_type>& A,
+ const Teuchos::RCP<Teuchos::FancyOStream>& out) :
+ A_ (A),
+ out_ (out),
+ isInitialized_ (false),
+ isComputed_ (false),
+ numInitialize_ (0),
+ numCompute_ (0),
+ numApply_ (0),
+ initializeTime_ (0.0),
+ computeTime_ (0.0),
+ applyTime_ (0.0)
+{
+ if (! out_.is_null ()) {
+ *out_ << ">>> DEBUG Ifpack2::LocalSparseTriangularSolver constructor"
+ << std::endl;
+ }
+ typedef typename Tpetra::CrsMatrix<scalar_type, local_ordinal_type,
+ global_ordinal_type, node_type> crs_matrix_type;
+ if (! A.is_null ()) {
+ Teuchos::RCP<const crs_matrix_type> A_crs =
+ Teuchos::rcp_dynamic_cast<const crs_matrix_type> (A);
+ TEUCHOS_TEST_FOR_EXCEPTION
+ (A_crs.is_null (), std::invalid_argument,
+ "Ifpack2::LocalSparseTriangularSolver constructor: "
+ "The input matrix A is not a Tpetra::CrsMatrix.");
+ A_crs_ = A_crs;
+ }
+}
+
 template<class MatrixType>
 LocalSparseTriangularSolver<MatrixType>::
 ~LocalSparseTriangularSolver ()
@@ -91,6 +123,9 @@ LocalSparseTriangularSolver<MatrixType>::
 initialize ()
 {
  const char prefix[] = "Ifpack2::LocalSparseTriangularSolver::initialize: ";
+ if (! out_.is_null ()) {
+ *out_ << ">>> DEBUG " << prefix << std::endl;
+ }
 
  TEUCHOS_TEST_FOR_EXCEPTION
  (A_.is_null (), std::runtime_error, prefix << "You must call "
@@ -128,6 +163,9 @@ LocalSparseTriangularSolver<MatrixType>::
 compute ()
 {
  const char prefix[] = "Ifpack2::LocalSparseTriangularSolver::compute: ";
+ if (! out_.is_null ()) {
+ *out_ << ">>> DEBUG " << prefix << std::endl;
+ }
 
  TEUCHOS_TEST_FOR_EXCEPTION
  (A_.is_null (), std::runtime_error, prefix << "You must call "
@@ -169,6 +207,23 @@ apply (const Tpetra::MultiVector<scalar_type, local_ordinal_type,
  typedef scalar_type ST;
  typedef Teuchos::ScalarTraits<ST> STS;
  const char prefix[] = "Ifpack2::LocalSparseTriangularSolver::apply: ";
+ if (! out_.is_null ()) {
+ *out_ << ">>> DEBUG " << prefix;
+ if (A_crs_.is_null ()) {
+ *out_ << "A_crs_ is null!" << std::endl;
+ }
+ else {
+ const std::string uplo = A_crs_->isUpperTriangular () ? "U" :
+ (A_crs_->isLowerTriangular () ? "L" : "N");
+ const std::string trans = (mode == Teuchos::CONJ_TRANS) ? "C" :
+ (mode == Teuchos::TRANS ? "T" : "N");
+ const std::string diag =
+ (A_crs_->getNodeNumDiags () < A_crs_->getNodeNumRows ()) ? "U" : "N";
+ *out_ << "uplo=\"" << uplo
+ << "\", trans=\"" << trans
+ << "\", diag=\"" << diag << "\"" << std::endl;
+ }
+ }
 
  TEUCHOS_TEST_FOR_EXCEPTION
  (! isComputed (), std::runtime_error, prefix << "If compute() has not yet "