Skip to content

Commit

Permalink
Merge pull request #8584 from kliegeois/hessian2
Browse files Browse the repository at this point in the history 
Piro: add HessianScaledThyraVector to ROL adapters and use it in Piro_PerformAnalysis
  • Loading branch information
@kliegeois
kliegeois authored Jan 25, 2021
2 parents 5e36d3d + a93322d commit e76f001
Show file tree
Hide file tree
Showing 11 changed files with 1,219 additions and 50 deletions.
60 changes: 56 additions & 4 deletions packages/piro/src/Piro_PerformAnalysis.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -63,6 +63,7 @@
#ifdef HAVE_PIRO_ROL
#include "ROL_ThyraVector.hpp"
#include "ROL_ScaledThyraVector.hpp"
#include "ROL_HessianScaledThyraVector.hpp"
#include "ROL_Thyra_BoundConstraint.hpp"
#include "ROL_ThyraME_Objective.hpp"
#include "ROL_ThyraProductME_Objective.hpp"
Expand All @@ -78,8 +79,11 @@
#include "Thyra_VectorDefaultBase.hpp"
#include "Thyra_DefaultProductVectorSpace.hpp"
#include "Thyra_DefaultProductVector.hpp"
#include "Thyra_DefaultBlockedLinearOp.hpp"
#endif

#include "Teko_InverseLibrary.hpp"
#include "Teko_PreconditionerFactory.hpp"

using std::cout; using std::endl; using std::string;
using Teuchos::RCP; using Teuchos::rcp; using Teuchos::ParameterList;
Expand Down Expand Up @@ -484,6 +488,8 @@ 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 @@ -499,15 +505,61 @@ Piro::PerformROLAnalysis(
ROL::Ptr<ROL::Algorithm<double> > algo;
algo = ROL::makePtr<ROL::Algorithm<double>>(step, status,false);

Teko::LinearOp H, invH;
if (useHessianDotProduct) {
*out << "\nStart the computation of H_pp" << std::endl;
Teko::BlockedLinearOp bH = Teko::createBlockedOp();
obj.hessian_22(bH, rol_x, rol_p);
*out << "End of the computation of H_pp" << std::endl;

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

Teuchos::ParameterList defaultParamList;
string defaultSolverType = "Belos";
defaultParamList.set("Linear Solver Type", "Belos");
Teuchos::ParameterList& belosList = defaultParamList.sublist("Linear Solver Types").sublist("Belos");
belosList.set("Solver Type", "Pseudo Block CG");
belosList.sublist("Solver Types").sublist("Pseudo Block CG").set<int>("Maximum Iterations", 1000);
belosList.sublist("Solver Types").sublist("Pseudo Block CG").set<double>("Convergence Tolerance", 1e-4);
belosList.sublist("Solver Types").sublist("Pseudo Block CG").set<int>("Num Blocks", 1000);
belosList.sublist("Solver Types").sublist("Pseudo Block CG").set("Verbosity", 0x7f);
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");

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);
diag[i] = Teko::buildInverse(*Teko::invFactoryFromParamList(pl, solverType), Teko::getBlock(i, i, bH));
}

H = Teko::toLinearOp(bH);
invH = Teko::createBlockUpperTriInverseOp(bH, diag);
}
else {
H = Teuchos::null;
invH = Teuchos::null;
}


//this is for testing the PrimalScaledThyraVector. At the moment the scaling is set to 1, so it is not changing the dot product
Teuchos::RCP<Thyra::VectorBase<double> > scaling_vector_p = p->clone_v();
Teuchos::RCP<Thyra::VectorBase<double> > scaling_vector_x = x->clone_v();
::Thyra::put_scalar<double>( 1.0, scaling_vector_p.ptr());
::Thyra::put_scalar<double>( 1.0, scaling_vector_x.ptr());
//::Thyra::randomize<double>( 0.5, 2.0, scaling_vector_p.ptr());
//::Thyra::randomize<double>( 0.5, 2.0, scaling_vector_x.ptr());
ROL::PrimalScaledThyraVector<double> rol_x_primal(x, scaling_vector_x);
ROL::PrimalScaledThyraVector<double> rol_p_primal(p, scaling_vector_p);
ROL::PrimalHessianScaledThyraVector<double> rol_p_primal(p, H, invH, removeMeanOfTheRHS);

