51 #include "Teuchos_oblackholestream.hpp" 52 #include "Teuchos_GlobalMPISession.hpp" 53 #include "Teuchos_XMLParameterListHelpers.hpp" 54 #include "Teuchos_LAPACK.hpp" 77 return std::sqrt(this->
dot(r,r));
80 Real
dot(
const std::vector<Real> &x,
const std::vector<Real> &y) {
82 Real c = (((int)x.size()==this->
nx_) ? 4.0 : 2.0);
83 for (
unsigned i=0; i<x.size(); i++) {
85 ip += this->dx_/6.0*(c*x[i] + x[i+1])*y[i];
87 else if ( i == x.size()-1 ) {
88 ip += this->dx_/6.0*(x[i-1] + c*x[i])*y[i];
91 ip += this->dx_/6.0*(x[i-1] + 4.0*x[i] + x[i+1])*y[i];
97 void update(std::vector<Real> &u,
const std::vector<Real> &s,
const Real alpha=1.0) {
98 for (
unsigned i=0; i<u.size(); i++) {
103 void scale(std::vector<Real> &u,
const Real alpha=0.0) {
104 for (
unsigned i=0; i<u.size(); i++) {
110 const std::vector<Real> &z) {
112 r.resize(this->nx_,0.0);
113 for (
int i=0; i<this->
nx_; i++) {
116 r[i] = this->nu_/this->dx_*(2.0*u[i]-u[i+1]);
118 else if (i==this->nx_-1) {
119 r[i] = this->nu_/this->dx_*(2.0*u[i]-u[i-1]);
122 r[i] = this->nu_/this->dx_*(2.0*u[i]-u[i-1]-u[i+1]);
126 r[i] += u[i+1]*(u[i]+u[i+1])/6.0;
129 r[i] -= u[i-1]*(u[i-1]+u[i])/6.0;
132 r[i] -= this->dx_/6.0*(z[i]+4.0*z[i+1]+z[i+2]);
134 r[i] -= this->dx_*this->
f_;
137 r[0] -= this->u0_*u[ 0]/6.0 + this->u0_*this->u0_/6.0 + this->nu_*this->u0_/this->
dx_;
138 r[this->nx_-1] += this->u1_*u[this->nx_-1]/6.0 + this->u1_*this->u1_/6.0 - this->nu_*this->u1_/this->
dx_;
142 const std::vector<Real> &u) {
145 d.resize(this->nx_,this->nu_*2.0/this->dx_);
147 dl.resize(this->nx_-1,-this->nu_/this->dx_);
149 du.resize(this->nx_-1,-this->nu_/this->dx_);
151 for (
int i=0; i<this->
nx_; i++) {
153 dl[i] += (-2.0*u[i]-u[i+1])/6.0;
158 du[i-1] += (u[i-1]+2.0*u[i])/6.0;
162 d[0] -= this->u0_/6.0;
163 d[this->nx_-1] += this->u1_/6.0;
166 void linear_solve(std::vector<Real> &u, std::vector<Real> &dl, std::vector<Real> &d, std::vector<Real> &du,
167 const std::vector<Real> &r,
const bool transpose =
false) {
168 u.assign(r.begin(),r.end());
170 Teuchos::LAPACK<int,Real> lp;
171 std::vector<Real> du2(this->nx_-2,0.0);
172 std::vector<int> ipiv(this->nx_,0);
176 lp.GTTRF(this->nx_,&dl[0],&d[0],&du[0],&du2[0],&ipiv[0],&info);
181 lp.GTTRS(trans,this->nx_,nhrs,&dl[0],&d[0],&du[0],&du2[0],&ipiv[0],&u[0],ldb,&info);
187 : nx_(nx), nu_(nu), u0_(u0), u1_(u1), f_(f) {
188 dx_ = 1.0/((Real)nx+1.0);
193 Teuchos::RCP<std::vector<Real> > cp =
194 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<
ROL::StdVector<Real> >(c)).getVector());
195 Teuchos::RCP<const std::vector<Real> > up =
197 Teuchos::RCP<const std::vector<Real> > zp =
204 Teuchos::RCP<std::vector<Real> > jvp =
205 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<
ROL::StdVector<Real> >(jv)).getVector());
206 Teuchos::RCP<const std::vector<Real> > vp =
208 Teuchos::RCP<const std::vector<Real> > up =
210 Teuchos::RCP<const std::vector<Real> > zp =
213 for (
int i = 0; i < this->
nx_; i++) {
214 (*jvp)[i] = this->nu_/this->dx_*2.0*(*vp)[i];
216 (*jvp)[i] += -this->nu_/this->dx_*(*vp)[i-1]
217 -(*up)[i-1]/6.0*(*vp)[i]
218 -((*up)[i]+2.0*(*up)[i-1])/6.0*(*vp)[i-1];
220 if ( i < this->nx_-1 ) {
221 (*jvp)[i] += -this->nu_/this->dx_*(*vp)[i+1]
222 +(*up)[i+1]/6.0*(*vp)[i]
223 +((*up)[i]+2.0*(*up)[i+1])/6.0*(*vp)[i+1];
226 (*jvp)[0] -= this->u0_/6.0*(*vp)[0];
227 (*jvp)[this->nx_-1] += this->u1_/6.0*(*vp)[this->nx_-1];
232 Teuchos::RCP<std::vector<Real> > jvp =
233 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<
ROL::StdVector<Real> >(jv)).getVector());
234 Teuchos::RCP<const std::vector<Real> > vp =
236 Teuchos::RCP<const std::vector<Real> > up =
238 Teuchos::RCP<const std::vector<Real> > zp =
240 for (
int i=0; i<this->
nx_; i++) {
242 (*jvp)[i] = -this->dx_/6.0*((*vp)[i]+4.0*(*vp)[i+1]+(*vp)[i+2]);
248 Teuchos::RCP<std::vector<Real> > ijvp =
249 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<
ROL::StdVector<Real> >(ijv)).getVector());
250 Teuchos::RCP<const std::vector<Real> > vp =
252 Teuchos::RCP<const std::vector<Real> > up =
254 Teuchos::RCP<const std::vector<Real> > zp =
257 std::vector<Real> d(this->nx_,0.0);
258 std::vector<Real> dl(this->nx_-1,0.0);
259 std::vector<Real> du(this->nx_-1,0.0);
267 Teuchos::RCP<std::vector<Real> > jvp =
268 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<
ROL::StdVector<Real> >(ajv)).getVector());
269 Teuchos::RCP<const std::vector<Real> > vp =
271 Teuchos::RCP<const std::vector<Real> > up =
273 Teuchos::RCP<const std::vector<Real> > zp =
276 for (
int i = 0; i < this->
nx_; i++) {
277 (*jvp)[i] = this->nu_/this->dx_*2.0*(*vp)[i];
279 (*jvp)[i] += -this->nu_/this->dx_*(*vp)[i-1]
280 -(*up)[i-1]/6.0*(*vp)[i]
281 +((*up)[i-1]+2.0*(*up)[i])/6.0*(*vp)[i-1];
283 if ( i < this->nx_-1 ) {
284 (*jvp)[i] += -this->nu_/this->dx_*(*vp)[i+1]
285 +(*up)[i+1]/6.0*(*vp)[i]
286 -((*up)[i+1]+2.0*(*up)[i])/6.0*(*vp)[i+1];
289 (*jvp)[0] -= this->u0_/6.0*(*vp)[0];
290 (*jvp)[this->nx_-1] += this->u1_/6.0*(*vp)[this->nx_-1];
295 Teuchos::RCP<std::vector<Real> > jvp =
296 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<
ROL::StdVector<Real> >(jv)).getVector());
297 Teuchos::RCP<const std::vector<Real> > vp =
299 Teuchos::RCP<const std::vector<Real> > up =
301 Teuchos::RCP<const std::vector<Real> > zp =
303 for (
int i=0; i<this->nx_+2; i++) {
305 (*jvp)[i] = -this->dx_/6.0*(*vp)[i];
308 (*jvp)[i] = -this->dx_/6.0*(4.0*(*vp)[i-1]+(*vp)[i]);
310 else if ( i == this->nx_ ) {
311 (*jvp)[i] = -this->dx_/6.0*(4.0*(*vp)[i-1]+(*vp)[i-2]);
313 else if ( i == this->nx_+1 ) {
314 (*jvp)[i] = -this->dx_/6.0*(*vp)[i-2];
317 (*jvp)[i] = -this->dx_/6.0*((*vp)[i-2]+4.0*(*vp)[i-1]+(*vp)[i]);
324 Teuchos::RCP<std::vector<Real> > iajvp =
325 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<
ROL::StdVector<Real> >(iajv)).getVector());
326 Teuchos::RCP<const std::vector<Real> > vp =
328 Teuchos::RCP<const std::vector<Real> > up =
331 std::vector<Real> d(this->nx_,0.0);
332 std::vector<Real> du(this->nx_-1,0.0);
333 std::vector<Real> dl(this->nx_-1,0.0);
341 Teuchos::RCP<std::vector<Real> > ahwvp =
342 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<
ROL::StdVector<Real> >(ahwv)).getVector());
343 Teuchos::RCP<const std::vector<Real> > wp =
345 Teuchos::RCP<const std::vector<Real> > vp =
347 Teuchos::RCP<const std::vector<Real> > up =
349 Teuchos::RCP<const std::vector<Real> > zp =
351 for (
int i=0; i<this->
nx_; i++) {
355 (*ahwvp)[i] += ((*wp)[i]*(*vp)[i+1] - (*wp)[i+1]*(2.0*(*vp)[i]+(*vp)[i+1]))/6.0;
358 (*ahwvp)[i] += ((*wp)[i-1]*((*vp)[i-1]+2.0*(*vp)[i]) - (*wp)[i]*(*vp)[i-1])/6.0;
395 case 1: val = ((x<0.5) ? 1.0 : 0.0);
break;
396 case 2: val = 1.0;
break;
397 case 3: val = std::abs(std::sin(8.0*M_PI*x));
break;
398 case 4: val = std::exp(-0.5*(x-0.5)*(x-0.5));
break;
403 Real
dot(
const std::vector<Real> &x,
const std::vector<Real> &y) {
405 Real c = (((int)x.size()==this->
nx_) ? 4.0 : 2.0);
406 for (
unsigned i=0; i<x.size(); i++) {
408 ip += this->dx_/6.0*(c*x[i] + x[i+1])*y[i];
410 else if ( i == x.size()-1 ) {
411 ip += this->dx_/6.0*(x[i-1] + c*x[i])*y[i];
414 ip += this->dx_/6.0*(x[i-1] + 4.0*x[i] + x[i+1])*y[i];
420 void apply_mass(std::vector<Real> &Mu,
const std::vector<Real> &u ) {
421 Mu.resize(u.size(),0.0);
422 Real c = (((int)u.size()==this->
nx_) ? 4.0 : 2.0);
423 for (
unsigned i=0; i<u.size(); i++) {
425 Mu[i] = this->dx_/6.0*(c*u[i] + u[i+1]);
427 else if ( i == u.size()-1 ) {
428 Mu[i] = this->dx_/6.0*(u[i-1] + c*u[i]);
431 Mu[i] = this->dx_/6.0*(u[i-1] + 4.0*u[i] + u[i+1]);
442 dx_ = 1.0/((Real)nx+1.0);
446 Teuchos::RCP<const std::vector<Real> > up =
448 Teuchos::RCP<const std::vector<Real> > zp =
451 Real res1 = 0.0, res2 = 0.0, res3 = 0.0;
452 Real valu = 0.0, valz = this->
dot(*zp,*zp);
453 for (
int i=0; i<this->
nx_; i++) {
455 res1 = (*up)[i]-evaluate_target((Real)(i+1)*this->dx_);
456 res2 = (*up)[i+1]-evaluate_target((Real)(i+2)*this->dx_);
457 valu += this->dx_/6.0*(4.0*res1 + res2)*res1;
459 else if ( i == this->nx_-1 ) {
460 res1 = (*up)[i-1]-evaluate_target((Real)i*this->dx_);
461 res2 = (*up)[i]-evaluate_target((Real)(i+1)*this->dx_);
462 valu += this->dx_/6.0*(res1 + 4.0*res2)*res2;
465 res1 = (*up)[i-1]-evaluate_target((Real)i*this->dx_);
466 res2 = (*up)[i]-evaluate_target((Real)(i+1)*this->dx_);
467 res3 = (*up)[i+1]-evaluate_target((Real)(i+2)*this->dx_);
468 valu += this->dx_/6.0*(res1 + 4.0*res2 + res3)*res2;
471 return 0.5*(valu + this->alpha_*valz);
476 Teuchos::RCP<std::vector<Real> > gup = Teuchos::rcp_const_cast<std::vector<Real> >(
479 Teuchos::RCP<const std::vector<Real> > up =
481 Teuchos::RCP<const std::vector<Real> > zp =
484 std::vector<Real> diff(this->nx_,0.0);
485 for (
int i=0; i<this->
nx_; i++) {
486 diff[i] = ((*up)[i]-this->evaluate_target((Real)(i+1)*this->dx_));
488 this->apply_mass(*gup,diff);
493 Teuchos::RCP<std::vector<Real> > gzp = Teuchos::rcp_const_cast<std::vector<Real> >(
496 Teuchos::RCP<const std::vector<Real> > up =
498 Teuchos::RCP<const std::vector<Real> > zp =
501 for (
int i=0; i<this->nx_+2; i++) {
503 (*gzp)[i] = this->alpha_*this->dx_/6.0*(2.0*(*zp)[i]+(*zp)[i+1]);
505 else if (i==this->nx_+1) {
506 (*gzp)[i] = this->alpha_*this->dx_/6.0*(2.0*(*zp)[i]+(*zp)[i-1]);
509 (*gzp)[i] = this->alpha_*this->dx_/6.0*((*zp)[i-1]+4.0*(*zp)[i]+(*zp)[i+1]);
516 Teuchos::RCP<std::vector<Real> > hvup = Teuchos::rcp_const_cast<std::vector<Real> >(
519 Teuchos::RCP<const std::vector<Real> > vup =
522 this->apply_mass(*hvup,*vup);
537 Teuchos::RCP<std::vector<Real> > hvzp = Teuchos::rcp_const_cast<std::vector<Real> >(
540 Teuchos::RCP<const std::vector<Real> > vzp =
543 for (
int i=0; i<this->nx_+2; i++) {
545 (*hvzp)[i] = this->alpha_*this->dx_/6.0*(2.0*(*vzp)[i]+(*vzp)[i+1]);
547 else if (i==this->nx_+1) {
548 (*hvzp)[i] = this->alpha_*this->dx_/6.0*(2.0*(*vzp)[i]+(*vzp)[i-1]);
551 (*hvzp)[i] = this->alpha_*this->dx_/6.0*((*vzp)[i-1]+4.0*(*vzp)[i]+(*vzp)[i+1]);
Provides the interface to evaluate simulation-based objective functions.
void compute_pde_jacobian(std::vector< Real > &dl, std::vector< Real > &d, std::vector< Real > &du, const std::vector< Real > &u)
void applyAdjointHessian_21(ROL::Vector< Real > &ahwv, const ROL::Vector< Real > &w, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Apply the adjoint of the partial constraint Hessian at , , to vector in direction ...
void linear_solve(std::vector< Real > &u, std::vector< Real > &dl, std::vector< Real > &d, std::vector< Real > &du, const std::vector< Real > &r, const bool transpose=false)
Real evaluate_target(Real x)
void applyAdjointHessian_22(ROL::Vector< Real > &ahwv, const ROL::Vector< Real > &w, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Apply the adjoint of the partial constraint Hessian at , , to vector in direction ...
void applyInverseJacobian_1(ROL::Vector< Real > &ijv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Apply the inverse partial constraint Jacobian at , , to the vector .
Real dot(const std::vector< Real > &x, const std::vector< Real > &y)
Contains definitions of custom data types in ROL.
EqualityConstraint_BurgersControl(int nx=128, Real nu=1.e-2, Real u0=1.0, Real u1=0.0, Real f=0.0)
Real compute_norm(const std::vector< Real > &r)
virtual void zero()
Set to zero vector.
Defines the linear algebra or vector space interface.
void compute_residual(std::vector< Real > &r, const std::vector< Real > &u, const std::vector< Real > &z)
void applyAdjointHessian_11(ROL::Vector< Real > &ahwv, const ROL::Vector< Real > &w, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Apply the adjoint of the partial constraint Hessian at , , to vector in direction ...
Defines the equality constraint operator interface for simulation-based optimization.
void hessVec_22(ROL::Vector< Real > &hv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Real value(const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Compute value.
void gradient_1(ROL::Vector< Real > &g, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Compute gradient with respect to first component.
void applyInverseAdjointJacobian_1(ROL::Vector< Real > &iajv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Apply the inverse of the adjoint of the partial constraint Jacobian at , , to the vector ...
void applyAdjointJacobian_1(ROL::Vector< Real > &ajv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Apply the adjoint of the partial constraint Jacobian at , , to the vector . This is the primary inter...
void applyAdjointJacobian_2(ROL::Vector< Real > &jv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Apply the adjoint of the partial constraint Jacobian at , , to vector . This is the primary interface...
void applyJacobian_1(ROL::Vector< Real > &jv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Apply the partial constraint Jacobian at , , to the vector .
Real dot(const std::vector< Real > &x, const std::vector< Real > &y)
Objective_BurgersControl(Real alpha=1.e-4, int nx=128)
void hessVec_21(ROL::Vector< Real > &hv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
void hessVec_12(ROL::Vector< Real > &hv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
void gradient_2(ROL::Vector< Real > &g, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Compute gradient with respect to second component.
void scale(std::vector< Real > &u, const Real alpha=0.0)
void applyJacobian_2(ROL::Vector< Real > &jv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Apply the partial constraint Jacobian at , , to the vector .
void hessVec_11(ROL::Vector< Real > &hv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Apply Hessian approximation to vector.
void value(ROL::Vector< Real > &c, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Evaluate the constraint operator at .
void apply_mass(std::vector< Real > &Mu, const std::vector< Real > &u)
void update(std::vector< Real > &u, const std::vector< Real > &s, const Real alpha=1.0)
void applyAdjointHessian_12(ROL::Vector< Real > &ahwv, const ROL::Vector< Real > &w, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Apply the adjoint of the partial constraint Hessian at , , to vector in direction ...