Merge Pull Request trilinos#12769 from cgcgcg/Trilinos/tpetraAddTest

Automatically Merged using Trilinos Pull Request AutoTester
PR Title: b'Tpetra, Teko: SpAdd issues'
PR Author: cgcgcg

Author: trilinos-autotester

packages/teko/src/Teko_Utilities.cpp

@@ -2195,34 +2195,43 @@ const ModifiableLinearOp explicitAdd(const LinearOp &opl_in, const LinearOp &opr
     // Build output operator
     RCP<Thyra::LinearOpBase<ST>> explicitOp;
     RCP<Tpetra::CrsMatrix<ST, LO, GO, NT>> explicitCrsOp;
-    if (destOp != Teuchos::null) {
+    if (!destOp.is_null()) {
       explicitOp = destOp;
-      try {  // try to reuse matrix sparsity with Add. If it fails, build new operator with add
-        RCP<Thyra::TpetraLinearOp<ST, LO, GO, NT>> tOp =
-            rcp_dynamic_cast<Thyra::TpetraLinearOp<ST, LO, GO, NT>>(destOp);
+      RCP<Thyra::TpetraLinearOp<ST, LO, GO, NT>> tOp =
+          rcp_dynamic_cast<Thyra::TpetraLinearOp<ST, LO, GO, NT>>(destOp);
+      if (!tOp.is_null())
         explicitCrsOp =
-            rcp_dynamic_cast<Tpetra::CrsMatrix<ST, LO, GO, NT>>(tOp->getTpetraOperator(), true);
-        Tpetra::MatrixMatrix::Add<ST, LO, GO, NT>(*tCrsOpl, transpl, scalarl, *tCrsOpr, transpr,
-                                                  scalarr, explicitCrsOp);
-      } catch (std::logic_error &e) {
-        RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
-        *out << "*** THROWN EXCEPTION ***\n";
-        *out << e.what() << std::endl;
-        *out << "************************\n";
-        *out << "Teko: explicitAdd unable to reuse existing operator. Creating new operator.\n"
-             << std::endl;
+            rcp_dynamic_cast<Tpetra::CrsMatrix<ST, LO, GO, NT>>(tOp->getTpetraOperator());
+      bool needNewTpetraMatrix =
+          (explicitCrsOp.is_null()) || (tCrsOpl == explicitCrsOp) || (tCrsOpr == explicitCrsOp);
+      if (!needNewTpetraMatrix) {
+        try {
+          // try to reuse matrix sparsity with Add. If it fails, build new operator with add
+          Tpetra::MatrixMatrix::Add<ST, LO, GO, NT>(*tCrsOpl, transpl, scalarl, *tCrsOpr, transpr,
+                                                    scalarr, explicitCrsOp);
+        } catch (std::logic_error &e) {
+          RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
+          *out << "*** THROWN EXCEPTION ***\n";
+          *out << e.what() << std::endl;
+          *out << "************************\n";
+          *out << "Teko: explicitAdd unable to reuse existing operator. Creating new operator.\n"
+               << std::endl;
+          needNewTpetraMatrix = true;
+        }
+      }
+      if (needNewTpetraMatrix)
+        // Do explicit matrix-matrix add
         explicitCrsOp = Tpetra::MatrixMatrix::add<ST, LO, GO, NT>(scalarl, transpl, *tCrsOpl,
                                                                   scalarr, transpr, *tCrsOpr);
-      }
     } else {
-      explicitOp    = rcp(new Thyra::TpetraLinearOp<ST, LO, GO, NT>());
+      explicitOp = rcp(new Thyra::TpetraLinearOp<ST, LO, GO, NT>());
+      // Do explicit matrix-matrix add
       explicitCrsOp = Tpetra::MatrixMatrix::add<ST, LO, GO, NT>(scalarl, transpl, *tCrsOpl, scalarr,
                                                                 transpr, *tCrsOpr);
     }
     RCP<Thyra::TpetraLinearOp<ST, LO, GO, NT>> tExplicitOp =
         rcp_dynamic_cast<Thyra::TpetraLinearOp<ST, LO, GO, NT>>(explicitOp);
 
-    // Do explicit matrix-matrix add
     tExplicitOp->initialize(Thyra::tpetraVectorSpace<ST, LO, GO, NT>(explicitCrsOp->getRangeMap()),
                             Thyra::tpetraVectorSpace<ST, LO, GO, NT>(explicitCrsOp->getDomainMap()),
                             explicitCrsOp);

packages/teko/tests/src/tExplicitOps_tpetra.cpp

@@ -463,6 +463,23 @@ bool tExplicitOps_tpetra::test_add_mod(int verbosity, std::ostream& os) {
       Teuchos::rcp_dynamic_cast<const Thyra::TpetraLinearOp<ST, LO, GO, NT> >(expOp2, true);
   RCP<const Tpetra::Operator<ST, LO, GO, NT> > eop3 = tOp3->getConstTpetraOperator();
 
+  // Create a target linear operator, with a known sparsity pattern (diagonal) for a second time
+  expOp2 = Teko::explicitAdd(H_, H_, expOp2);
+  RCP<const Thyra::TpetraLinearOp<ST, LO, GO, NT> > tOp3b =
+      Teuchos::rcp_dynamic_cast<const Thyra::TpetraLinearOp<ST, LO, GO, NT> >(expOp2, true);
+  RCP<const Tpetra::Operator<ST, LO, GO, NT> > eop3b = tOp3b->getConstTpetraOperator();
+
+  // check that underlying pointers are the same
+  {
+    std::stringstream ss;
+    Teuchos::FancyOStream fos(rcpFromRef(ss), "      |||");
+    TEST_ASSERT(eop3.getRawPtr() == eop3b.getRawPtr(),
+                std::endl
+                    << " tExplicitOps_tpetra::test_add_mod"
+                    << ": Testing matrix addition returns new pointer");
+    if (not(eop3.getRawPtr() == eop3b.getRawPtr()) || verbosity >= 10) os << ss.str();
+  }
+
   // Try to add matricies with sparsity patterns that differ from the target operator
   expOp2 = Teko::explicitAdd(Teko::scale(-4.0, F_), Teko::adjoint(G_), expOp2);
   RCP<const Thyra::TpetraLinearOp<ST, LO, GO, NT> > tOp4 =
@@ -493,6 +510,52 @@ bool tExplicitOps_tpetra::test_add_mod(int verbosity, std::ostream& os) {
     if (not result || verbosity >= 10) os << ss.str();
   }
 
+  // Create a target linear operator, with a known sparsity pattern (diagonal) and aliased args
+  expOp2 = Teko::explicitAdd(H_, H_, expOp2);
+  tOp3   = Teuchos::rcp_dynamic_cast<const Thyra::TpetraLinearOp<ST, LO, GO, NT> >(expOp2, true);
+  eop3   = tOp3->getConstTpetraOperator();
+  auto expOp3 = Teko::explicitAdd(H_, expOp2, expOp2);
+  RCP<const Thyra::TpetraLinearOp<ST, LO, GO, NT> > tOp5 =
+      Teuchos::rcp_dynamic_cast<const Thyra::TpetraLinearOp<ST, LO, GO, NT> >(expOp3, true);
+  RCP<const Tpetra::Operator<ST, LO, GO, NT> > eop5 = tOp5->getConstTpetraOperator();
+
+  // check that underlying pointers are different
+  // since the sparsity pattern of expOp entering the explicitAdd
+  // does not match the sparsity pattern returned by the explicitAdd
+  // a new tpetra operator will be created
+  {
+    std::stringstream ss;
+    Teuchos::FancyOStream fos(rcpFromRef(ss), "      |||");
+    TEST_ASSERT(eop3.getRawPtr() != eop5.getRawPtr(),
+                std::endl
+                    << " tExplicitOps_tpetra::test_add_mod2"
+                    << ": Testing matrix addition returns new pointer");
+    if (not(eop3.getRawPtr() == eop5.getRawPtr()) || verbosity >= 10) os << ss.str();
+  }
+
+  // Create a target linear operator, with a known sparsity pattern (diagonal) and aliased args
+  expOp2 = Teko::explicitAdd(H_, H_, expOp2);
+  tOp3   = Teuchos::rcp_dynamic_cast<const Thyra::TpetraLinearOp<ST, LO, GO, NT> >(expOp2, true);
+  eop3   = tOp3->getConstTpetraOperator();
+  auto expOp4 = Teko::explicitAdd(expOp2, H_, expOp2);
+  RCP<const Thyra::TpetraLinearOp<ST, LO, GO, NT> > tOp6 =
+      Teuchos::rcp_dynamic_cast<const Thyra::TpetraLinearOp<ST, LO, GO, NT> >(expOp4, true);
+  RCP<const Tpetra::Operator<ST, LO, GO, NT> > eop6 = tOp6->getConstTpetraOperator();
+
+  // check that underlying pointers are different
+  // since the sparsity pattern of expOp entering the explicitAdd
+  // does not match the sparsity pattern returned by the explicitAdd
+  // a new tpetra operator will be created
+  {
+    std::stringstream ss;
+    Teuchos::FancyOStream fos(rcpFromRef(ss), "      |||");
+    TEST_ASSERT(eop3.getRawPtr() != eop6.getRawPtr(),
+                std::endl
+                    << " tExplicitOps_tpetra::test_add_mod3"
+                    << ": Testing matrix addition returns new pointer");
+    if (not(eop3.getRawPtr() == eop6.getRawPtr()) || verbosity >= 10) os << ss.str();
+  }
+
   return allPassed;
 }
 

packages/tpetra/core/test/MatrixMatrix/MatrixMatrix_UnitTests.cpp

@@ -44,6 +44,7 @@
 #include <Teuchos_UnitTestHarness.hpp>
 
 #include "TpetraExt_MatrixMatrix.hpp"
+#include "Teuchos_Assert.hpp"
 #include "TpetraExt_TripleMatrixMultiply.hpp"
 #include "Tpetra_MatrixIO.hpp"
 #include "Tpetra_Core.hpp"
@@ -436,6 +437,56 @@ add_test_results add_into_test(
   return toReturn;
 }
 
+
+template<class Matrix_t>
+add_test_results reuse_add_test(
+    const std::string& name,
+    RCP<Matrix_t > A,
+    RCP<Matrix_t > B,
+    bool AT,
+    bool BT,
+    RCP<Matrix_t > C,
+    RCP<const Comm<int> > comm)
+{
+  typedef typename Matrix_t::scalar_type SC;
+  typedef typename Matrix_t::local_ordinal_type LO;
+  typedef typename Matrix_t::global_ordinal_type GO;
+  typedef typename Matrix_t::node_type NT;
+  typedef Map<LO,GO,NT> Map_t;
+
+  add_test_results toReturn;
+  toReturn.correctNorm = C->getFrobeniusNorm ();
+
+  RCP<const Map_t > rowmap = AT ? A->getDomainMap() : A->getRowMap();
+  size_t estSize = A->getGlobalMaxNumRowEntries() + B->getGlobalMaxNumRowEntries();
+         // estSize is upper bound for A, B; estimate only for AT, BT.
+  RCP<Matrix_t> computedC = rcp( new Matrix_t(rowmap, estSize));
+
+  SC one = Teuchos::ScalarTraits<SC>::one();
+  Tpetra::MatrixMatrix::Add(*A, AT, one, *B, BT, one, computedC);
+  computedC->fillComplete(A->getDomainMap(), A->getRangeMap());
+
+  // Call Add a second time
+  Tpetra::MatrixMatrix::Add(*A, AT, one, *B, BT, one, computedC);
+
+  TEUCHOS_ASSERT(A->getDomainMap()->isSameAs(*computedC->getDomainMap()));
+  TEUCHOS_ASSERT(A->getRangeMap()->isSameAs(*computedC->getRangeMap()));
+
+  toReturn.computedNorm = computedC->getFrobeniusNorm ();
+  toReturn.epsilon = fabs(toReturn.correctNorm - toReturn.computedNorm);
+
+#if 0
+  Tpetra::MatrixMarket::Writer<Matrix_t>::writeSparseFile(
+    name+"_calculated.mtx", computedC);
+  Tpetra::MatrixMarket::Writer<Matrix_t>::writeSparseFile(
+    name+"_real.mtx", C);
+#endif
+
+
+  return toReturn;
+}
+
+
 template<class Matrix_t>
 mult_test_results multiply_test_manualfc(
   const std::string& name,
@@ -1194,6 +1245,18 @@ TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(Tpetra_MatMat, operations_test,SC,LO, GO, NT)
       newOut << "\tComputed norm: " << results.computedNorm << endl;
       newOut << "\tEpsilon: " << results.epsilon << endl;
 
+      if (verbose)
+        newOut << "Running 3-argument add reuse test (nonnull C on input) for "
+               << currentSystem.name() << endl;
+
+      results = reuse_add_test(name+"_add",A, B, AT, BT, C, comm);
+
+      TEST_COMPARE(results.epsilon, <, epsilon);
+      newOut << "Reuse Add Test Results: " << endl;
+      newOut << "\tCorrect Norm: " << results.correctNorm << endl;
+      newOut << "\tComputed norm: " << results.computedNorm << endl;
+      newOut << "\tEpsilon: " << results.epsilon << endl;
+
       // FIXME (mfh 09 May 2013) This test doesn't currently pass.  I
       // don't think anyone ever exercised the case where C is null on
       // input before.  I'm disabling this test for now until I have a