56 #include "Teuchos_LAPACK.hpp" 84 void update(std::vector<Real> &u,
const std::vector<Real> &s,
const Real alpha=1.0)
const {
85 for (
unsigned i=0; i<u.size(); i++) {
90 void axpy(std::vector<Real> &out,
const Real a,
const std::vector<Real> &x,
const std::vector<Real> &y)
const {
91 for (
unsigned i=0; i < x.size(); i++) {
92 out[i] = a*x[i] + y[i];
96 void scale(std::vector<Real> &u,
const Real alpha=0.0)
const {
97 for (
unsigned i=0; i<u.size(); i++) {
102 void linear_solve(std::vector<Real> &u, std::vector<Real> &dl, std::vector<Real> &d, std::vector<Real> &du,
103 const std::vector<Real> &r,
const bool transpose =
false)
const {
104 if ( r.size() == 1 ) {
105 u.resize(1,r[0]/d[0]);
107 else if ( r.size() == 2 ) {
109 Real det = d[0]*d[1] - dl[0]*du[0];
110 u[0] = (d[1]*r[0] - du[0]*r[1])/det;
111 u[1] = (d[0]*r[1] - dl[0]*r[0])/det;
114 u.assign(r.begin(),r.end());
116 Teuchos::LAPACK<int,Real> lp;
117 std::vector<Real> du2(r.size()-2,0.0);
118 std::vector<int> ipiv(r.size(),0);
123 lp.GTTRF(dim,&dl[0],&d[0],&du[0],&du2[0],&ipiv[0],&info);
128 lp.GTTRS(trans,dim,nhrs,&dl[0],&d[0],&du[0],&du2[0],&ipiv[0],&u[0],ldb,&info);
133 BurgersFEM(
int nx = 128, Real nl = 1.0, Real cH1 = 1.0, Real cL2 = 1.0)
134 : nx_(nx), dx_(1.0/((Real)nx+1.0)), nl_(nl), cH1_(cH1), cL2_(cL2) {}
137 nu_ = std::pow(10.0,nu-2.0);
155 Real
compute_L2_dot(
const std::vector<Real> &x,
const std::vector<Real> &y)
const {
157 Real c = (((int)x.size()==
nx_) ? 4.0 : 2.0);
158 for (
unsigned i=0; i<x.size(); i++) {
160 ip += dx_/6.0*(c*x[i] + x[i+1])*y[i];
162 else if ( i == x.size()-1 ) {
163 ip += dx_/6.0*(x[i-1] + c*x[i])*y[i];
166 ip += dx_/6.0*(x[i-1] + 4.0*x[i] + x[i+1])*y[i];
178 void apply_mass(std::vector<Real> &Mu,
const std::vector<Real> &u )
const {
179 Mu.resize(u.size(),0.0);
180 Real c = (((int)u.size()==
nx_) ? 4.0 : 2.0);
181 for (
unsigned i=0; i<u.size(); i++) {
183 Mu[i] = dx_/6.0*(c*u[i] + u[i+1]);
185 else if ( i == u.size()-1 ) {
186 Mu[i] = dx_/6.0*(u[i-1] + c*u[i]);
189 Mu[i] = dx_/6.0*(u[i-1] + 4.0*u[i] + u[i+1]);
196 unsigned nx = u.size();
198 std::vector<Real> dl(nx-1,dx_/6.0);
199 std::vector<Real> d(nx,2.0*dx_/3.0);
200 std::vector<Real> du(nx-1,dx_/6.0);
201 if ( (
int)nx != nx_ ) {
210 std::vector<Real> u(nx_,0.0), Mu(nx_,0.0), iMMu(nx_,0.0), diff(nx_,0.0);
211 for (
int i = 0; i <
nx_; i++) {
212 u[i] = 2.0*(Real)rand()/(Real)RAND_MAX - 1.0;
216 axpy(diff,-1.0,iMMu,u);
219 outStream <<
"Test Inverse State Mass Matrix\n";
220 outStream <<
" ||u - inv(M)Mu|| = " << error <<
"\n";
221 outStream <<
" ||u|| = " << normu <<
"\n";
222 outStream <<
" Relative Error = " << error/normu <<
"\n";
225 u.resize(nx_+2,0.0); Mu.resize(nx_+2,0.0); iMMu.resize(nx_+2,0.0); diff.resize(nx_+2,0.0);
226 for (
int i = 0; i < nx_+2; i++) {
227 u[i] = 2.0*(Real)rand()/(Real)RAND_MAX - 1.0;
231 axpy(diff,-1.0,iMMu,u);
234 outStream <<
"Test Inverse Control Mass Matrix\n";
235 outStream <<
" ||z - inv(M)Mz|| = " << error <<
"\n";
236 outStream <<
" ||z|| = " << normu <<
"\n";
237 outStream <<
" Relative Error = " << error/normu <<
"\n";
245 Real
compute_H1_dot(
const std::vector<Real> &x,
const std::vector<Real> &y)
const {
247 for (
int i=0; i<
nx_; i++) {
249 ip += cL2_*dx_/6.0*(4.0*x[i] + x[i+1])*y[i];
250 ip += cH1_*(2.0*x[i] - x[i+1])/dx_*y[i];
252 else if ( i == nx_-1 ) {
253 ip += cL2_*dx_/6.0*(x[i-1] + 4.0*x[i])*y[i];
254 ip += cH1_*(2.0*x[i] - x[i-1])/dx_*y[i];
257 ip += cL2_*dx_/6.0*(x[i-1] + 4.0*x[i] + x[i+1])*y[i];
258 ip += cH1_*(2.0*x[i] - x[i-1] - x[i+1])/dx_*y[i];
270 void apply_H1(std::vector<Real> &Mu,
const std::vector<Real> &u )
const {
272 for (
int i=0; i<
nx_; i++) {
274 Mu[i] = cL2_*dx_/6.0*(4.0*u[i] + u[i+1])
275 + cH1_*(2.0*u[i] - u[i+1])/
dx_;
277 else if ( i == nx_-1 ) {
278 Mu[i] = cL2_*dx_/6.0*(u[i-1] + 4.0*u[i])
279 + cH1_*(2.0*u[i] - u[i-1])/
dx_;
282 Mu[i] = cL2_*dx_/6.0*(u[i-1] + 4.0*u[i] + u[i+1])
283 + cH1_*(2.0*u[i] - u[i-1] - u[i+1])/
dx_;
291 std::vector<Real> dl(nx_-1,cL2_*dx_/6.0 - cH1_/dx_);
292 std::vector<Real> d(nx_,2.0*(cL2_*dx_/3.0 + cH1_/dx_));
293 std::vector<Real> du(nx_-1,cL2_*dx_/6.0 - cH1_/dx_);
298 std::vector<Real> u(nx_,0.0), Mu(nx_,0.0), iMMu(nx_,0.0), diff(nx_,0.0);
299 for (
int i = 0; i <
nx_; i++) {
300 u[i] = 2.0*(Real)rand()/(Real)RAND_MAX - 1.0;
304 axpy(diff,-1.0,iMMu,u);
307 outStream <<
"Test Inverse State H1 Matrix\n";
308 outStream <<
" ||u - inv(M)Mu|| = " << error <<
"\n";
309 outStream <<
" ||u|| = " << normu <<
"\n";
310 outStream <<
" Relative Error = " << error/normu <<
"\n";
319 const std::vector<Real> &z)
const {
322 for (
int i=0; i<
nx_; i++) {
325 r[i] = nu_/dx_*(2.0*u[i]-u[i+1]);
328 r[i] = nu_/dx_*(2.0*u[i]-u[i-1]);
331 r[i] = nu_/dx_*(2.0*u[i]-u[i-1]-u[i+1]);
335 r[i] += nl_*u[i+1]*(u[i]+u[i+1])/6.0;
338 r[i] -= nl_*u[i-1]*(u[i-1]+u[i])/6.0;
341 r[i] -= dx_/6.0*(z[i]+4.0*z[i+1]+z[i+2]);
346 r[0] -= nl_*(u0_*u[ 0]/6.0 + u0_*u0_/6.0) + nu_*u0_/
dx_;
347 r[nx_-1] += nl_*(u1_*u[nx_-1]/6.0 + u1_*u1_/6.0) - nu_*u1_/
dx_;
355 std::vector<Real> &d,
356 std::vector<Real> &du,
357 const std::vector<Real> &u)
const {
360 d.resize(nx_,nu_*2.0/dx_);
362 dl.resize(nx_-1,-nu_/dx_);
364 du.resize(nx_-1,-nu_/dx_);
366 for (
int i=0; i<
nx_; i++) {
368 dl[i] += nl_*(-2.0*u[i]-u[i+1])/6.0;
369 d[i] += nl_*u[i+1]/6.0;
372 d[i] -= nl_*u[i-1]/6.0;
373 du[i-1] += nl_*(u[i-1]+2.0*u[i])/6.0;
377 d[ 0] -= nl_*u0_/6.0;
378 d[nx_-1] += nl_*u1_/6.0;
383 const std::vector<Real> &v,
384 const std::vector<Real> &u,
385 const std::vector<Real> &z)
const {
387 for (
int i = 0; i <
nx_; i++) {
388 jv[i] = nu_/dx_*2.0*v[i];
390 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]);
393 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]);
396 jv[ 0] -= nl_*u0_/6.0*v[0];
397 jv[nx_-1] += nl_*u1_/6.0*v[nx_-1];
402 const std::vector<Real> &v,
403 const std::vector<Real> &u,
404 const std::vector<Real> &z)
const {
406 std::vector<Real> d(nx_,0.0);
407 std::vector<Real> dl(nx_-1,0.0);
408 std::vector<Real> du(nx_-1,0.0);
416 const std::vector<Real> &v,
417 const std::vector<Real> &u,
418 const std::vector<Real> &z)
const {
420 for (
int i = 0; i <
nx_; i++) {
421 ajv[i] = nu_/dx_*2.0*v[i];
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[i] += -nu_/dx_*v[i+1]+nl_*(u[i+1]/6.0*v[i]
428 -(u[i+1]+2.0*u[i])/6.0*v[i+1]);
431 ajv[ 0] -= nl_*u0_/6.0*v[0];
432 ajv[nx_-1] += nl_*u1_/6.0*v[nx_-1];
437 const std::vector<Real> &v,
438 const std::vector<Real> &u,
439 const std::vector<Real> &z)
const {
441 std::vector<Real> d(nx_,0.0);
442 std::vector<Real> du(nx_-1,0.0);
443 std::vector<Real> dl(nx_-1,0.0);
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++) {
459 jv[i] = -dx_/6.0*(v[i]+4.0*v[i+1]+v[i+2]);
465 const std::vector<Real> &v,
466 const std::vector<Real> &u,
467 const std::vector<Real> &z)
const {
468 for (
int i=0; i<nx_+2; i++) {
470 jv[i] = -dx_/6.0*v[i];
473 jv[i] = -dx_/6.0*(4.0*v[i-1]+v[i]);
475 else if ( i == nx_ ) {
476 jv[i] = -dx_/6.0*(4.0*v[i-1]+v[i-2]);
478 else if ( i == nx_+1 ) {
479 jv[i] = -dx_/6.0*v[i-2];
482 jv[i] = -dx_/6.0*(v[i-2]+4.0*v[i-1]+v[i]);
491 const std::vector<Real> &w,
492 const std::vector<Real> &v,
493 const std::vector<Real> &u,
494 const std::vector<Real> &z)
const {
495 for (
int i=0; i<
nx_; i++) {
499 ahwv[i] += (w[i]*v[i+1] - w[i+1]*(2.0*v[i]+v[i+1]))/6.0;
502 ahwv[i] += (w[i-1]*(v[i-1]+2.0*v[i]) - w[i]*v[i-1])/6.0;
508 const std::vector<Real> &w,
509 const std::vector<Real> &v,
510 const std::vector<Real> &u,
511 const std::vector<Real> &z) {
512 ahwv.assign(u.size(),0.0);
515 const std::vector<Real> &w,
516 const std::vector<Real> &v,
517 const std::vector<Real> &u,
518 const std::vector<Real> &z) {
519 ahwv.assign(z.size(),0.0);
522 const std::vector<Real> &w,
523 const std::vector<Real> &v,
524 const std::vector<Real> &u,
525 const std::vector<Real> &z) {
526 ahwv.assign(z.size(),0.0);
533 Teuchos::RCP<std::vector<Real> > vec_;
534 Teuchos::RCP<BurgersFEM<Real> > fem_;
536 mutable Teuchos::RCP<L2VectorDual<Real> > dual_vec_;
541 : vec_(vec), fem_(fem), dual_vec_(Teuchos::null) {}
545 const std::vector<Real>& xval = *ex.
getVector();
546 std::copy(xval.begin(),xval.end(),vec_->begin());
551 const std::vector<Real>& xval = *ex.
getVector();
552 unsigned dimension = vec_->size();
553 for (
unsigned i=0; i<dimension; i++) {
554 (*vec_)[i] += xval[i];
559 unsigned dimension = vec_->size();
560 for (
unsigned i=0; i<dimension; i++) {
567 const std::vector<Real>& xval = *ex.
getVector();
568 return fem_->compute_L2_dot(xval,*vec_);
573 val = std::sqrt( dot(*
this) );
577 Teuchos::RCP<ROL::Vector<Real> >
clone()
const {
578 return Teuchos::rcp(
new L2VectorPrimal( Teuchos::rcp(
new std::vector<Real>(vec_->size(),0.0)),fem_));
581 Teuchos::RCP<const std::vector<Real> >
getVector()
const {
589 Teuchos::RCP<ROL::Vector<Real> >
basis(
const int i )
const {
590 Teuchos::RCP<L2VectorPrimal> e
591 = Teuchos::rcp(
new L2VectorPrimal( Teuchos::rcp(
new std::vector<Real>(vec_->size(),0.0)),fem_));
592 (*e->getVector())[i] = 1.0;
602 Teuchos::rcp(
new std::vector<Real>(*vec_)),fem_));
604 fem_->apply_mass(*(Teuchos::rcp_const_cast<std::vector<Real> >(dual_vec_->getVector())),*vec_);
613 Teuchos::RCP<std::vector<Real> > vec_;
614 Teuchos::RCP<BurgersFEM<Real> > fem_;
616 mutable Teuchos::RCP<L2VectorPrimal<Real> > dual_vec_;
621 : vec_(vec), fem_(fem), dual_vec_(Teuchos::null) {}
625 const std::vector<Real>& xval = *ex.
getVector();
626 std::copy(xval.begin(),xval.end(),vec_->begin());
631 const std::vector<Real>& xval = *ex.
getVector();
632 unsigned dimension = vec_->size();
633 for (
unsigned i=0; i<dimension; i++) {
634 (*vec_)[i] += xval[i];
639 unsigned dimension = vec_->size();
640 for (
unsigned i=0; i<dimension; i++) {
647 const std::vector<Real>& xval = *ex.
getVector();
648 unsigned dimension = vec_->size();
649 std::vector<Real> Mx(dimension,0.0);
650 fem_->apply_inverse_mass(Mx,xval);
652 for (
unsigned i = 0; i < dimension; i++) {
653 val += Mx[i]*(*vec_)[i];
660 val = std::sqrt( dot(*
this) );
664 Teuchos::RCP<ROL::Vector<Real> >
clone()
const {
665 return Teuchos::rcp(
new L2VectorDual( Teuchos::rcp(
new std::vector<Real>(vec_->size(),0.0)),fem_));
668 Teuchos::RCP<const std::vector<Real> >
getVector()
const {
676 Teuchos::RCP<ROL::Vector<Real> >
basis(
const int i )
const {
677 Teuchos::RCP<L2VectorDual> e
678 = Teuchos::rcp(
new L2VectorDual( Teuchos::rcp(
new std::vector<Real>(vec_->size(),0.0)),fem_));
679 (*e->getVector())[i] = 1.0;
689 Teuchos::rcp(
new std::vector<Real>(*vec_)),fem_));
691 fem_->apply_inverse_mass(*(Teuchos::rcp_const_cast<std::vector<Real> >(dual_vec_->getVector())),*vec_);
700 Teuchos::RCP<std::vector<Real> > vec_;
701 Teuchos::RCP<BurgersFEM<Real> > fem_;
703 mutable Teuchos::RCP<H1VectorDual<Real> > dual_vec_;
708 : vec_(vec), fem_(fem), dual_vec_(Teuchos::null) {}
712 const std::vector<Real>& xval = *ex.
getVector();
713 std::copy(xval.begin(),xval.end(),vec_->begin());
718 const std::vector<Real>& xval = *ex.
getVector();
719 unsigned dimension = vec_->size();
720 for (
unsigned i=0; i<dimension; i++) {
721 (*vec_)[i] += xval[i];
726 unsigned dimension = vec_->size();
727 for (
unsigned i=0; i<dimension; i++) {
734 const std::vector<Real>& xval = *ex.
getVector();
735 return fem_->compute_H1_dot(xval,*vec_);
740 val = std::sqrt( dot(*
this) );
744 Teuchos::RCP<ROL::Vector<Real> >
clone()
const {
745 return Teuchos::rcp(
new H1VectorPrimal( Teuchos::rcp(
new std::vector<Real>(vec_->size(),0.0)),fem_));
748 Teuchos::RCP<const std::vector<Real> >
getVector()
const {
756 Teuchos::RCP<ROL::Vector<Real> >
basis(
const int i )
const {
757 Teuchos::RCP<H1VectorPrimal> e
758 = Teuchos::rcp(
new H1VectorPrimal( Teuchos::rcp(
new std::vector<Real>(vec_->size(),0.0)),fem_));
759 (*e->getVector())[i] = 1.0;
769 Teuchos::rcp(
new std::vector<Real>(*vec_)),fem_));
771 fem_->apply_H1(*(Teuchos::rcp_const_cast<std::vector<Real> >(dual_vec_->getVector())),*vec_);
780 Teuchos::RCP<std::vector<Real> > vec_;
781 Teuchos::RCP<BurgersFEM<Real> > fem_;
783 mutable Teuchos::RCP<H1VectorPrimal<Real> > dual_vec_;
788 : vec_(vec), fem_(fem), dual_vec_(Teuchos::null) {}
792 const std::vector<Real>& xval = *ex.
getVector();
793 std::copy(xval.begin(),xval.end(),vec_->begin());
798 const std::vector<Real>& xval = *ex.
getVector();
799 unsigned dimension = vec_->size();
800 for (
unsigned i=0; i<dimension; i++) {
801 (*vec_)[i] += xval[i];
806 unsigned dimension = vec_->size();
807 for (
unsigned i=0; i<dimension; i++) {
814 const std::vector<Real>& xval = *ex.
getVector();
815 unsigned dimension = vec_->size();
816 std::vector<Real> Mx(dimension,0.0);
817 fem_->apply_inverse_H1(Mx,xval);
819 for (
unsigned i = 0; i < dimension; i++) {
820 val += Mx[i]*(*vec_)[i];
827 val = std::sqrt( dot(*
this) );
831 Teuchos::RCP<ROL::Vector<Real> >
clone()
const {
832 return Teuchos::rcp(
new H1VectorDual( Teuchos::rcp(
new std::vector<Real>(vec_->size(),0.0)),fem_));
835 Teuchos::RCP<const std::vector<Real> >
getVector()
const {
843 Teuchos::RCP<ROL::Vector<Real> >
basis(
const int i )
const {
844 Teuchos::RCP<H1VectorDual> e
845 = Teuchos::rcp(
new H1VectorDual( Teuchos::rcp(
new std::vector<Real>(vec_->size(),0.0)),fem_));
846 (*e->getVector())[i] = 1.0;
856 Teuchos::rcp(
new std::vector<Real>(*vec_)),fem_));
858 fem_->apply_inverse_H1(*(Teuchos::rcp_const_cast<std::vector<Real> >(dual_vec_->getVector())),*vec_);
877 Teuchos::RCP<BurgersFEM<Real> > fem_;
882 : fem_(fem), useHessian_(useHessian) {}
886 Teuchos::RCP<std::vector<Real> > cp =
887 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<PrimalConstraintVector>(c)).getVector());
888 Teuchos::RCP<const std::vector<Real> > up =
889 (Teuchos::dyn_cast<PrimalStateVector>(
const_cast<ROL::Vector<Real> &
>(u))).getVector();
890 Teuchos::RCP<const std::vector<Real> > zp =
891 (Teuchos::dyn_cast<PrimalControlVector>(
const_cast<ROL::Vector<Real> &
>(z))).getVector();
893 const std::vector<Real> param
895 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
897 fem_->compute_residual(*cp,*up,*zp);
902 Teuchos::RCP<std::vector<Real> > jvp =
903 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<PrimalConstraintVector>(jv)).getVector());
904 Teuchos::RCP<const std::vector<Real> > vp =
905 (Teuchos::dyn_cast<PrimalStateVector>(
const_cast<ROL::Vector<Real> &
>(v))).getVector();
906 Teuchos::RCP<const std::vector<Real> > up =
907 (Teuchos::dyn_cast<PrimalStateVector>(
const_cast<ROL::Vector<Real> &
>(u))).getVector();
908 Teuchos::RCP<const std::vector<Real> > zp =
909 (Teuchos::dyn_cast<PrimalControlVector>(
const_cast<ROL::Vector<Real> &
>(z))).getVector();
911 const std::vector<Real> param
913 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
915 fem_->apply_pde_jacobian(*jvp,*vp,*up,*zp);
920 Teuchos::RCP<std::vector<Real> > jvp =
921 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<PrimalConstraintVector>(jv)).getVector());
922 Teuchos::RCP<const std::vector<Real> > vp =
923 (Teuchos::dyn_cast<PrimalControlVector>(
const_cast<ROL::Vector<Real> &
>(v))).getVector();
924 Teuchos::RCP<const std::vector<Real> > up =
925 (Teuchos::dyn_cast<PrimalStateVector>(
const_cast<ROL::Vector<Real> &
>(u))).getVector();
926 Teuchos::RCP<const std::vector<Real> > zp =
927 (Teuchos::dyn_cast<PrimalControlVector>(
const_cast<ROL::Vector<Real> &
>(z))).getVector();
929 const std::vector<Real> param
931 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
933 fem_->apply_control_jacobian(*jvp,*vp,*up,*zp);
938 Teuchos::RCP<std::vector<Real> > ijvp =
939 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<PrimalStateVector>(ijv)).getVector());
940 Teuchos::RCP<const std::vector<Real> > vp =
941 (Teuchos::dyn_cast<PrimalConstraintVector>(
const_cast<ROL::Vector<Real> &
>(v))).getVector();
942 Teuchos::RCP<const std::vector<Real> > up =
943 (Teuchos::dyn_cast<PrimalStateVector>(
const_cast<ROL::Vector<Real> &
>(u))).getVector();
944 Teuchos::RCP<const std::vector<Real> > zp =
945 (Teuchos::dyn_cast<PrimalControlVector>(
const_cast<ROL::Vector<Real> &
>(z))).getVector();
947 const std::vector<Real> param
949 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
951 fem_->apply_inverse_pde_jacobian(*ijvp,*vp,*up,*zp);
956 Teuchos::RCP<std::vector<Real> > jvp =
957 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<DualStateVector>(ajv)).getVector());
958 Teuchos::RCP<const std::vector<Real> > vp =
959 (Teuchos::dyn_cast<DualConstraintVector>(
const_cast<ROL::Vector<Real> &
>(v))).getVector();
960 Teuchos::RCP<const std::vector<Real> > up =
961 (Teuchos::dyn_cast<PrimalStateVector>(
const_cast<ROL::Vector<Real> &
>(u))).getVector();
962 Teuchos::RCP<const std::vector<Real> > zp =
963 (Teuchos::dyn_cast<PrimalControlVector>(
const_cast<ROL::Vector<Real> &
>(z))).getVector();
965 const std::vector<Real> param
967 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
969 fem_->apply_adjoint_pde_jacobian(*jvp,*vp,*up,*zp);
974 Teuchos::RCP<std::vector<Real> > jvp =
975 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<DualControlVector>(jv)).getVector());
976 Teuchos::RCP<const std::vector<Real> > vp =
977 (Teuchos::dyn_cast<DualConstraintVector>(
const_cast<ROL::Vector<Real> &
>(v))).getVector();
978 Teuchos::RCP<const std::vector<Real> > up =
979 (Teuchos::dyn_cast<PrimalStateVector>(
const_cast<ROL::Vector<Real> &
>(u))).getVector();
980 Teuchos::RCP<const std::vector<Real> > zp =
981 (Teuchos::dyn_cast<PrimalControlVector>(
const_cast<ROL::Vector<Real> &
>(z))).getVector();
983 const std::vector<Real> param
985 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
987 fem_->apply_adjoint_control_jacobian(*jvp,*vp,*up,*zp);
992 Teuchos::RCP<std::vector<Real> > iajvp =
993 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<DualConstraintVector>(iajv)).getVector());
994 Teuchos::RCP<const std::vector<Real> > vp =
995 (Teuchos::dyn_cast<DualStateVector>(
const_cast<ROL::Vector<Real> &
>(v))).getVector();
996 Teuchos::RCP<const std::vector<Real> > up =
997 (Teuchos::dyn_cast<PrimalStateVector>(
const_cast<ROL::Vector<Real> &
>(u))).getVector();
998 Teuchos::RCP<const std::vector<Real> > zp =
999 (Teuchos::dyn_cast<PrimalControlVector>(
const_cast<ROL::Vector<Real> &
>(z))).getVector();
1001 const std::vector<Real> param
1003 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
1005 fem_->apply_inverse_adjoint_pde_jacobian(*iajvp,*vp,*up,*zp);
1010 if ( useHessian_ ) {
1011 Teuchos::RCP<std::vector<Real> > ahwvp =
1012 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<DualStateVector>(ahwv)).getVector());
1013 Teuchos::RCP<const std::vector<Real> > wp =
1014 (Teuchos::dyn_cast<DualConstraintVector>(
const_cast<ROL::Vector<Real> &
>(w))).getVector();
1015 Teuchos::RCP<const std::vector<Real> > vp =
1016 (Teuchos::dyn_cast<PrimalStateVector>(
const_cast<ROL::Vector<Real> &
>(v))).getVector();
1017 Teuchos::RCP<const std::vector<Real> > up =
1018 (Teuchos::dyn_cast<PrimalStateVector>(
const_cast<ROL::Vector<Real> &
>(u))).getVector();
1019 Teuchos::RCP<const std::vector<Real> > zp =
1020 (Teuchos::dyn_cast<PrimalControlVector>(
const_cast<ROL::Vector<Real> &
>(z))).getVector();
1022 const std::vector<Real> param
1024 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
1026 fem_->apply_adjoint_pde_hessian(*ahwvp,*wp,*vp,*up,*zp);
1035 if ( useHessian_ ) {
1036 Teuchos::RCP<std::vector<Real> > ahwvp =
1037 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<DualControlVector>(ahwv)).getVector());
1038 Teuchos::RCP<const std::vector<Real> > wp =
1039 (Teuchos::dyn_cast<DualConstraintVector>(
const_cast<ROL::Vector<Real> &
>(w))).getVector();
1040 Teuchos::RCP<const std::vector<Real> > vp =
1041 (Teuchos::dyn_cast<PrimalStateVector>(
const_cast<ROL::Vector<Real> &
>(v))).getVector();
1042 Teuchos::RCP<const std::vector<Real> > up =
1043 (Teuchos::dyn_cast<PrimalStateVector>(
const_cast<ROL::Vector<Real> &
>(u))).getVector();
1044 Teuchos::RCP<const std::vector<Real> > zp =
1045 (Teuchos::dyn_cast<PrimalControlVector>(
const_cast<ROL::Vector<Real> &
>(z))).getVector();
1047 const std::vector<Real> param
1049 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
1051 fem_->apply_adjoint_control_pde_hessian(*ahwvp,*wp,*vp,*up,*zp);
1059 if ( useHessian_ ) {
1060 Teuchos::RCP<std::vector<Real> > ahwvp =
1061 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<DualStateVector>(ahwv)).getVector());
1062 Teuchos::RCP<const std::vector<Real> > wp =
1063 (Teuchos::dyn_cast<DualConstraintVector>(
const_cast<ROL::Vector<Real> &
>(w))).getVector();
1064 Teuchos::RCP<const std::vector<Real> > vp =
1065 (Teuchos::dyn_cast<PrimalControlVector>(
const_cast<ROL::Vector<Real> &
>(v))).getVector();
1066 Teuchos::RCP<const std::vector<Real> > up =
1067 (Teuchos::dyn_cast<PrimalStateVector>(
const_cast<ROL::Vector<Real> &
>(u))).getVector();
1068 Teuchos::RCP<const std::vector<Real> > zp =
1069 (Teuchos::dyn_cast<PrimalControlVector>(
const_cast<ROL::Vector<Real> &
>(z))).getVector();
1071 const std::vector<Real> param
1073 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
1075 fem_->apply_adjoint_pde_control_hessian(*ahwvp,*wp,*vp,*up,*zp);
1083 if ( useHessian_ ) {
1084 Teuchos::RCP<std::vector<Real> > ahwvp =
1085 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<DualControlVector>(ahwv)).getVector());
1086 Teuchos::RCP<const std::vector<Real> > wp =
1087 (Teuchos::dyn_cast<DualConstraintVector>(
const_cast<ROL::Vector<Real> &
>(w))).getVector();
1088 Teuchos::RCP<const std::vector<Real> > vp =
1089 (Teuchos::dyn_cast<PrimalControlVector>(
const_cast<ROL::Vector<Real> &
>(v))).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> > zp =
1093 (Teuchos::dyn_cast<PrimalControlVector>(
const_cast<ROL::Vector<Real> &
>(z))).getVector();
1095 const std::vector<Real> param
1097 fem_->set_problem_data(param[0],param[1],param[2],param[3]);
1099 fem_->apply_adjoint_control_hessian(*ahwvp,*wp,*vp,*up,*zp);
1107 template<
class Real>
1118 Teuchos::RCP<BurgersFEM<Real> > fem_;
1119 Teuchos::RCP<ROL::Vector<Real> > ud_;
1120 Teuchos::RCP<ROL::Vector<Real> > diff_;
1125 Real alpha = 1.e-4) : alpha_(alpha), fem_(fem), ud_(ud) {
1126 diff_ = ud_->clone();
1130 Teuchos::RCP<const std::vector<Real> > up =
1131 (Teuchos::dyn_cast<PrimalStateVector>(
const_cast<ROL::Vector<Real> &
>(u))).getVector();
1132 Teuchos::RCP<const std::vector<Real> > zp =
1133 (Teuchos::dyn_cast<PrimalControlVector>(
const_cast<ROL::Vector<Real> &
>(z))).getVector();
1134 Teuchos::RCP<const std::vector<Real> > udp =
1137 std::vector<Real> diff(udp->size(),0.0);
1138 for (
unsigned i = 0; i < udp->size(); i++) {
1139 diff[i] = (*up)[i] - (*udp)[i];
1141 return 0.5*(fem_->compute_L2_dot(diff,diff) + alpha_*fem_->compute_L2_dot(*zp,*zp));
1145 Teuchos::RCP<std::vector<Real> > gp =
1146 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<DualStateVector>(g)).getVector());
1147 Teuchos::RCP<const std::vector<Real> > up =
1148 (Teuchos::dyn_cast<PrimalStateVector>(
const_cast<ROL::Vector<Real> &
>(u))).getVector();
1149 Teuchos::RCP<const std::vector<Real> > udp =
1152 std::vector<Real> diff(udp->size(),0.0);
1153 for (
unsigned i = 0; i < udp->size(); i++) {
1154 diff[i] = (*up)[i] - (*udp)[i];
1156 fem_->apply_mass(*gp,diff);
1160 Teuchos::RCP<std::vector<Real> > gp =
1161 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<DualControlVector>(g)).getVector());
1162 Teuchos::RCP<const std::vector<Real> > zp =
1163 (Teuchos::dyn_cast<PrimalControlVector>(
const_cast<ROL::Vector<Real> &
>(z))).getVector();
1165 fem_->apply_mass(*gp,*zp);
1166 for (
unsigned i = 0; i < zp->size(); i++) {
1173 Teuchos::RCP<std::vector<Real> > hvp =
1174 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<DualStateVector>(hv)).getVector());
1175 Teuchos::RCP<const std::vector<Real> > vp =
1176 (Teuchos::dyn_cast<PrimalStateVector>(
const_cast<ROL::Vector<Real> &
>(v))).getVector();
1178 fem_->apply_mass(*hvp,*vp);
1193 Teuchos::RCP<std::vector<Real> > hvp =
1194 Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<DualControlVector>(hv)).getVector());
1195 Teuchos::RCP<const std::vector<Real> > vp =
1196 (Teuchos::dyn_cast<PrimalControlVector>(
const_cast<ROL::Vector<Real> &
>(v))).getVector();
1198 fem_->apply_mass(*hvp,*vp);
1199 for (
unsigned i = 0; i < vp->size(); i++) {
1200 (*hvp)[i] *= alpha_;
1205 template<
class Real>
1209 std::vector<Real> x_lo_;
1210 std::vector<Real> x_up_;
1213 Teuchos::RCP<BurgersFEM<Real> > fem_;
1218 xvec = Teuchos::rcp_const_cast<std::vector<Real> >(
1221 catch (std::exception &e) {
1222 xvec = Teuchos::rcp_const_cast<std::vector<Real> >(
1233 catch (std::exception &e) {
1239 void axpy(std::vector<Real> &out,
const Real a,
1240 const std::vector<Real> &x,
const std::vector<Real> &y)
const{
1241 out.resize(dim_,0.0);
1242 for (
unsigned i = 0; i < dim_; i++) {
1243 out[i] = a*x[i] + y[i];
1248 for (
int i = 0; i < dim_; i++ ) {
1249 x[i] = std::max(x_lo_[i],std::min(x_up_[i],x[i]));
1256 : x_lo_(l), x_up_(u), scale_(
scale), fem_(fem) {
1257 dim_ = x_lo_.size();
1258 for (
int i = 0; i < dim_; i++ ) {
1260 min_diff_ = x_up_[i] - x_lo_[i];
1263 min_diff_ = ( (min_diff_ < (x_up_[i] - x_lo_[i])) ? min_diff_ : (x_up_[i] - x_lo_[i]) );
1270 Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1273 for (
int i = 0; i < dim_; i++ ) {
1274 if ( (*ex)[i] >= x_lo_[i] && (*ex)[i] <= x_up_[i] ) { cnt *= 1; }
1277 if ( cnt == 0 ) { val =
false; }
1282 Teuchos::RCP<std::vector<Real> > ex; cast_vector(ex,x);
1287 Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1288 Teuchos::RCP<std::vector<Real> > ev; cast_vector(ev,v);
1289 Real epsn = std::min(scale_*eps,min_diff_);
1290 for (
int i = 0; i < dim_; i++ ) {
1291 if ( ((*ex)[i] <= x_lo_[i]+epsn) ) {
1298 Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1299 Teuchos::RCP<std::vector<Real> > ev; cast_vector(ev,v);
1300 Real epsn = std::min(scale_*eps,min_diff_);
1301 for (
int i = 0; i < dim_; i++ ) {
1302 if ( ((*ex)[i] >= x_up_[i]-epsn) ) {
1309 Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1310 Teuchos::RCP<std::vector<Real> > ev; cast_vector(ev,v);
1311 Real epsn = std::min(scale_*eps,min_diff_);
1312 for (
int i = 0; i < dim_; i++ ) {
1313 if ( ((*ex)[i] <= x_lo_[i]+epsn) ||
1314 ((*ex)[i] >= x_up_[i]-epsn) ) {
1321 Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1322 Teuchos::RCP<const std::vector<Real> > eg; cast_const_vector(eg,g);
1323 Teuchos::RCP<std::vector<Real> > ev; cast_vector(ev,v);
1324 Real epsn = std::min(scale_*eps,min_diff_);
1325 for (
int i = 0; i < dim_; i++ ) {
1326 if ( ((*ex)[i] <= x_lo_[i]+epsn && (*eg)[i] > 0.0) ) {
1333 Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1334 Teuchos::RCP<const std::vector<Real> > eg; cast_const_vector(eg,g);
1335 Teuchos::RCP<std::vector<Real> > ev; cast_vector(ev,v);
1336 Real epsn = std::min(scale_*eps,min_diff_);
1337 for (
int i = 0; i < dim_; i++ ) {
1338 if ( ((*ex)[i] >= x_up_[i]-epsn && (*eg)[i] < 0.0) ) {
1345 Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1346 Teuchos::RCP<const std::vector<Real> > eg; cast_const_vector(eg,g);
1347 Teuchos::RCP<std::vector<Real> > ev; cast_vector(ev,v);
1348 Real epsn = std::min(scale_*eps,min_diff_);
1349 for (
int i = 0; i < dim_; i++ ) {
1350 if ( ((*ex)[i] <= x_lo_[i]+epsn && (*eg)[i] > 0.0) ||
1351 ((*ex)[i] >= x_up_[i]-epsn && (*eg)[i] < 0.0) ) {
1358 Teuchos::RCP<std::vector<Real> > us = Teuchos::rcp(
new std::vector<Real>(dim_,0.0) );
1359 us->assign(x_up_.begin(),x_up_.end());
1365 Teuchos::RCP<std::vector<Real> > ls = Teuchos::rcp(
new std::vector<Real>(dim_,0.0) );
1366 ls->assign(x_lo_.begin(),x_lo_.end());
1372 template<
class Real>
1376 std::vector<Real> x_lo_;
1377 std::vector<Real> x_up_;
1380 Teuchos::RCP<BurgersFEM<Real> > fem_;
1385 xvec = Teuchos::rcp_const_cast<std::vector<Real> >(
1388 catch (std::exception &e) {
1389 xvec = Teuchos::rcp_const_cast<std::vector<Real> >(
1400 catch (std::exception &e) {
1406 void axpy(std::vector<Real> &out,
const Real a,
1407 const std::vector<Real> &x,
const std::vector<Real> &y)
const{
1408 out.resize(dim_,0.0);
1409 for (
unsigned i = 0; i < dim_; i++) {
1410 out[i] = a*x[i] + y[i];
1415 for (
int i = 0; i < dim_; i++ ) {
1416 x[i] = std::max(x_lo_[i],std::min(x_up_[i],x[i]));
1423 : x_lo_(l), x_up_(u), scale_(
scale), fem_(fem) {
1424 dim_ = x_lo_.size();
1425 for (
int i = 0; i < dim_; i++ ) {
1427 min_diff_ = x_up_[i] - x_lo_[i];
1430 min_diff_ = ( (min_diff_ < (x_up_[i] - x_lo_[i])) ? min_diff_ : (x_up_[i] - x_lo_[i]) );
1437 Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1440 for (
int i = 0; i < dim_; i++ ) {
1441 if ( (*ex)[i] >= x_lo_[i] && (*ex)[i] <= x_up_[i] ) { cnt *= 1; }
1444 if ( cnt == 0 ) { val =
false; }
1449 Teuchos::RCP<std::vector<Real> > ex; cast_vector(ex,x);
1454 Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1455 Teuchos::RCP<std::vector<Real> > ev; cast_vector(ev,v);
1456 Real epsn = std::min(scale_*eps,min_diff_);
1457 for (
int i = 0; i < dim_; i++ ) {
1458 if ( ((*ex)[i] <= x_lo_[i]+epsn) ) {
1465 Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1466 Teuchos::RCP<std::vector<Real> > ev; cast_vector(ev,v);
1467 Real epsn = std::min(scale_*eps,min_diff_);
1468 for (
int i = 0; i < dim_; i++ ) {
1469 if ( ((*ex)[i] >= x_up_[i]-epsn) ) {
1476 Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1477 Teuchos::RCP<std::vector<Real> > ev; cast_vector(ev,v);
1478 Real epsn = std::min(scale_*eps,min_diff_);
1479 for (
int i = 0; i < dim_; i++ ) {
1480 if ( ((*ex)[i] <= x_lo_[i]+epsn) ||
1481 ((*ex)[i] >= x_up_[i]-epsn) ) {
1488 Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1489 Teuchos::RCP<const std::vector<Real> > eg; cast_const_vector(eg,g);
1490 Teuchos::RCP<std::vector<Real> > ev; cast_vector(ev,v);
1491 Real epsn = std::min(scale_*eps,min_diff_);
1492 for (
int i = 0; i < dim_; i++ ) {
1493 if ( ((*ex)[i] <= x_lo_[i]+epsn && (*eg)[i] > 0.0) ) {
1500 Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1501 Teuchos::RCP<const std::vector<Real> > eg; cast_const_vector(eg,g);
1502 Teuchos::RCP<std::vector<Real> > ev; cast_vector(ev,v);
1503 Real epsn = std::min(scale_*eps,min_diff_);
1504 for (
int i = 0; i < dim_; i++ ) {
1505 if ( ((*ex)[i] >= x_up_[i]-epsn && (*eg)[i] < 0.0) ) {
1512 Teuchos::RCP<const std::vector<Real> > ex; cast_const_vector(ex,x);
1513 Teuchos::RCP<const std::vector<Real> > eg; cast_const_vector(eg,g);
1514 Teuchos::RCP<std::vector<Real> > ev; cast_vector(ev,v);
1515 Real epsn = std::min(scale_*eps,min_diff_);
1516 for (
int i = 0; i < dim_; i++ ) {
1517 if ( ((*ex)[i] <= x_lo_[i]+epsn && (*eg)[i] > 0.0) ||
1518 ((*ex)[i] >= x_up_[i]-epsn && (*eg)[i] < 0.0) ) {
1525 Teuchos::RCP<std::vector<Real> > us = Teuchos::rcp(
new std::vector<Real>(dim_,0.0) );
1526 us->assign(x_up_.begin(),x_up_.end());
1532 Teuchos::RCP<std::vector<Real> > ls = Teuchos::rcp(
new std::vector<Real>(dim_,0.0) );
1533 ls->assign(x_lo_.begin(),x_lo_.end());
1539 template<
class Real,
class Ordinal>
1545 xvec = Teuchos::rcp_const_cast<std::vector<Real> >(
1548 catch (std::exception &e) {
1549 xvec = Teuchos::rcp_const_cast<std::vector<Real> >(
1556 :
ROL::TeuchosBatchManager<Real,Ordinal>(comm) {}
1558 Teuchos::RCP<std::vector<Real> > input_ptr;
1559 cast_vector(input_ptr,input);
1560 int dim_i = input_ptr->size();
1561 Teuchos::RCP<std::vector<Real> > output_ptr;
1562 cast_vector(output_ptr,output);
1563 int dim_o = output_ptr->size();
1564 if ( dim_i != dim_o ) {
1565 std::cout <<
"L2VectorBatchManager: DIMENSION MISMATCH ON RANK " 1574 template<
class Real,
class Ordinal>
1580 xvec = Teuchos::rcp_const_cast<std::vector<Real> >(
1583 catch (std::exception &e) {
1584 xvec = Teuchos::rcp_const_cast<std::vector<Real> >(
1591 :
ROL::TeuchosBatchManager<Real,Ordinal>(comm) {}
1593 Teuchos::RCP<std::vector<Real> > input_ptr;
1594 cast_vector(input_ptr,input);
1595 int dim_i = input_ptr->size();
1596 Teuchos::RCP<std::vector<Real> > output_ptr;
1597 cast_vector(output_ptr,output);
1598 int dim_o = output_ptr->size();
1599 if ( dim_i != dim_o ) {
1600 std::cout <<
"H1VectorBatchManager: DIMENSION MISMATCH ON RANK " 1609 template<
class Real>
1610 Real
random(
const Teuchos::RCP<
const Teuchos::Comm<int> > &comm) {
1612 if ( Teuchos::rank<int>(*comm)==0 ) {
1613 val = (Real)rand()/(Real)RAND_MAX;
1615 Teuchos::broadcast<int,Real>(*comm,0,1,&val);
BurgersFEM(int nx=128, Real nl=1.0, Real cH1=1.0, Real cL2=1.0)
L2VectorBatchManager(const Teuchos::RCP< const Teuchos::Comm< Ordinal > > &comm)
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
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.
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)
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 .
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.
Real random(const Teuchos::RCP< const Teuchos::Comm< int > > &comm)
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 .
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)
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
void cast_vector(Teuchos::RCP< std::vector< Real > > &xvec, ROL::Vector< Real > &x) const
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
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 .
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)
const std::vector< Real > getParameter(void) const
Provides the interface to apply upper and lower bound constraints.
void sumAll(Real *input, Real *output, int dim)
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 set_problem_data(const Real nu, const Real f, const Real u0, const Real u1)
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)
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
void cast_vector(Teuchos::RCP< std::vector< Real > > &xvec, ROL::Vector< Real > &x) 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.
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 sumAll(ROL::Vector< Real > &input, ROL::Vector< Real > &output)
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)
void sumAll(ROL::Vector< Real > &input, ROL::Vector< Real > &output)
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)
H1VectorBatchManager(const Teuchos::RCP< const Teuchos::Comm< Ordinal > > &comm)
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.