Merge 'trilinos/Trilinos:develop' (2e8c2c0) into 'tcad-charon/Trilino…

…s:develop' (02d09c3). * trilinos-develop: Piro: enable "Explicit Matrix Transpose" option when computing senstivities w/out ROL (trilinos#9427) Tpetra: Fixing deprecated code
tcad-charon · Jul 15, 2021 · 5db8a9a · 5db8a9a
2 parents 02d09c3 + 2e8c2c0
commit 5db8a9a
Show file tree

Hide file tree

Showing 2 changed files with 38 additions and 13 deletions.
diff --git a/packages/piro/src/Piro_SteadyStateSolver_Def.hpp b/packages/piro/src/Piro_SteadyStateSolver_Def.hpp
@@ -541,6 +541,7 @@ void Piro::SteadyStateSolver<Scalar>::evalConvergedModelResponsesAndSensitivitie
       }
     }
   } else if(computeAdjointSensitivities) {
+
     // Compute adjoint layouts of df/dp, dg/dx depending on
     for (int i=0; i<num_p_; i++) {
       // p
@@ -900,6 +901,13 @@ void Piro::SteadyStateSolver<Scalar>::evalConvergedModelResponsesAndSensitivitie
     }
   } else if (computeAdjointSensitivities) {
 
+    bool explicitlyTransposeMatrix = appParams.isParameter("Enable Explicit Matrix Transpose") ?
+        appParams.get<bool>("Enable Explicit Matrix Transpose") :
+        false;
+
+    if(explicitlyTransposeMatrix)
+      appParams.sublist("Optimization Status").set("Compute Transposed Jacobian", true);
+
     // Create implicitly transpose Jacobian and preconditioner
     Teuchos::RCP< const Thyra::LinearOpWithSolveFactoryBase<double> > lows_factory = model_->get_W_factory();
     TEUCHOS_ASSERT(Teuchos::nonnull(lows_factory));
@@ -937,15 +945,24 @@ void Piro::SteadyStateSolver<Scalar>::evalConvergedModelResponsesAndSensitivitie
       model_->evalModel(in_args, out_args);
     }
 
-    if(Teuchos::nonnull(prec))
-      Thyra::initializePreconditionedOp<double>(*lows_factory,
-          Thyra::transpose<double>(lop),
-          Thyra::unspecifiedPrec<double>(::Thyra::transpose<double>(prec->getUnspecifiedPrecOp())),
-          jacobian.ptr());
-    else
-      Thyra::initializeOp<double>(*lows_factory,
-          Thyra::transpose<double>(lop),
+    if(Teuchos::nonnull(prec)) {
+      if(explicitlyTransposeMatrix) {
+        Thyra::initializePreconditionedOp<double>(*lows_factory,
+          lop,
+          prec,
           jacobian.ptr());
+      } else {
+        Thyra::initializePreconditionedOp<double>(*lows_factory,
+            Thyra::transpose<double>(lop),
+            Thyra::unspecifiedPrec<double>(::Thyra::transpose<double>(prec->getUnspecifiedPrecOp())),
+            jacobian.ptr());
+      }
+    } else {
+      if(explicitlyTransposeMatrix)
+        Thyra::initializeOp<double>(*lows_factory, lop, jacobian.ptr());
+      else
+        Thyra::initializeOp<double>(*lows_factory, Thyra::transpose<double>(lop), jacobian.ptr());
+    }
 
 
     for (int j=0; j<num_g_; j++) {
@@ -1138,6 +1155,8 @@ void Piro::SteadyStateSolver<Scalar>::evalConvergedModelResponsesAndSensitivitie
         }
       }
     }
+    if(explicitlyTransposeMatrix)
+      appParams.sublist("Optimization Status").set("Compute Transposed Jacobian", false);
   }
 }
 

diff --git a/packages/tpetra/core/test/FECrsMatrix/FECrsMatrix_UnitTests.cpp b/packages/tpetra/core/test/FECrsMatrix/FECrsMatrix_UnitTests.cpp
@@ -310,6 +310,9 @@ TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL( FECrsMatrix, Assemble1D_2, LO, GO, Scalar, No
   using FEMAT = typename Tpetra::FECrsMatrix<Scalar,LO,GO,Node>;
   using CMAT = typename Tpetra::CrsMatrix<Scalar,LO,GO,Node>;
   using FEG = typename Tpetra::FECrsGraph<LO,GO,Node>;
+  using cols_type = typename CMAT::nonconst_local_inds_host_view_type;
+  using vals_type = typename CMAT::nonconst_values_host_view_type;
+
 
   // get a comm
   RCP<const Comm<int> > comm = getDefaultComm();
@@ -367,14 +370,17 @@ TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL( FECrsMatrix, Assemble1D_2, LO, GO, Scalar, No
   if (fe_matrix.getRowMap()->isNodeLocalElement(local_row))
   {
     auto num_entries = fe_matrix.getNodeNumEntries();
-    Teuchos::Array<LO> cols(num_entries);
+    cols_type cols("cols",num_entries);
     Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
     Scalar one = Teuchos::ScalarTraits<Scalar>::one();
-    Teuchos::Array<Scalar> vals(num_entries, zero);
-    fe_matrix.getLocalRowCopy(local_row, cols(), vals(), num_entries);
-    std::fill(vals.begin(), vals.end(), zero);
+    vals_type vals("vals",num_entries);
+    fe_matrix.getLocalRowCopy(local_row, cols, vals, num_entries);
+    Kokkos::deep_copy(vals,zero);
     vals[0] = one;
-    fe_matrix.replaceLocalValues(local_row, cols(), vals());
+    auto cols_sub = Kokkos::subview(cols,Kokkos::make_pair((size_t)0,num_entries));
+    auto vals_sub = Kokkos::subview(vals,Kokkos::make_pair((size_t)0,num_entries));
+
+    fe_matrix.replaceLocalValues(local_row, cols_sub, vals_sub);
   }
   fe_matrix.endModify();
   TPETRA_GLOBAL_SUCCESS_CHECK(out,comm,success)