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)),
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_;
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]);
351 std::vector<Real> &d,
352 std::vector<Real> &du,
353 const std::vector<Real> &u)
const {
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;
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++) {
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];
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];
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 ROL::Ptr<std::vector<Real> >
vec_;
530 ROL::Ptr<BurgersFEM<Real> >
fem_;
532 mutable ROL::Ptr<L2VectorDual<Real> >
dual_vec_;
540 const L2VectorPrimal &ex = dynamic_cast<const L2VectorPrimal&>(x);
541 const std::vector<Real>& xval = *ex.
getVector();
542 std::copy(xval.begin(),xval.end(),
vec_->begin());
546 const L2VectorPrimal &ex = dynamic_cast<const L2VectorPrimal&>(x);
547 const std::vector<Real>& xval = *ex.
getVector();
550 (*vec_)[i] += xval[i];
562 const L2VectorPrimal & ex = dynamic_cast<const L2VectorPrimal&>(x);
563 const std::vector<Real>& xval = *ex.
getVector();
564 return fem_->compute_L2_dot(xval,*
vec_);
569 val = std::sqrt(
dot(*
this) );
573 ROL::Ptr<ROL::Vector<Real> >
clone()
const {
574 return ROL::makePtr<L2VectorPrimal>( ROL::makePtr<std::vector<Real>>(
vec_->size(),0.0),
fem_);
585 ROL::Ptr<ROL::Vector<Real> >
basis(
const int i )
const {
586 ROL::Ptr<L2VectorPrimal> e
587 = ROL::makePtr<L2VectorPrimal>( ROL::makePtr<std::vector<Real>>(
vec_->size(),0.0),
fem_);
588 (*e->getVector())[i] = 1.0;
597 dual_vec_ = ROL::makePtr<L2VectorDual<Real>>(
598 ROL::makePtr<std::vector<Real>>(*vec_),
fem_);
600 fem_->apply_mass(*(ROL::constPtrCast<std::vector<Real> >(
dual_vec_->getVector())),*
vec_);
609 ROL::Ptr<std::vector<Real> >
vec_;
610 ROL::Ptr<BurgersFEM<Real> >
fem_;
612 mutable ROL::Ptr<L2VectorPrimal<Real> >
dual_vec_;
620 const L2VectorDual &ex = dynamic_cast<const L2VectorDual&>(x);
621 const std::vector<Real>& xval = *ex.
getVector();
622 std::copy(xval.begin(),xval.end(),
vec_->begin());
626 const L2VectorDual &ex = dynamic_cast<const L2VectorDual&>(x);
627 const std::vector<Real>& xval = *ex.
getVector();
630 (*vec_)[i] += xval[i];
642 const L2VectorDual & ex = dynamic_cast<const L2VectorDual&>(x);
643 const std::vector<Real>& xval = *ex.
getVector();
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 ROL::Ptr<ROL::Vector<Real> >
clone()
const {
661 return ROL::makePtr<L2VectorDual>( ROL::makePtr<std::vector<Real>>(
vec_->size(),0.0),
fem_);
672 ROL::Ptr<ROL::Vector<Real> >
basis(
const int i )
const {
673 ROL::Ptr<L2VectorDual> e
674 = ROL::makePtr<L2VectorDual>( ROL::makePtr<std::vector<Real>>(
vec_->size(),0.0),
fem_);
675 (*e->getVector())[i] = 1.0;
684 dual_vec_ = ROL::makePtr<L2VectorPrimal<Real>>(
685 ROL::makePtr<std::vector<Real>>(*vec_),
fem_);
687 fem_->apply_inverse_mass(*(ROL::constPtrCast<std::vector<Real> >(
dual_vec_->getVector())),*
vec_);
696 ROL::Ptr<std::vector<Real> >
vec_;
697 ROL::Ptr<BurgersFEM<Real> >
fem_;
699 mutable ROL::Ptr<H1VectorDual<Real> >
dual_vec_;
707 const H1VectorPrimal &ex = dynamic_cast<const H1VectorPrimal&>(x);
708 const std::vector<Real>& xval = *ex.
getVector();
709 std::copy(xval.begin(),xval.end(),
vec_->begin());
713 const H1VectorPrimal &ex = dynamic_cast<const H1VectorPrimal&>(x);
714 const std::vector<Real>& xval = *ex.
getVector();
717 (*vec_)[i] += xval[i];
729 const H1VectorPrimal & ex = dynamic_cast<const H1VectorPrimal&>(x);
730 const std::vector<Real>& xval = *ex.
getVector();
731 return fem_->compute_H1_dot(xval,*
vec_);
736 val = std::sqrt(
dot(*
this) );
740 ROL::Ptr<ROL::Vector<Real> >
clone()
const {
741 return ROL::makePtr<H1VectorPrimal>( ROL::makePtr<std::vector<Real>>(
vec_->size(),0.0),
fem_);
752 ROL::Ptr<ROL::Vector<Real> >
basis(
const int i )
const {
753 ROL::Ptr<H1VectorPrimal> e
754 = ROL::makePtr<H1VectorPrimal>( ROL::makePtr<std::vector<Real>>(
vec_->size(),0.0),
fem_);
755 (*e->getVector())[i] = 1.0;
764 dual_vec_ = ROL::makePtr<H1VectorDual<Real>>(
765 ROL::makePtr<std::vector<Real>>(*vec_),
fem_);
767 fem_->apply_H1(*(ROL::constPtrCast<std::vector<Real> >(
dual_vec_->getVector())),*
vec_);
776 ROL::Ptr<std::vector<Real> >
vec_;
777 ROL::Ptr<BurgersFEM<Real> >
fem_;
779 mutable ROL::Ptr<H1VectorPrimal<Real> >
dual_vec_;
787 const H1VectorDual &ex = dynamic_cast<const H1VectorDual&>(x);
788 const std::vector<Real>& xval = *ex.
getVector();
789 std::copy(xval.begin(),xval.end(),
vec_->begin());
793 const H1VectorDual &ex = dynamic_cast<const H1VectorDual&>(x);
794 const std::vector<Real>& xval = *ex.
getVector();
797 (*vec_)[i] += xval[i];
809 const H1VectorDual & ex = dynamic_cast<const H1VectorDual&>(x);
810 const std::vector<Real>& xval = *ex.
getVector();
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 ROL::Ptr<ROL::Vector<Real> >
clone()
const {
828 return ROL::makePtr<H1VectorDual>( ROL::makePtr<std::vector<Real>>(
vec_->size(),0.0),
fem_);
839 ROL::Ptr<ROL::Vector<Real> >
basis(
const int i )
const {
840 ROL::Ptr<H1VectorDual> e
841 = ROL::makePtr<H1VectorDual>( ROL::makePtr<std::vector<Real>>(
vec_->size(),0.0),
fem_);
842 (*e->getVector())[i] = 1.0;
851 dual_vec_ = ROL::makePtr<H1VectorPrimal<Real>>(
852 ROL::makePtr<std::vector<Real>>(*vec_),
fem_);
854 fem_->apply_inverse_H1(*(ROL::constPtrCast<std::vector<Real> >(
dual_vec_->getVector())),*
vec_);
873 ROL::Ptr<BurgersFEM<Real> >
fem_;
882 ROL::Ptr<std::vector<Real> > cp =
883 dynamic_cast<PrimalConstraintVector&>(c).getVector();
884 ROL::Ptr<const std::vector<Real> > up =
885 dynamic_cast<const PrimalStateVector&>(u).getVector();
886 ROL::Ptr<const std::vector<Real> > zp =
887 dynamic_cast<const PrimalControlVector&>(z).getVector();
888 fem_->compute_residual(*cp,*up,*zp);
895 ROL::Ptr<std::vector<Real> > jvp =
896 dynamic_cast<PrimalConstraintVector&>(jv).getVector();
897 ROL::Ptr<const std::vector<Real> > vp =
898 dynamic_cast<const PrimalStateVector&>(v).getVector();
899 ROL::Ptr<const std::vector<Real> > up =
900 dynamic_cast<const PrimalStateVector&>(u).getVector();
901 ROL::Ptr<const std::vector<Real> > zp =
902 dynamic_cast<const PrimalControlVector&>(z).getVector();
903 fem_->apply_pde_jacobian(*jvp,*vp,*up,*zp);
908 ROL::Ptr<std::vector<Real> > jvp =
909 dynamic_cast<PrimalConstraintVector&>(jv).getVector();
910 ROL::Ptr<const std::vector<Real> > vp =
911 dynamic_cast<const PrimalControlVector&>(v).getVector();
912 ROL::Ptr<const std::vector<Real> > up =
913 dynamic_cast<const PrimalStateVector&>(u).getVector();
914 ROL::Ptr<const std::vector<Real> > zp =
915 dynamic_cast<const PrimalControlVector&>(z).getVector();
916 fem_->apply_control_jacobian(*jvp,*vp,*up,*zp);
921 ROL::Ptr<std::vector<Real> > ijvp =
922 dynamic_cast<PrimalStateVector&>(ijv).getVector();
923 ROL::Ptr<const std::vector<Real> > vp =
924 dynamic_cast<const PrimalConstraintVector&>(v).getVector();
925 ROL::Ptr<const std::vector<Real> > up =
926 dynamic_cast<const PrimalStateVector&>(u).getVector();
927 ROL::Ptr<const std::vector<Real> > zp =
928 dynamic_cast<const PrimalControlVector&>(z).getVector();
929 fem_->apply_inverse_pde_jacobian(*ijvp,*vp,*up,*zp);
934 ROL::Ptr<std::vector<Real> > jvp =
935 dynamic_cast<DualStateVector&>(ajv).getVector();
936 ROL::Ptr<const std::vector<Real> > vp =
937 dynamic_cast<const DualConstraintVector&>(v).getVector();
938 ROL::Ptr<const std::vector<Real> > up =
939 dynamic_cast<const PrimalStateVector&>(u).getVector();
940 ROL::Ptr<const std::vector<Real> > zp =
941 dynamic_cast<const PrimalControlVector&>(z).getVector();
942 fem_->apply_adjoint_pde_jacobian(*jvp,*vp,*up,*zp);
947 ROL::Ptr<std::vector<Real> > jvp =
948 dynamic_cast<DualControlVector&>(jv).getVector();
949 ROL::Ptr<const std::vector<Real> > vp =
950 dynamic_cast<const DualConstraintVector&>(v).getVector();
951 ROL::Ptr<const std::vector<Real> > up =
952 dynamic_cast<const PrimalStateVector&>(u).getVector();
953 ROL::Ptr<const std::vector<Real> > zp =
954 dynamic_cast<const PrimalControlVector&>(z).getVector();
955 fem_->apply_adjoint_control_jacobian(*jvp,*vp,*up,*zp);
960 ROL::Ptr<std::vector<Real> > iajvp =
961 dynamic_cast<DualConstraintVector&>(iajv).getVector();
962 ROL::Ptr<const std::vector<Real> > vp =
963 dynamic_cast<const DualStateVector&>(v).getVector();
964 ROL::Ptr<const std::vector<Real> > up =
965 dynamic_cast<const PrimalStateVector&>(u).getVector();
966 ROL::Ptr<const std::vector<Real> > zp =
967 dynamic_cast<const PrimalControlVector&>(z).getVector();
968 fem_->apply_inverse_adjoint_pde_jacobian(*iajvp,*vp,*up,*zp);
974 ROL::Ptr<std::vector<Real> > ahwvp =
975 dynamic_cast<DualStateVector&>(ahwv).getVector();
976 ROL::Ptr<const std::vector<Real> > wp =
977 dynamic_cast<const DualConstraintVector&>(w).getVector();
978 ROL::Ptr<const std::vector<Real> > vp =
979 dynamic_cast<const PrimalStateVector&>(v).getVector();
980 ROL::Ptr<const std::vector<Real> > up =
981 dynamic_cast<const PrimalStateVector&>(u).getVector();
982 ROL::Ptr<const std::vector<Real> > zp =
983 dynamic_cast<const PrimalControlVector&>(z).getVector();
984 fem_->apply_adjoint_pde_hessian(*ahwvp,*wp,*vp,*up,*zp);
994 ROL::Ptr<std::vector<Real> > ahwvp =
995 dynamic_cast<DualControlVector&>(ahwv).getVector();
996 ROL::Ptr<const std::vector<Real> > wp =
997 dynamic_cast<const DualConstraintVector&>(w).getVector();
998 ROL::Ptr<const std::vector<Real> > vp =
999 dynamic_cast<const PrimalStateVector&>(v).getVector();
1000 ROL::Ptr<const std::vector<Real> > up =
1001 dynamic_cast<const PrimalStateVector&>(u).getVector();
1002 ROL::Ptr<const std::vector<Real> > zp =
1003 dynamic_cast<const PrimalControlVector&>(z).getVector();
1004 fem_->apply_adjoint_control_pde_hessian(*ahwvp,*wp,*vp,*up,*zp);
1013 ROL::Ptr<std::vector<Real> > ahwvp =
1014 dynamic_cast<DualStateVector&>(ahwv).getVector();
1015 ROL::Ptr<const std::vector<Real> > wp =
1016 dynamic_cast<const DualConstraintVector&>(w).getVector();
1017 ROL::Ptr<const std::vector<Real> > vp =
1018 dynamic_cast<const PrimalControlVector&>(v).getVector();
1019 ROL::Ptr<const std::vector<Real> > up =
1020 dynamic_cast<const PrimalStateVector&>(u).getVector();
1021 ROL::Ptr<const std::vector<Real> > zp =
1022 dynamic_cast<const PrimalControlVector&>(z).getVector();
1023 fem_->apply_adjoint_pde_control_hessian(*ahwvp,*wp,*vp,*up,*zp);
1032 ROL::Ptr<std::vector<Real> > ahwvp =
1033 dynamic_cast<DualControlVector&>(ahwv).getVector();
1034 ROL::Ptr<const std::vector<Real> > wp =
1035 dynamic_cast<const DualConstraintVector&>(w).getVector();
1036 ROL::Ptr<const std::vector<Real> > vp =
1037 dynamic_cast<const PrimalControlVector&>(v).getVector();
1038 ROL::Ptr<const std::vector<Real> > up =
1039 dynamic_cast<const PrimalStateVector&>(u).getVector();
1040 ROL::Ptr<const std::vector<Real> > zp =
1041 dynamic_cast<const PrimalControlVector&>(z).getVector();
1042 fem_->apply_adjoint_control_hessian(*ahwvp,*wp,*vp,*up,*zp);
1050 template<
class Real>
1062 ROL::Ptr<ROL::Vector<Real> >
ud_;
1073 ROL::Ptr<const std::vector<Real> > up =
1074 dynamic_cast<const PrimalStateVector&>(u).getVector();
1075 ROL::Ptr<const std::vector<Real> > zp =
1076 dynamic_cast<const PrimalControlVector&>(z).getVector();
1077 ROL::Ptr<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 ROL::Ptr<std::vector<Real> > gp =
1089 dynamic_cast<DualStateVector&>(g).getVector();
1090 ROL::Ptr<const std::vector<Real> > up =
1091 dynamic_cast<const PrimalStateVector&>(u).getVector();
1092 ROL::Ptr<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 ROL::Ptr<std::vector<Real> > gp =
1104 dynamic_cast<DualControlVector&>(g).getVector();
1105 ROL::Ptr<const std::vector<Real> > zp =
1106 dynamic_cast<const PrimalControlVector&>(z).getVector();
1108 fem_->apply_mass(*gp,*zp);
1109 for (
unsigned i = 0; i < zp->size(); i++) {
1116 ROL::Ptr<std::vector<Real> > hvp =
1117 dynamic_cast<DualStateVector&>(hv).getVector();
1118 ROL::Ptr<const std::vector<Real> > vp =
1119 dynamic_cast<const PrimalStateVector&>(v).getVector();
1121 fem_->apply_mass(*hvp,*vp);
1136 ROL::Ptr<std::vector<Real> > hvp =
1137 dynamic_cast<DualControlVector&>(hv).getVector();
1138 ROL::Ptr<const std::vector<Real> > vp =
1139 dynamic_cast<const PrimalControlVector&>(v).getVector();
1141 fem_->apply_mass(*hvp,*vp);
1142 for (
unsigned i = 0; i < vp->size(); i++) {
1148 template<
class Real>
1157 ROL::Ptr<ROL::Vector<Real> >
l_;
1158 ROL::Ptr<ROL::Vector<Real> >
u_;
1165 catch (std::exception &e) {
1175 catch (std::exception &e) {
1180 void axpy(std::vector<Real> &out,
const Real a,
1181 const std::vector<Real> &x,
const std::vector<Real> &y)
const{
1182 out.resize(
dim_,0.0);
1183 for (
unsigned i = 0; i <
dim_; i++) {
1184 out[i] = a*x[i] + y[i];
1189 for (
int i = 0; i <
dim_; i++ ) {
1190 x[i] = std::max(
x_lo_[i],std::min(
x_up_[i],x[i]));
1199 for (
int i = 0; i <
dim_; i++ ) {
1208 l_ = ROL::makePtr<L2VectorPrimal<Real>>(
1209 ROL::makePtr<std::vector<Real>>(l), fem);
1210 u_ = ROL::makePtr<L2VectorPrimal<Real>>(
1211 ROL::makePtr<std::vector<Real>>(u), fem);
1218 for (
int i = 0; i <
dim_; i++ ) {
1219 if ( (*ex)[i] >=
x_lo_[i] && (*ex)[i] <=
x_up_[i] ) { cnt *= 1; }
1222 if ( cnt == 0 ) { val =
false; }
1227 ROL::Ptr<std::vector<Real> > ex;
cast_vector(ex,x);
1233 ROL::Ptr<std::vector<Real> > ev;
cast_vector(ev,v);
1235 for (
int i = 0; i <
dim_; i++ ) {
1236 if ( ((*ex)[i] <=
x_lo_[i]+epsn) ) {
1244 ROL::Ptr<std::vector<Real> > ev;
cast_vector(ev,v);
1246 for (
int i = 0; i <
dim_; i++ ) {
1247 if ( ((*ex)[i] >=
x_up_[i]-epsn) ) {
1255 ROL::Ptr<std::vector<Real> > ev;
cast_vector(ev,v);
1257 for (
int i = 0; i <
dim_; i++ ) {
1258 if ( ((*ex)[i] <=
x_lo_[i]+epsn) ||
1259 ((*ex)[i] >=
x_up_[i]-epsn) ) {
1268 ROL::Ptr<std::vector<Real> > ev;
cast_vector(ev,v);
1270 for (
int i = 0; i <
dim_; i++ ) {
1271 if ( ((*ex)[i] <=
x_lo_[i]+epsn && (*eg)[i] > 0.0) ) {
1280 ROL::Ptr<std::vector<Real> > ev;
cast_vector(ev,v);
1282 for (
int i = 0; i <
dim_; i++ ) {
1283 if ( ((*ex)[i] >=
x_up_[i]-epsn && (*eg)[i] < 0.0) ) {
1292 ROL::Ptr<std::vector<Real> > ev;
cast_vector(ev,v);
1294 for (
int i = 0; i <
dim_; i++ ) {
1295 if ( ((*ex)[i] <=
x_lo_[i]+epsn && (*eg)[i] > 0.0) ||
1296 ((*ex)[i] >=
x_up_[i]-epsn && (*eg)[i] < 0.0) ) {
1311 template<
class Real>
1320 ROL::Ptr<ROL::Vector<Real> >
l_;
1321 ROL::Ptr<ROL::Vector<Real> >
u_;
1328 catch (std::exception &e) {
1338 catch (std::exception &e) {
1343 void axpy(std::vector<Real> &out,
const Real a,
1344 const std::vector<Real> &x,
const std::vector<Real> &y)
const{
1345 out.resize(
dim_,0.0);
1346 for (
unsigned i = 0; i <
dim_; i++) {
1347 out[i] = a*x[i] + y[i];
1352 for (
int i = 0; i <
dim_; i++ ) {
1353 x[i] = std::max(
x_lo_[i],std::min(
x_up_[i],x[i]));
1362 for (
int i = 0; i <
dim_; i++ ) {
1371 l_ = ROL::makePtr<H1VectorPrimal<Real>>(
1372 ROL::makePtr<std::vector<Real>>(l), fem);
1373 u_ = ROL::makePtr<H1VectorPrimal<Real>>(
1374 ROL::makePtr<std::vector<Real>>(u), fem);
1381 for (
int i = 0; i <
dim_; i++ ) {
1382 if ( (*ex)[i] >=
x_lo_[i] && (*ex)[i] <=
x_up_[i] ) { cnt *= 1; }
1385 if ( cnt == 0 ) { val =
false; }
1390 ROL::Ptr<std::vector<Real> > ex;
cast_vector(ex,x);
1396 ROL::Ptr<std::vector<Real> > ev;
cast_vector(ev,v);
1398 for (
int i = 0; i <
dim_; i++ ) {
1399 if ( ((*ex)[i] <=
x_lo_[i]+epsn) ) {
1407 ROL::Ptr<std::vector<Real> > ev;
cast_vector(ev,v);
1409 for (
int i = 0; i <
dim_; i++ ) {
1410 if ( ((*ex)[i] >=
x_up_[i]-epsn) ) {
1418 ROL::Ptr<std::vector<Real> > ev;
cast_vector(ev,v);
1420 for (
int i = 0; i <
dim_; i++ ) {
1421 if ( ((*ex)[i] <=
x_lo_[i]+epsn) ||
1422 ((*ex)[i] >=
x_up_[i]-epsn) ) {
1431 ROL::Ptr<std::vector<Real> > ev;
cast_vector(ev,v);
1433 for (
int i = 0; i <
dim_; i++ ) {
1434 if ( ((*ex)[i] <=
x_lo_[i]+epsn && (*eg)[i] > 0.0) ) {
1443 ROL::Ptr<std::vector<Real> > ev;
cast_vector(ev,v);
1445 for (
int i = 0; i <
dim_; i++ ) {
1446 if ( ((*ex)[i] >=
x_up_[i]-epsn && (*eg)[i] < 0.0) ) {
1455 ROL::Ptr<std::vector<Real> > ev;
cast_vector(ev,v);
1457 for (
int i = 0; i <
dim_; i++ ) {
1458 if ( ((*ex)[i] <=
x_lo_[i]+epsn && (*eg)[i] > 0.0) ||
1459 ((*ex)[i] >=
x_up_[i]-epsn && (*eg)[i] < 0.0) ) {