57 #include "Teuchos_oblackholestream.hpp" 58 #include "Teuchos_GlobalMPISession.hpp" 70 typedef typename vector::size_type
uint;
86 Teuchos::SerialDenseMatrix<int, Real>
H_;
88 Teuchos::RCP<const vector>
getVector(
const V& x ) {
89 using Teuchos::dyn_cast;
94 using Teuchos::dyn_cast;
102 : nz_(nz), u0_(1.0), u1_(0.0), alpha_(alpha), useCorrection_(false) {
104 hu_ = 1.0/((Real)nu_+1.0);
108 void apply_mass(std::vector<Real> &Mz,
const std::vector<Real> &z ) {
110 for (uint i=0; i<
nu_; i++) {
112 Mz[i] = hu_/6.0*(2.0*z[i] + z[i+1]);
114 else if ( i == nu_-1 ) {
115 Mz[i] = hu_/6.0*(z[i-1] + 2.0*z[i]);
118 Mz[i] = hu_/6.0*(z[i-1] + 4.0*z[i] + z[i+1]);
124 return (x <= 0.5) ? 1.0 : 0.0;
128 const std::vector<Real> &d,
const std::vector<Real> &u,
129 bool addBC =
true ) {
132 for (uint i = 0; i <
nu_; i++) {
134 Bd[i] = 1.0/hu_*( u[i]*d[i] + (u[i]-u[i+1])*d[i+1] );
136 else if ( i == nu_-1 ) {
137 Bd[i] = 1.0/hu_*( (u[i]-u[i-1])*d[i] + u[i]*d[i+1] );
140 Bd[i] = 1.0/hu_*( (u[i]-u[i-1])*d[i] + (u[i]-u[i+1])*d[i+1] );
144 Bd[ 0] -= u0_*d[ 0]/
hu_;
145 Bd[nu_-1] -= u1_*d[nz_-1]/
hu_;
150 const std::vector<Real> &d,
const std::vector<Real> &u,
151 bool addBC =
true ) {
154 for (uint i = 0; i <
nz_; i++) {
156 Bd[i] = 1.0/hu_*u[i]*d[i];
158 else if ( i == nz_-1 ) {
159 Bd[i] = 1.0/hu_*u[i-1]*d[i-1];
162 Bd[i] = 1.0/hu_*( (u[i]-u[i-1])*(d[i]-d[i-1]) );
166 Bd[ 0] -= u0_*d[ 0]/
hu_;
167 Bd[nz_-1] -= u1_*d[nu_-1]/
hu_;
174 std::vector<Real> d(nu_,1.0);
175 std::vector<Real> o(nu_-1,1.0);
176 for ( uint i = 0; i <
nu_; i++ ) {
177 d[i] = (z[i] + z[i+1])/hu_;
185 u[ 0] = z[ 0]/hu_ *
u0_;
186 u[nu_-1] = z[nz_-1]/hu_ *
u1_;
188 Teuchos::LAPACK<int,Real> lp;
192 lp.PTTRF(nu_,&d[0],&o[0],&info);
193 lp.PTTRS(nu_,nhrs,&d[0],&o[0],&u[0],ldb,&info);
200 for ( uint i = 0; i <
nu_; i++ ) {
201 d[i] = (z[i] + z[i+1])/hu_;
208 for (uint i = 0; i <
nu_; i++) {
215 Teuchos::LAPACK<int,Real> lp;
219 lp.PTTRF(nu_,&d[0],&o[0],&info);
220 lp.PTTRS(nu_,nhrs,&d[0],&o[0],&p[0],ldb,&info);
224 const std::vector<Real> &u,
const std::vector<Real> &z) {
228 for ( uint i = 0; i <
nu_; i++ ) {
229 d[i] = (z[i] + z[i+1])/hu_;
239 Teuchos::LAPACK<int,Real> lp;
243 lp.PTTRF(nu_,&d[0],&o[0],&info);
244 lp.PTTRS(nu_,nhrs,&d[0],&o[0],&w[0],ldb,&info);
248 const std::vector<Real> &v,
const std::vector<Real> &p,
249 const std::vector<Real> &u,
const std::vector<Real> &z) {
253 for ( uint i = 0; i <
nu_; i++ ) {
254 d[i] = (z[i] + z[i+1])/hu_;
263 std::vector<Real> res(nu_,0.0);
265 for (uint i = 0; i <
nu_; i++) {
269 Teuchos::LAPACK<int,Real> lp;
273 lp.PTTRF(nu_,&d[0],&o[0],&info);
274 lp.PTTRS(nu_,nhrs,&d[0],&o[0],&q[0],ldb,&info);
281 if ( flag && useCorrection_ ) {
284 RCP<V> e = z.
clone();
285 RCP<V> h = z.
clone();
286 for ( uint i = 0; i <
nz_; i++ ) {
289 for ( uint j = 0; j <
nz_; j++ ) {
291 (
H_)(j,i) = e->dot(*h);
294 std::vector<vector> eigenvals = ROL::computeEigenvalues<Real>(
H_);
295 std::sort((eigenvals[0]).begin(), (eigenvals[0]).end());
296 Real inertia = (eigenvals[0])[0];
297 Real correction = 0.0;
298 if ( inertia <= 0.0 ) {
300 if ( inertia == 0.0 ) {
302 while ( eigenvals[0][cnt] == 0.0 ) {
306 if ( cnt == nz_-1 ) {
310 for ( uint i = 0; i <
nz_; i++ ) {
311 (
H_)(i,i) += correction;
329 for (uint i=0; i<
nz_;i++) {
330 val += hz_*alpha_*0.5*(*zp)[i]*(*zp)[i];
336 for (uint i=0; i<
nu_; i++) {
340 res = hu_/6.0*(2.0*res1 + res2)*res1;
342 else if ( i == nu_-1 ) {
345 res = hu_/6.0*(res1 + 2.0*res2)*res2;
351 res = hu_/6.0*(res1 + 4.0*res2 + res3)*res2;
376 for ( uint i = 0; i <
nz_; i++ ) {
377 (*gp)[i] += hz_*alpha_*(*zp)[i];
382 if ( useCorrection_ ) {
391 useCorrection_ =
true;
395 useCorrection_ =
false;
425 std::vector<Real> tmp(nz_,0.0);
427 for (uint i=0; i <
nz_; i++) {
431 for (uint i=0; i <
nz_; i++) {
432 (*hvp)[i] += hz_*alpha_*(*vp)[i];
439 using Teuchos::rcp_const_cast;
444 RCP<vector> vp = rcp_const_cast<vector>(
getVector(v));
446 Teuchos::SerialDenseVector<int, Real> hv_teuchos(Teuchos::View, &((*hvp)[0]), static_cast<int>(nz_));
447 Teuchos::SerialDenseVector<int, Real> v_teuchos(Teuchos::View, &(( *vp)[0]), static_cast<int>(nz_));
448 hv_teuchos.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, H_, v_teuchos, 0.0);
457 int main(
int argc,
char *argv[]) {
459 typedef std::vector<RealT>
vector;
463 typedef typename vector::size_type
uint;
465 using Teuchos::RCP;
using Teuchos::rcp;
467 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
470 int iprint = argc - 1;
471 Teuchos::RCP<std::ostream> outStream;
472 Teuchos::oblackholestream bhs;
474 outStream = Teuchos::rcp(&std::cout,
false);
476 outStream = Teuchos::rcp(&bhs,
false);
489 RCP<vector> x_rcp = rcp(
new vector(dim, 0.0) );
490 RCP<vector> y_rcp = rcp(
new vector(dim, 0.0) );
493 for (uint i=0; i<dim; i++) {
494 (*x_rcp)[i] = (
RealT)rand()/(
RealT)RAND_MAX + 1.e2;
495 (*y_rcp)[i] = (
RealT)rand()/(
RealT)RAND_MAX + 1.e2;
504 RCP<vector> l_rcp = rcp(
new vector(dim,1.0) );
505 RCP<vector> u_rcp = rcp(
new vector(dim,10.0) );
507 RCP<V> lo = rcp(
new SV(l_rcp) );
508 RCP<V> up = rcp(
new SV(u_rcp) );
512 Teuchos::ParameterList parlist;
515 parlist.sublist(
"General").sublist(
"Krylov").set(
"Absolute Tolerance",1.e-8);
516 parlist.sublist(
"General").sublist(
"Krylov").set(
"Relative Tolerance",1.e-4);
517 parlist.sublist(
"General").sublist(
"Krylov").set(
"Iteration Limit",static_cast<int>(dim));
519 parlist.sublist(
"Step").sublist(
"Primal Dual Active Set").set(
"Relative Step Tolerance",1.e-8);
520 parlist.sublist(
"Step").sublist(
"Primal Dual Active Set").set(
"Relative Gradient Tolerance",1.e-6);
521 parlist.sublist(
"Step").sublist(
"Primal Dual Active Set").set(
"Iteration Limit", 10);
522 parlist.sublist(
"Step").sublist(
"Primal Dual Active Set").set(
"Dual Scaling",(alpha>0.0)?alpha:1.e-4);
523 parlist.sublist(
"General").sublist(
"Secant").set(
"Use as Hessian",
true);
525 parlist.sublist(
"Status Test").set(
"Gradient Tolerance",1.e-12);
526 parlist.sublist(
"Status Test").set(
"Step Tolerance",1.e-14);
527 parlist.sublist(
"Status Test").set(
"Iteration Limit",100);
534 algo.
run(x,obj,icon,
true,*outStream);
538 file.open(
"control_PDAS.txt");
539 for ( uint i = 0; i < dim; i++ ) {
540 file << (*x_rcp)[i] <<
"\n";
546 parlist.sublist(
"General").sublist(
"Secant").set(
"Use as Hessian",
false);
547 parlist.sublist(
"Step").sublist(
"Trust Region").set(
"Subproblem Solver",
"Truncated CG");
552 algo_tr.
run(y,obj,icon,
true,*outStream);
554 std::ofstream file_tr;
555 file_tr.open(
"control_TR.txt");
556 for ( uint i = 0; i < dim; i++ ) {
557 file_tr << (*y_rcp)[i] <<
"\n";
561 std::vector<RealT> u;
563 std::ofstream file_u;
564 file_u.open(
"state.txt");
565 for ( uint i = 0; i < (dim-1); i++ ) {
566 file_u << u[i] <<
"\n";
570 RCP<V> diff = x.clone();
573 RealT error = diff->norm()/std::sqrt((
RealT)dim-1.0);
574 std::cout <<
"\nError between PDAS solution and TR solution is " << error <<
"\n";
578 catch (std::logic_error err) {
579 *outStream << err.what() <<
"\n";
584 std::cout <<
"End Result: TEST FAILED\n";
586 std::cout <<
"End Result: TEST PASSED\n";
Provides the interface to evaluate objective functions.
std::vector< Real > 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 Teuchos::RCP< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
Defines the linear algebra or vector space interface.
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)
Teuchos::RCP< const vector > getVector(const V &x)
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.
int main(int argc, char *argv[])
void hessVec_true(ROL::Vector< Real > &hv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &z, Real &tol)
ROL::StdVector< Real > SV
Provides the interface to apply upper and lower bound constraints.
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)
void hessVec(ROL::Vector< Real > &hv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &z, Real &tol)
Apply Hessian approximation to vector.
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.
virtual std::vector< std::string > run(Vector< Real > &x, Objective< Real > &obj, bool print=false, std::ostream &outStream=std::cout)
Run algorithm on unconstrained problems (Type-U). This is the primary Type-U interface.
virtual Teuchos::RCP< Vector > basis(const int i) const
Return i-th basis vector.
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)
Teuchos::RCP< vector > getVector(V &x)
Teuchos::SerialDenseMatrix< int, Real > H_
static const double ROL_EPSILON
Platform-dependent machine epsilon.