61 #include "Teuchos_GlobalMPISession.hpp"
89 Teuchos::SerialDenseMatrix<int, Real>
H_;
107 hu_ = 1.0/((Real)
nu_+1.0);
111 void apply_mass(std::vector<Real> &Mz,
const std::vector<Real> &z ) {
115 Mz[i] =
hu_/6.0*(2.0*z[i] + z[i+1]);
117 else if ( i == nu_-1 ) {
118 Mz[i] =
hu_/6.0*(z[i-1] + 2.0*z[i]);
121 Mz[i] =
hu_/6.0*(z[i-1] + 4.0*z[i] + z[i+1]);
127 return (x <= 0.5) ? 1.0 : 0.0;
131 const std::vector<Real> &d,
const std::vector<Real> &u,
132 bool addBC =
true ) {
135 for (
uint i = 0; i <
nu_; i++) {
137 Bd[i] = 1.0/
hu_*( u[i]*d[i] + (u[i]-u[i+1])*d[i+1] );
139 else if ( i == nu_-1 ) {
140 Bd[i] = 1.0/
hu_*( (u[i]-u[i-1])*d[i] + u[i]*d[i+1] );
143 Bd[i] = 1.0/
hu_*( (u[i]-u[i-1])*d[i] + (u[i]-u[i+1])*d[i+1] );
153 const std::vector<Real> &d,
const std::vector<Real> &u,
154 bool addBC =
true ) {
157 for (
uint i = 0; i <
nz_; i++) {
159 Bd[i] = 1.0/
hu_*u[i]*d[i];
161 else if ( i == nz_-1 ) {
162 Bd[i] = 1.0/
hu_*u[i-1]*d[i-1];
165 Bd[i] = 1.0/
hu_*( (u[i]-u[i-1])*(d[i]-d[i-1]) );
177 std::vector<Real> d(
nu_,1.0);
178 std::vector<Real> o(
nu_-1,1.0);
179 for (
uint i = 0; i <
nu_; i++ ) {
180 d[i] = (z[i] + z[i+1])/
hu_;
191 Teuchos::LAPACK<int,Real> lp;
195 lp.PTTRF(nu_,&d[0],&o[0],&info);
196 lp.PTTRS(nu_,nhrs,&d[0],&o[0],&u[0],ldb,&info);
203 for (
uint i = 0; i <
nu_; i++ ) {
204 d[i] = (z[i] + z[i+1])/
hu_;
211 for (
uint i = 0; i <
nu_; i++) {
218 Teuchos::LAPACK<int,Real> lp;
222 lp.PTTRF(nu_,&d[0],&o[0],&info);
223 lp.PTTRS(nu_,nhrs,&d[0],&o[0],&p[0],ldb,&info);
227 const std::vector<Real> &u,
const std::vector<Real> &z) {
231 for (
uint i = 0; i <
nu_; i++ ) {
232 d[i] = (z[i] + z[i+1])/
hu_;
242 Teuchos::LAPACK<int,Real> lp;
246 lp.PTTRF(nu_,&d[0],&o[0],&info);
247 lp.PTTRS(nu_,nhrs,&d[0],&o[0],&w[0],ldb,&info);
251 const std::vector<Real> &v,
const std::vector<Real> &p,
252 const std::vector<Real> &u,
const std::vector<Real> &z) {
256 for (
uint i = 0; i <
nu_; i++ ) {
257 d[i] = (z[i] + z[i+1])/
hu_;
266 std::vector<Real> res(nu_,0.0);
268 for (
uint i = 0; i <
nu_; i++) {
272 Teuchos::LAPACK<int,Real> lp;
276 lp.PTTRF(nu_,&d[0],&o[0],&info);
277 lp.PTTRS(nu_,nhrs,&d[0],&o[0],&q[0],ldb,&info);
285 Real tol = std::sqrt(ROL::ROL_EPSILON<Real>());
287 ROL::Ptr<V> e = z.
clone();
288 ROL::Ptr<V> h = z.
clone();
289 for (
uint i = 0; i <
nz_; i++ ) {
292 for (
uint j = 0; j <
nz_; j++ ) {
294 (
H_)(j,i) = e->dot(*h);
297 std::vector<vector> eigenvals = ROL::computeEigenvalues<Real>(
H_);
298 std::sort((eigenvals[0]).begin(), (eigenvals[0]).end());
299 Real inertia = (eigenvals[0])[0];
300 Real correction = 0.0;
301 if ( inertia <= 0.0 ) {
302 correction = (1.0+std::sqrt(ROL::ROL_EPSILON<Real>()))*std::abs(inertia);
303 if ( inertia == 0.0 ) {
305 while ( eigenvals[0][cnt] == 0.0 ) {
308 correction = std::sqrt(ROL::ROL_EPSILON<Real>())*eigenvals[0][cnt];
309 if ( cnt == nz_-1 ) {
313 for (
uint i = 0; i <
nz_; i++ ) {
314 (
H_)(i,i) += correction;
324 ROL::Ptr<const vector> zp =
getVector(z);
343 res =
hu_/6.0*(2.0*res1 + res2)*res1;
345 else if ( i == nu_-1 ) {
348 res = hu_/6.0*(res1 + 2.0*res2)*res2;
354 res = hu_/6.0*(res1 + 4.0*res2 + res3)*res2;
365 ROL::Ptr<const vector> zp =
getVector(z);
379 for (
uint i = 0; i <
nz_; i++ ) {
405 ROL::Ptr<const vector> vp =
getVector(v);
406 ROL::Ptr<const vector> zp =
getVector(z);
428 std::vector<Real> tmp(
nz_,0.0);
442 using ROL::constPtrCast;
447 ROL::Ptr<vector> vp = ROL::constPtrCast<vector>(
getVector(v));
449 Teuchos::SerialDenseVector<int, Real> hv_teuchos(Teuchos::View, &((*hvp)[0]), static_cast<int>(
nz_));
450 Teuchos::SerialDenseVector<int, Real> v_teuchos(Teuchos::View, &(( *vp)[0]), static_cast<int>(
nz_));
451 hv_teuchos.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0,
H_, v_teuchos, 0.0);
460 int main(
int argc,
char *argv[]) {
462 typedef std::vector<RealT> vector;
470 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
473 int iprint = argc - 1;
474 ROL::Ptr<std::ostream> outStream;
477 outStream = ROL::makePtrFromRef(std::cout);
479 outStream = ROL::makePtrFromRef(bhs);
492 ROL::Ptr<vector> x_ptr = ROL::makePtr<vector>(dim, 0.0);
493 ROL::Ptr<vector> y_ptr = ROL::makePtr<vector>(dim, 0.0);
496 for (uint i=0; i<dim; i++) {
497 (*x_ptr)[i] = (
RealT)rand()/(
RealT)RAND_MAX + 1.e2;
498 (*y_ptr)[i] = (
RealT)rand()/(
RealT)RAND_MAX + 1.e2;
507 ROL::Ptr<vector> l_ptr = ROL::makePtr<vector>(dim,1.0);
508 ROL::Ptr<vector> u_ptr = ROL::makePtr<vector>(dim,10.0);
510 ROL::Ptr<V> lo = ROL::makePtr<SV>(l_ptr);
511 ROL::Ptr<V> up = ROL::makePtr<SV>(u_ptr);
515 ROL::ParameterList parlist;
518 parlist.sublist(
"General").sublist(
"Krylov").set(
"Absolute Tolerance",1.e-8);
519 parlist.sublist(
"General").sublist(
"Krylov").set(
"Relative Tolerance",1.e-4);
520 parlist.sublist(
"General").sublist(
"Krylov").set(
"Iteration Limit",static_cast<int>(dim));
522 parlist.sublist(
"Step").sublist(
"Primal Dual Active Set").set(
"Relative Step Tolerance",1.e-8);
523 parlist.sublist(
"Step").sublist(
"Primal Dual Active Set").set(
"Relative Gradient Tolerance",1.e-6);
524 parlist.sublist(
"Step").sublist(
"Primal Dual Active Set").set(
"Iteration Limit", 10);
525 parlist.sublist(
"Step").sublist(
"Primal Dual Active Set").set(
"Dual Scaling",(alpha>0.0)?alpha:1.e-4);
526 parlist.sublist(
"General").sublist(
"Secant").set(
"Use as Hessian",
true);
528 parlist.sublist(
"Status Test").set(
"Gradient Tolerance",1.e-12);
529 parlist.sublist(
"Status Test").set(
"Step Tolerance",1.e-14);
530 parlist.sublist(
"Status Test").set(
"Iteration Limit",100);
537 algo.
run(x,obj,icon,
true,*outStream);
541 file.open(
"control_PDAS.txt");
542 for ( uint i = 0; i < dim; i++ ) {
543 file << (*x_ptr)[i] <<
"\n";
549 parlist.sublist(
"General").sublist(
"Secant").set(
"Use as Hessian",
false);
550 parlist.sublist(
"Step").sublist(
"Trust Region").set(
"Subproblem Solver",
"Truncated CG");
551 parlist.sublist(
"Step").sublist(
"Trust Region").set(
"Initial Radius", 1e3);
552 parlist.sublist(
"Step").sublist(
"Trust Region").set(
"Maximum Radius", 1e8);
557 algo_tr.
run(y,obj,icon,
true,*outStream);
559 std::ofstream file_tr;
560 file_tr.open(
"control_TR.txt");
561 for ( uint i = 0; i < dim; i++ ) {
562 file_tr << (*y_ptr)[i] <<
"\n";
566 std::vector<RealT> u;
568 std::ofstream file_u;
569 file_u.open(
"state.txt");
570 for ( uint i = 0; i < (dim-1); i++ ) {
571 file_u << u[i] <<
"\n";
575 ROL::Ptr<V> diff = x.clone();
578 RealT error = diff->norm()/std::sqrt((
RealT)dim-1.0);
579 std::cout <<
"\nError between PDAS solution and TR solution is " << error <<
"\n";
580 errorFlag = ((error > 1.e2*std::sqrt(ROL::ROL_EPSILON<RealT>())) ? 1 : 0);
583 catch (std::logic_error err) {
584 *outStream << err.what() <<
"\n";
589 std::cout <<
"End Result: TEST FAILED\n";
591 std::cout <<
"End Result: TEST PASSED\n";
Provides the interface to evaluate objective functions.
typename PV< Real >::size_type size_type
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
std::vector< Real > vector
virtual ROL::Ptr< Vector > basis(const int i) const
Return i-th basis vector.
void solve_state_equation(std::vector< Real > &u, const std::vector< Real > &z)
void apply_transposed_linearized_control_operator(std::vector< Real > &Bd, const std::vector< Real > &z, const std::vector< Real > &d, const std::vector< Real > &u, bool addBC=true)
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.
Contains definitions for helper functions in ROL.
Defines the linear algebra or vector space interface.
Defines a no-output stream class ROL::NullStream and a function makeStreamPtr which either wraps a re...
void solve_adjoint_sensitivity_equation(std::vector< Real > &q, const std::vector< Real > &w, const std::vector< Real > &v, const std::vector< Real > &p, const std::vector< Real > &u, const std::vector< Real > &z)
virtual std::vector< std::vector< Real > > checkGradient(const Vector< Real > &x, const Vector< Real > &d, const bool printToStream=true, std::ostream &outStream=std::cout, const int numSteps=ROL_NUM_CHECKDERIV_STEPS, const int order=1)
Finite-difference gradient check.
void apply_mass(std::vector< Real > &Mz, const std::vector< Real > &z)
Real value(const ROL::Vector< Real > &z, Real &tol)
Compute value.
void update(const ROL::Vector< Real > &z, bool flag, int iter)
Update objective function.
void deactivateInertia(void)
Provides an interface to run optimization algorithms.
Provides the elementwise interface to apply upper and lower bound constraints.
void hessVec_true(ROL::Vector< Real > &hv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &z, Real &tol)
ROL::StdVector< Real > SV
void solve_adjoint_equation(std::vector< Real > &p, const std::vector< Real > &u, const std::vector< Real > &z)
void solve_state_sensitivity_equation(std::vector< Real > &w, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z)
basic_nullstream< char, char_traits< char >> nullstream
void hessVec(ROL::Vector< Real > &hv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &z, Real &tol)
Apply Hessian approximation to vector.
int main(int argc, char *argv[])
virtual std::vector< std::vector< Real > > checkHessVec(const Vector< Real > &x, const Vector< Real > &v, const bool printToStream=true, std::ostream &outStream=std::cout, const int numSteps=ROL_NUM_CHECKDERIV_STEPS, const int order=1)
Finite-difference Hessian-applied-to-vector check.
void apply_linearized_control_operator(std::vector< Real > &Bd, const std::vector< Real > &z, const std::vector< Real > &d, const std::vector< Real > &u, bool addBC=true)
void activateInertia(void)
void hessVec_inertia(ROL::Vector< Real > &hv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &z, Real &tol)
Objective_PoissonInversion(int nz=32, Real alpha=1.e-4)
void gradient(ROL::Vector< Real > &g, const ROL::Vector< Real > &z, Real &tol)
Compute gradient.
Real evaluate_target(Real x)
ROL::Ptr< const vector > getVector(const V &x)
Teuchos::SerialDenseMatrix< int, Real > H_
ROL::Ptr< vector > getVector(V &x)