// Run Algorithm
std::vector<std::string> output;
Expand Down
2 changes: 2 additions & 0 deletions packages/piro/test/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -116,6 +116,8 @@ IF (Piro_ENABLE_NOX)
Main_AnalysisDriver_Tpetra.cpp
MockModelEval_A_Tpetra.cpp
MockModelEval_A_Tpetra.hpp
MockModelEval_B_Tpetra.cpp
MockModelEval_B_Tpetra.hpp
NUM_MPI_PROCS 1-4
ARGS -v
PASS_REGULAR_EXPRESSION "TEST PASSED"
Expand Down
99 changes: 58 additions & 41 deletions packages/piro/test/Main_AnalysisDriver_Tpetra.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,6 +44,7 @@
#include <string>

#include "MockModelEval_A_Tpetra.hpp"
#include "MockModelEval_B_Tpetra.hpp"
//#include "ObserveSolution_Epetra.hpp"

#include "Piro_SolverFactory.hpp"
Expand Down Expand Up @@ -122,61 +123,77 @@ int main(int argc, char *argv[]) {

try {

// Create (1) a Model Evaluator and (2) a ParameterList
const RCP<Thyra::ModelEvaluator<double>> Model = rcp(new MockModelEval_A_Tpetra(appComm));
std::vector<std::string> mockModels = {"MockModelEval_A_Tpetra", "MockModelEval_B_Tpetra"};
for (auto mockModel : mockModels) {

// BEGIN Builder
const RCP<Teuchos::ParameterList> appParams = rcp(new Teuchos::ParameterList("Application Parameters"));
Teuchos::updateParametersFromXmlFile(inputFile, Teuchos::ptr(appParams.get()));
// Create (1) a Model Evaluator and (2) a ParameterList
RCP<Thyra::ModelEvaluator<double>> Model;
if (mockModel=="MockModelEval_A_Tpetra")
Model = rcp(new MockModelEval_A_Tpetra(appComm));
if (mockModel=="MockModelEval_B_Tpetra")
Model = rcp(new MockModelEval_B_Tpetra(appComm));

const RCP<Teuchos::ParameterList> piroParams = Teuchos::sublist(appParams,"Piro");
Teuchos::ParameterList& analysisParams = piroParams->sublist("Analysis");
// BEGIN Builder
const RCP<Teuchos::ParameterList> appParams = rcp(new Teuchos::ParameterList("Application Parameters"));
Teuchos::updateParametersFromXmlFile(inputFile, Teuchos::ptr(appParams.get()));

const RCP<Teuchos::ParameterList> piroParams = Teuchos::sublist(appParams,"Piro");
Teuchos::ParameterList& analysisParams = piroParams->sublist("Analysis");

Stratimikos::DefaultLinearSolverBuilder linearSolverBuilder;

#ifdef HAVE_PIRO_IFPACK2
typedef Thyra::PreconditionerFactoryBase<double> Base;
typedef Thyra::Ifpack2PreconditionerFactory<Tpetra_CrsMatrix> Impl;
linearSolverBuilder.setPreconditioningStrategyFactory(
Teuchos::abstractFactoryStd<Base, Impl>(), "Ifpack2");
#endif
Stratimikos::DefaultLinearSolverBuilder linearSolverBuilder;

#ifdef HAVE_PIRO_MUELU
using local_ordinal_type = Tpetra::Map<>::local_ordinal_type;
using global_ordinal_type = Tpetra::Map<>::global_ordinal_type;
Stratimikos::enableMueLu<local_ordinal_type, global_ordinal_type>(linearSolverBuilder);
#endif
#ifdef HAVE_PIRO_IFPACK2
typedef Thyra::PreconditionerFactoryBase<double> Base;
typedef Thyra::Ifpack2PreconditionerFactory<Tpetra_CrsMatrix> Impl;
linearSolverBuilder.setPreconditioningStrategyFactory(
Teuchos::abstractFactoryStd<Base, Impl>(), "Ifpack2");
#endif

const Teuchos::RCP<Teuchos::ParameterList> stratList = Piro::extractStratimikosParams(piroParams);
linearSolverBuilder.setParameterList(stratList);
#ifdef HAVE_PIRO_MUELU
using local_ordinal_type = Tpetra::Map<>::local_ordinal_type;
using global_ordinal_type = Tpetra::Map<>::global_ordinal_type;
Stratimikos::enableMueLu<local_ordinal_type, global_ordinal_type>(linearSolverBuilder);
#endif

const Teuchos::RCP<Teuchos::ParameterList> stratList = Piro::extractStratimikosParams(piroParams);
linearSolverBuilder.setParameterList(stratList);

const RCP<Thyra::LinearOpWithSolveFactoryBase<double>> lowsFactory =
createLinearSolveStrategy(linearSolverBuilder);

RCP<Thyra::ModelEvaluator<double>> ModelWithSolve = rcp(new Thyra::DefaultModelEvaluatorWithSolveFactory<double>(
Model, lowsFactory));
const RCP<Thyra::LinearOpWithSolveFactoryBase<double>> lowsFactory =
createLinearSolveStrategy(linearSolverBuilder);

const RCP<Thyra::ModelEvaluatorDefaultBase<double>> piro = solverFactory.createSolver(piroParams, ModelWithSolve);
RCP<Thyra::ModelEvaluator<double>> ModelWithSolve = rcp(new Thyra::DefaultModelEvaluatorWithSolveFactory<double>(
Model, lowsFactory));

// Call the analysis routine
RCP<Thyra::VectorBase<double>> p;
status = Piro::PerformAnalysis(*piro, analysisParams, p);
const RCP<Thyra::ModelEvaluatorDefaultBase<double>> piro = solverFactory.createSolver(piroParams, ModelWithSolve);

if (Teuchos::nonnull(p)) { //p might be null if the package ROL is not enabled
Thyra::DetachedVectorView<double> p_view(p);
double p_exact[2] = {1,3};
double tol = 1e-5;
// Call the analysis routine
RCP<Thyra::VectorBase<double>> p;
status = Piro::PerformAnalysis(*piro, analysisParams, p);

double l2_diff = std::sqrt(std::pow(p_view(0)-p_exact[0],2) + std::pow(p_view(1)-p_exact[1],2));
if(l2_diff > tol) {
status+=100;
if (Proc==0) {
std::cout << "\nPiro_AnalysisDrvier: Optimum parameter values are: {"
<< p_exact[0] << ", " << p_exact[1] << "}, but computed values are: {"
<< p_view(0) << ", " << p_view(1) << "}." <<
"\n Difference in l2 norm: " << l2_diff << " > tol: " << tol << std::endl;
if (Teuchos::nonnull(p)) { //p might be null if the package ROL is not enabled
Thyra::DetachedVectorView<double> p_view(p);
double p_exact[2];
if (mockModel=="MockModelEval_A_Tpetra") {
p_exact[0] = 1;
p_exact[1] = 3;
}
if (mockModel=="MockModelEval_B_Tpetra") {
p_exact[0] = 6;
p_exact[1] = 4;
}
double tol = 1e-5;

double l2_diff = std::sqrt(std::pow(p_view(0)-p_exact[0],2) + std::pow(p_view(1)-p_exact[1],2));
if(l2_diff > tol) {
status+=100;
if (Proc==0) {
std::cout << "\nPiro_AnalysisDrvier: Optimum parameter values are: {"
<< p_exact[0] << ", " << p_exact[1] << "}, but computed values are: {"
<< p_view(0) << ", " << p_view(1) << "}." <<
"\n Difference in l2 norm: " << l2_diff << " > tol: " << tol << std::endl;
}
}
}
}
Expand Down
53 changes: 48 additions & 5 deletions packages/piro/test/MockModelEval_A_Tpetra.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -88,13 +88,27 @@ MockModelEval_A_Tpetra::MockModelEval_A_Tpetra(const Teuchos::RCP<const Teuchos:

//set up jacobian graph
crs_graph = rcp(new Tpetra_CrsGraph(x_map, vecLength));
std::vector<typename Tpetra_CrsGraph::global_ordinal_type> indices(vecLength);
for (int i=0; i<vecLength; i++) indices[i]=i;
const int nodeNumElements = x_map->getNodeNumElements();
for (int i=0; i<nodeNumElements; i++)
crs_graph->insertGlobalIndices(x_map->getGlobalElement(i), vecLength, &indices[0]);
{
std::vector<typename Tpetra_CrsGraph::global_ordinal_type> indices(vecLength);
for (int i=0; i<vecLength; i++) indices[i]=i;
const int nodeNumElements = x_map->getNodeNumElements();
for (int i=0; i<nodeNumElements; i++)
crs_graph->insertGlobalIndices(x_map->getGlobalElement(i), vecLength, &indices[0]);
}
crs_graph->fillComplete();

//set up hessian graph
hess_crs_graph = rcp(new Tpetra_CrsGraph(p_map, numParameters));
if (comm->getRank() == 0)
{
std::vector<typename Tpetra_CrsGraph::global_ordinal_type> indices(numParameters);
for (int i=0; i<numParameters; i++) indices[i]=i;
const int nodeNumElements = p_map->getNodeNumElements();
for (int i=0; i<nodeNumElements; i++)
hess_crs_graph->insertGlobalIndices(p_map->getGlobalElement(i), numParameters, &indices[0]);
}
hess_crs_graph->fillComplete();

// Setup nominal values, lower and upper bounds
nominalValues = this->createInArgsImpl();
lowerBounds = this->createInArgsImpl();
Expand Down Expand Up @@ -196,6 +210,14 @@ MockModelEval_A_Tpetra::get_W_factory() const
return Teuchos::null;
}

Teuchos::RCP<Thyra::LinearOpBase<double>>
MockModelEval_A_Tpetra::create_hess_g_pp( int j, int l1, int l2 ) const
{
const Teuchos::RCP<Tpetra_Operator> H =
Teuchos::rcp(new Tpetra_CrsMatrix(hess_crs_graph));
return Thyra::createLinearOp(H);
}

Thyra::ModelEvaluatorBase::InArgs<double>
MockModelEval_A_Tpetra::getNominalValues() const
{
Expand Down Expand Up @@ -441,6 +463,27 @@ void MockModelEval_A_Tpetra::evalModelImpl(
W_out_crs->fillComplete();
}

const Teuchos::RCP<Tpetra_Operator> H_pp_out =
Teuchos::nonnull(outArgs.get_hess_g_pp(0,0,0)) ?
ConverterT::getTpetraOperator(outArgs.get_hess_g_pp(0,0,0)) :
Teuchos::null;

if (H_pp_out != Teuchos::null) {
Teuchos::RCP<Tpetra_CrsMatrix> H_pp_out_crs =
Teuchos::rcp_dynamic_cast<Tpetra_CrsMatrix>(H_pp_out, true);
H_pp_out_crs->resumeFill();
H_pp_out_crs->setAllToScalar(0.0);

if (comm->getRank() == 0) {
std::vector<double> vals = {2, 1};
std::vector<typename Tpetra_CrsGraph::global_ordinal_type> indices = {0, 1};
H_pp_out_crs->replaceGlobalValues(0, 2, &vals[0], &indices[0]);
vals[0] = 1;
H_pp_out_crs->replaceGlobalValues(1, 2, &vals[0], &indices[0]);
}
H_pp_out_crs->fillComplete();
}

if (Teuchos::nonnull(dfdp_out)) {
dfdp_out->putScalar(0.0);
auto dfdp_out_data_0 = dfdp_out->getVectorNonConst(0)->getDataNonConst();
Expand Down
5 changes: 5 additions & 0 deletions packages/piro/test/MockModelEval_A_Tpetra.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -113,6 +113,10 @@ class MockModelEval_A_Tpetra
Teuchos::RCP<const Thyra::LinearOpWithSolveFactoryBase<double>>
get_W_factory() const;

/** \brief . */
Teuchos::RCP<Thyra::LinearOpBase<double>>
create_hess_g_pp( int j, int l1, int l2 ) const;

/** \brief . */
Thyra::ModelEvaluatorBase::InArgs<double>
createInArgs() const;
Expand Down Expand Up @@ -166,6 +170,7 @@ class MockModelEval_A_Tpetra
Teuchos::RCP<const Tpetra_Map> p_map;
Teuchos::RCP<const Tpetra_Map> g_map;
Teuchos::RCP<Tpetra_CrsGraph> crs_graph;
Teuchos::RCP<Tpetra_CrsGraph> hess_crs_graph;
Teuchos::RCP<const Teuchos::Comm<int> > comm;

Teuchos::RCP<Tpetra_Vector> p_vec;
Expand Down
Loading

0 comments on commit e76f001

Please sign in to comment.