Skip to content

Commit

Permalink
Piro: work on the Hessain based implementation of ROL dot product (tr…
Browse files Browse the repository at this point in the history 
…ilinos#9501)

Generalization of  the parameter list for defining Hesian based dot product.
Now it looks like:

  Matrix Based Dot Product:
    Matrix Type: Hessian Of Response
    Matrix Types:
      Hessian Of Response:
        Response Index: 0
          Remove Mean Of The Right-hand Side: false
          Block Diagonal Solver:
            Block 0:
             Linear Solver Type: Belos

At the moment Matrix Type can either be the Identity (usual l2 dot product) or the Hessian of a Response.
One also has to provide the Idex of the Reponse used to compute the Hessian dot product

Modified Piro test input file accordingly.
  • Loading branch information
@mperego
mperego authored Jul 29, 2021
1 parent 4a9c5bc commit 19be9ac
Show file tree
Hide file tree
Showing 3 changed files with 74 additions and 24 deletions.
48 changes: 32 additions & 16 deletions packages/piro/src/Piro_PerformAnalysis.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -493,8 +493,6 @@ Piro::PerformROLAnalysis(
}

bool useFullSpace = rolParams.get("Full Space",false);
bool useHessianDotProduct = rolParams.get("Hessian Dot Product",false);
bool removeMeanOfTheRHS = rolParams.get("Remove Mean Of The Right-hand Side",false);

*out << "\nROL options:" << std::endl;
rolParams.sublist("ROL Options").print(*out);
Expand All @@ -510,16 +508,37 @@ Piro::PerformROLAnalysis(
ROL::Ptr<ROL::Algorithm<double> > algo;
algo = ROL::makePtr<ROL::Algorithm<double>>(step, status, true);

#ifdef HAVE_PIRO_TEKO
bool useHessianDotProduct = false;
Teuchos::ParameterList hessianDotProductList;
if(rolParams.isSublist("Matrix Based Dot Product")) {
const Teuchos::ParameterList& matrixDotProductList = rolParams.sublist("Matrix Based Dot Product");
auto matrixType = matrixDotProductList.get<std::string>("Matrix Type");
if(matrixType == "Hessian Of Response") {
useHessianDotProduct = true;
hessianDotProductList = matrixDotProductList.sublist("Matrix Types").sublist("Hessian Of Response");
}
else if (matrixType == "Identity")
useHessianDotProduct = false;
else {
TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter,
std::endl << "Error in Piro::PerformROLAnalysis: " <<
"Matrix Type not recognized. Available options are: \n" <<
"\"Identity\" and \"Hessian Of Response\""<<std::endl);
}
}

#ifdef HAVE_PIRO_TEKO
Teko::LinearOp H, invH;
if (useHessianDotProduct) {
int hessianResponseIndex = hessianDotProductList.get<int>("Response Index");
*out << "\nStart the computation of H_pp" << std::endl;
Teko::BlockedLinearOp bH = Teko::createBlockedOp();
obj.hessian_22(bH, rol_x, rol_p);
obj.hessian_22(bH, rol_x, rol_p, hessianResponseIndex);
*out << "End of the computation of H_pp" << std::endl;

int numBlocks = bH->productRange()->numBlocks();

/* Not using defaults to increase user awareness
Teuchos::ParameterList defaultParamList;
string defaultSolverType = "Belos";
defaultParamList.set("Linear Solver Type", "Belos");
Expand All @@ -532,22 +551,14 @@ Piro::PerformROLAnalysis(
belosList.sublist("Solver Types").sublist("Pseudo Block CG").set("Output Frequency", 100);
belosList.sublist("VerboseObject").set("Verbosity Level", "medium");
defaultParamList.set("Preconditioner Type", "None");
*/

Teuchos::ParameterList dHess;
if(rolParams.isSublist("Hessian Diagonal Inverse"))
dHess = rolParams.sublist("Hessian Diagonal Inverse");

Teuchos::ParameterList dHess = hessianDotProductList.sublist("Block Diagonal Solver");
std::vector<Teko::LinearOp> diag(numBlocks);

for (int i=0; i<numBlocks; ++i) {
string blockName = "Block "+std::to_string(i);
string blockSolverType = "Block solver type "+std::to_string(i);
Teuchos::ParameterList pl;
if(dHess.isSublist(blockName))
pl = dHess.sublist(blockName);
else
pl = defaultParamList;
string solverType = dHess.get<string>(blockSolverType, defaultSolverType);
Teuchos::ParameterList pl = dHess.sublist(blockName);
std::string solverType = pl.get<string>("Linear Solver Type");
diag[i] = Teko::buildInverse(*Teko::invFactoryFromParamList(pl, solverType), Teko::getBlock(i, i, bH));
}

Expand All @@ -558,6 +569,10 @@ Piro::PerformROLAnalysis(
H = Teuchos::null;
invH = Teuchos::null;
}
#else
TEUCHOS_TEST_FOR_EXCEPTION(useHessianDotProduct, Teuchos::Exceptions::InvalidParameter,
std::endl << "Error in Piro::PerformROLAnalysis: " <<
"Teko is required for computing the Hessian based dot Product"<<std::endl);
#endif


Expand All @@ -567,6 +582,7 @@ Piro::PerformROLAnalysis(
//::Thyra::randomize<double>( 0.5, 2.0, scaling_vector_x.ptr());
ROL::PrimalScaledThyraVector<double> rol_x_primal(x, scaling_vector_x);
#ifdef HAVE_PIRO_TEKO
bool removeMeanOfTheRHS = hessianDotProductList.get("Remove Mean Of The Right-hand Side",false);
ROL::PrimalHessianScaledThyraVector<double> rol_p_primal(p, H, invH, removeMeanOfTheRHS);
#else
Teuchos::RCP<Thyra::VectorBase<double> > scaling_vector_p = p->clone_v();
Expand Down
15 changes: 8 additions & 7 deletions packages/piro/src/Piro_ThyraProductME_Objective_SimOpt.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -268,14 +268,15 @@ class ThyraProductME_Objective_SimOpt : public ROL::Objective_SimOpt<Real> {

void hessian_22(const Teuchos::RCP<Thyra::PhysicallyBlockedLinearOpBase<Real>> H,
const ROL::Vector<Real> &u,
const ROL::Vector<Real> &z) {
const ROL::Vector<Real> &z,
const int g_idx) {
if(verbosityLevel >= Teuchos::VERB_MEDIUM)
*out << "ROL::ThyraProductME_Objective_SimOpt::hessian_22" << std::endl;

Thyra::ModelEvaluatorBase::OutArgs<Real> outArgs = thyra_model.createOutArgs();
bool supports_deriv = true;
for(std::size_t i=0; i<p_indices.size(); ++i)
supports_deriv = supports_deriv && outArgs.supports(Thyra::ModelEvaluatorBase::OUT_ARG_hess_g_pp, g_index, p_indices[i], p_indices[i]);
supports_deriv = supports_deriv && outArgs.supports(Thyra::ModelEvaluatorBase::OUT_ARG_hess_g_pp, g_idx, p_indices[i], p_indices[i]);

if(supports_deriv) { //use derivatives computed by model evaluator
const ROL::ThyraVector<Real> & thyra_p = dynamic_cast<const ROL::ThyraVector<Real>&>(z);
Expand All @@ -294,18 +295,18 @@ class ThyraProductME_Objective_SimOpt : public ROL::Objective_SimOpt<Real> {
}
inArgs.set_x(thyra_x.getVector());

Teuchos::RCP< Thyra::VectorBase<Real> > multiplier_g = Thyra::createMember<Real>(thyra_model.get_g_multiplier_space(g_index));
Teuchos::RCP< Thyra::VectorBase<Real> > multiplier_g = Thyra::createMember<Real>(thyra_model.get_g_multiplier_space(g_idx));
Thyra::put_scalar(1.0, multiplier_g.ptr());
inArgs.set_g_multiplier(g_index, multiplier_g);
inArgs.set_g_multiplier(g_idx, multiplier_g);

Thyra::ModelEvaluatorBase::OutArgs<Real> outArgs = thyra_model.createOutArgs();

for(std::size_t i=0; i<p_indices.size(); ++i) {
bool supports_deriv = outArgs.supports(Thyra::ModelEvaluatorBase::OUT_ARG_hess_g_pp, g_index, p_indices[i], p_indices[i]);
bool supports_deriv = outArgs.supports(Thyra::ModelEvaluatorBase::OUT_ARG_hess_g_pp, g_idx, p_indices[i], p_indices[i]);
ROL_TEST_FOR_EXCEPTION( !supports_deriv, std::logic_error, "ROL::ThyraProductME_Objective_SimOpt: H_pp is not supported");

Teuchos::RCP<Thyra::LinearOpBase<Real>> hess_g_pp = thyra_model.create_hess_g_pp(g_index, p_indices[i], p_indices[i]);
outArgs.set_hess_g_pp(g_index, p_indices[i], p_indices[i], hess_g_pp);
Teuchos::RCP<Thyra::LinearOpBase<Real>> hess_g_pp = thyra_model.create_hess_g_pp(g_idx, p_indices[i], p_indices[i]);
outArgs.set_hess_g_pp(g_idx, p_indices[i], p_indices[i], hess_g_pp);
H->setBlock(p_indices[i], p_indices[i], hess_g_pp);
}
H->endBlockFill();
Expand Down
35 changes: 34 additions & 1 deletion packages/piro/test/_input_Analysis_ROL_Tpetra.xml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -181,7 +181,40 @@

<Parameter name="Step Method" type="string" value="Line Search" />
<Parameter name="Full Space" type="bool" value="false" />
<Parameter name="Hessian Dot Product" type="bool" value="true" />

<ParameterList name="Matrix Based Dot Product">
<Parameter name="Matrix Type" type="string" value="Hessian Of Response" />
<ParameterList name="Matrix Types">
<ParameterList name="Hessian Of Response">
<Parameter name="Response Index" type="int" value="0" />
<ParameterList name="Block Diagonal Solver">
<ParameterList name="Block 0">
<Parameter name="Linear Solver Type" type="string" value="Belos" />
<ParameterList name="Linear Solver Types">
<ParameterList name="Belos">
<Parameter name="Solver Type" type="string" value="Block GMRES" />
<ParameterList name="Solver Types">
<ParameterList name="Block GMRES">
<Parameter name="Maximum Iterations" type="int" value="200" />
<Parameter name="Num Blocks" type="int" value="200" />
<Parameter name="Verbosity" type="int" value="33" />
<Parameter name="Output Frequency" type="int" value="20" />
<Parameter name="Output Style" type="int" value="1" />
<Parameter name="Convergence Tolerance" type="double" value="1e-7" />
</ParameterList>
</ParameterList>
<ParameterList name="VerboseObject">
<Parameter name="Verbosity Level" type="string" value="medium" />
</ParameterList>
</ParameterList>
</ParameterList>
<Parameter name="Preconditioner Type" type="string" value="None" />
</ParameterList>
</ParameterList>
</ParameterList>
</ParameterList>
</ParameterList>

<Parameter name="Use NOX Solver" type="bool" value="false" />

<!-- =========== BEGIN ROL INPUT PARAMETER SUBLIST =========== -->
Expand Down

0 comments on commit 19be9ac

Please sign in to comment.