diff --git a/packages/piro/src/Piro_PerformAnalysis.cpp b/packages/piro/src/Piro_PerformAnalysis.cpp
index e20f99fdf37f..a3db2d7d0ade 100644
--- a/packages/piro/src/Piro_PerformAnalysis.cpp
+++ b/packages/piro/src/Piro_PerformAnalysis.cpp
@@ -210,8 +210,19 @@ Piro::PerformROLAnalysis(
     RCP< Thyra::VectorBase<double> >& p)
 {
   auto rolParams = analysisParams.sublist("ROL");
+
 #ifdef HAVE_PIRO_ROL
 
+  int verbose = rolParams.get<int>("Verbosity Level", 3);
+  Teuchos::EVerbosityLevel verbosityLevel;
+  switch(verbose) {
+    case 1: verbosityLevel= Teuchos::VERB_LOW; break;
+    case 2: verbosityLevel= Teuchos::VERB_MEDIUM; break;
+    case 3: verbosityLevel= Teuchos::VERB_HIGH; break;
+    case 4: verbosityLevel= Teuchos::VERB_EXTREME; break;
+    default: verbosityLevel= Teuchos::VERB_NONE;
+  }
+
   if(rolParams.isParameter("Use Old Reduced Space Interface") && rolParams.get<bool>("Use Old Reduced Space Interface")) {
 
     using std::string;
@@ -244,7 +255,7 @@ Piro::PerformROLAnalysis(
     ROL::ThyraVector<double> rol_p(p_prod);
 
 
-    ROL::ThyraProductME_Objective<double> obj(piroModel, g_index, p_indices, Teuchos::rcp(&analysisParams.sublist("Optimization Status"),false));
+    ROL::ThyraProductME_Objective<double> obj(piroModel, g_index, p_indices, Teuchos::rcp(&analysisParams.sublist("Optimization Status"),false),verbosityLevel);
 
 
     bool print = rolParams.get<bool>("Print Output", false);
@@ -450,9 +461,8 @@ Piro::PerformROLAnalysis(
     Teuchos::RCP<Thyra::VectorBase<double>> lambda_vec = Thyra::createMember(x_space);
     ROL::ThyraVector<double> rol_lambda(lambda_vec);
 
-    bool always_recompute = true;
-    ThyraProductME_Objective_SimOpt<double> obj(*model, g_index, p_indices, Teuchos::rcp(&analysisParams.sublist("Optimization Status"),false));
-    ThyraProductME_Constraint_SimOpt<double> constr(*model, g_index, p_indices, Teuchos::rcp(&analysisParams.sublist("Optimization Status"),false),always_recompute);
+    ThyraProductME_Objective_SimOpt<double> obj(*model, g_index, p_indices, Teuchos::rcp(&analysisParams.sublist("Optimization Status"),false),verbosityLevel);
+    ThyraProductME_Constraint_SimOpt<double> constr(*model, g_index, p_indices, Teuchos::rcp(&analysisParams.sublist("Optimization Status"),false),verbosityLevel);
 
     constr.setSolveParameters(rolParams.sublist("ROL Options"));
 
diff --git a/packages/rol/adapters/thyra/src/function/ROL_ThyraProductME_Constraint_SimOpt.hpp b/packages/rol/adapters/thyra/src/function/ROL_ThyraProductME_Constraint_SimOpt.hpp
index da6a5e0ae984..087daa2f8012 100644
--- a/packages/rol/adapters/thyra/src/function/ROL_ThyraProductME_Constraint_SimOpt.hpp
+++ b/packages/rol/adapters/thyra/src/function/ROL_ThyraProductME_Constraint_SimOpt.hpp
@@ -52,6 +52,7 @@
 #include "ROL_Constraint_SimOpt.hpp"
 
 #include "Tpetra_CrsMatrix.hpp"
+#include "Teuchos_VerbosityLevel.hpp"
 
 using namespace ROL;
 
@@ -61,16 +62,22 @@ class ThyraProductME_Constraint_SimOpt : public Constraint_SimOpt<Real> {
 
 public:
 
-  ThyraProductME_Constraint_SimOpt(Thyra::ModelEvaluatorDefaultBase<double>& thyra_model_, int g_index_, const std::vector<int>& p_indices_,Teuchos::RCP<Teuchos::ParameterList> params_ = Teuchos::null, bool recompute = false) :
-    thyra_model(thyra_model_), g_index(g_index_), p_indices(p_indices_), params(params_) {
+  ThyraProductME_Constraint_SimOpt(Thyra::ModelEvaluatorDefaultBase<double>& thyra_model_, int g_index_, const std::vector<int>& p_indices_,
+    Teuchos::RCP<Teuchos::ParameterList> params_ = Teuchos::null, Teuchos::EVerbosityLevel verbLevel= Teuchos::VERB_HIGH) :
+        thyra_model(thyra_model_), g_index(g_index_), p_indices(p_indices_), params(params_),
+        out(Teuchos::VerboseObjectBase::getDefaultOStream()),
+        verbosityLevel(verbLevel){
     thyra_solver = Teuchos::null;
-    updateValue = updateJacobian1 = true;
-    value_ = 0;
+    computeValue = computeJacobian1 = solveConstraint = true;
     num_responses = -1;
-    x_ptr = Teuchos::null;
-    grad_ptr = Teuchos::null;
-    if(params != Teuchos::null)
+    value_ptr_ = Teuchos::null;
+    rol_u_ptr = Teuchos::null;
+    rol_z_ptr = Teuchos::null;
+    jac1 = Teuchos::null;
+    if(params != Teuchos::null) {
       params->set<int>("Optimizer Iteration Number", -1);
+      params->set<Teuchos::RCP<Vector<Real> > >("Optimization Variable", Teuchos::null);
+    }
   };
 
   void setExternalSolver(Teuchos::RCP<Thyra::ModelEvaluatorDefaultBase<double>> thyra_solver_) {
@@ -82,55 +89,83 @@ class ThyraProductME_Constraint_SimOpt : public Constraint_SimOpt<Real> {
   }
 
   void value(Vector<Real> &c, const Vector<Real> &u, const Vector<Real> &z, Real &tol) {
+#ifdef  HAVE_ROL_DEBUG
+    //u and z should be updated in the update functions before calling value
+    TEUCHOS_ASSERT(!u_hasChanged(u));
+    TEUCHOS_ASSERT(!z_hasChanged(z));
+#endif
 
-    const ThyraVector<Real>  & thyra_p = dynamic_cast<const ThyraVector<Real>&>(z);
-    const ThyraVector<Real>  & thyra_x = dynamic_cast<const ThyraVector<Real>&>(u);
-    ThyraVector<Real>  & thyra_f = dynamic_cast<ThyraVector<Real>&>(c);
-    Teuchos::RCP<const Thyra::ProductVectorBase<Real> > thyra_prodvec_p = Teuchos::rcp_dynamic_cast<const Thyra::ProductVectorBase<Real>>(thyra_p.getVector());
-
-    Thyra::ModelEvaluatorBase::InArgs<Real> inArgs = thyra_model.createInArgs();
-    Thyra::ModelEvaluatorBase::OutArgs<Real> outArgs = thyra_model.createOutArgs();
-
-    outArgs.set_f(thyra_f.getVector());
-    for(std::size_t i=0; i<p_indices.size(); ++i)
-      inArgs.set_p(p_indices[i], thyra_prodvec_p->getVectorBlock(i));
-    inArgs.set_x(thyra_x.getVector());
-
-    thyra_model.evalModel(inArgs, outArgs);
-
+    if(verbosityLevel >= Teuchos::VERB_MEDIUM)
+      *out << "ROL::ThyraProductME_Constraint_SimOpt::value" << std::endl;
 
-    /*  {
-          const ThyraVector<Real>  & thyra_p = dynamic_cast<const ThyraVector<Real>&>(z);
-          const ThyraVector<Real>  & thyra_x = dynamic_cast<const ThyraVector<Real>&>(u);
-          const ThyraVector<Real>  & thyra_f = dynamic_cast<const ThyraVector<Real>&>(c);
+    if(!computeValue) {
+      if(verbosityLevel >= Teuchos::VERB_HIGH)
+        *out << "ROL::ThyraProductME_Constraint_SimOpt::value, Skipping Value Computation" << std::endl;
 
+      TEUCHOS_ASSERT(Teuchos::nonnull(value_ptr_));
+      c.set(*value_ptr_);
+    }
+    else {
+      const ThyraVector<Real>  & thyra_p = dynamic_cast<const ThyraVector<Real>&>(z);
+      const ThyraVector<Real>  & thyra_x = dynamic_cast<const ThyraVector<Real>&>(u);
+      ThyraVector<Real>  & thyra_f = dynamic_cast<ThyraVector<Real>&>(c);
+      Teuchos::RCP<const Thyra::ProductVectorBase<Real> > thyra_prodvec_p = Teuchos::rcp_dynamic_cast<const Thyra::ProductVectorBase<Real>>(thyra_p.getVector());
 
-                Thyra::ConstDetachedVectorView<Real> x_view(thyra_x.getVector());
-                    Thyra::ConstDetachedVectorView<Real> p_view(thyra_p.getVector());
-                    Thyra::ConstDetachedVectorView<Real> f_view(thyra_f.getVector());
+      Thyra::ModelEvaluatorBase::InArgs<Real> inArgs = thyra_model.createInArgs();
+      Thyra::ModelEvaluatorBase::OutArgs<Real> outArgs = thyra_model.createOutArgs();
 
-                    std::cout << "\nEnd of value... x:" << " ";
-                    for (std::size_t i=0; i<x_view.subDim(); ++i)
-                      std::cout << x_view(i) << " ";
-                    std::cout << "\np:" << " ";
-                    for (std::size_t i=0; i<p_view.subDim(); ++i)
-                      std::cout << p_view(i) << " ";
-                    std::cout << "\nf:" << " ";
-                    for (std::size_t i=0; i<f_view.subDim(); ++i)
-                      std::cout << f_view(i) << " ";
+      outArgs.set_f(thyra_f.getVector());
+      for(std::size_t i=0; i<p_indices.size(); ++i)
+        inArgs.set_p(p_indices[i], thyra_prodvec_p->getVectorBlock(i));
+      inArgs.set_x(thyra_x.getVector());
 
-                    std::cout << "Norm: " << c.norm() <<std::endl;;
-          }*/
+      thyra_model.evalModel(inArgs, outArgs);
 
+      /*
+      {
+        const ThyraVector<Real>  & thyra_p = dynamic_cast<const ThyraVector<Real>&>(z);
+        const ThyraVector<Real>  & thyra_x = dynamic_cast<const ThyraVector<Real>&>(u);
+        const ThyraVector<Real>  & thyra_f = dynamic_cast<const ThyraVector<Real>&>(c);
+
+        Thyra::ConstDetachedVectorView<Real> x_view(thyra_x.getVector());
+        Thyra::ConstDetachedVectorView<Real> p_view(thyra_p.getVector());
+        Thyra::ConstDetachedVectorView<Real> f_view(thyra_f.getVector());
+
+        std::cout << "\nEnd of value... x:" << " ";
+        for (std::size_t i=0; i<x_view.subDim(); ++i)
+          std::cout << x_view(i) << " ";
+        std::cout << "\np:" << " ";
+        for (std::size_t i=0; i<p_view.subDim(); ++i)
+          std::cout << p_view(i) << " ";
+        std::cout << "\nf:" << " ";
+        for (std::size_t i=0; i<f_view.subDim(); ++i)
+          std::cout << f_view(i) << " ";
+
+        std::cout << "Norm: " << c.norm() <<std::endl;
+      }
+      */
 
+      if (Teuchos::is_null(value_ptr_))
+        value_ptr_ = c.clone();
+      value_ptr_->set(c);
 
-    updateValue = false;
+      computeValue = false;
+    }
   }
 
   void applyJacobian_1(Vector<Real> &jv, const Vector<Real> &v, const Vector<Real> &u,
       const Vector<Real> &z, Real &tol) {
 
-    if(updateJacobian1) {
+    if(verbosityLevel >= Teuchos::VERB_MEDIUM)
+      *out << "ROL::ThyraProductME_Constraint_SimOpt::applyJacobian_1" << std::endl;
+
+#ifdef  HAVE_ROL_DEBUG
+    //u and z should be updated in the update functions before calling applyJacobian_1
+    TEUCHOS_ASSERT(!u_hasChanged(u));
+    TEUCHOS_ASSERT(!z_hasChanged(z));
+#endif
+
+    if(computeJacobian1) {
       // Create Jacobian
       const ThyraVector<Real>  & thyra_x = dynamic_cast<const ThyraVector<Real>&>(u);
       const ThyraVector<Real>  & thyra_p = dynamic_cast<const ThyraVector<Real>&>(z);
@@ -150,7 +185,11 @@ class ThyraProductME_Constraint_SimOpt : public Constraint_SimOpt<Real> {
       outArgs.set_W_op(lop);
       thyra_model.evalModel(inArgs, outArgs);
       jac1 = lop;
-      updateJacobian1 = false;
+
+      computeJacobian1 = false;
+    } else {
+      if(verbosityLevel >= Teuchos::VERB_HIGH)
+        *out << "ROL::ThyraProductME_Constraint_SimOpt::applyJacobian_1, Skipping Jacobian Computation" << std::endl;
     }
 
     const ThyraVector<Real>  & thyra_v = dynamic_cast<const ThyraVector<Real>&>(v);
@@ -161,7 +200,15 @@ class ThyraProductME_Constraint_SimOpt : public Constraint_SimOpt<Real> {
 
   void applyJacobian_2(Vector<Real> &jv, const Vector<Real> &v, const Vector<Real> &u,
       const Vector<Real> &z, Real &tol) {
-    //std::cout <<  "Jacobian 2: " << tol <<std::endl;
+
+    if(verbosityLevel >= Teuchos::VERB_MEDIUM)
+      *out << "ROL::ThyraProductME_Constraint_SimOpt::applyJacobian_2" << std::endl;
+
+#ifdef  HAVE_ROL_DEBUG
+    //u and z should be updated in the update functions before calling applyJacobian_1
+    TEUCHOS_ASSERT(!u_hasChanged(u));
+    TEUCHOS_ASSERT(!z_hasChanged(z));
+#endif
 
     const ThyraVector<Real>  & thyra_p = dynamic_cast<const ThyraVector<Real>&>(z);
     const ThyraVector<Real>  & thyra_x = dynamic_cast<const ThyraVector<Real>&>(u);
@@ -181,59 +228,56 @@ class ThyraProductME_Constraint_SimOpt : public Constraint_SimOpt<Real> {
 
       auto p_space = thyra_model.get_p_space(i);
       auto f_space = thyra_model.get_f_space();
-      int num_params = p_space->dim();
-      int num_resids = f_space->dim();
       auto p_space_plus = Teuchos::rcp_dynamic_cast<const Thyra::SpmdVectorSpaceDefaultBase<Real>>(p_space);
       auto f_space_plus = Teuchos::rcp_dynamic_cast<const Thyra::SpmdVectorSpaceDefaultBase<Real>>(f_space);
-      bool p_dist = !p_space_plus->isLocallyReplicated();//p_space->DistributedGlobal();
-      bool f_dist = !f_space_plus->isLocallyReplicated();//f_space->DistributedGlobal();
       Thyra::ModelEvaluatorBase::DerivativeSupport ds =  outArgs.supports(Thyra::ModelEvaluatorBase::OUT_ARG_DfDp,i);
       // Determine which layout to use for df/dp.  Ideally one would look
-      // at num_params, num_resids, what is supported by the underlying
+      // at the parameter and residual dimensions, what is supported by the underlying
       // model evaluator, and the sensitivity method, and make the best
       // choice to minimze the number of solves.  However this choice depends
       // also on what layout of dg/dx is supported (e.g., if only the operator
       // form is supported for forward sensitivities, then df/dp must be
       // DERIV_MV_BY_COL).  For simplicity, we order the conditional tests
       // to get the right layout in most situations.
-      enum DerivativeLayout { OP, COL, ROW } dfdp_layout;
-      { // if (sensitivity_method == "Adjoint")
-        if (ds.supports(Thyra::ModelEvaluatorBase::DERIV_LINEAR_OP))
-          dfdp_layout = OP;
-        else if (ds.supports(Thyra::ModelEvaluatorBase::DERIV_MV_GRADIENT_FORM) && !f_dist)
-          dfdp_layout = ROW;
-        else if (ds.supports(Thyra::ModelEvaluatorBase::DERIV_MV_JACOBIAN_FORM) && !p_dist)
-          dfdp_layout = COL;
+
+      if (ds.supports(Thyra::ModelEvaluatorBase::DERIV_LINEAR_OP)) {
+        auto dfdp_op = thyra_model.create_DfDp_op(i);
+        TEUCHOS_TEST_FOR_EXCEPTION(
+            dfdp_op == Teuchos::null, std::logic_error,
+            std::endl << "ROL::ThyraProductME_Constraint_SimOpt::applyJacobian_2():  " <<
+            "Needed df/dp operator (" << i << ") is null!" << std::endl);
+        outArgs.set_DfDp(i,dfdp_op);
+      } else {
+        TEUCHOS_TEST_FOR_EXCEPTION(!ds.supports(Thyra::ModelEvaluatorBase::DERIV_LINEAR_OP),
+          std::logic_error,
+          std::endl <<
+          "ROL::ThyraProductME_Constraint_SimOpt::applyJacobian_2():  " <<
+          "The code related to df/dp multivector has been commented out because never tested.  " <<
+          std::endl);
+
+          /*
+          if (ds.supports(Thyra::ModelEvaluatorBase::DERIV_MV_GRADIENT_FORM) && f_space_plus->isLocallyReplicated()) {
+          auto dfdp = Thyra::createMembers(p_space, f_space->dim());
+
+          Thyra::ModelEvaluatorBase::DerivativeMultiVector<Real>
+          dmv_dfdp(dfdp, Thyra::ModelEvaluatorBase::DERIV_MV_GRADIENT_FORM);
+          outArgs.set_DfDp(i,dmv_dfdp);
+        } else if (ds.supports(Thyra::ModelEvaluatorBase::DERIV_MV_JACOBIAN_FORM) && p_space_plus->isLocallyReplicated()) {
+          auto dfdp = Thyra::createMembers(f_space, p_space->dim());
+          Thyra::ModelEvaluatorBase::DerivativeMultiVector<Real>
+          dmv_dfdp(dfdp, Thyra::ModelEvaluatorBase::DERIV_MV_JACOBIAN_FORM);
+          outArgs.set_DfDp(i,dmv_dfdp);
+        }
         else
           TEUCHOS_TEST_FOR_EXCEPTION(
               true, std::logic_error,
-              std::endl << "Piro::NOXSolver::evalModel():  " <<
+              std::endl << "ROL::ThyraProductME_Constraint_SimOpt::applyJacobian_2():  " <<
               "For df/dp(" << i <<") with adjoint sensitivities, " <<
               "underlying ModelEvaluator must support DERIV_LINEAR_OP, " <<
               "DERIV_MV_BY_COL with p not distributed, or "
               "DERIV_TRANS_MV_BY_ROW with f not distributed." <<
               std::endl);
-      }
-      if (dfdp_layout == COL) {
-        auto dfdp = Thyra::createMembers(f_space, num_params);
-        Thyra::ModelEvaluatorBase::DerivativeMultiVector<Real>
-        dmv_dfdp(dfdp, Thyra::ModelEvaluatorBase::DERIV_MV_JACOBIAN_FORM);
-        outArgs.set_DfDp(i,dmv_dfdp);
-      }
-      else if (dfdp_layout == ROW) {
-        auto dfdp = Thyra::createMembers(p_space, num_resids);
-
-        Thyra::ModelEvaluatorBase::DerivativeMultiVector<Real>
-        dmv_dfdp(dfdp, Thyra::ModelEvaluatorBase::DERIV_MV_GRADIENT_FORM);
-        outArgs.set_DfDp(i,dmv_dfdp);
-      }
-      else if (dfdp_layout == OP) {
-        auto dfdp_op = thyra_model.create_DfDp_op(i);
-        TEUCHOS_TEST_FOR_EXCEPTION(
-            dfdp_op == Teuchos::null, std::logic_error,
-            std::endl << "Piro::NOXSolver::evalModel():  " <<
-            "Needed df/dp operator (" << i << ") is null!" << std::endl);
-        outArgs.set_DfDp(i,dfdp_op);
+              */
       }
     }
 
@@ -251,6 +295,14 @@ class ThyraProductME_Constraint_SimOpt : public Constraint_SimOpt<Real> {
         dfdp_op->apply(Thyra::NOTRANS,*thyra_prodvec_v->getVectorBlock(i), temp_jv_ptr->getVector().ptr(),1.0, 0.0);
         thyra_jv.axpy(1.0, *temp_jv_ptr);
       } else {
+        TEUCHOS_TEST_FOR_EXCEPTION(
+          dfdp_op == Teuchos::null,
+          std::logic_error,
+          std::endl <<
+          "ROL::ThyraProductME_Constraint_SimOpt::applyJacobian_2():  " <<
+          "The code related to df/dp multivector has been commented out because never tested.  " <<
+          std::endl);
+        /*
         Thyra::ModelEvaluatorBase::EDerivativeMultiVectorOrientation dfdp_orient =
             outArgs.get_DfDp(i).getMultiVectorOrientation();
         Thyra::ModelEvaluatorBase::EDerivativeMultiVectorOrientation dgdx_orient =
@@ -284,16 +336,27 @@ class ThyraProductME_Constraint_SimOpt : public Constraint_SimOpt<Real> {
                 jv_view(ix) += v_view(0,ip)*dfdp_view(ix,ip);
           }
         }
+        */
       }
     }
   }
 
 
-  virtual void applyInverseJacobian_1(Vector<Real> &ijv,
+  void applyInverseJacobian_1(Vector<Real> &ijv,
       const Vector<Real> &v,
       const Vector<Real> &u,
       const Vector<Real> &z,
       Real &tol) {
+
+    if(verbosityLevel >= Teuchos::VERB_MEDIUM)
+      *out << "ROL::ThyraProductME_Constraint_SimOpt::applyInverseJacobian_1" << std::endl;
+
+#ifdef  HAVE_ROL_DEBUG
+    //u and z should be updated in the update functions before calling applyJacobian_2
+    TEUCHOS_ASSERT(!u_hasChanged(u));
+    TEUCHOS_ASSERT(!z_hasChanged(z));
+#endif
+
     const ThyraVector<Real>  & thyra_p = dynamic_cast<const ThyraVector<Real>&>(z);
     const ThyraVector<Real>  & thyra_x = dynamic_cast<const ThyraVector<Real>&>(u);
     const ThyraVector<Real>  & thyra_v = dynamic_cast<const ThyraVector<Real>&>(v);
@@ -314,10 +377,13 @@ class ThyraProductME_Constraint_SimOpt : public Constraint_SimOpt<Real> {
     Teuchos::RCP< const Thyra::LinearOpWithSolveFactoryBase<double> > lows_factory = thyra_model.get_W_factory();
     TEUCHOS_ASSERT(Teuchos::nonnull(lows_factory));
     Teuchos::RCP< Thyra::LinearOpBase<double> > lop;
-    if(updateJacobian1)
+    if(computeJacobian1)
       lop = thyra_model.create_W_op();
-    else
+    else {
+      if(verbosityLevel >= Teuchos::VERB_HIGH)
+        *out << "ROL::ThyraProductME_Constraint_SimOpt::applyInverseJacobian_1, Skipping Jacobian Computation" << std::endl;
       lop = jac1;
+    }
 
 
     Teuchos::RCP< const ::Thyra::DefaultLinearOpSource<double> > losb = Teuchos::rcp(new ::Thyra::DefaultLinearOpSource<double>(lop));
@@ -331,14 +397,15 @@ class ThyraProductME_Constraint_SimOpt : public Constraint_SimOpt<Real> {
     }
     const Teuchos::RCP<Thyra::LinearOpWithSolveBase<Real> > jacobian = lows_factory->createOp();
 
-    if(updateJacobian1)
+    if(computeJacobian1)
     {
       outArgs.set_W_op(lop);
       thyra_model.evalModel(inArgs, outArgs);
       outArgs.set_W_op(Teuchos::null);
       inArgs.set_x(Teuchos::null);
       jac1 = lop;
-      updateJacobian1 = false;
+
+      computeJacobian1 = false;
     }
 
     if (Teuchos::nonnull(prec_factory))
@@ -372,8 +439,18 @@ class ThyraProductME_Constraint_SimOpt : public Constraint_SimOpt<Real> {
   void applyAdjointJacobian_1(Vector<Real> &ajv, const Vector<Real> &v, const Vector<Real> &u,
       const Vector<Real> &z, Real &tol) {
 
+#ifdef  HAVE_ROL_DEBUG
+    //u and z should be updated in the update functions before calling applyAdjointJacobian_1
+    TEUCHOS_ASSERT(!u_hasChanged(u));
+    TEUCHOS_ASSERT(!z_hasChanged(z));
+#endif
+
+    if(verbosityLevel >= Teuchos::VERB_MEDIUM)
+      *out << "ROL::ThyraProductME_Constraint_SimOpt::applyAdjointJacobian_1" << std::endl;
+
+
     Teuchos::RCP< Thyra::LinearOpBase<double> > lop;
-    if(updateJacobian1){
+    if(computeJacobian1){
       const ThyraVector<Real>  & thyra_p = dynamic_cast<const ThyraVector<Real>&>(z);
       const ThyraVector<Real>  & thyra_x = dynamic_cast<const ThyraVector<Real>&>(u);
 
@@ -395,10 +472,15 @@ class ThyraProductME_Constraint_SimOpt : public Constraint_SimOpt<Real> {
       outArgs.set_W_op(lop);
       thyra_model.evalModel(inArgs, outArgs);
       jac1 = lop;
-      updateJacobian1 = false;
+
+      computeJacobian1 = false;
     }
-    else
+    else {
+      if(verbosityLevel >= Teuchos::VERB_HIGH)
+        *out << "ROL::ThyraProductME_Constraint_SimOpt::applyAdjointJacobian_1, Skipping Jacobian Computation" << std::endl;
       lop = jac1;
+    }
+
 
     const ThyraVector<Real>  & thyra_v = dynamic_cast<const ThyraVector<Real>&>(v);
     ThyraVector<Real>  & thyra_ajv = dynamic_cast<ThyraVector<Real>&>(ajv);
@@ -406,12 +488,21 @@ class ThyraProductME_Constraint_SimOpt : public Constraint_SimOpt<Real> {
 
   }
 
-  virtual void applyInverseAdjointJacobian_1(Vector<Real> &iajv,
+  void applyInverseAdjointJacobian_1(Vector<Real> &iajv,
       const Vector<Real> &v,
       const Vector<Real> &u,
       const Vector<Real> &z,
       Real &tol) {
 
+    if(verbosityLevel >= Teuchos::VERB_MEDIUM)
+      *out << "ROL::ThyraProductME_Constraint_SimOpt::applyInverseAdjointJacobian_1" << std::endl;
+
+#ifdef  HAVE_ROL_DEBUG
+    //u and z should be updated in the update functions before calling applyInverseAdjointJacobian_1
+    TEUCHOS_ASSERT(!u_hasChanged(u));
+    TEUCHOS_ASSERT(!z_hasChanged(z));
+#endif
+
     const ThyraVector<Real>  & thyra_p = dynamic_cast<const ThyraVector<Real>&>(z);
     const ThyraVector<Real>  & thyra_x = dynamic_cast<const ThyraVector<Real>&>(u);
     const ThyraVector<Real>  & thyra_v = dynamic_cast<const ThyraVector<Real>&>(v);
@@ -431,10 +522,13 @@ class ThyraProductME_Constraint_SimOpt : public Constraint_SimOpt<Real> {
     TEUCHOS_ASSERT(Teuchos::nonnull(lows_factory));
     Teuchos::RCP< Thyra::LinearOpBase<double> > lop;
 
-    if(updateJacobian1)
+    if(computeJacobian1)
       lop = thyra_model.create_W_op();
-    else
+    else {
+      if(verbosityLevel >= Teuchos::VERB_HIGH)
+        *out << "ROL::ThyraProductME_Constraint_SimOpt::applyInverseAdjointJacobian_1, Skipping Jacobian Computation" << std::endl;
       lop = jac1;
+    }
 
 
     Teuchos::RCP< const ::Thyra::DefaultLinearOpSource<double> > losb = Teuchos::rcp(new ::Thyra::DefaultLinearOpSource<double>(lop));
@@ -448,13 +542,14 @@ class ThyraProductME_Constraint_SimOpt : public Constraint_SimOpt<Real> {
     }
     const Teuchos::RCP<Thyra::LinearOpWithSolveBase<Real> > jacobian = lows_factory->createOp();
 
-    if(updateJacobian1)
+    if(computeJacobian1)
     {
       outArgs.set_W_op(lop);
       thyra_model.evalModel(inArgs, outArgs);
       outArgs.set_W_op(Teuchos::null);
       jac1 = lop;
-      updateJacobian1 = false;
+
+      computeJacobian1 = false;
     }
 
     if (Teuchos::nonnull(prec_factory))
@@ -487,13 +582,21 @@ class ThyraProductME_Constraint_SimOpt : public Constraint_SimOpt<Real> {
 
   void applyAdjointJacobian_2(Vector<Real> &ajv, const Vector<Real> &v, const Vector<Real> &u,
       const Vector<Real> &z, Real &tol) {
-    // std::cout << "Adjoint Jacobian 2" <<std::endl;
+
+    if(verbosityLevel >= Teuchos::VERB_MEDIUM)
+      *out << "ROL::ThyraProductME_Constraint_SimOpt::applyAdjointJacobian_2" << std::endl;
+
+#ifdef  HAVE_ROL_DEBUG
+    //u and z should be updated in the update functions before calling applyAdjointJacobian_2
+    TEUCHOS_ASSERT(!u_hasChanged(u));
+    TEUCHOS_ASSERT(!z_hasChanged(z));
+#endif
 
     const ThyraVector<Real>  & thyra_p = dynamic_cast<const ThyraVector<Real>&>(z);
     const ThyraVector<Real>  & thyra_x = dynamic_cast<const ThyraVector<Real>&>(u);
     const ThyraVector<Real>  & thyra_v = dynamic_cast<const ThyraVector<Real>&>(v);
+
     ThyraVector<Real>  & thyra_ajv = dynamic_cast<ThyraVector<Real>&>(ajv);
-    //Teuchos::RCP< Thyra::VectorBase<Real> > thyra_f = Thyra::createMember<Real>(thyra_model.get_f_space());
     Teuchos::RCP<const Thyra::ProductVectorBase<Real> > thyra_prodvec_p = Teuchos::rcp_dynamic_cast<const Thyra::ProductVectorBase<Real>>(thyra_p.getVector());
     Teuchos::RCP<Thyra::ProductVectorBase<Real> > thyra_prodvec_ajv = Teuchos::rcp_dynamic_cast<Thyra::ProductVectorBase<Real>>(thyra_ajv.getVector());
 
@@ -508,59 +611,50 @@ class ThyraProductME_Constraint_SimOpt : public Constraint_SimOpt<Real> {
 
       auto p_space = thyra_model.get_p_space(i);
       auto f_space = thyra_model.get_f_space();
-      int num_params = p_space->dim();
-      int num_resids = f_space->dim();
       auto p_space_plus = Teuchos::rcp_dynamic_cast<const Thyra::SpmdVectorSpaceDefaultBase<Real>>(p_space);
       auto f_space_plus = Teuchos::rcp_dynamic_cast<const Thyra::SpmdVectorSpaceDefaultBase<Real>>(f_space);
-      bool p_dist = !p_space_plus->isLocallyReplicated();//p_space->DistributedGlobal();
-      bool f_dist = !f_space_plus->isLocallyReplicated();//f_space->DistributedGlobal();
       Thyra::ModelEvaluatorBase::DerivativeSupport ds =  outArgs.supports(Thyra::ModelEvaluatorBase::OUT_ARG_DfDp,i);
-      // Determine which layout to use for df/dp.  Ideally one would look
-      // at num_params, num_resids, what is supported by the underlying
-      // model evaluator, and the sensitivity method, and make the best
-      // choice to minimze the number of solves.  However this choice depends
-      // also on what layout of dg/dx is supported (e.g., if only the operator
-      // form is supported for forward sensitivities, then df/dp must be
-      // DERIV_MV_BY_COL).  For simplicity, we order the conditional tests
-      // to get the right layout in most situations.
-      enum DerivativeLayout { OP, COL, ROW } dfdp_layout;
-      { // if (sensitivity_method == "Adjoint")
-        if (ds.supports(Thyra::ModelEvaluatorBase::DERIV_LINEAR_OP))
-          dfdp_layout = OP;
-        else if (ds.supports(Thyra::ModelEvaluatorBase::DERIV_MV_GRADIENT_FORM) && !f_dist)
-          dfdp_layout = ROW;
-        else if (ds.supports(Thyra::ModelEvaluatorBase::DERIV_MV_JACOBIAN_FORM) && !p_dist)
-          dfdp_layout = COL;
+      // Determine which layout to use for df/dp.
+
+      if (ds.supports(Thyra::ModelEvaluatorBase::DERIV_LINEAR_OP)) {
+        auto dfdp_op = thyra_model.create_DfDp_op(i);
+        TEUCHOS_TEST_FOR_EXCEPTION(
+            dfdp_op == Teuchos::null, std::logic_error,
+            std::endl << "ROL::ThyraProductME_Constraint_SimOpt::applyAdjointJacobian_2:  " <<
+            "Needed df/dp operator (" << i << ") is null!" << std::endl);
+        outArgs.set_DfDp(i,dfdp_op);
+      } else {
+        TEUCHOS_TEST_FOR_EXCEPTION(
+          !ds.supports(Thyra::ModelEvaluatorBase::DERIV_LINEAR_OP),
+          std::logic_error,
+          std::endl <<
+          "ROL::ThyraProductME_Constraint_SimOpt::applyAdjointJacobian_2():  " <<
+          "The code related to df/dp multivector has been commented out because never tested.  " <<
+          std::endl);
+
+          /*
+          if (ds.supports(Thyra::ModelEvaluatorBase::DERIV_MV_GRADIENT_FORM) && f_space_plus->isLocallyReplicated()) {
+          auto dfdp = Thyra::createMembers(p_space, f_space->dim());
+
+          Thyra::ModelEvaluatorBase::DerivativeMultiVector<Real>
+          dmv_dfdp(dfdp, Thyra::ModelEvaluatorBase::DERIV_MV_GRADIENT_FORM);
+          outArgs.set_DfDp(i,dmv_dfdp);
+        } else if (ds.supports(Thyra::ModelEvaluatorBase::DERIV_MV_JACOBIAN_FORM) && p_space_plus->isLocallyReplicated()) {
+          auto dfdp = Thyra::createMembers(f_space, p_space->dim());
+          Thyra::ModelEvaluatorBase::DerivativeMultiVector<Real>
+          dmv_dfdp(dfdp, Thyra::ModelEvaluatorBase::DERIV_MV_JACOBIAN_FORM);
+          outArgs.set_DfDp(i,dmv_dfdp);
+        }
         else
           TEUCHOS_TEST_FOR_EXCEPTION(
               true, std::logic_error,
-              std::endl << "Piro::NOXSolver::evalModel():  " <<
+              std::endl << "ROL::ThyraProductME_Constraint_SimOpt::applyAdjointJacobian_2():  " <<
               "For df/dp(" << i <<") with adjoint sensitivities, " <<
               "underlying ModelEvaluator must support DERIV_LINEAR_OP, " <<
               "DERIV_MV_BY_COL with p not distributed, or "
               "DERIV_TRANS_MV_BY_ROW with f not distributed." <<
               std::endl);
-      }
-      if (dfdp_layout == COL) {
-        auto dfdp = Thyra::createMembers(f_space, num_params);
-        Thyra::ModelEvaluatorBase::DerivativeMultiVector<Real>
-        dmv_dfdp(dfdp, Thyra::ModelEvaluatorBase::DERIV_MV_JACOBIAN_FORM);
-        outArgs.set_DfDp(i,dmv_dfdp);
-      }
-      else if (dfdp_layout == ROW) {
-        auto dfdp = Thyra::createMembers(p_space, num_resids);
-
-        Thyra::ModelEvaluatorBase::DerivativeMultiVector<Real>
-        dmv_dfdp(dfdp, Thyra::ModelEvaluatorBase::DERIV_MV_GRADIENT_FORM);
-        outArgs.set_DfDp(i,dmv_dfdp);
-      }
-      else if (dfdp_layout == OP) {
-        auto dfdp_op = thyra_model.create_DfDp_op(i);
-        TEUCHOS_TEST_FOR_EXCEPTION(
-            dfdp_op == Teuchos::null, std::logic_error,
-            std::endl << "Piro::NOXSolver::evalModel():  " <<
-            "Needed df/dp operator (" << i << ") is null!" << std::endl);
-        outArgs.set_DfDp(i,dfdp_op);
+              */
       }
     }
 
@@ -576,6 +670,14 @@ class ThyraProductME_Constraint_SimOpt : public Constraint_SimOpt<Real> {
         dfdp_op->apply(Thyra::TRANS,*thyra_v.getVector(), thyra_prodvec_ajv->getNonconstVectorBlock(i).ptr(),1.0, 0.0);
         // Thyra::update(1.0,  *tmp, thyra_ajv.getMultiVector().ptr());
       } else {
+        TEUCHOS_TEST_FOR_EXCEPTION(
+          dfdp_op == Teuchos::null,
+          std::logic_error,
+          std::endl <<
+          "ROL::ThyraProductME_Constraint_SimOpt::applyAdjointJacobian_2():  " <<
+          "The code related to df/dp multivector has been commented out because never tested.  " <<
+          std::endl);
+        /*
         Thyra::ModelEvaluatorBase::EDerivativeMultiVectorOrientation dfdp_orient =
             outArgs.get_DfDp(i).getMultiVectorOrientation();
         Thyra::ModelEvaluatorBase::EDerivativeMultiVectorOrientation dgdx_orient =
@@ -608,6 +710,7 @@ class ThyraProductME_Constraint_SimOpt : public Constraint_SimOpt<Real> {
                 ajv_view(ip) += v_view(0,ix)*dfdp_view(ix,ip);
           }
         }
+        */
       }
     }
   }
@@ -617,6 +720,17 @@ class ThyraProductME_Constraint_SimOpt : public Constraint_SimOpt<Real> {
       const Vector<Real> &z,
       Real &tol) {
 
+    if(verbosityLevel >= Teuchos::VERB_MEDIUM)
+      *out << "ROL::ThyraProductME_Constraint_SimOpt::solve" << std::endl;
+
+
+    if(!solveConstraint) {
+      TEUCHOS_ASSERT(Teuchos::nonnull(rol_u_ptr));
+      u.set(*rol_u_ptr);
+      value(c, u, z, tol);
+      return;
+    }
+
     if(thyra_solver.is_null())
       Constraint_SimOpt<Real>::solve(c,u,z,tol);
     else {
@@ -654,43 +768,91 @@ class ThyraProductME_Constraint_SimOpt : public Constraint_SimOpt<Real> {
       this->update_1(u);
     }
 
-    updateValue = false;
+    if (Teuchos::is_null(value_ptr_))
+      value_ptr_ = c.clone();
+    value_ptr_->set(c);
+
+    computeValue = solveConstraint = false;
   }
 
 
   /** \brief Update constraint functions with respect to Sim variable.
                 x is the optimization variable,
-                flag = true if optimization variable is changed,
+                flag = ??,
                 iter is the outer algorithm iterations count.
    */
-  void update_1( const Vector<Real> &u, bool flag = true, int iter = -1 ) {
-    if (flag == true) {
-      updateValue = true;
-      updateJacobian1 = true;
+  void update_1( const Vector<Real> &u, bool /*flag*/ = true, int iter = -1 ) {
+    if(u_hasChanged(u)) {
+      if(verbosityLevel >= Teuchos::VERB_HIGH)
+        *out << "ROL::ThyraProductME_Constraint_SimOpt::update_1, The State Changed" << std::endl;
+      computeValue = computeJacobian1 = true;
+
+      if (Teuchos::is_null(rol_u_ptr))
+        rol_u_ptr = u.clone();
+      rol_u_ptr->set(u);
     }
+
     if(params != Teuchos::null)
       params->set<int>("Optimizer Iteration Number", iter);
   }
 
   /** \brief Update constraint functions with respect to Opt variable.
                 x is the optimization variable,
-                flag = true if optimization variable is changed,
+                flag = ??,
                 iter is the outer algorithm iterations count.
    */
-  void update_2( const Vector<Real> &z, bool flag = true, int iter = -1 ) {
-    if (flag == true) {
-      updateValue = true;
-      updateJacobian1 = true;
+  void update_2( const Vector<Real> &z, bool /*flag*/ = true, int iter = -1 ) {
+    if(z_hasChanged(z)) {
+      if(verbosityLevel >= Teuchos::VERB_HIGH)
+        *out << "ROL::ThyraProductME_Constraint_SimOpt::update_2, The Parameter Changed" << std::endl;
+      computeValue = computeJacobian1 = solveConstraint = true;
+
+      if (Teuchos::is_null(rol_z_ptr))
+        rol_z_ptr = z.clone();
+      rol_z_ptr->set(z);
     }
+
     if(Teuchos::nonnull(params)) {
+      auto& z_stored_ptr = params->get<Teuchos::RCP<Vector<Real> > >("Optimization Variable");
+      if(Teuchos::is_null(z_stored_ptr) || z_hasChanged(*z_stored_ptr)) {
+        if(verbosityLevel >= Teuchos::VERB_HIGH)
+          *out << "ROL::ThyraProductME_Constraint_SimOpt::update_2, Signaling That Parameter Changed" << std::endl;
+        params->set<bool>("Optimization Variables Changed", true);
+        if(Teuchos::is_null(z_stored_ptr))
+          z_stored_ptr = z.clone();
+        z_stored_ptr->set(z);
+      }
+
       params->set<int>("Optimizer Iteration Number", iter);
-      if(flag == true)
-        params->set<bool>("Optimization Variables Changed",true);
     }
   }
 
+  bool z_hasChanged(const Vector<Real> &rol_z) const {
+    bool changed = true;
+    if (Teuchos::nonnull(rol_z_ptr)) {
+      auto diff = rol_z.clone();
+      diff->set(*rol_z_ptr);
+      diff->axpy( -1.0, rol_z );
+      Real norm = diff->norm();
+      changed = (norm == 0) ? false : true;
+    }
+    return changed;
+  }
+
+  bool u_hasChanged(const Vector<Real> &rol_u) const {
+    bool changed = true;
+    if (Teuchos::nonnull(rol_u_ptr)) {
+      auto diff = rol_u.clone();
+      diff->set(*rol_u_ptr);
+      diff->axpy( -1.0, rol_u );
+      Real norm = diff->norm();
+      changed = (norm == 0) ? false : true;
+    }
+    return changed;
+  }
+
 public:
-  bool updateValue, updateJacobian1;
+  bool computeValue, computeJacobian1, solveConstraint;
 
 private:
   Teuchos::RCP<Thyra::ModelEvaluatorDefaultBase<double>> thyra_solver;
@@ -698,10 +860,12 @@ class ThyraProductME_Constraint_SimOpt : public Constraint_SimOpt<Real> {
   const int g_index;
   const std::vector<int> p_indices;
   int num_responses;
-  Real value_;
-  Teuchos::RCP<Vector<Real> > x_ptr, grad_ptr;
+  Teuchos::RCP<Vector<Real> > value_ptr_;
+  Teuchos::RCP<Vector<Real> > rol_u_ptr, rol_z_ptr;
   Teuchos::RCP<Teuchos::ParameterList> params;
+  Teuchos::RCP<Teuchos::FancyOStream> out;
   Teuchos::RCP< Thyra::LinearOpBase<double> > jac1;
+  Teuchos::EVerbosityLevel verbosityLevel;
 
 };
 
diff --git a/packages/rol/adapters/thyra/src/function/ROL_ThyraProductME_Objective.hpp b/packages/rol/adapters/thyra/src/function/ROL_ThyraProductME_Objective.hpp
index a5745b8a8b2a..0bb75042d7fd 100644
--- a/packages/rol/adapters/thyra/src/function/ROL_ThyraProductME_Objective.hpp
+++ b/packages/rol/adapters/thyra/src/function/ROL_ThyraProductME_Objective.hpp
@@ -50,7 +50,7 @@
 //#include "Thyra_DefaultProductVectorSpace.hpp"
 #include "Thyra_ProductVectorBase.hpp"
 #include <vector>
-#include <iostream>
+#include "Teuchos_VerbosityLevel.hpp"
 
 /** \class ROL::ThyraProductME_Objective
     \brief Implements the ROL::Objective interface for a Thyra Model Evaluator Objective.
@@ -62,11 +62,14 @@ template <class Real>
 class ThyraProductME_Objective : public Objective<Real> {
 public:
 
-  ThyraProductME_Objective(Thyra::ModelEvaluatorDefaultBase<double>& thyra_model_, int g_index_, const std::vector<int>& p_indices_,Teuchos::RCP<Teuchos::ParameterList> params_ = Teuchos::null) :
-    thyra_model(thyra_model_), g_index(g_index_), p_indices(p_indices_), params(params_) {
-    computeValue = true;
+  ThyraProductME_Objective(Thyra::ModelEvaluatorDefaultBase<double>& thyra_model_, int g_index_, const std::vector<int>& p_indices_,
+      Teuchos::RCP<Teuchos::ParameterList> params_ = Teuchos::null, Teuchos::EVerbosityLevel verbLevel= Teuchos::VERB_HIGH) :
+    thyra_model(thyra_model_), g_index(g_index_), p_indices(p_indices_), params(params_),
+    out(Teuchos::VerboseObjectBase::getDefaultOStream()),
+    verbosityLevel(verbLevel) {
+    computeValue = computeGradient = true;
     value_ = 0;
-    if(params != Teuchos::null) {
+    if(Teuchos::nonnull(params)) {
       params->set<int>("Optimizer Iteration Number", -1);
       params->set<bool>("Compute State", true);
     }
@@ -80,8 +83,19 @@ class ThyraProductME_Objective : public Objective<Real> {
   */
   Real value( const Vector<Real> &rol_x, Real &tol ) {
 
-    if(!computeValue)
-       return value_;
+#ifdef  HAVE_ROL_DEBUG
+    //x should be updated in the update function before calling value
+    TEUCHOS_ASSERT(!x_hasChanged(rol_x));
+#endif
+
+    if(verbosityLevel >= Teuchos::VERB_MEDIUM)
+      *out << "ROL::ThyraProductME_Objective::value" << std::endl;
+
+    if(!computeValue) {
+      if(verbosityLevel >= Teuchos::VERB_HIGH)
+      *out << "ROL::ThyraProductME_Objective::value, Skipping Value Computation" << std::endl;
+      return value_;
+    }
 
     // Real norm = rol_x.norm();
     // std::cout << "Value norm: " << norm << std::endl;
@@ -119,6 +133,20 @@ class ThyraProductME_Objective : public Objective<Real> {
   */
   void gradient( Vector<Real> &rol_g, const Vector<Real> &rol_x, Real &tol ) {
 
+#ifdef  HAVE_ROL_DEBUG
+    //x should be updated in the update function before calling value
+    TEUCHOS_ASSERT(!x_hasChanged(rol_x));
+#endif
+
+    if(verbosityLevel >= Teuchos::VERB_MEDIUM)
+      *out << "ROL::ThyraProductME_Objective::gradient" << std::endl;
+
+    if(!computeGradient) {
+      if(verbosityLevel >= Teuchos::VERB_HIGH)
+        *out << "ROL::ThyraProductME_Objective::gradient, Skipping Gradient Computation" << std::endl;
+      return rol_g.set(*grad_ptr_);
+    }
+
     // Real norm = rol_x.norm();
     // std::cout << "In Gradient, Value norm: " << norm << std::endl;
 
@@ -139,6 +167,8 @@ class ThyraProductME_Objective : public Objective<Real> {
     Teuchos::RCP< Thyra::VectorBase<Real> > g;
 
     if(computeValue) {
+      if(verbosityLevel >= Teuchos::VERB_HIGH)
+        *out << "ROL::ThyraProductME_Objective::gradient, Computing Value" << std::endl;
       g = Thyra::createMember<Real>(thyra_model.get_g_space(g_index));
       outArgs.set_g(g_index, g);
     }
@@ -164,28 +194,58 @@ class ThyraProductME_Objective : public Objective<Real> {
       value_ = ::Thyra::get_ele(*g,0);
       computeValue = false;
     }
+    
+    if (grad_ptr_ == Teuchos::null)
+      grad_ptr_ = rol_g.clone();
+    grad_ptr_->set(rol_g);
+    
+    computeGradient = false;
   };
 
-  void update( const Vector<Real> & /*x*/, bool flag = true, int iter = -1 ) {
-    computeValue = flag;
+  void update( const Vector<Real> & x, bool flag = true, int iter = -1 ) {
     if(Teuchos::nonnull(params)) {
       params->set<int>("Optimizer Iteration Number", iter);
-      if(flag == true) {
+    }
+
+    if(x_hasChanged(x)) {
+
+      if(verbosityLevel >= Teuchos::VERB_HIGH)
+        *out << "ROL::ThyraProductME_Objective::update, The Parameter Changed" << std::endl;
+      computeValue = computeGradient = true;
+
+      if(Teuchos::nonnull(params)) {
         params->set<bool>("Compute State", true);
         params->set<bool>("Optimization Variables Changed", true);
       }
     }
   }
+  
+  bool x_hasChanged(const Vector<Real> &rol_x) {
+    bool changed = true;
+    if (Teuchos::nonnull(rol_x_ptr)) {
+      rol_x_ptr->axpy( -1.0, rol_x );
+      Real norm = rol_x_ptr->norm();
+      changed = (norm == 0) ? false : true;
+    } else {
+      rol_x_ptr = rol_x.clone();
+    }
+    rol_x_ptr->set(rol_x);
+    return changed;
+  }  
 
 public:
-  bool computeValue;
+  bool computeValue, computeGradient;
 
 private:
   Thyra::ModelEvaluatorDefaultBase<Real>& thyra_model;
   const int g_index;
   const std::vector<int> p_indices;
   Real value_;
+  Teuchos::RCP<Vector<Real> > grad_ptr_;
+  Teuchos::RCP<Vector<Real> > rol_x_ptr;
   Teuchos::RCP<Teuchos::ParameterList> params;
+  Teuchos::RCP<Teuchos::FancyOStream> out;
+  Teuchos::EVerbosityLevel verbosityLevel;
 
 }; // class Objective
 
diff --git a/packages/rol/adapters/thyra/src/function/ROL_ThyraProductME_Objective_SimOpt.hpp b/packages/rol/adapters/thyra/src/function/ROL_ThyraProductME_Objective_SimOpt.hpp
index 3faec1f2b38d..fae187499e71 100644
--- a/packages/rol/adapters/thyra/src/function/ROL_ThyraProductME_Objective_SimOpt.hpp
+++ b/packages/rol/adapters/thyra/src/function/ROL_ThyraProductME_Objective_SimOpt.hpp
@@ -48,6 +48,7 @@
 #include "ROL_StdVector.hpp"
 #include "ROL_Objective_SimOpt.hpp"
 #include "ROL_Types.hpp"
+#include "Teuchos_VerbosityLevel.hpp"
 
 using namespace ROL;
 
@@ -57,44 +58,77 @@ class ThyraProductME_Objective_SimOpt : public Objective_SimOpt<Real> {
 public:
 
 
-  ThyraProductME_Objective_SimOpt(Thyra::ModelEvaluatorDefaultBase<double>& thyra_model_, int g_index_, const std::vector<int>& p_indices_,Teuchos::RCP<Teuchos::ParameterList> params_ = Teuchos::null) :
-    thyra_model(thyra_model_), g_index(g_index_), p_indices(p_indices_), params(params_) {
-    updateValue = true;
+  ThyraProductME_Objective_SimOpt(Thyra::ModelEvaluatorDefaultBase<double>& thyra_model_, int g_index_, const std::vector<int>& p_indices_,
+      Teuchos::RCP<Teuchos::ParameterList> params_ = Teuchos::null, Teuchos::EVerbosityLevel verbLevel= Teuchos::VERB_HIGH) :
+    thyra_model(thyra_model_), g_index(g_index_), p_indices(p_indices_), params(params_),
+    out(Teuchos::VerboseObjectBase::getDefaultOStream()),
+    verbosityLevel(verbLevel) {
+    computeValue = computeGradient1 = computeGradient2 = true;
     value_ = 0;
-    x_ptr = Teuchos::null;
+    rol_u_ptr = rol_z_ptr = Teuchos::null;
     if(params != Teuchos::null) {
       params->set<int>("Optimizer Iteration Number", -1);
+      params->set<Teuchos::RCP<Vector<Real> > >("Optimization Variable", Teuchos::null);
     }
   };
 
 
   Real value(const Vector<Real> &u, const Vector<Real> &z, Real &tol ) {
 
-    if(updateValue) {
-      const ThyraVector<Real>  & thyra_p = dynamic_cast<const ThyraVector<Real>&>(z);
-      const ThyraVector<Real>  & thyra_x = dynamic_cast<const ThyraVector<Real>&>(u);
-      Teuchos::RCP< Thyra::VectorBase<Real> > g = Thyra::createMember<Real>(thyra_model.get_g_space(g_index));
-      Teuchos::RCP<const Thyra::ProductVectorBase<Real> > thyra_prodvec_p = Teuchos::rcp_dynamic_cast<const Thyra::ProductVectorBase<Real>>(thyra_p.getVector());
+#ifdef  HAVE_ROL_DEBUG
+    //u and z should be updated in the update functions before calling applyAdjointJacobian_2
+    TEUCHOS_ASSERT(!u_hasChanged(u));
+    TEUCHOS_ASSERT(!z_hasChanged(z));
+#endif
 
-      Thyra::ModelEvaluatorBase::InArgs<Real> inArgs = thyra_model.createInArgs();
-      Thyra::ModelEvaluatorBase::OutArgs<Real> outArgs = thyra_model.createOutArgs();
+    if(verbosityLevel >= Teuchos::VERB_MEDIUM)
+      *out << "ROL::ThyraProductME_Objective_SimOpt::value" << std::endl;
 
-      outArgs.set_g(g_index, g);
-      for(std::size_t i=0; i<p_indices.size(); ++i)
-        inArgs.set_p(p_indices[i], thyra_prodvec_p->getVectorBlock(i));
-      inArgs.set_x(thyra_x.getVector());
+    if(!computeValue) {
+      if(verbosityLevel >= Teuchos::VERB_HIGH)
+        *out << "ROL::ThyraProductME_Objective_SimOpt::value, Skipping Computation of Value" << std::endl;
+      return value_;
+    }
 
-      thyra_model.evalModel(inArgs, outArgs);
+    const ThyraVector<Real>  & thyra_p = dynamic_cast<const ThyraVector<Real>&>(z);
+    const ThyraVector<Real>  & thyra_x = dynamic_cast<const ThyraVector<Real>&>(u);
+    Teuchos::RCP< Thyra::VectorBase<Real> > g = Thyra::createMember<Real>(thyra_model.get_g_space(g_index));
+    Teuchos::RCP<const Thyra::ProductVectorBase<Real> > thyra_prodvec_p = Teuchos::rcp_dynamic_cast<const Thyra::ProductVectorBase<Real>>(thyra_p.getVector());
 
-      value_ = ::Thyra::get_ele(*g,0);
+    Thyra::ModelEvaluatorBase::InArgs<Real> inArgs = thyra_model.createInArgs();
+    Thyra::ModelEvaluatorBase::OutArgs<Real> outArgs = thyra_model.createOutArgs();
+
+    outArgs.set_g(g_index, g);
+    for(std::size_t i=0; i<p_indices.size(); ++i)
+      inArgs.set_p(p_indices[i], thyra_prodvec_p->getVectorBlock(i));
+    inArgs.set_x(thyra_x.getVector());
+
+    thyra_model.evalModel(inArgs, outArgs);
+
+    value_ = ::Thyra::get_ele(*g,0);
+
+    computeValue = false;
 
-      updateValue = false;
-    }
     return value_;
   }
 
   void gradient_1(Vector<Real> &g, const Vector<Real> &u, const Vector<Real> &z, Real &tol ) {
 
+#ifdef  HAVE_ROL_DEBUG
+    //u and z should be updated in the update functions before calling gradient_1
+    TEUCHOS_ASSERT(!u_hasChanged(u));
+    TEUCHOS_ASSERT(!z_hasChanged(z));
+#endif
+
+    if(verbosityLevel >= Teuchos::VERB_MEDIUM)
+      *out << "ROL::ThyraProductME_Objective_SimOpt::gradient_1" << std::endl;
+
+    if(!computeGradient1) {
+      if(verbosityLevel >= Teuchos::VERB_HIGH)
+        *out << "ROL::ThyraProductME_Objective_SimOpt::gradient_1, Skipping Computation of Gradient 1" << std::endl;
+      return g.set(*grad1_ptr_);
+    }
+
     const ThyraVector<Real>  & thyra_p = dynamic_cast<const ThyraVector<Real>&>(z);
     const ThyraVector<Real>  & thyra_x = dynamic_cast<const ThyraVector<Real>&>(u);
 
@@ -111,7 +145,7 @@ class ThyraProductME_Objective_SimOpt : public Objective_SimOpt<Real> {
 
     Teuchos::RCP< Thyra::VectorBase<Real> > thyra_g;
 
-    if(updateValue) {
+    if(computeValue) {
       thyra_g = Thyra::createMember<Real>(thyra_model.get_g_space(g_index));
       outArgs.set_g(g_index, thyra_g);
     }
@@ -133,14 +167,38 @@ class ThyraProductME_Objective_SimOpt : public Objective_SimOpt<Real> {
     }
     thyra_model.evalModel(inArgs, outArgs);
 
-    if(updateValue) {
+    if(computeValue) {
+      if(verbosityLevel >= Teuchos::VERB_HIGH)
+      *out << "ROL::ThyraProductME_Objective_SimOpt::gradient_1, Computing Value" << std::endl;
       value_ = ::Thyra::get_ele(*thyra_g,0);
-      updateValue = false;
+      computeValue = false;
     }
+
+    if (Teuchos::is_null(grad1_ptr_))
+      grad1_ptr_ = g.clone();
+    grad1_ptr_->set(g);
+
+    computeGradient1 = false;
   }
 
   void gradient_2(Vector<Real> &g, const Vector<Real> &u, const Vector<Real> &z, Real &tol ) {
 
+#ifdef  HAVE_ROL_DEBUG
+    //u and z should be updated in the update functions before calling gradient_2
+    TEUCHOS_ASSERT(!u_hasChanged(u));
+    TEUCHOS_ASSERT(!z_hasChanged(z));
+#endif
+
+    if(verbosityLevel >= Teuchos::VERB_MEDIUM)
+      *out << "ROL::ThyraProductME_Objective_SimOpt::gradient_2" << std::endl;
+
+
+    if(!computeGradient2) {
+      if(verbosityLevel >= Teuchos::VERB_HIGH)
+        *out << "ROL::ThyraProductME_Objective_SimOpt::gradient_2, Skipping Computation of Gradient 2" << std::endl;
+      return g.set(*grad2_ptr_);
+    }
+
     const ThyraVector<Real>  & thyra_p = dynamic_cast<const ThyraVector<Real>&>(z);
     const ThyraVector<Real>  & thyra_x = dynamic_cast<const ThyraVector<Real>&>(u);
 
@@ -160,7 +218,9 @@ class ThyraProductME_Objective_SimOpt : public Objective_SimOpt<Real> {
 
     Teuchos::RCP< Thyra::VectorBase<Real> > thyra_g;
 
-    if(updateValue) {
+    if(computeValue) {
+      if(verbosityLevel >= Teuchos::VERB_HIGH)
+        *out << "ROL::ThyraProductME_Objective_SimOpt::gradient_2, Computing Value" << std::endl;
       thyra_g = Thyra::createMember<Real>(thyra_model.get_g_space(g_index));
       outArgs.set_g(g_index, thyra_g);
     }
@@ -182,31 +242,86 @@ class ThyraProductME_Objective_SimOpt : public Objective_SimOpt<Real> {
     }
     thyra_model.evalModel(inArgs, outArgs);
 
-    if(updateValue) {
+    if(computeValue) {
       value_ = ::Thyra::get_ele(*thyra_g,0);
-      updateValue = false;
+      computeValue = false;
     }
+
+    if (grad2_ptr_ == Teuchos::null)
+      grad2_ptr_ = g.clone();
+    grad2_ptr_->set(g);
+
+    computeGradient2 = false;
   }
 
-  void update( const Vector<Real> &/*u*/, const Vector<Real> &/*z*/, bool flag = true, int iter = -1) {
-    updateValue = flag;
+  void update( const Vector<Real> &u, const Vector<Real> &z, bool /*flag*/ = true, int iter = -1) {
+    if(z_hasChanged(z) || u_hasChanged(u)) {
+      if(verbosityLevel >= Teuchos::VERB_HIGH)
+        *out << "ROL::ThyraProductME_Objective_SimOpt::update, Either The State Or The Parameters Changed" << std::endl;
+      computeValue = computeGradient1 = computeGradient2 = true;
+
+      if (Teuchos::is_null(rol_z_ptr))
+        rol_z_ptr = z.clone();
+      rol_z_ptr->set(z);
+
+      if (Teuchos::is_null(rol_u_ptr))
+        rol_u_ptr = u.clone();
+      rol_u_ptr->set(u);
+    }
+
     if(params != Teuchos::null) {
+      auto& z_stored_ptr = params->get<Teuchos::RCP<Vector<Real> > >("Optimization Variable");
+      if(Teuchos::is_null(z_stored_ptr) || z_hasChanged(*z_stored_ptr)) {
+        if(verbosityLevel >= Teuchos::VERB_HIGH)
+          *out << "ROL::ThyraProductME_Objective_SimOpt::update, Signaling That Parameter Changed" << std::endl;
+        params->set<bool>("Optimization Variables Changed", true);
+        if(Teuchos::is_null(z_stored_ptr))
+          z_stored_ptr = z.clone();
+        z_stored_ptr->set(z);
+      }
       params->set<int>("Optimizer Iteration Number", iter);
-      if(flag == true)
-        params->set<bool>("Optimization Variables Changed",true);
     }
   }
 
+  bool z_hasChanged(const Vector<Real> &rol_z) const {
+    bool changed = true;
+    if (Teuchos::nonnull(rol_z_ptr)) {
+      auto diff = rol_z.clone();
+      diff->set(*rol_z_ptr);
+      diff->axpy( -1.0, rol_z );
+      Real norm = diff->norm();
+      changed = (norm == 0) ? false : true;
+    }
+    return changed;
+  }
+
+  bool u_hasChanged(const Vector<Real> &rol_u) const {
+    bool changed = true;
+    if (Teuchos::nonnull(rol_u_ptr)) {
+      auto diff = rol_u.clone();
+      diff->set(*rol_u_ptr);
+      diff->axpy( -1.0, rol_u );
+      Real norm = diff->norm();
+      changed = (norm == 0) ? false : true;
+    }
+    return changed;
+  }
+
 public:
-  bool updateValue;
+  bool computeValue, computeGradient1, computeGradient2;
 
 private:
   Thyra::ModelEvaluatorDefaultBase<Real>& thyra_model;
   const int g_index;
   const std::vector<int> p_indices;
   Real value_;
-  Teuchos::RCP<Vector<Real> > x_ptr;
+  Teuchos::RCP<Vector<Real> > grad1_ptr_;
+  Teuchos::RCP<Vector<Real> > grad2_ptr_;
+  Teuchos::RCP<Vector<Real> > rol_z_ptr;
+  Teuchos::RCP<Vector<Real> > rol_u_ptr;
   Teuchos::RCP<Teuchos::ParameterList> params;
+  Teuchos::RCP<Teuchos::FancyOStream> out;
+  Teuchos::EVerbosityLevel verbosityLevel;
 
 };
 
diff --git a/packages/rol/src/function/constraint/ROL_ConstraintDef.hpp b/packages/rol/src/function/constraint/ROL_ConstraintDef.hpp
index 9e9645f22257..22e75d8a810b 100644
--- a/packages/rol/src/function/constraint/ROL_ConstraintDef.hpp
+++ b/packages/rol/src/function/constraint/ROL_ConstraintDef.hpp
@@ -592,6 +592,7 @@ Real Constraint<Real>::checkAdjointConsistencyJacobian(const Vector<Real> &w,
   ROL::Ptr<Vector<Real> > Jv = dualw.clone();
   ROL::Ptr<Vector<Real> > Jw = dualv.clone();
   
+  this->update(x);
   applyJacobian(*Jv,v,x,tol);
   applyAdjointJacobian(*Jw,w,x,tol);
 
diff --git a/packages/rol/src/step/ROL_AugmentedLagrangianStep.hpp b/packages/rol/src/step/ROL_AugmentedLagrangianStep.hpp
index b32976ff1817..9d29feb3ae6a 100644
--- a/packages/rol/src/step/ROL_AugmentedLagrangianStep.hpp
+++ b/packages/rol/src/step/ROL_AugmentedLagrangianStep.hpp
@@ -450,6 +450,7 @@ class AugmentedLagrangianStep : public Step<Real> {
     algo_state.snorm = s.norm();
     algo_state.iter++;
     // Update objective function value
+    obj.update(x);
     algo_state.value = augLag.getObjectiveValue(x);
     // Update constraint value
     augLag.getConstraintVec(*(state->constraintVec),x);
diff --git a/packages/rol/src/step/augmentedlagrangian/ROL_QuadraticPenalty.hpp b/packages/rol/src/step/augmentedlagrangian/ROL_QuadraticPenalty.hpp
index 29dd5f426013..b9504bc34ebf 100644
--- a/packages/rol/src/step/augmentedlagrangian/ROL_QuadraticPenalty.hpp
+++ b/packages/rol/src/step/augmentedlagrangian/ROL_QuadraticPenalty.hpp
@@ -180,6 +180,7 @@ class QuadraticPenalty : public Objective<Real> {
   virtual void hessVec( Vector<Real> &hv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
     // Apply objective Hessian to a vector
     if (HessianApprox_ < 3) {
+      con_->update(x);
       con_->applyJacobian(*primalConVector_,v,x,tol);
       con_->applyAdjointJacobian(hv,primalConVector_->dual(),x,tol);
       if (!useScaling_) {