55 #include "Teuchos_LAPACK.hpp" 83 void update(std::vector<Real> &u,
const std::vector<Real> &s,
const Real alpha=1.0)
const {
84 for (
unsigned i=0; i<u.size(); i++) {
89 void axpy(std::vector<Real> &out,
const Real a,
const std::vector<Real> &x,
const std::vector<Real> &y)
const {
90 for (
unsigned i=0; i < x.size(); i++) {
91 out[i] = a*x[i] + y[i];
95 void scale(std::vector<Real> &u,
const Real alpha=0.0)
const {
96 for (
unsigned i=0; i<u.size(); i++) {
101 void linear_solve(std::vector<Real> &u, std::vector<Real> &dl, std::vector<Real> &d, std::vector<Real> &du,
102 const std::vector<Real> &r,
const bool transpose =
false)
const {
103 if ( r.size() == 1 ) {
104 u.resize(1,r[0]/d[0]);
106 else if ( r.size() == 2 ) {
108 Real det = d[0]*d[1] - dl[0]*du[0];
109 u[0] = (d[1]*r[0] - du[0]*r[1])/det;
110 u[1] = (d[0]*r[1] - dl[0]*r[0])/det;
113 u.assign(r.begin(),r.end());
115 Teuchos::LAPACK<int,Real> lp;
116 std::vector<Real> du2(r.size()-2,0.0);
117 std::vector<int> ipiv(r.size(),0);
122 lp.GTTRF(dim,&dl[0],&d[0],&du[0],&du2[0],&ipiv[0],&info);
127 lp.GTTRS(trans,dim,nhrs,&dl[0],&d[0],&du[0],&du2[0],&ipiv[0],&u[0],ldb,&info);
133 Real u0 = 1.0, Real u1 = 0.0, Real f = 0.0,
134 Real cH1 = 1.0, Real cL2 = 1.0)
135 : nx_(nx), dx_(1.0/((Real)nx+1.0)),
136 nu_(nu), nl_(nl), u0_(u0), u1_(u1), f_(f),
137 cH1_(cH1), cL2_(cL2) {}
151 Real
compute_L2_dot(
const std::vector<Real> &x,
const std::vector<Real> &y)
const {
153 Real c = (((int)x.size()==
nx_) ? 4.0 : 2.0);
154 for (
unsigned i=0; i<x.size(); i++) {
156 ip += dx_/6.0*(c*x[i] + x[i+1])*y[i];
158 else if ( i == x.size()-1 ) {
159 ip += dx_/6.0*(x[i-1] + c*x[i])*y[i];
162 ip += dx_/6.0*(x[i-1] + 4.0*x[i] + x[i+1])*y[i];
174 void apply_mass(std::vector<Real> &Mu,
const std::vector<Real> &u )
const {
175 Mu.resize(u.size(),0.0);
176 Real c = (((int)u.size()==
nx_) ? 4.0 : 2.0);
177 for (
unsigned i=0; i<u.size(); i++) {
179 Mu[i] = dx_/6.0*(c*u[i] + u[i+1]);
181 else if ( i == u.size()-1 ) {
182 Mu[i] = dx_/6.0*(u[i-1] + c*u[i]);
185 Mu[i] = dx_/6.0*(u[i-1] + 4.0*u[i] + u[i+1]);
192 unsigned nx = u.size();
194 std::vector<Real> dl(nx-1,dx_/6.0);
195 std::vector<Real> d(nx,2.0*dx_/3.0);
196 std::vector<Real> du(nx-1,dx_/6.0);
197 if ( (
int)nx != nx_ ) {
206 std::vector<Real> u(nx_,0.0), Mu(nx_,0.0), iMMu(nx_,0.0), diff(nx_,0.0);
207 for (
int i = 0; i <
nx_; i++) {
208 u[i] = 2.0*(Real)rand()/(Real)RAND_MAX - 1.0;
212 axpy(diff,-1.0,iMMu,u);
215 outStream <<
"Test Inverse State Mass Matrix\n";
216 outStream <<
" ||u - inv(M)Mu|| = " << error <<
"\n";
217 outStream <<
" ||u|| = " << normu <<
"\n";
218 outStream <<
" Relative Error = " << error/normu <<
"\n";
221 u.resize(nx_+2,0.0); Mu.resize(nx_+2,0.0); iMMu.resize(nx_+2,0.0); diff.resize(nx_+2,0.0);
222 for (
int i = 0; i < nx_+2; i++) {
223 u[i] = 2.0*(Real)rand()/(Real)RAND_MAX - 1.0;
227 axpy(diff,-1.0,iMMu,u);
230 outStream <<
"Test Inverse Control Mass Matrix\n";
231 outStream <<
" ||z - inv(M)Mz|| = " << error <<
"\n";
232 outStream <<
" ||z|| = " << normu <<
"\n";
233 outStream <<
" Relative Error = " << error/normu <<
"\n";
241 Real
compute_H1_dot(
const std::vector<Real> &x,
const std::vector<Real> &y)
const {
243 for (
int i=0; i<
nx_; i++) {
245 ip += cL2_*dx_/6.0*(4.0*x[i] + x[i+1])*y[i];
246 ip += cH1_*(2.0*x[i] - x[i+1])/dx_*y[i];
248 else if ( i == nx_-1 ) {
249 ip += cL2_*dx_/6.0*(x[i-1] + 4.0*x[i])*y[i];
250 ip += cH1_*(2.0*x[i] - x[i-1])/dx_*y[i];
253 ip += cL2_*dx_/6.0*(x[i-1] + 4.0*x[i] + x[i+1])*y[i];
254 ip += cH1_*(2.0*x[i] - x[i-1] - x[i+1])/dx_*y[i];
266 void apply_H1(std::vector<Real> &Mu,
const std::vector<Real> &u )
const {
268 for (
int i=0; i<
nx_; i++) {
270 Mu[i] = cL2_*dx_/6.0*(4.0*u[i] + u[i+1])
271 + cH1_*(2.0*u[i] - u[i+1])/
dx_;
273 else if ( i == nx_-1 ) {
274 Mu[i] = cL2_*dx_/6.0*(u[i-1] + 4.0*u[i])
275 + cH1_*(2.0*u[i] - u[i-1])/
dx_;
278 Mu[i] = cL2_*dx_/6.0*(u[i-1] + 4.0*u[i] + u[i+1])
279 + cH1_*(2.0*u[i] - u[i-1] - u[i+1])/
dx_;
287 std::vector<Real> dl(nx_-1,cL2_*dx_/6.0 - cH1_/dx_);
288 std::vector<Real> d(nx_,2.0*(cL2_*dx_/3.0 + cH1_/dx_));
289 std::vector<Real> du(nx_-1,cL2_*dx_/6.0 - cH1_/dx_);
294 std::vector<Real> u(nx_,0.0), Mu(nx_,0.0), iMMu(nx_,0.0), diff(nx_,0.0);
295 for (
int i = 0; i <
nx_; i++) {
296 u[i] = 2.0*(Real)rand()/(Real)RAND_MAX - 1.0;
300 axpy(diff,-1.0,iMMu,u);
303 outStream <<
"Test Inverse State H1 Matrix\n";
304 outStream <<
" ||u - inv(M)Mu|| = " << error <<
"\n";
305 outStream <<
" ||u|| = " << normu <<
"\n";
306 outStream <<
" Relative Error = " << error/normu <<
"\n";
315 const std::vector<Real> &z)
const {
318 for (
int i=0; i<
nx_; i++) {
321 r[i] = nu_/dx_*(2.0*u[i]-u[i+1]);
324 r[i] = nu_/dx_*(2.0*u[i]-u[i-1]);
327 r[i] = nu_/dx_*(2.0*u[i]-u[i-1]-u[i+1]);
331 r[i] += nl_*u[i+1]*(u[i]+u[i+1])/6.0;
334 r[i] -= nl_*u[i-1]*(u[i-1]+u[i])/6.0;
337 r[i] -= dx_/6.0*(z[i]+4.0*z[i+1]+z[i+2]);
342 r[0] -= nl_*(u0_*u[ 0]/6.0 + u0_*u0_/6.0) + nu_*u0_/
dx_;
343 r[nx_-1] += nl_*(u1_*u[nx_-1]/6.0 + u1_*u1_/6.0) - nu_*u1_/
dx_;
351 std::vector<Real> &d,
352 std::vector<Real> &du,
353 const std::vector<Real> &u)
const {
356 d.resize(nx_,nu_*2.0/dx_);
358 dl.resize(nx_-1,-nu_/dx_);
360 du.resize(nx_-1,-nu_/dx_);
362 for (
int i=0; i<
nx_; i++) {
364 dl[i] += nl_*(-2.0*u[i]-u[i+1])/6.0;
365 d[i] += nl_*u[i+1]/6.0;
368 d[i] -= nl_*u[i-1]/6.0;
369 du[i-1] += nl_*(u[i-1]+2.0*u[i])/6.0;
373 d[ 0] -= nl_*u0_/6.0;
374 d[nx_-1] += nl_*u1_/6.0;
379 const std::vector<Real> &v,
380 const std::vector<Real> &u,
381 const std::vector<Real> &z)
const {
383 for (
int i = 0; i <
nx_; i++) {
384 jv[i] = nu_/dx_*2.0*v[i];
386 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]);
389 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]);
392 jv[ 0] -= nl_*u0_/6.0*v[0];
393 jv[nx_-1] += nl_*u1_/6.0*v[nx_-1];
398 const std::vector<Real> &v,
399 const std::vector<Real> &u,
400 const std::vector<Real> &z)
const {
402 std::vector<Real> d(nx_,0.0);
403 std::vector<Real> dl(nx_-1,0.0);
404 std::vector<Real> du(nx_-1,0.0);
412 const std::vector<Real> &v,
413 const std::vector<Real> &u,
414 const std::vector<Real> &z)
const {
416 for (
int i = 0; i <
nx_; i++) {
417 ajv[i] = nu_/dx_*2.0*v[i];
419 ajv[i] += -nu_/dx_*v[i-1]-nl_*(u[i-1]/6.0*v[i]
420 -(u[i-1]+2.0*u[i])/6.0*v[i-1]);
423 ajv[i] += -nu_/dx_*v[i+1]+nl_*(u[i+1]/6.0*v[i]
424 -(u[i+1]+2.0*u[i])/6.0*v[i+1]);
427 ajv[ 0] -= nl_*u0_/6.0*v[0];
428 ajv[nx_-1] += nl_*u1_/6.0*v[nx_-1];
433 const std::vector<Real> &v,
434 const std::vector<Real> &u,
435 const std::vector<Real> &z)
const {
437 std::vector<Real> d(nx_,0.0);
438 std::vector<Real> du(nx_-1,0.0);
439 std::vector<Real> dl(nx_-1,0.0);
450 const std::vector<Real> &v,
451 const std::vector<Real> &u,
452 const std::vector<Real> &z)
const {
453 for (
int i=0; i<
nx_; i++) {
455 jv[i] = -dx_/6.0*(v[i]+4.0*v[i+1]+v[i+2]);
461 const std::vector<Real> &v,
462 const std::vector<Real> &u,
463 const std::vector<Real> &z)
const {
464 for (
int i=0; i<nx_+2; i++) {
466 jv[i] = -dx_/6.0*v[i];
469 jv[i] = -dx_/6.0*(4.0*v[i-1]+v[i]);
471 else if ( i == nx_ ) {
472 jv[i] = -dx_/6.0*(4.0*v[i-1]+v[i-2]);
474 else if ( i == nx_+1 ) {
475 jv[i] = -dx_/6.0*v[i-2];
478 jv[i] = -dx_/6.0*(v[i-2]+4.0*v[i-1]+v[i]);
487 const std::vector<Real> &w,
488 const std::vector<Real> &v,
489 const std::vector<Real> &u,
490 const std::vector<Real> &z)
const {
491 for (
int i=0; i<
nx_; i++) {
495 ahwv[i] += (w[i]*v[i+1] - w[i+1]*(2.0*v[i]+v[i+1]))/6.0;
498 ahwv[i] += (w[i-1]*(v[i-1]+2.0*v[i]) - w[i]*v[i-1])/6.0;
504 const std::vector<Real> &w,
505 const std::vector<Real> &v,
506 const std::vector<Real> &u,
507 const std::vector<Real> &z) {
508 ahwv.assign(u.size(),0.0);
511 const std::vector<Real> &w,
512 const std::vector<Real> &v,
513 const std::vector<Real> &u,
514 const std::vector<Real> &z) {
515 ahwv.assign(z.size(),0.0);
518 const std::vector<Real> &w,
519 const std::vector<Real> &v,
520 const std::vector<Real> &u,
521 const std::vector<Real> &z) {
522 ahwv.assign(z.size(),0.0);
529 Teuchos::RCP<std::vector<Real> > vec_;
530 Teuchos::RCP<BurgersFEM<Real> > fem_;
532 mutable Teuchos::RCP<L2VectorDual<Real> > dual_vec_;
537 : vec_(vec), fem_(fem), dual_vec_(Teuchos::null) {}
541 const std::vector<Real>& xval = *ex.
getVector();
542 std::copy(xval.begin(),xval.end(),vec_->begin());
547 const std::vector<Real>& xval = *ex.
getVector();
548 unsigned dimension = vec_->size();
549 for (
unsigned i=0; i<dimension; i++) {
550 (*vec_)[i] += xval[i];
555 unsigned dimension = vec_->size();
556 for (
unsigned i=0; i<dimension; i++) {
563 const std::vector<Real>& xval = *ex.
getVector();
564 return fem_->compute_L2_dot(xval,*vec_);
569 val = std::sqrt( dot(*
this) );
573 Teuchos::RCP<ROL::Vector<Real> >
clone()
const {
574 return Teuchos::rcp(
new L2VectorPrimal( Teuchos::rcp(
new std::vector<Real>(vec_->size(),0.0)),fem_));
577 Teuchos::RCP<const std::vector<Real> >
getVector()
const {
585 Teuchos::RCP<ROL::Vector<Real> >
basis(
const int i )
const {
586 Teuchos::RCP<L2VectorPrimal> e
587 = Teuchos::rcp(
new L2VectorPrimal( Teuchos::rcp(
new std::vector<Real>(vec_->size(),0.0)),fem_));
588 (*e->getVector())[i] = 1.0;
598 Teuchos::rcp(
new std::vector<Real>(*vec_)),fem_));
600 fem_->apply_mass(*(Teuchos::rcp_const_cast<std::vector<Real> >(dual_vec_->getVector())),*vec_);
609 Teuchos::RCP<std::vector<Real> > vec_;
610 Teuchos::RCP<BurgersFEM<Real> > fem_;
612 mutable Teuchos::RCP<L2VectorPrimal<Real> > dual_vec_;
617 : vec_(vec), fem_(fem), dual_vec_(Teuchos::null) {}
621 const std::vector<Real>& xval = *ex.
getVector();
622 std::copy(xval.begin(),xval.end(),vec_->begin());
627 const std::vector<Real>& xval = *ex.
getVector();
628 unsigned dimension = vec_->size();
629 for (
unsigned i=0; i<dimension; i++) {
630 (*vec_)[i] += xval[i];
635 unsigned dimension = vec_->size();
636 for (
unsigned i=0; i<dimension; i++) {
643 const std::vector<Real>& xval = *ex.
getVector();
644 unsigned dimension = vec_->size();
645 std::vector<Real> Mx(dimension,0.0);
646 fem_->apply_inverse_mass(Mx,xval);
648 for (
unsigned i = 0; i < dimension; i++) {
649 val += Mx[i]*(*vec_)[i];
656 val = std::sqrt( dot(*
this) );
660 Teuchos::RCP<ROL::Vector<Real> >
clone()
const {
661 return Teuchos::rcp(
new L2VectorDual( Teuchos::rcp(
new std::vector<Real>(vec_->size(),0.0)),fem_));
664 Teuchos::RCP<const std::vector<Real> >
getVector()
const {
672 Teuchos::RCP<ROL::Vector<Real> >
basis(
const int i )
const {
673 Teuchos::RCP<L2VectorDual> e
674 = Teuchos::rcp(
new L2VectorDual( Teuchos::rcp(
new std::vector<Real>(vec_->size(),0.0)),fem_));
675 (*e->getVector())[i] = 1.0;
685 Teuchos::rcp(
new std::vector<Real>(*vec_)),fem_));
687 fem_->apply_inverse_mass(*(Teuchos::rcp_const_cast<std::vector<Real> >(dual_vec_->getVector())),*vec_);
696 Teuchos::RCP<std::vector<Real> > vec_;
697 Teuchos::RCP<BurgersFEM<Real> > fem_;
699 mutable Teuchos::RCP<H1VectorDual<Real> > dual_vec_;
704 : vec_(vec), fem_(fem), dual_vec_(Teuchos::null) {}
708 const std::vector<Real>& xval = *ex.
getVector();
709 std::copy(xval.begin(),xval.end(),vec_->begin());
714 const std::vector<Real>& xval = *ex.
getVector();
715 unsigned dimension = vec_->size();
716 for (
unsigned i=0; i<dimension; i++) {
717 (*vec_)[i] += xval[i];
722 unsigned dimension = vec_->size();
723 for (
unsigned i=0; i<dimension; i++) {
730 const std::vector<Real>& xval = *ex.
getVector();
731 return fem_->compute_H1_dot(xval,*vec_);
736 val = std::sqrt( dot(*
this) );
740 Teuchos::RCP<ROL::Vector<Real> >
clone()
const {
741 return Teuchos::rcp(
new H1VectorPrimal( Teuchos::rcp(
new std::vector<Real>(vec_->size(),0.0)),fem_));
744 Teuchos::RCP<const std::vector<Real> >
getVector()
const {
752 Teuchos::RCP<ROL::Vector<Real> >
basis(
const int i )
const {
753 Teuchos::RCP<H1VectorPrimal> e
754 = Teuchos::rcp(
new H1VectorPrimal( Teuchos::rcp(
new std::vector<Real>(vec_->size(),0.0)),fem_));
755 (*e->getVector())[i] = 1.0;
765 Teuchos::rcp(
new std::vector<Real>(*vec_)),fem_));
767 fem_->apply_H1(*(Teuchos::rcp_const_cast<std::vector<Real> >(dual_vec_->getVector())),*vec_);
776 Teuchos::RCP<std::vector<Real> > vec_;
777 Teuchos::RCP<BurgersFEM<Real> > fem_;
779 mutable Teuchos::RCP<H1VectorPrimal<Real> > dual_vec_;
784 : vec_(vec), fem_(fem), dual_vec_(Teuchos::null) {}
788 const std::vector<Real>& xval = *ex.
getVector();
789 std::copy(xval.begin(),xval.end(),vec_->begin());
794 const std::vector<Real>& xval = *ex.
getVector();
795 unsigned dimension = vec_->size();
796 for (
unsigned i=0; i<dimension; i++) {
797 (*vec_)[i] += xval[i];
802 unsigned dimension = vec_->size();
803 for (
unsigned i=0; i<dimension; i++) {
810 const std::vector<Real>& xval = *ex.
getVector();
811 unsigned dimension = vec_->size();
812 std::vector<Real> Mx(dimension,0.0);
813 fem_->apply_inverse_H1(Mx,xval);
815 for (
unsigned i = 0; i < dimension; i++) {
816 val += Mx[i]*(*vec_)[i];
823 val = std::sqrt( dot(*
this) );
827 Teuchos::RCP<ROL::Vector<Real> >
clone()
const {
828 return Teuchos::rcp(
new H1VectorDual( Teuchos::rcp(
new std::vector<Real>(vec_->size(),0.0)),fem_));
831 Teuchos::RCP<const std::vector<Real> >
getVector()
const {
839 Teuchos::RCP<ROL::Vector<Real> >
basis(
const int i )
const {
840 Teuchos::RCP<H1VectorDual> e
841 = Teuchos::rcp(
new H1VectorDual( Teuchos::rcp(
new std::vector<Real>(vec_->size(),0.0)),fem_));
842 (*e->getVector())[i] = 1.0;
852 Teuchos::rcp(
new std::vector<Real>(*vec_)),fem_));
854 fem_->apply_inverse_H1(*(Teuchos::rcp_const_cast<std::vector<Real> >(dual_vec_->getVector())),*vec_);
873 Teuchos::RCP<BurgersFEM<Real> > fem_;
878 : fem_(fem), useHessian_(useHessian) {}
882 Teuchos::RCP<std::vector<Real> > cp =
883 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<PrimalConstraintVector>(c)).getVector());
884 Teuchos::RCP<const std::vector<Real> > up =
885 (Teuchos::dyn_cast<PrimalStateVector>(
const_cast<ROL::Vector<Real> &
>(u))).getVector();
886 Teuchos::RCP<const std::vector<Real> > zp =
887 (Teuchos::dyn_cast<PrimalControlVector>(
const_cast<ROL::Vector<Real> &
>(z))).getVector();
888 fem_->compute_residual(*cp,*up,*zp);
895 Teuchos::RCP<std::vector<Real> > jvp =
896 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<PrimalConstraintVector>(jv)).getVector());
897 Teuchos::RCP<const std::vector<Real> > vp =
898 (Teuchos::dyn_cast<PrimalStateVector>(
const_cast<ROL::Vector<Real> &
>(v))).getVector();
899 Teuchos::RCP<const std::vector<Real> > up =
900 (Teuchos::dyn_cast<PrimalStateVector>(
const_cast<ROL::Vector<Real> &
>(u))).getVector();
901 Teuchos::RCP<const std::vector<Real> > zp =
902 (Teuchos::dyn_cast<PrimalControlVector>(
const_cast<ROL::Vector<Real> &
>(z))).getVector();
903 fem_->apply_pde_jacobian(*jvp,*vp,*up,*zp);
908 Teuchos::RCP<std::vector<Real> > jvp =
909 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<PrimalConstraintVector>(jv)).getVector());
910 Teuchos::RCP<const std::vector<Real> > vp =
911 (Teuchos::dyn_cast<PrimalControlVector>(
const_cast<ROL::Vector<Real> &
>(v))).getVector();
912 Teuchos::RCP<const std::vector<Real> > up =
913 (Teuchos::dyn_cast<PrimalStateVector>(
const_cast<ROL::Vector<Real> &
>(u))).getVector();
914 Teuchos::RCP<const std::vector<Real> > zp =
915 (Teuchos::dyn_cast<PrimalControlVector>(
const_cast<ROL::Vector<Real> &
>(z))).getVector();
916 fem_->apply_control_jacobian(*jvp,*vp,*up,*zp);
921 Teuchos::RCP<std::vector<Real> > ijvp =
922 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<PrimalStateVector>(ijv)).getVector());
923 Teuchos::RCP<const std::vector<Real> > vp =
924 (Teuchos::dyn_cast<PrimalConstraintVector>(
const_cast<ROL::Vector<Real> &
>(v))).getVector();
925 Teuchos::RCP<const std::vector<Real> > up =
926 (Teuchos::dyn_cast<PrimalStateVector>(
const_cast<ROL::Vector<Real> &
>(u))).getVector();
927 Teuchos::RCP<const std::vector<Real> > zp =
928 (Teuchos::dyn_cast<PrimalControlVector>(
const_cast<ROL::Vector<Real> &
>(z))).getVector();
929 fem_->apply_inverse_pde_jacobian(*ijvp,*vp,*up,*zp);
934 Teuchos::RCP<std::vector<Real> > jvp =
935 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<DualStateVector>(ajv)).getVector());
936 Teuchos::RCP<const std::vector<Real> > vp =
937 (Teuchos::dyn_cast<DualConstraintVector>(
const_cast<ROL::Vector<Real> &
>(v))).getVector();
938 Teuchos::RCP<const std::vector<Real> > up =
939 (Teuchos::dyn_cast<PrimalStateVector>(
const_cast<ROL::Vector<Real> &
>(u))).getVector();
940 Teuchos::RCP<const std::vector<Real> > zp =
941 (Teuchos::dyn_cast<PrimalControlVector>(
const_cast<ROL::Vector<Real> &
>(z))).getVector();
942 fem_->apply_adjoint_pde_jacobian(*jvp,*vp,*up,*zp);
947 Teuchos::RCP<std::vector<Real> > jvp =
948 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<DualControlVector>(jv)).getVector());
949 Teuchos::RCP<const std::vector<Real> > vp =
950 (Teuchos::dyn_cast<DualConstraintVector>(
const_cast<ROL::Vector<Real> &
>(v))).getVector();
951 Teuchos::RCP<const std::vector<Real> > up =
952 (Teuchos::dyn_cast<PrimalStateVector>(
const_cast<ROL::Vector<Real> &
>(u))).getVector();
953 Teuchos::RCP<const std::vector<Real> > zp =
954 (Teuchos::dyn_cast<PrimalControlVector>(
const_cast<ROL::Vector<Real> &
>(z))).getVector();
955 fem_->apply_adjoint_control_jacobian(*jvp,*vp,*up,*zp);
960 Teuchos::RCP<std::vector<Real> > iajvp =
961 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<DualConstraintVector>(iajv)).getVector());
962 Teuchos::RCP<const std::vector<Real> > vp =
963 (Teuchos::dyn_cast<DualStateVector>(
const_cast<ROL::Vector<Real> &
>(v))).getVector();
964 Teuchos::RCP<const std::vector<Real> > up =
965 (Teuchos::dyn_cast<PrimalStateVector>(
const_cast<ROL::Vector<Real> &
>(u))).getVector();
966 Teuchos::RCP<const std::vector<Real> > zp =
967 (Teuchos::dyn_cast<PrimalControlVector>(
const_cast<ROL::Vector<Real> &
>(z))).getVector();
968 fem_->apply_inverse_adjoint_pde_jacobian(*iajvp,*vp,*up,*zp);
974 Teuchos::RCP<std::vector<Real> > ahwvp =
975 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<DualStateVector>(ahwv)).getVector());
976 Teuchos::RCP<const std::vector<Real> > wp =
977 (Teuchos::dyn_cast<DualConstraintVector>(
const_cast<ROL::Vector<Real> &
>(w))).getVector();
978 Teuchos::RCP<const std::vector<Real> > vp =
979 (Teuchos::dyn_cast<PrimalStateVector>(
const_cast<ROL::Vector<Real> &
>(v))).getVector();
980 Teuchos::RCP<const std::vector<Real> > up =
981 (Teuchos::dyn_cast<PrimalStateVector>(
const_cast<ROL::Vector<Real> &
>(u))).getVector();
982 Teuchos::RCP<const std::vector<Real> > zp =
983 (Teuchos::dyn_cast<PrimalControlVector>(
const_cast<ROL::Vector<Real> &
>(z))).getVector();
984 fem_->apply_adjoint_pde_hessian(*ahwvp,*wp,*vp,*up,*zp);
994 Teuchos::RCP<std::vector<Real> > ahwvp =
995 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<DualControlVector>(ahwv)).getVector());
996 Teuchos::RCP<const std::vector<Real> > wp =
997 (Teuchos::dyn_cast<DualConstraintVector>(
const_cast<ROL::Vector<Real> &
>(w))).getVector();
998 Teuchos::RCP<const std::vector<Real> > vp =
999 (Teuchos::dyn_cast<PrimalStateVector>(
const_cast<ROL::Vector<Real> &
>(v))).getVector();
1000 Teuchos::RCP<const std::vector<Real> > up =
1001 (Teuchos::dyn_cast<PrimalStateVector>(
const_cast<ROL::Vector<Real> &
>(u))).getVector();
1002 Teuchos::RCP<const std::vector<Real> > zp =
1003 (Teuchos::dyn_cast<PrimalControlVector>(
const_cast<ROL::Vector<Real> &
>(z))).getVector();
1004 fem_->apply_adjoint_control_pde_hessian(*ahwvp,*wp,*vp,*up,*zp);
1012 if ( useHessian_ ) {
1013 Teuchos::RCP<std::vector<Real> > ahwvp =
1014 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<DualStateVector>(ahwv)).getVector());
1015 Teuchos::RCP<const std::vector<Real> > wp =
1016 (Teuchos::dyn_cast<DualConstraintVector>(
const_cast<ROL::Vector<Real> &
>(w))).getVector();
1017 Teuchos::RCP<const std::vector<Real> > vp =
1018 (Teuchos::dyn_cast<PrimalControlVector>(
const_cast<ROL::Vector<Real> &
>(v))).getVector();
1019 Teuchos::RCP<const std::vector<Real> > up =
1020 (Teuchos::dyn_cast<PrimalStateVector>(
const_cast<ROL::Vector<Real> &
>(u))).getVector();
1021 Teuchos::RCP<const std::vector<Real> > zp =
1022 (Teuchos::dyn_cast<PrimalControlVector>(
const_cast<ROL::Vector<Real> &
>(z))).getVector();
1023 fem_->apply_adjoint_pde_control_hessian(*ahwvp,*wp,*vp,*up,*zp);
1031 if ( useHessian_ ) {
1032 Teuchos::RCP<std::vector<Real> > ahwvp =
1033 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<DualControlVector>(ahwv)).getVector());
1034 Teuchos::RCP<const std::vector<Real> > wp =
1035 (Teuchos::dyn_cast<DualConstraintVector>(
const_cast<ROL::Vector<Real> &
>(w))).getVector();
1036 Teuchos::RCP<const std::vector<Real> > vp =
1037 (Teuchos::dyn_cast<PrimalControlVector>(
const_cast<ROL::Vector<Real> &
>(v))).getVector();
1038 Teuchos::RCP<const std::vector<Real> > up =
1039 (Teuchos::dyn_cast<PrimalStateVector>(
const_cast<ROL::Vector<Real> &
>(u))).getVector();
1040 Teuchos::RCP<const std::vector<Real> > zp =
1041 (Teuchos::dyn_cast<PrimalControlVector>(
const_cast<ROL::Vector<Real> &
>(z))).getVector();
1042 fem_->apply_adjoint_control_hessian(*ahwvp,*wp,*vp,*up,*zp);
1050 template<
class Real>
1061 Teuchos::RCP<BurgersFEM<Real> >
fem_;
1062 Teuchos::RCP<ROL::Vector<Real> >
ud_;
1068 Real alpha = 1.e-4) : alpha_(alpha), fem_(fem), ud_(ud) {
1069 diff_ = ud_->clone();
1073 Teuchos::RCP<const std::vector<Real> > up =
1074 (Teuchos::dyn_cast<PrimalStateVector>(
const_cast<ROL::Vector<Real> &
>(u))).getVector();
1075 Teuchos::RCP<const std::vector<Real> > zp =
1076 (Teuchos::dyn_cast<PrimalControlVector>(
const_cast<ROL::Vector<Real> &
>(z))).getVector();
1077 Teuchos::RCP<const std::vector<Real> > udp =
1080 std::vector<Real> diff(udp->size(),0.0);
1081 for (
unsigned i = 0; i < udp->size(); i++) {
1082 diff[i] = (*up)[i] - (*udp)[i];
1084 return 0.5*(fem_->compute_L2_dot(diff,diff) + alpha_*fem_->compute_L2_dot(*zp,*zp));
1088 Teuchos::RCP<std::vector<Real> > gp =
1089 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<DualStateVector>(g)).getVector());
1090 Teuchos::RCP<const std::vector<Real> > up =
1091 (Teuchos::dyn_cast<PrimalStateVector>(
const_cast<ROL::Vector<Real> &
>(u))).getVector();
1092 Teuchos::RCP<const std::vector<Real> > udp =
1095 std::vector<Real> diff(udp->size(),0.0);
1096 for (
unsigned i = 0; i < udp->size(); i++) {
1097 diff[i] = (*up)[i] - (*udp)[i];
1099 fem_->apply_mass(*gp,diff);
1103 Teuchos::RCP<std::vector<Real> > gp =
1104 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<DualControlVector>(g)).getVector());
1105 Teuchos::RCP<const std::vector<Real> > zp =
1106 (Teuchos::dyn_cast<PrimalControlVector>(
const_cast<ROL::Vector<Real> &
>(z))).getVector();
1108 fem_->apply_mass(*gp,*zp);
1109 for (
unsigned i = 0; i < zp->size(); i++) {
1116 Teuchos::RCP<std::vector<Real> > hvp =
1117 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<DualStateVector>(hv)).getVector());
1118 Teuchos::RCP<const std::vector<Real> > vp =
1119 (Teuchos::dyn_cast<PrimalStateVector>(
const_cast<ROL::Vector<Real> &
>(v))).getVector();
1121 fem_->apply_mass(*hvp,*vp);
1136 Teuchos::RCP<std::vector<Real> > hvp =
1137 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<DualControlVector>(hv)).getVector());
1138 Teuchos::RCP<const std::vector<Real> > vp =
1139 (Teuchos::dyn_cast<PrimalControlVector>(
const_cast<ROL::Vector<Real> &
>(v))).getVector();
1141 fem_->apply_mass(*hvp,*vp);
1142 for (
unsigned i = 0; i < vp->size(); i++) {
1143 (*hvp)[i] *= alpha_;
1148 template<
class Real>
1156 Teuchos::RCP<BurgersFEM<Real> >
fem_;
1161 xvec = Teuchos::rcp_const_cast<std::vector<Real> >(
1164 catch (std::exception &e) {
1165 xvec = Teuchos::rcp_const_cast<std::vector<Real> >(
1176 catch (std::exception &e) {
1182 void axpy(std::vector<Real> &out,
const Real a,
1183 const std::vector<Real> &x,
const std::vector<Real> &y)
const{
1184 out.resize(dim_,0.0);
1185 for (
unsigned i = 0; i < dim_; i++) {
1186 out[i] = a*x[i] + y[i];
1191 for (
int i = 0; i < dim_; i++ ) {
1192 x[i] = std::max(x_lo_[i],std::min(x_up_[i],x[i]));
1199 : x_lo_(l), x_up_(u), scale_(
scale), fem_(fem) {
1200 dim_ = x_lo_.size();
1201 for (
int i = 0; i < dim_; i++ ) {
1203 min_diff_ = x_up_[i] - x_lo_[i];
1206 min_diff_ = ( (min_diff_ < (x_up_[i] - x_lo_[i])) ? min_diff_ : (x_up_[i] - x_lo_[i]) );
1213 Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1216 for (
int i = 0; i < dim_; i++ ) {
1217 if ( (*ex)[i] >= x_lo_[i] && (*ex)[i] <= x_up_[i] ) { cnt *= 1; }
1220 if ( cnt == 0 ) { val =
false; }
1225 Teuchos::RCP<std::vector<Real> > ex; cast_vector(ex,x);
1230 Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1231 Teuchos::RCP<std::vector<Real> > ev; cast_vector(ev,v);
1232 Real epsn = std::min(scale_*eps,min_diff_);
1233 for (
int i = 0; i < dim_; i++ ) {
1234 if ( ((*ex)[i] <= x_lo_[i]+epsn) ) {
1241 Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1242 Teuchos::RCP<std::vector<Real> > ev; cast_vector(ev,v);
1243 Real epsn = std::min(scale_*eps,min_diff_);
1244 for (
int i = 0; i < dim_; i++ ) {
1245 if ( ((*ex)[i] >= x_up_[i]-epsn) ) {
1252 Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1253 Teuchos::RCP<std::vector<Real> > ev; cast_vector(ev,v);
1254 Real epsn = std::min(scale_*eps,min_diff_);
1255 for (
int i = 0; i < dim_; i++ ) {
1256 if ( ((*ex)[i] <= x_lo_[i]+epsn) ||
1257 ((*ex)[i] >= x_up_[i]-epsn) ) {
1264 Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1265 Teuchos::RCP<const std::vector<Real> > eg; cast_const_vector(eg,g);
1266 Teuchos::RCP<std::vector<Real> > ev; cast_vector(ev,v);
1267 Real epsn = std::min(scale_*eps,min_diff_);
1268 for (
int i = 0; i < dim_; i++ ) {
1269 if ( ((*ex)[i] <= x_lo_[i]+epsn && (*eg)[i] > 0.0) ) {
1276 Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1277 Teuchos::RCP<const std::vector<Real> > eg; cast_const_vector(eg,g);
1278 Teuchos::RCP<std::vector<Real> > ev; cast_vector(ev,v);
1279 Real epsn = std::min(scale_*eps,min_diff_);
1280 for (
int i = 0; i < dim_; i++ ) {
1281 if ( ((*ex)[i] >= x_up_[i]-epsn && (*eg)[i] < 0.0) ) {
1288 Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1289 Teuchos::RCP<const std::vector<Real> > eg; cast_const_vector(eg,g);
1290 Teuchos::RCP<std::vector<Real> > ev; cast_vector(ev,v);
1291 Real epsn = std::min(scale_*eps,min_diff_);
1292 for (
int i = 0; i < dim_; i++ ) {
1293 if ( ((*ex)[i] <= x_lo_[i]+epsn && (*eg)[i] > 0.0) ||
1294 ((*ex)[i] >= x_up_[i]-epsn && (*eg)[i] < 0.0) ) {
1301 Teuchos::RCP<std::vector<Real> > us = Teuchos::rcp(
new std::vector<Real>(dim_,0.0) );
1302 us->assign(x_up_.begin(),x_up_.end());
1308 Teuchos::RCP<std::vector<Real> > ls = Teuchos::rcp(
new std::vector<Real>(dim_,0.0) );
1309 ls->assign(x_lo_.begin(),x_lo_.end());
1315 template<
class Real>
1323 Teuchos::RCP<BurgersFEM<Real> >
fem_;
1328 xvec = Teuchos::rcp_const_cast<std::vector<Real> >(
1331 catch (std::exception &e) {
1332 xvec = Teuchos::rcp_const_cast<std::vector<Real> >(
1343 catch (std::exception &e) {
1349 void axpy(std::vector<Real> &out,
const Real a,
1350 const std::vector<Real> &x,
const std::vector<Real> &y)
const{
1351 out.resize(dim_,0.0);
1352 for (
unsigned i = 0; i < dim_; i++) {
1353 out[i] = a*x[i] + y[i];
1358 for (
int i = 0; i < dim_; i++ ) {
1359 x[i] = std::max(x_lo_[i],std::min(x_up_[i],x[i]));
1366 : x_lo_(l), x_up_(u), scale_(
scale), fem_(fem) {
1367 dim_ = x_lo_.size();
1368 for (
int i = 0; i < dim_; i++ ) {
1370 min_diff_ = x_up_[i] - x_lo_[i];
1373 min_diff_ = ( (min_diff_ < (x_up_[i] - x_lo_[i])) ? min_diff_ : (x_up_[i] - x_lo_[i]) );
1380 Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1383 for (
int i = 0; i < dim_; i++ ) {
1384 if ( (*ex)[i] >= x_lo_[i] && (*ex)[i] <= x_up_[i] ) { cnt *= 1; }
1387 if ( cnt == 0 ) { val =
false; }
1392 Teuchos::RCP<std::vector<Real> > ex; cast_vector(ex,x);
1397 Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1398 Teuchos::RCP<std::vector<Real> > ev; cast_vector(ev,v);
1399 Real epsn = std::min(scale_*eps,min_diff_);
1400 for (
int i = 0; i < dim_; i++ ) {
1401 if ( ((*ex)[i] <= x_lo_[i]+epsn) ) {
1408 Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1409 Teuchos::RCP<std::vector<Real> > ev; cast_vector(ev,v);
1410 Real epsn = std::min(scale_*eps,min_diff_);
1411 for (
int i = 0; i < dim_; i++ ) {
1412 if ( ((*ex)[i] >= x_up_[i]-epsn) ) {
1419 Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1420 Teuchos::RCP<std::vector<Real> > ev; cast_vector(ev,v);
1421 Real epsn = std::min(scale_*eps,min_diff_);
1422 for (
int i = 0; i < dim_; i++ ) {
1423 if ( ((*ex)[i] <= x_lo_[i]+epsn) ||
1424 ((*ex)[i] >= x_up_[i]-epsn) ) {
1431 Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1432 Teuchos::RCP<const std::vector<Real> > eg; cast_const_vector(eg,g);
1433 Teuchos::RCP<std::vector<Real> > ev; cast_vector(ev,v);
1434 Real epsn = std::min(scale_*eps,min_diff_);
1435 for (
int i = 0; i < dim_; i++ ) {
1436 if ( ((*ex)[i] <= x_lo_[i]+epsn && (*eg)[i] > 0.0) ) {
1443 Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1444 Teuchos::RCP<const std::vector<Real> > eg; cast_const_vector(eg,g);
1445 Teuchos::RCP<std::vector<Real> > ev; cast_vector(ev,v);
1446 Real epsn = std::min(scale_*eps,min_diff_);
1447 for (
int i = 0; i < dim_; i++ ) {
1448 if ( ((*ex)[i] >= x_up_[i]-epsn && (*eg)[i] < 0.0) ) {
1455 Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1456 Teuchos::RCP<const std::vector<Real> > eg; cast_const_vector(eg,g);
1457 Teuchos::RCP<std::vector<Real> > ev; cast_vector(ev,v);
1458 Real epsn = std::min(scale_*eps,min_diff_);
1459 for (
int i = 0; i < dim_; i++ ) {
1460 if ( ((*ex)[i] <= x_lo_[i]+epsn && (*eg)[i] > 0.0) ||
1461 ((*ex)[i] >= x_up_[i]-epsn && (*eg)[i] < 0.0) ) {
1468 Teuchos::RCP<std::vector<Real> > us = Teuchos::rcp(
new std::vector<Real>(dim_,0.0) );
1469 us->assign(x_up_.begin(),x_up_.end());
1475 Teuchos::RCP<std::vector<Real> > ls = Teuchos::rcp(
new std::vector<Real>(dim_,0.0) );
1476 ls->assign(x_lo_.begin(),x_lo_.end());
Provides the interface to evaluate simulation-based objective functions.
void pruneUpperActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &x, Real eps)
Set variables to zero if they correspond to the upper -active set.
void cast_vector(Teuchos::RCP< std::vector< Real > > &xvec, ROL::Vector< Real > &x) const
std::vector< Real > x_up_
L2VectorPrimal< Real > PrimalControlVector
Real compute_H1_dot(const std::vector< Real > &x, const std::vector< Real > &y) const
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 ...
Real norm() const
Returns where .
void cast_const_vector(Teuchos::RCP< const std::vector< Real > > &xvec, const ROL::Vector< Real > &x) const
void apply_mass(std::vector< Real > &Mu, const std::vector< Real > &u) const
Teuchos::RCP< const std::vector< Real > > getVector() const
bool isFeasible(const ROL::Vector< Real > &x)
Check if the vector, v, is feasible.
void axpy(std::vector< Real > &out, const Real a, const std::vector< Real > &x, const std::vector< Real > &y) const
Teuchos::RCP< std::vector< Real > > getVector()
Teuchos::RCP< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
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)
Teuchos::RCP< const std::vector< Real > > getVector() const
int dimension() const
Return dimension of the vector space.
Teuchos::RCP< ROL::Vector< Real > > ud_
void plus(const ROL::Vector< Real > &x)
Compute , where .
Objective_BurgersControl(const Teuchos::RCP< BurgersFEM< Real > > &fem, const Teuchos::RCP< ROL::Vector< Real > > &ud, Real alpha=1.e-4)
std::vector< Real > x_up_
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 ...
Teuchos::RCP< BurgersFEM< Real > > fem_
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 .
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
L2VectorDual(const Teuchos::RCP< std::vector< Real > > &vec, const Teuchos::RCP< BurgersFEM< Real > > &fem)
Teuchos::RCP< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
Real norm() const
Returns where .
void pruneActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &g, const ROL::Vector< Real > &x, Real eps)
Set variables to zero if they correspond to the -binding set.
Contains definitions of custom data types in ROL.
Real dot(const ROL::Vector< Real > &x) const
Compute where .
void plus(const ROL::Vector< Real > &x)
Compute , where .
void pruneLowerActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &x, Real eps)
Set variables to zero if they correspond to the lower -active set.
Teuchos::RCP< const std::vector< Real > > getVector() const
Real compute_H1_norm(const std::vector< Real > &r) const
Real norm() const
Returns where .
void pruneUpperActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &g, const ROL::Vector< Real > &x, Real eps)
Set variables to zero if they correspond to the upper -binding set.
Teuchos::RCP< std::vector< Real > > getVector()
void project(ROL::Vector< Real > &x)
Project optimization variables onto the bounds.
Teuchos::RCP< std::vector< Real > > getVector()
virtual void zero()
Set to zero vector.
Defines the linear algebra or vector space interface.
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 ...
H1VectorDual< Real > DualStateVector
void scale(const Real alpha)
Compute where .
Defines the equality constraint operator interface for simulation-based optimization.
L2VectorDual< Real > DualControlVector
Real dot(const ROL::Vector< Real > &x) const
Compute where .
void hessVec_22(ROL::Vector< Real > &hv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Teuchos::RCP< ROL::Vector< Real > > diff_
Real value(const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Compute value.
void pruneActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &x, Real eps)
Set variables to zero if they correspond to the -active set.
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
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 ...
bool isFeasible(const ROL::Vector< Real > &x)
Check if the vector, v, is feasible.
int dimension() const
Return dimension of the vector space.
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
H1BoundConstraint(std::vector< Real > &l, std::vector< Real > &u, const Teuchos::RCP< BurgersFEM< Real > > &fem, Real scale=1.0)
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 pruneActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &x, Real eps)
Set variables to zero if they correspond to the -active set.
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 pruneLowerActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &g, const ROL::Vector< Real > &x, Real eps)
Set variables to zero if they correspond to the lower -binding set.
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 .
Teuchos::RCP< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
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
H1VectorPrimal< Real > PrimalStateVector
Teuchos::RCP< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
Real compute_L2_dot(const std::vector< Real > &x, const std::vector< Real > &y) const
void pruneLowerActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &x, Real eps)
Set variables to zero if they correspond to the lower -active set.
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_
void apply_inverse_mass(std::vector< Real > &Mu, const std::vector< Real > &u) const
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)
void setVectorToUpperBound(ROL::Vector< Real > &u)
Set the input vector to the upper bound.
H1VectorDual< Real > PrimalConstraintVector
Teuchos::RCP< BurgersFEM< Real > > fem_
void update(std::vector< Real > &u, const std::vector< Real > &s, const Real alpha=1.0) const
void pruneLowerActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &g, const ROL::Vector< Real > &x, Real eps)
Set variables to zero if they correspond to the lower -binding set.
void scale(const Real alpha)
Compute where .
Teuchos::RCP< BurgersFEM< Real > > fem_
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 cast_vector(Teuchos::RCP< std::vector< Real > > &xvec, ROL::Vector< Real > &x) const
H1VectorPrimal< Real > DualConstraintVector
void setVectorToLowerBound(ROL::Vector< Real > &l)
Set the input vector to the lower bound.
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)
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(Teuchos::RCP< const std::vector< Real > > &xvec, const ROL::Vector< Real > &x) const
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 .
H1VectorDual< Real > DualStateVector
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 apply_control_jacobian(std::vector< Real > &jv, const std::vector< Real > &v, const std::vector< Real > &u, const std::vector< Real > &z) const
void test_inverse_H1(std::ostream &outStream=std::cout)
void value(ROL::Vector< Real > &c, const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, Real &tol)
Evaluate the constraint operator at .
const ROL::Vector< Real > & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
H1VectorDual(const Teuchos::RCP< std::vector< Real > > &vec, const Teuchos::RCP< BurgersFEM< Real > > &fem)
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 pruneUpperActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &g, const ROL::Vector< Real > &x, Real eps)
Set variables to zero if they correspond to the upper -binding set.
void pruneActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &g, const ROL::Vector< Real > &x, Real eps)
Set variables to zero if they correspond to the -binding set.
H1VectorPrimal< Real > PrimalStateVector
L2VectorPrimal(const Teuchos::RCP< std::vector< Real > > &vec, const Teuchos::RCP< BurgersFEM< Real > > &fem)
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
EqualityConstraint_BurgersControl(Teuchos::RCP< BurgersFEM< Real > > &fem, bool useHessian=true)
Teuchos::RCP< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
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 apply_H1(std::vector< Real > &Mu, const std::vector< Real > &u) const
void plus(const ROL::Vector< Real > &x)
Compute , where .
virtual void set(const Vector &x)
Set where .
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 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)
int dimension() const
Return dimension of the vector space.
L2VectorDual< Real > DualControlVector
const ROL::Vector< Real > & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
Teuchos::RCP< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
H1VectorPrimal(const Teuchos::RCP< std::vector< Real > > &vec, const Teuchos::RCP< BurgersFEM< Real > > &fem)
Real norm() const
Returns where .
void scale(std::vector< Real > &u, const Real alpha=0.0) const
Teuchos::RCP< std::vector< Real > > getVector()
const ROL::Vector< Real > & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
Teuchos::RCP< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
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 ...
L2BoundConstraint(std::vector< Real > &l, std::vector< Real > &u, const Teuchos::RCP< BurgersFEM< Real > > &fem, Real scale=1.0)
void pruneUpperActive(ROL::Vector< Real > &v, const ROL::Vector< Real > &x, Real eps)
Set variables to zero if they correspond to the upper -active set.
Teuchos::RCP< const std::vector< Real > > getVector() const
void setVectorToLowerBound(ROL::Vector< Real > &l)
Set the input vector to the lower bound.
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
Teuchos::RCP< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
void setVectorToUpperBound(ROL::Vector< Real > &u)
Set the input vector to the upper bound.