58 #include "Teuchos_GlobalMPISession.hpp"
65 int main(
int argc,
char *argv[]) {
67 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
70 int iprint = argc - 1;
71 ROL::Ptr<std::ostream> outStream;
74 outStream = ROL::makePtrFromRef(std::cout);
76 outStream = ROL::makePtrFromRef(bhs);
88 ROL::ParameterList parlist;
89 std::string stepname =
"Trust Region";
90 parlist.sublist(
"Step").sublist(stepname).set(
"Subproblem Solver",
"Truncated CG");
91 parlist.sublist(
"General").sublist(
"Krylov").set(
"Iteration Limit",50);
92 parlist.sublist(
"General").sublist(
"Krylov").set(
"Relative Tolerance",1e-2);
93 parlist.sublist(
"General").sublist(
"Krylov").set(
"Absolute Tolerance",1e-4);
94 parlist.sublist(
"Status Test").set(
"Gradient Tolerance",1.e-12);
95 parlist.sublist(
"Status Test").set(
"Step Tolerance",1.e-14);
96 parlist.sublist(
"Status Test").set(
"Iteration Limit",100);
100 ROL::Ptr<std::vector<RealT> > x_ptr = ROL::makePtr<std::vector<RealT>>(dim, 0.0);
102 for (
int i=0; i<dim; i++) {
108 algo.
run(x, obj,
true, *outStream);
112 H = ROL::computeDenseHessian<RealT>(obj, x);
116 std::vector<std::vector<RealT> > eigenvals = ROL::computeEigenvalues<RealT>(H);
118 *outStream <<
"\nEigenvalues:\n";
119 for (
unsigned i=0; i<(eigenvals[0]).size(); i++) {
121 *outStream << std::right
122 << std::setw(28) <<
"Real"
123 << std::setw(28) <<
"Imag"
126 *outStream << std::scientific << std::setprecision(16) << std::right
127 << std::setw(28) << (eigenvals[0])[i]
128 << std::setw(28) << (eigenvals[1])[i]
135 std::vector<std::vector<RealT> > genEigenvals = ROL::computeGenEigenvalues<RealT>(H, M);
137 *outStream <<
"\nGeneralized eigenvalues:\n";
138 for (
unsigned i=0; i<(genEigenvals[0]).size(); i++) {
140 *outStream << std::right
141 << std::setw(28) <<
"Real"
142 << std::setw(28) <<
"Imag"
145 *outStream << std::scientific << std::setprecision(16) << std::right
146 << std::setw(28) << (genEigenvals[0])[i]
147 << std::setw(28) << (genEigenvals[1])[i]
152 std::sort((eigenvals[0]).begin(), (eigenvals[0]).end());
153 std::sort((eigenvals[1]).begin(), (eigenvals[1]).end());
154 std::sort((genEigenvals[0]).begin(), (genEigenvals[0]).end());
155 std::sort((genEigenvals[1]).begin(), (genEigenvals[1]).end());
157 RealT errtol = std::sqrt(ROL::ROL_EPSILON<RealT>());
158 for (
unsigned i=0; i<(eigenvals[0]).size(); i++) {
159 if ( std::abs( (genEigenvals[0])[i] - (eigenvals[0])[i] ) > errtol*((eigenvals[0])[i]+ROL::ROL_THRESHOLD<RealT>()) ) {
161 *outStream << std::scientific << std::setprecision(20) <<
"Real genEigenvals - eigenvals (" << i <<
") = " << std::abs( (genEigenvals[0])[i] - (eigenvals[0])[i] ) <<
" > " << errtol*((eigenvals[0])[i]+1e4*ROL::ROL_THRESHOLD<RealT>()) <<
"\n";
163 if ( std::abs( (genEigenvals[1])[i] - (eigenvals[1])[i] ) > errtol*((eigenvals[1])[i]+ROL::ROL_THRESHOLD<RealT>()) ) {
165 *outStream << std::scientific << std::setprecision(20) <<
"Imag genEigenvals - eigenvals (" << i <<
") = " << std::abs( (genEigenvals[1])[i] - (eigenvals[1])[i] ) <<
" > " << errtol*((eigenvals[1])[i]+ROL::ROL_THRESHOLD<RealT>()) <<
"\n";
170 Teuchos::SerialDenseMatrix<int, RealT> invH = ROL::computeInverse<RealT>(H);
171 Teuchos::SerialDenseMatrix<int, RealT> HinvH(H);
174 HinvH.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, H, invH, 0.0);
177 if (HinvH.normOne() > errtol) {
179 *outStream << std::scientific << std::setprecision(20) <<
"1-norm of H*inv(H) - I = " << HinvH.normOne() <<
" > " << errtol <<
"\n";
183 stepname =
"Line Search";
184 parlist.sublist(
"Step").sublist(stepname).sublist(
"Descent Method").set(
"Type",
"Newton-Krylov");
188 for (
int i=0; i<dim; i++) {
193 newton_algo.
run(x, obj,
true, *outStream);
195 ROL::Ptr<const ROL::AlgorithmState<RealT> > new_state = newton_algo.
getState();
196 ROL::Ptr<const ROL::AlgorithmState<RealT> > old_state = algo.
getState();
197 *outStream <<
"old_optimal_value = " << old_state->value << std::endl;
198 *outStream <<
"new_optimal_value = " << new_state->value << std::endl;
199 if ( std::abs(new_state->value - old_state->value) / std::abs(old_state->value) > errtol ) {
201 *outStream << std::scientific << std::setprecision(20) <<
"\nabs(new_optimal_value - old_optimal_value) / abs(old_optimal_value) = " << std::abs(new_state->value - old_state->value) / std::abs(old_state->value) <<
" > " << errtol <<
"\n";
205 catch (std::logic_error err) {
206 *outStream << err.what() <<
"\n";
211 std::cout <<
"End Result: TEST FAILED\n";
213 std::cout <<
"End Result: TEST PASSED\n";
Contains definitions of custom data types in ROL.
virtual std::vector< std::string > run(Vector< Real > &x, Objective< Real > &obj, bool print=false, std::ostream &outStream=std::cout, bool printVectors=false, std::ostream &vectorStream=std::cout)
Run algorithm on unconstrained problems (Type-U). This is the primary Type-U interface.
Defines a no-output stream class ROL::NullStream and a function makeStreamPtr which either wraps a re...
Contains definitions for Poisson material inversion.
Provides the ROL::Vector interface for scalar values, to be used, for example, with scalar constraint...
Provides an interface to run optimization algorithms.
ROL::Ptr< const AlgorithmState< Real > > getState(void) const
int dimension() const
Return dimension of the vector space.
Teuchos::SerialDenseMatrix< int, Real > computeDotMatrix(const Vector< Real > &x)
basic_nullstream< char, char_traits< char >> nullstream
int main(int argc, char *argv[])
Poisson material inversion.