57 #include "Teuchos_oblackholestream.hpp" 58 #include "Teuchos_GlobalMPISession.hpp" 70 typedef typename vector::size_type
uint;
86 Teuchos::SerialDenseMatrix<int, Real>
H_;
89 using Teuchos::dyn_cast;
94 using Teuchos::dyn_cast;
104 hu_ = 1.0/((Real)
nu_+1.0);
108 void apply_mass(std::vector<Real> &Mz,
const std::vector<Real> &z ) {
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] );
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]) );
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_;
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);
282 Real tol = std::sqrt(ROL::ROL_EPSILON<Real>());
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 ) {
299 correction = (1.0+std::sqrt(ROL::ROL_EPSILON<Real>()))*std::abs(inertia);
300 if ( inertia == 0.0 ) {
302 while ( eigenvals[0][cnt] == 0.0 ) {
305 correction = std::sqrt(ROL::ROL_EPSILON<Real>())*eigenvals[0][cnt];
306 if ( cnt ==
nz_-1 ) {
310 for (
uint i = 0; i <
nz_; i++ ) {
311 (
H_)(i,i) += correction;
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++ ) {
425 std::vector<Real> tmp(
nz_,0.0);
439 using Teuchos::rcp_const_cast;
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");
548 parlist.sublist(
"Step").sublist(
"Trust Region").set(
"Initial Radius", 1e3);
549 parlist.sublist(
"Step").sublist(
"Trust Region").set(
"Maximum Radius", 1e8);
554 algo_tr.
run(y,obj,icon,
true,*outStream);
556 std::ofstream file_tr;
557 file_tr.open(
"control_TR.txt");
558 for ( uint i = 0; i < dim; i++ ) {
559 file_tr << (*y_rcp)[i] <<
"\n";
563 std::vector<RealT> u;
565 std::ofstream file_u;
566 file_u.open(
"state.txt");
567 for ( uint i = 0; i < (dim-1); i++ ) {
568 file_u << u[i] <<
"\n";
572 RCP<V> diff = x.clone();
575 RealT error = diff->norm()/std::sqrt((
RealT)dim-1.0);
576 std::cout <<
"\nError between PDAS solution and TR solution is " << error <<
"\n";
577 errorFlag = ((error > 1.e2*std::sqrt(ROL::ROL_EPSILON<RealT>())) ? 1 : 0);
580 catch (std::logic_error err) {
581 *outStream << err.what() <<
"\n";
586 std::cout <<
"End Result: TEST FAILED\n";
588 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.
virtual Teuchos::RCP< Vector > basis(const int i) const
Return i-th basis vector.
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.
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_