21 #include "Teuchos_LAPACK.hpp"
49 void update(std::vector<Real> &u,
const std::vector<Real> &s,
const Real alpha=1.0)
const {
50 for (
unsigned i=0; i<u.size(); i++) {
55 void axpy(std::vector<Real> &out,
const Real a,
const std::vector<Real> &x,
const std::vector<Real> &y)
const {
56 for (
unsigned i=0; i < x.size(); i++) {
57 out[i] = a*x[i] + y[i];
61 void scale(std::vector<Real> &u,
const Real alpha=0.0)
const {
62 for (
unsigned i=0; i<u.size(); i++) {
67 void linear_solve(std::vector<Real> &u, std::vector<Real> &dl, std::vector<Real> &d, std::vector<Real> &du,
68 const std::vector<Real> &r,
const bool transpose =
false)
const {
69 if ( r.size() == 1 ) {
70 u.resize(1,r[0]/d[0]);
72 else if ( r.size() == 2 ) {
74 Real det = d[0]*d[1] - dl[0]*du[0];
75 u[0] = (d[1]*r[0] - du[0]*r[1])/det;
76 u[1] = (d[0]*r[1] - dl[0]*r[0])/det;
79 u.assign(r.begin(),r.end());
81 Teuchos::LAPACK<int,Real> lp;
82 std::vector<Real> du2(r.size()-2,0.0);
83 std::vector<int> ipiv(r.size(),0);
88 lp.GTTRF(dim,&dl[0],&d[0],&du[0],&du2[0],&ipiv[0],&info);
93 lp.GTTRS(trans,dim,nhrs,&dl[0],&d[0],&du[0],&du2[0],&ipiv[0],&u[0],ldb,&info);
99 Real u0 = 1.0, Real u1 = 0.0, Real f = 0.0,
100 Real cH1 = 1.0, Real cL2 = 1.0)
101 :
nx_(nx),
dx_(1.0/((Real)nx+1.0)),
117 Real
compute_L2_dot(
const std::vector<Real> &x,
const std::vector<Real> &y)
const {
119 Real c = (((int)x.size()==
nx_) ? 4.0 : 2.0);
120 for (
unsigned i=0; i<x.size(); i++) {
122 ip +=
dx_/6.0*(c*x[i] + x[i+1])*y[i];
124 else if ( i == x.size()-1 ) {
125 ip +=
dx_/6.0*(x[i-1] + c*x[i])*y[i];
128 ip +=
dx_/6.0*(x[i-1] + 4.0*x[i] + x[i+1])*y[i];
140 void apply_mass(std::vector<Real> &Mu,
const std::vector<Real> &u )
const {
141 Mu.resize(u.size(),0.0);
142 Real c = (((int)u.size()==
nx_) ? 4.0 : 2.0);
143 for (
unsigned i=0; i<u.size(); i++) {
145 Mu[i] =
dx_/6.0*(c*u[i] + u[i+1]);
147 else if ( i == u.size()-1 ) {
148 Mu[i] =
dx_/6.0*(u[i-1] + c*u[i]);
151 Mu[i] =
dx_/6.0*(u[i-1] + 4.0*u[i] + u[i+1]);
158 unsigned nx = u.size();
160 std::vector<Real> dl(nx-1,
dx_/6.0);
161 std::vector<Real> d(nx,2.0*
dx_/3.0);
162 std::vector<Real> du(nx-1,
dx_/6.0);
163 if ( (
int)nx !=
nx_ ) {
172 std::vector<Real> u(
nx_,0.0), Mu(
nx_,0.0), iMMu(
nx_,0.0), diff(
nx_,0.0);
173 for (
int i = 0; i <
nx_; i++) {
174 u[i] = 2.0*(Real)rand()/(Real)RAND_MAX - 1.0;
178 axpy(diff,-1.0,iMMu,u);
181 outStream <<
"Test Inverse State Mass Matrix\n";
182 outStream <<
" ||u - inv(M)Mu|| = " << error <<
"\n";
183 outStream <<
" ||u|| = " << normu <<
"\n";
184 outStream <<
" Relative Error = " << error/normu <<
"\n";
187 u.resize(nx_+2,0.0); Mu.resize(nx_+2,0.0); iMMu.resize(nx_+2,0.0); diff.resize(nx_+2,0.0);
188 for (
int i = 0; i < nx_+2; i++) {
189 u[i] = 2.0*(Real)rand()/(Real)RAND_MAX - 1.0;
193 axpy(diff,-1.0,iMMu,u);
196 outStream <<
"Test Inverse Control Mass Matrix\n";
197 outStream <<
" ||z - inv(M)Mz|| = " << error <<
"\n";
198 outStream <<
" ||z|| = " << normu <<
"\n";
199 outStream <<
" Relative Error = " << error/normu <<
"\n";
207 Real
compute_H1_dot(
const std::vector<Real> &x,
const std::vector<Real> &y)
const {
209 for (
int i=0; i<
nx_; i++) {
211 ip +=
cL2_*
dx_/6.0*(4.0*x[i] + x[i+1])*y[i];
212 ip +=
cH1_*(2.0*x[i] - x[i+1])/
dx_*y[i];
214 else if ( i == nx_-1 ) {
215 ip +=
cL2_*
dx_/6.0*(x[i-1] + 4.0*x[i])*y[i];
216 ip +=
cH1_*(2.0*x[i] - x[i-1])/
dx_*y[i];
219 ip +=
cL2_*
dx_/6.0*(x[i-1] + 4.0*x[i] + x[i+1])*y[i];
220 ip +=
cH1_*(2.0*x[i] - x[i-1] - x[i+1])/
dx_*y[i];
232 void apply_H1(std::vector<Real> &Mu,
const std::vector<Real> &u )
const {
234 for (
int i=0; i<
nx_; i++) {
236 Mu[i] =
cL2_*
dx_/6.0*(4.0*u[i] + u[i+1])
237 +
cH1_*(2.0*u[i] - u[i+1])/
dx_;
239 else if ( i == nx_-1 ) {
240 Mu[i] =
cL2_*
dx_/6.0*(u[i-1] + 4.0*u[i])
241 +
cH1_*(2.0*u[i] - u[i-1])/
dx_;
244 Mu[i] =
cL2_*
dx_/6.0*(u[i-1] + 4.0*u[i] + u[i+1])
245 +
cH1_*(2.0*u[i] - u[i-1] - u[i+1])/
dx_;
260 std::vector<Real> u(
nx_,0.0), Mu(
nx_,0.0), iMMu(
nx_,0.0), diff(
nx_,0.0);
261 for (
int i = 0; i <
nx_; i++) {
262 u[i] = 2.0*(Real)rand()/(Real)RAND_MAX - 1.0;
266 axpy(diff,-1.0,iMMu,u);
269 outStream <<
"Test Inverse State H1 Matrix\n";
270 outStream <<
" ||u - inv(M)Mu|| = " << error <<
"\n";
271 outStream <<
" ||u|| = " << normu <<
"\n";
272 outStream <<
" Relative Error = " << error/normu <<
"\n";
281 const std::vector<Real> &z)
const {
284 for (
int i=0; i<
nx_; i++) {
287 r[i] =
nu_/
dx_*(2.0*u[i]-u[i+1]);
290 r[i] =
nu_/
dx_*(2.0*u[i]-u[i-1]);
293 r[i] =
nu_/
dx_*(2.0*u[i]-u[i-1]-u[i+1]);
297 r[i] +=
nl_*u[i+1]*(u[i]+u[i+1])/6.0;
300 r[i] -=
nl_*u[i-1]*(u[i-1]+u[i])/6.0;
303 r[i] -=
dx_/6.0*(z[i]+4.0*z[i+1]+z[i+2]);
317 std::vector<Real> &d,
318 std::vector<Real> &du,
319 const std::vector<Real> &u)
const {
328 for (
int i=0; i<
nx_; i++) {
330 dl[i] +=
nl_*(-2.0*u[i]-u[i+1])/6.0;
331 d[i] +=
nl_*u[i+1]/6.0;
334 d[i] -=
nl_*u[i-1]/6.0;
335 du[i-1] +=
nl_*(u[i-1]+2.0*u[i])/6.0;
345 const std::vector<Real> &v,
346 const std::vector<Real> &u,
347 const std::vector<Real> &z)
const {
349 for (
int i = 0; i <
nx_; i++) {
352 jv[i] += -
nu_/
dx_*v[i-1]-
nl_*(u[i-1]/6.0*v[i]+(u[i]+2.0*u[i-1])/6.0*v[i-1]);
355 jv[i] += -
nu_/
dx_*v[i+1]+
nl_*(u[i+1]/6.0*v[i]+(u[i]+2.0*u[i+1])/6.0*v[i+1]);
358 jv[ 0] -=
nl_*
u0_/6.0*v[0];
359 jv[nx_-1] +=
nl_*
u1_/6.0*v[nx_-1];
364 const std::vector<Real> &v,
365 const std::vector<Real> &u,
366 const std::vector<Real> &z)
const {
368 std::vector<Real> d(
nx_,0.0);
369 std::vector<Real> dl(
nx_-1,0.0);
370 std::vector<Real> du(
nx_-1,0.0);
378 const std::vector<Real> &v,
379 const std::vector<Real> &u,
380 const std::vector<Real> &z)
const {
382 for (
int i = 0; i <
nx_; i++) {
383 ajv[i] =
nu_/
dx_*2.0*v[i];
385 ajv[i] += -
nu_/
dx_*v[i-1]-
nl_*(u[i-1]/6.0*v[i]
386 -(u[i-1]+2.0*u[i])/6.0*v[i-1]);
389 ajv[i] += -
nu_/
dx_*v[i+1]+
nl_*(u[i+1]/6.0*v[i]
390 -(u[i+1]+2.0*u[i])/6.0*v[i+1]);
393 ajv[ 0] -=
nl_*
u0_/6.0*v[0];
394 ajv[nx_-1] +=
nl_*
u1_/6.0*v[nx_-1];
399 const std::vector<Real> &v,
400 const std::vector<Real> &u,
401 const std::vector<Real> &z)
const {
403 std::vector<Real> d(
nx_,0.0);
404 std::vector<Real> du(
nx_-1,0.0);
405 std::vector<Real> dl(
nx_-1,0.0);
416 const std::vector<Real> &v,
417 const std::vector<Real> &u,
418 const std::vector<Real> &z)
const {
419 for (
int i=0; i<
nx_; i++) {
421 jv[i] = -
dx_/6.0*(v[i]+4.0*v[i+1]+v[i+2]);
427 const std::vector<Real> &v,
428 const std::vector<Real> &u,
429 const std::vector<Real> &z)
const {
430 for (
int i=0; i<
nx_+2; i++) {
432 jv[i] = -
dx_/6.0*v[i];
435 jv[i] = -
dx_/6.0*(4.0*v[i-1]+v[i]);
437 else if ( i == nx_ ) {
438 jv[i] = -
dx_/6.0*(4.0*v[i-1]+v[i-2]);
440 else if ( i == nx_+1 ) {
441 jv[i] = -
dx_/6.0*v[i-2];
444 jv[i] = -
dx_/6.0*(v[i-2]+4.0*v[i-1]+v[i]);
453 const std::vector<Real> &w,
454 const std::vector<Real> &v,
455 const std::vector<Real> &u,
456 const std::vector<Real> &z)
const {
457 for (
int i=0; i<
nx_; i++) {
461 ahwv[i] += (w[i]*v[i+1] - w[i+1]*(2.0*v[i]+v[i+1]))/6.0;
464 ahwv[i] += (w[i-1]*(v[i-1]+2.0*v[i]) - w[i]*v[i-1])/6.0;
470 const std::vector<Real> &w,
471 const std::vector<Real> &v,
472 const std::vector<Real> &u,
473 const std::vector<Real> &z) {
474 ahwv.assign(u.size(),0.0);
477 const std::vector<Real> &w,
478 const std::vector<Real> &v,
479 const std::vector<Real> &u,
480 const std::vector<Real> &z) {
481 ahwv.assign(z.size(),0.0);
484 const std::vector<Real> &w,
485 const std::vector<Real> &v,
486 const std::vector<Real> &u,
487 const std::vector<Real> &z) {
488 ahwv.assign(z.size(),0.0);
495 ROL::Ptr<std::vector<Real> >
vec_;
496 ROL::Ptr<BurgersFEM<Real> >
fem_;
498 mutable ROL::Ptr<L2VectorDual<Real> >
dual_vec_;
507 const std::vector<Real>& xval = *ex.
getVector();
508 std::copy(xval.begin(),xval.end(),
vec_->begin());
513 const std::vector<Real>& xval = *ex.
getVector();
516 (*vec_)[i] += xval[i];
529 const std::vector<Real>& xval = *ex.
getVector();
530 return fem_->compute_L2_dot(xval,*
vec_);
535 val = std::sqrt(
dot(*
this) );
539 ROL::Ptr<ROL::Vector<Real> >
clone()
const {
540 return ROL::makePtr<L2VectorPrimal>( ROL::makePtr<std::vector<Real>>(
vec_->size(),0.0),
fem_);
551 ROL::Ptr<ROL::Vector<Real> >
basis(
const int i )
const {
552 ROL::Ptr<L2VectorPrimal> e
553 = ROL::makePtr<L2VectorPrimal>( ROL::makePtr<std::vector<Real>>(
vec_->size(),0.0),
fem_);
554 (*e->getVector())[i] = 1.0;
563 dual_vec_ = ROL::makePtr<L2VectorDual<Real>>(
564 ROL::makePtr<std::vector<Real>>(*vec_),
fem_);
566 fem_->apply_mass(*(ROL::constPtrCast<std::vector<Real> >(
dual_vec_->getVector())),*
vec_);
572 const std::vector<Real>& xval = *ex.
getVector();
573 return std::inner_product(
vec_->begin(),
vec_->end(), xval.begin(), Real(0));
581 ROL::Ptr<std::vector<Real> >
vec_;
582 ROL::Ptr<BurgersFEM<Real> >
fem_;
584 mutable ROL::Ptr<L2VectorPrimal<Real> >
dual_vec_;
593 const std::vector<Real>& xval = *ex.
getVector();
594 std::copy(xval.begin(),xval.end(),
vec_->begin());
599 const std::vector<Real>& xval = *ex.
getVector();
602 (*vec_)[i] += xval[i];
615 const std::vector<Real>& xval = *ex.
getVector();
617 std::vector<Real> Mx(dimension,0.0);
618 fem_->apply_inverse_mass(Mx,xval);
620 for (
unsigned i = 0; i <
dimension; i++) {
621 val += Mx[i]*(*vec_)[i];
628 val = std::sqrt(
dot(*
this) );
632 ROL::Ptr<ROL::Vector<Real> >
clone()
const {
633 return ROL::makePtr<L2VectorDual>( ROL::makePtr<std::vector<Real>>(
vec_->size(),0.0),
fem_);
644 ROL::Ptr<ROL::Vector<Real> >
basis(
const int i )
const {
645 ROL::Ptr<L2VectorDual> e
646 = ROL::makePtr<L2VectorDual>( ROL::makePtr<std::vector<Real>>(
vec_->size(),0.0),
fem_);
647 (*e->getVector())[i] = 1.0;
656 dual_vec_ = ROL::makePtr<L2VectorPrimal<Real>>(
657 ROL::makePtr<std::vector<Real>>(*vec_),
fem_);
659 fem_->apply_inverse_mass(*(ROL::constPtrCast<std::vector<Real> >(
dual_vec_->getVector())),*
vec_);
665 const std::vector<Real>& xval = *ex.
getVector();
666 return std::inner_product(
vec_->begin(),
vec_->end(), xval.begin(), Real(0));
674 ROL::Ptr<std::vector<Real> >
vec_;
675 ROL::Ptr<BurgersFEM<Real> >
fem_;
677 mutable ROL::Ptr<H1VectorDual<Real> >
dual_vec_;
686 const std::vector<Real>& xval = *ex.
getVector();
687 std::copy(xval.begin(),xval.end(),
vec_->begin());
692 const std::vector<Real>& xval = *ex.
getVector();
695 (*vec_)[i] += xval[i];
708 const std::vector<Real>& xval = *ex.
getVector();
709 return fem_->compute_H1_dot(xval,*
vec_);
714 val = std::sqrt(
dot(*
this) );
718 ROL::Ptr<ROL::Vector<Real> >
clone()
const {
719 return ROL::makePtr<H1VectorPrimal>( ROL::makePtr<std::vector<Real>>(
vec_->size(),0.0),
fem_);
730 ROL::Ptr<ROL::Vector<Real> >
basis(
const int i )
const {
731 ROL::Ptr<H1VectorPrimal> e
732 = ROL::makePtr<H1VectorPrimal>( ROL::makePtr<std::vector<Real>>(
vec_->size(),0.0),
fem_);
733 (*e->getVector())[i] = 1.0;
742 dual_vec_ = ROL::makePtr<H1VectorDual<Real>>(
743 ROL::makePtr<std::vector<Real>>(*vec_),
fem_);
745 fem_->apply_H1(*(ROL::constPtrCast<std::vector<Real> >(
dual_vec_->getVector())),*
vec_);
751 const std::vector<Real>& xval = *ex.
getVector();
752 return std::inner_product(
vec_->begin(),
vec_->end(), xval.begin(), Real(0));
760 ROL::Ptr<std::vector<Real> >
vec_;
761 ROL::Ptr<BurgersFEM<Real> >
fem_;
763 mutable ROL::Ptr<H1VectorPrimal<Real> >
dual_vec_;
772 const std::vector<Real>& xval = *ex.
getVector();
773 std::copy(xval.begin(),xval.end(),
vec_->begin());
778 const std::vector<Real>& xval = *ex.
getVector();
781 (*vec_)[i] += xval[i];
794 const std::vector<Real>& xval = *ex.
getVector();
796 std::vector<Real> Mx(dimension,0.0);
797 fem_->apply_inverse_H1(Mx,xval);
799 for (
unsigned i = 0; i <
dimension; i++) {
800 val += Mx[i]*(*vec_)[i];
807 val = std::sqrt(
dot(*
this) );
811 ROL::Ptr<ROL::Vector<Real> >
clone()
const {
812 return ROL::makePtr<H1VectorDual>( ROL::makePtr<std::vector<Real>>(
vec_->size(),0.0),
fem_);
823 ROL::Ptr<ROL::Vector<Real> >
basis(
const int i )
const {
824 ROL::Ptr<H1VectorDual> e
825 = ROL::makePtr<H1VectorDual>( ROL::makePtr<std::vector<Real>>(
vec_->size(),0.0),
fem_);
826 (*e->getVector())[i] = 1.0;
835 dual_vec_ = ROL::makePtr<H1VectorPrimal<Real>>(
836 ROL::makePtr<std::vector<Real>>(*vec_),
fem_);
838 fem_->apply_inverse_H1(*(ROL::constPtrCast<std::vector<Real> >(
dual_vec_->getVector())),*
vec_);
844 const std::vector<Real>& xval = *ex.
getVector();
845 return std::inner_product(
vec_->begin(),
vec_->end(), xval.begin(), Real(0));
863 ROL::Ptr<BurgersFEM<Real> >
fem_;
872 ROL::Ptr<std::vector<Real> > cp =
874 ROL::Ptr<const std::vector<Real> > up =
876 ROL::Ptr<const std::vector<Real> > zp =
878 fem_->compute_residual(*cp,*up,*zp);
885 ROL::Ptr<std::vector<Real> > jvp =
887 ROL::Ptr<const std::vector<Real> > vp =
889 ROL::Ptr<const std::vector<Real> > up =
891 ROL::Ptr<const std::vector<Real> > zp =
893 fem_->apply_pde_jacobian(*jvp,*vp,*up,*zp);
898 ROL::Ptr<std::vector<Real> > jvp =
900 ROL::Ptr<const std::vector<Real> > vp =
902 ROL::Ptr<const std::vector<Real> > up =
904 ROL::Ptr<const std::vector<Real> > zp =
906 fem_->apply_control_jacobian(*jvp,*vp,*up,*zp);
911 ROL::Ptr<std::vector<Real> > ijvp =
913 ROL::Ptr<const std::vector<Real> > vp =
915 ROL::Ptr<const std::vector<Real> > up =
917 ROL::Ptr<const std::vector<Real> > zp =
919 fem_->apply_inverse_pde_jacobian(*ijvp,*vp,*up,*zp);
924 ROL::Ptr<std::vector<Real> > jvp =
926 ROL::Ptr<const std::vector<Real> > vp =
928 ROL::Ptr<const std::vector<Real> > up =
930 ROL::Ptr<const std::vector<Real> > zp =
932 fem_->apply_adjoint_pde_jacobian(*jvp,*vp,*up,*zp);
937 ROL::Ptr<std::vector<Real> > jvp =
939 ROL::Ptr<const std::vector<Real> > vp =
941 ROL::Ptr<const std::vector<Real> > up =
943 ROL::Ptr<const std::vector<Real> > zp =
945 fem_->apply_adjoint_control_jacobian(*jvp,*vp,*up,*zp);
950 ROL::Ptr<std::vector<Real> > iajvp =
952 ROL::Ptr<const std::vector<Real> > vp =
954 ROL::Ptr<const std::vector<Real> > up =
956 ROL::Ptr<const std::vector<Real> > zp =
958 fem_->apply_inverse_adjoint_pde_jacobian(*iajvp,*vp,*up,*zp);
964 ROL::Ptr<std::vector<Real> > ahwvp =
966 ROL::Ptr<const std::vector<Real> > wp =
968 ROL::Ptr<const std::vector<Real> > vp =
970 ROL::Ptr<const std::vector<Real> > up =
972 ROL::Ptr<const std::vector<Real> > zp =
974 fem_->apply_adjoint_pde_hessian(*ahwvp,*wp,*vp,*up,*zp);
984 ROL::Ptr<std::vector<Real> > ahwvp =
986 ROL::Ptr<const std::vector<Real> > wp =
988 ROL::Ptr<const std::vector<Real> > vp =
990 ROL::Ptr<const std::vector<Real> > up =
992 ROL::Ptr<const std::vector<Real> > zp =
994 fem_->apply_adjoint_control_pde_hessian(*ahwvp,*wp,*vp,*up,*zp);
1003 ROL::Ptr<std::vector<Real> > ahwvp =
1005 ROL::Ptr<const std::vector<Real> > wp =
1007 ROL::Ptr<const std::vector<Real> > vp =
1009 ROL::Ptr<const std::vector<Real> > up =
1011 ROL::Ptr<const std::vector<Real> > zp =
1013 fem_->apply_adjoint_pde_control_hessian(*ahwvp,*wp,*vp,*up,*zp);
1022 ROL::Ptr<std::vector<Real> > ahwvp =
1024 ROL::Ptr<const std::vector<Real> > wp =
1026 ROL::Ptr<const std::vector<Real> > vp =
1028 ROL::Ptr<const std::vector<Real> > up =
1030 ROL::Ptr<const std::vector<Real> > zp =
1032 fem_->apply_adjoint_control_hessian(*ahwvp,*wp,*vp,*up,*zp);
1040 template<
class Real>
1052 ROL::Ptr<ROL::Vector<Real> >
ud_;
1063 ROL::Ptr<const std::vector<Real> > up =
1065 ROL::Ptr<const std::vector<Real> > zp =
1067 ROL::Ptr<const std::vector<Real> > udp =
1070 std::vector<Real> diff(udp->size(),0.0);
1071 for (
unsigned i = 0; i < udp->size(); i++) {
1072 diff[i] = (*up)[i] - (*udp)[i];
1074 return 0.5*(
fem_->compute_L2_dot(diff,diff) +
alpha_*
fem_->compute_L2_dot(*zp,*zp));
1078 ROL::Ptr<std::vector<Real> > gp =
1080 ROL::Ptr<const std::vector<Real> > up =
1082 ROL::Ptr<const std::vector<Real> > udp =
1085 std::vector<Real> diff(udp->size(),0.0);
1086 for (
unsigned i = 0; i < udp->size(); i++) {
1087 diff[i] = (*up)[i] - (*udp)[i];
1089 fem_->apply_mass(*gp,diff);
1093 ROL::Ptr<std::vector<Real> > gp =
1095 ROL::Ptr<const std::vector<Real> > zp =
1098 fem_->apply_mass(*gp,*zp);
1099 for (
unsigned i = 0; i < zp->size(); i++) {
1106 ROL::Ptr<std::vector<Real> > hvp =
1108 ROL::Ptr<const std::vector<Real> > vp =
1111 fem_->apply_mass(*hvp,*vp);
1126 ROL::Ptr<std::vector<Real> > hvp =
1128 ROL::Ptr<const std::vector<Real> > vp =
1131 fem_->apply_mass(*hvp,*vp);
1132 for (
unsigned i = 0; i < vp->size(); i++) {
1138 template<
class Real>
1147 ROL::Ptr<ROL::Vector<Real> >
l_;
1148 ROL::Ptr<ROL::Vector<Real> >
u_;
1155 catch (std::exception &e) {
1165 catch (std::exception &e) {
1170 void axpy(std::vector<Real> &out,
const Real a,
1171 const std::vector<Real> &x,
const std::vector<Real> &y)
const{
1172 out.resize(
dim_,0.0);
1173 for (
unsigned i = 0; i <
dim_; i++) {
1174 out[i] = a*x[i] + y[i];
1179 for (
int i = 0; i <
dim_; i++ ) {
1180 x[i] = std::max(
x_lo_[i],std::min(
x_up_[i],x[i]));
1189 for (
int i = 0; i <
dim_; i++ ) {
1198 l_ = ROL::makePtr<L2VectorPrimal<Real>>(
1199 ROL::makePtr<std::vector<Real>>(l), fem);
1200 u_ = ROL::makePtr<L2VectorPrimal<Real>>(
1201 ROL::makePtr<std::vector<Real>>(u), fem);
1208 for (
int i = 0; i <
dim_; i++ ) {
1209 if ( (*ex)[i] >=
x_lo_[i] && (*ex)[i] <=
x_up_[i] ) { cnt *= 1; }
1212 if ( cnt == 0 ) { val =
false; }
1217 ROL::Ptr<std::vector<Real> > ex;
cast_vector(ex,x);
1223 ROL::Ptr<std::vector<Real> > ev;
cast_vector(ev,v);
1225 for (
int i = 0; i <
dim_; i++ ) {
1226 if ( ((*ex)[i] <=
x_lo_[i]+epsn) ) {
1234 ROL::Ptr<std::vector<Real> > ev;
cast_vector(ev,v);
1236 for (
int i = 0; i <
dim_; i++ ) {
1237 if ( ((*ex)[i] >=
x_up_[i]-epsn) ) {
1245 ROL::Ptr<std::vector<Real> > ev;
cast_vector(ev,v);
1247 for (
int i = 0; i <
dim_; i++ ) {
1248 if ( ((*ex)[i] <=
x_lo_[i]+epsn) ||
1249 ((*ex)[i] >=
x_up_[i]-epsn) ) {
1258 ROL::Ptr<std::vector<Real> > ev;
cast_vector(ev,v);
1260 for (
int i = 0; i <
dim_; i++ ) {
1261 if ( ((*ex)[i] <=
x_lo_[i]+epsn && (*eg)[i] > geps) ) {
1270 ROL::Ptr<std::vector<Real> > ev;
cast_vector(ev,v);
1272 for (
int i = 0; i <
dim_; i++ ) {
1273 if ( ((*ex)[i] >=
x_up_[i]-epsn && (*eg)[i] < geps) ) {
1282 ROL::Ptr<std::vector<Real> > ev;
cast_vector(ev,v);
1284 for (
int i = 0; i <
dim_; i++ ) {
1285 if ( ((*ex)[i] <=
x_lo_[i]+epsn && (*eg)[i] > geps) ||
1286 ((*ex)[i] >=
x_up_[i]-epsn && (*eg)[i] < -geps) ) {
1301 template<
class Real>
1310 ROL::Ptr<ROL::Vector<Real> >
l_;
1311 ROL::Ptr<ROL::Vector<Real> >
u_;
1318 catch (std::exception &e) {
1328 catch (std::exception &e) {
1333 void axpy(std::vector<Real> &out,
const Real a,
1334 const std::vector<Real> &x,
const std::vector<Real> &y)
const{
1335 out.resize(
dim_,0.0);
1336 for (
unsigned i = 0; i <
dim_; i++) {
1337 out[i] = a*x[i] + y[i];
1342 for (
int i = 0; i <
dim_; i++ ) {
1343 x[i] = std::max(
x_lo_[i],std::min(
x_up_[i],x[i]));
1352 for (
int i = 0; i <
dim_; i++ ) {
1361 l_ = ROL::makePtr<H1VectorPrimal<Real>>(
1362 ROL::makePtr<std::vector<Real>>(l), fem);
1363 u_ = ROL::makePtr<H1VectorPrimal<Real>>(
1364 ROL::makePtr<std::vector<Real>>(u), fem);
1371 for (
int i = 0; i <
dim_; i++ ) {
1372 if ( (*ex)[i] >=
x_lo_[i] && (*ex)[i] <=
x_up_[i] ) { cnt *= 1; }
1375 if ( cnt == 0 ) { val =
false; }
1380 ROL::Ptr<std::vector<Real> > ex;
cast_vector(ex,x);
1386 ROL::Ptr<std::vector<Real> > ev;
cast_vector(ev,v);
1388 for (
int i = 0; i <
dim_; i++ ) {
1389 if ( ((*ex)[i] <=
x_lo_[i]+epsn) ) {
1397 ROL::Ptr<std::vector<Real> > ev;
cast_vector(ev,v);
1399 for (
int i = 0; i <
dim_; i++ ) {
1400 if ( ((*ex)[i] >=
x_up_[i]-epsn) ) {
1408 ROL::Ptr<std::vector<Real> > ev;
cast_vector(ev,v);
1410 for (
int i = 0; i <
dim_; i++ ) {
1411 if ( ((*ex)[i] <=
x_lo_[i]+epsn) ||
1412 ((*ex)[i] >=
x_up_[i]-epsn) ) {
1421 ROL::Ptr<std::vector<Real> > ev;
cast_vector(ev,v);
1423 for (
int i = 0; i <
dim_; i++ ) {
1424 if ( ((*ex)[i] <=
x_lo_[i]+epsn && (*eg)[i] > geps) ) {
1433 ROL::Ptr<std::vector<Real> > ev;
cast_vector(ev,v);
1435 for (
int i = 0; i <
dim_; i++ ) {
1436 if ( ((*ex)[i] >=
x_up_[i]-epsn && (*eg)[i] < -geps) ) {
1445 ROL::Ptr<std::vector<Real> > ev;
cast_vector(ev,v);
1447 for (
int i = 0; i <
dim_; i++ ) {
1448 if ( ((*ex)[i] <=
x_lo_[i]+epsn && (*eg)[i] > geps) ||
1449 ((*ex)[i] >=
x_up_[i]-epsn && (*eg)[i] < -geps) ) {
H1VectorPrimal< Real > DualConstraintVector
L2VectorDual(const ROL::Ptr< std::vector< Real > > &vec, const ROL::Ptr< BurgersFEM< Real > > &fem)
ROL::Ptr< std::vector< Real > > vec_
Provides the interface to evaluate simulation-based objective functions.
std::vector< Real > x_up_
ROL::Ptr< ROL::Vector< Real > > diff_
Real compute_H1_dot(const std::vector< Real > &x, const std::vector< Real > &y) const
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...
Real norm() const
Returns where .
void apply_mass(std::vector< Real > &Mu, const std::vector< Real > &u) const
bool isFeasible(const ROL::Vector< Real > &x)
Check if the vector, v, is feasible.
void pruneActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &g, const ROL::Vector< Real > &x, Real xeps=Real(0), Real geps=Real(0))
void axpy(std::vector< Real > &out, const Real a, const std::vector< Real > &x, const std::vector< Real > &y) const
ROL::Ptr< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
void cast_vector(ROL::Ptr< std::vector< Real > > &xvec, ROL::Vector< Real > &x) const
Real dot(const ROL::Vector< Real > &x) const override
Compute where .
Real dot(const ROL::Vector< Real > &x) const
Compute where .
void apply_adjoint_pde_control_hessian(std::vector< Real > &ahwv, const std::vector< Real > &w, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z)
ROL::Ptr< std::vector< Real > > vec_
int dimension() const
Return dimension of the vector space.
ROL::Ptr< L2VectorDual< Real > > dual_vec_
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 optimization-space derivative of the adjoint of the constraint simulation-space Jacobian at...
void plus(const ROL::Vector< Real > &x)
Compute , where .
std::vector< Real > x_up_
void pruneUpperActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &g, const ROL::Vector< Real > &x, Real xeps=Real(0), Real geps=Real(0))
Set variables to zero if they correspond to the upper -binding set.
void apply_adjoint_pde_hessian(std::vector< Real > &ahwv, const std::vector< Real > &w, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z) const
ROL::Ptr< BurgersFEM< Real > > fem_
H1VectorPrimal(const ROL::Ptr< std::vector< Real > > &vec, const ROL::Ptr< BurgersFEM< Real > > &fem)
Real norm() const
Returns where .
Contains definitions of custom data types in ROL.
Real dot(const ROL::Vector< Real > &x) const
Compute where .
void cast_vector(ROL::Ptr< std::vector< Real > > &xvec, ROL::Vector< Real > &x) const
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 simulation-space derivative of the adjoint of the constraint simulation-space Jacobian at ...
ROL::Ptr< const std::vector< Real > > getVector() const
const ROL::Ptr< const ROL::Vector< Real > > getUpperBound(void) const
Return the ref count pointer to the upper bound vector.
void plus(const ROL::Vector< Real > &x)
Compute , where .
ROL::Ptr< BurgersFEM< Real > > fem_
Real compute_H1_norm(const std::vector< Real > &r) const
L2BoundConstraint(std::vector< Real > &l, std::vector< Real > &u, const ROL::Ptr< BurgersFEM< Real > > &fem, Real scale=1.0)
Real norm() const
Returns where .
Real apply(const ROL::Vector< Real > &x) const
Apply to a dual vector. This is equivalent to the call .
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 ...
ROL::Ptr< std::vector< Real > > vec_
void pruneLowerActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &x, Real eps=Real(0))
Set variables to zero if they correspond to the lower -active set.
H1VectorDual(const ROL::Ptr< std::vector< Real > > &vec, const ROL::Ptr< BurgersFEM< Real > > &fem)
void project(ROL::Vector< Real > &x)
Project optimization variables onto the bounds.
virtual void zero()
Set to zero vector.
ROL::Ptr< std::vector< Real > > getVector()
Defines the linear algebra or vector space interface.
ROL::Ptr< BurgersFEM< Real > > fem_
H1VectorDual< Real > DualStateVector
H1VectorDual< Real > PrimalConstraintVector
void pruneLowerActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &x, Real eps=Real(0))
Set variables to zero if they correspond to the lower -active set.
Real apply(const ROL::Vector< Real > &x) const
Apply to a dual vector. This is equivalent to the call .
void scale(const Real alpha)
Compute where .
void pruneLowerActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &g, const ROL::Vector< Real > &x, Real xeps=Real(0), Real geps=Real(0))
Set variables to zero if they correspond to the -binding set.
L2VectorPrimal(const ROL::Ptr< std::vector< Real > > &vec, const ROL::Ptr< BurgersFEM< Real > > &fem)
Real dot(const ROL::Vector< Real > &x) const override
Compute where .
Real dot(const ROL::Vector< Real > &x) const
Compute where .
H1VectorPrimal< Real > PrimalStateVector
void hessVec_22(ROL::Vector< Real > &hv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Constraint_BurgersControl(ROL::Ptr< BurgersFEM< Real > > &fem, bool useHessian=true)
Real value(const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Compute value.
void plus(const ROL::Vector< Real > &x)
Compute , where .
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.
Real mesh_spacing(void) const
ROL::Ptr< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
bool isFeasible(const ROL::Vector< Real > &x)
Check if the vector, v, is feasible.
int dimension() const
Return dimension of the vector space.
ROL::Ptr< std::vector< Real > > vec_
void test_inverse_mass(std::ostream &outStream=std::cout)
void apply_pde_jacobian(std::vector< Real > &jv, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z) const
ROL::Ptr< BurgersFEM< Real > > fem_
const ROL::Ptr< const ROL::Vector< Real > > getLowerBound(void) const
Return the ref count pointer to the lower bound vector.
ROL::Ptr< ROL::Vector< Real > > ud_
ROL::Ptr< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
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 simulation-space derivative of the adjoint of the constraint optimization-space Jacobian at...
void pruneUpperActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &x, Real eps=Real(0))
Set variables to zero if they correspond to the upper -active set.
ROL::Ptr< const std::vector< Real > > getVector() const
ROL::Ptr< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
ROL::Ptr< const std::vector< Real > > getVector() const
ROL::Ptr< BurgersFEM< Real > > fem_
void set(const ROL::Vector< Real > &x)
Set where .
void scale(const Real alpha)
Compute where .
void compute_pde_jacobian(std::vector< Real > &dl, std::vector< Real > &d, std::vector< Real > &du, const std::vector< Real > &u) const
ROL::Ptr< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
H1VectorDual< Real > DualStateVector
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 compute_L2_dot(const std::vector< Real > &x, const std::vector< Real > &y) const
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 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) const
std::vector< Real > x_lo_
ROL::Ptr< std::vector< Real > > getVector()
L2VectorPrimal< Real > PrimalControlVector
ROL::Ptr< ROL::Vector< Real > > u_
void apply_inverse_mass(std::vector< Real > &Mu, const std::vector< Real > &u) const
void pruneActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &g, const ROL::Vector< Real > &x, Real xeps=Real(0), Real geps=Real(0))
void axpy(std::vector< Real > &out, const Real a, const std::vector< Real > &x, const std::vector< Real > &y) const
void hessVec_12(ROL::Vector< Real > &hv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
ROL::Ptr< std::vector< Real > > getVector()
ROL::Ptr< BurgersFEM< Real > > fem_
ROL::Ptr< const std::vector< Real > > getVector() const
ROL::Ptr< BurgersFEM< Real > > fem_
void update(std::vector< Real > &u, const std::vector< Real > &s, const Real alpha=1.0) const
void pruneUpperActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &x, Real eps=Real(0))
Set variables to zero if they correspond to the upper -active set.
void scale(const Real alpha)
Compute where .
void projection(std::vector< Real > &x)
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 pruneActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &x, Real eps=Real(0))
void set(const ROL::Vector< Real > &x)
Set where .
ROL::Ptr< ROL::Vector< Real > > l_
Real apply(const ROL::Vector< Real > &x) const
Apply to a dual vector. This is equivalent to the call .
void apply_adjoint_control_hessian(std::vector< Real > &ahwv, const std::vector< Real > &w, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z)
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 .
Provides the interface to apply upper and lower bound constraints.
void compute_residual(std::vector< Real > &r, const std::vector< Real > &u, const std::vector< Real > &z) const
void projection(std::vector< Real > &x)
void cast_const_vector(ROL::Ptr< const std::vector< Real > > &xvec, const ROL::Vector< Real > &x) const
ROL::Ptr< ROL::Vector< Real > > l_
ROL::Ptr< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) 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 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 apply_control_jacobian(std::vector< Real > &jv, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z) const
Real dot(const ROL::Vector< Real > &x) const override
Compute where .
void test_inverse_H1(std::ostream &outStream=std::cout)
const ROL::Vector< Real > & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
void apply_adjoint_control_jacobian(std::vector< Real > &jv, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z) const
Real dot(const ROL::Vector< Real > &x) const
Compute where .
L2VectorPrimal< Real > PrimalControlVector
void cast_const_vector(ROL::Ptr< const std::vector< Real > > &xvec, const ROL::Vector< Real > &x) const
H1VectorPrimal< Real > PrimalStateVector
std::vector< Real > x_lo_
void apply_inverse_pde_jacobian(std::vector< Real > &ijv, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z) const
void axpy(std::vector< Real > &out, const Real a, const std::vector< Real > &x, const std::vector< Real > &y) const
L2VectorDual< Real > DualControlVector
BurgersFEM(int nx=128, Real nu=1.e-2, Real nl=1.0, Real u0=1.0, Real u1=0.0, Real f=0.0, Real cH1=1.0, Real cL2=1.0)
void pruneLowerActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &g, const ROL::Vector< Real > &x, Real xeps=Real(0), Real geps=Real(0))
Set variables to zero if they correspond to the -binding set.
ROL::Ptr< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
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 optimization-space derivative of the adjoint of the constraint optimization-space Jacobian ...
void apply_H1(std::vector< Real > &Mu, const std::vector< Real > &u) const
void plus(const ROL::Vector< Real > &x)
Compute , where .
const ROL::Ptr< const ROL::Vector< Real > > getLowerBound(void) const
Return the ref count pointer to the lower bound vector.
void apply_inverse_H1(std::vector< Real > &Mu, const std::vector< Real > &u) const
void project(ROL::Vector< Real > &x)
Project optimization variables onto the bounds.
void apply_adjoint_pde_jacobian(std::vector< Real > &ajv, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z) const
void set(const ROL::Vector< Real > &x)
Set where .
ROL::Ptr< BurgersFEM< Real > > fem_
ROL::Ptr< ROL::Vector< Real > > u_
Real apply(const ROL::Vector< Real > &x) const
Apply to a dual vector. This is equivalent to the call .
ROL::Ptr< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
void apply_inverse_adjoint_pde_jacobian(std::vector< Real > &iajv, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z) const
void apply_adjoint_control_pde_hessian(std::vector< Real > &ahwv, const std::vector< Real > &w, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z)
void pruneUpperActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &g, const ROL::Vector< Real > &x, Real xeps=Real(0), Real geps=Real(0))
Set variables to zero if they correspond to the upper -binding set.
void pruneActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &x, Real eps=Real(0))
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...
int dimension() const
Return dimension of the vector space.
Defines the constraint operator interface for simulation-based optimization.
L2VectorDual< Real > DualControlVector
ROL::Ptr< H1VectorPrimal< Real > > dual_vec_
const ROL::Vector< Real > & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
void value(ROL::Vector< Real > &c, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Evaluate the constraint operator at .
void set(const ROL::Vector< Real > &x)
Set where .
Real norm() const
Returns where .
void scale(std::vector< Real > &u, const Real alpha=0.0) const
const ROL::Vector< Real > & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
ROL::Ptr< H1VectorDual< Real > > dual_vec_
Objective_BurgersControl(const ROL::Ptr< BurgersFEM< Real > > &fem, const ROL::Ptr< ROL::Vector< Real > > &ud, Real alpha=1.e-4)
const ROL::Ptr< const ROL::Vector< Real > > getUpperBound(void) const
Return the ref count pointer to the upper bound vector.
Real dot(const ROL::Vector< Real > &x) const override
Compute where .
ROL::Ptr< std::vector< Real > > getVector()
H1BoundConstraint(std::vector< Real > &l, std::vector< Real > &u, const ROL::Ptr< BurgersFEM< Real > > &fem, Real scale=1.0)
void scale(const Real alpha)
Compute where .
int dimension() const
Return dimension of the vector space.
const ROL::Vector< Real > & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
Real compute_L2_norm(const std::vector< Real > &r) const
ROL::Ptr< L2VectorPrimal< Real > > dual_vec_