84 #include "Teuchos_oblackholestream.hpp" 85 #include "Teuchos_GlobalMPISession.hpp" 86 #include "Teuchos_XMLParameterListHelpers.hpp" 100 template <
class Real,
class Element=Real>
103 template <
class Real,
class Element=Real>
106 template <
class Real,
class Element=Real>
109 template <
class Real,
class Element=Real>
116 template <
class Real,
class Element>
120 typedef typename vector::size_type
uint;
123 Teuchos::RCP<std::vector<Element> > std_vec_;
124 mutable Teuchos::RCP<OptDualStdVector<Real> > dual_vec_;
126 Teuchos::RCP<FiniteDifference<Real> >
fd_;
132 std_vec_(std_vec), dual_vec_(Teuchos::null), fd_(fd) {}
136 Teuchos::RCP<const vector> xvalptr = ex.
getVector();
137 uint dimension = std_vec_->size();
138 for (
uint i=0; i<dimension; i++) {
139 (*std_vec_)[i] += (*xvalptr)[i];
144 uint dimension = std_vec_->size();
145 for (
uint i=0; i<dimension; i++) {
146 (*std_vec_)[i] *= alpha;
155 Teuchos::RCP<const vector> xvalptr = ex.
getVector();
157 Teuchos::RCP<vector> kxvalptr = Teuchos::rcp(
new vector(std_vec_->size(), 0.0) );
159 fd_->apply(xvalptr,kxvalptr);
161 uint dimension = std_vec_->size();
162 for (
uint i=0; i<dimension; i++) {
163 val += (*std_vec_)[i]*(*kxvalptr)[i];
170 val = std::sqrt( dot(*
this) );
174 Teuchos::RCP<Vector<Real> >
clone()
const {
175 return Teuchos::rcp(
new OptStdVector( Teuchos::rcp(
new vector(std_vec_->size()) ),fd_ ) );
186 Teuchos::RCP<Vector<Real> >
basis(
const int i )
const {
187 Teuchos::RCP<vector> e_rcp = Teuchos::rcp(
new vector(std_vec_->size(),0.0) );
188 Teuchos::RCP<OptStdVector> e = Teuchos::rcp(
new OptStdVector( e_rcp, fd_ ) );
193 int dimension()
const {
return static_cast<int>(std_vec_->size());}
198 Teuchos::RCP<vector> dual_vecp = Teuchos::rcp(
new vector(*std_vec_));
200 fd_->apply(dual_vecp);
208 template <
class Real,
class Element>
212 typedef typename vector::size_type
uint;
215 Teuchos::RCP<std::vector<Element> > std_vec_;
216 mutable Teuchos::RCP<OptStdVector<Real> > dual_vec_;
217 Teuchos::RCP<FiniteDifference<Real> >
fd_;
222 std_vec_(std_vec), dual_vec_(Teuchos::null), fd_(fd) {}
226 Teuchos::RCP<const vector> xvalptr = ex.
getVector();
227 uint dimension = std_vec_->size();
228 for (
uint i=0; i<dimension; i++) {
229 (*std_vec_)[i] += (*xvalptr)[i];
234 uint dimension = std_vec_->size();
235 for (
uint i=0; i<dimension; i++) {
236 (*std_vec_)[i] *= alpha;
243 Teuchos::RCP<const vector> kxvalptr = ex.
getVector();
244 Teuchos::RCP<vector> xvalptr = Teuchos::rcp(
new vector(std_vec_->size(), 0.0) );
245 fd_->solve(kxvalptr,xvalptr);
247 uint dimension = std_vec_->size();
248 for (
unsigned i=0; i<dimension; i++) {
249 val += (*std_vec_)[i]*(*xvalptr)[i];
256 val = std::sqrt( dot(*
this) );
260 Teuchos::RCP<Vector<Real> >
clone()
const {
261 return Teuchos::rcp(
new OptDualStdVector( Teuchos::rcp(
new std::vector<Element>(std_vec_->size()) ), fd_ ) );
264 Teuchos::RCP<const std::vector<Element> >
getVector()
const {
272 Teuchos::RCP<Vector<Real> >
basis(
const int i )
const {
273 Teuchos::RCP<vector> e_rcp = Teuchos::rcp(
new vector(std_vec_->size(), 0.0 ) );
274 Teuchos::RCP<OptDualStdVector> e = Teuchos::rcp(
new OptDualStdVector( e_rcp,fd_ ) );
279 int dimension()
const {
return static_cast<int>(std_vec_->size());}
282 Teuchos::RCP<vector> dual_vecp = Teuchos::rcp(
new vector(*std_vec_));
285 fd_->solve(dual_vecp);
295 template <
class Real,
class Element>
299 typedef typename vector::size_type
uint;
302 Teuchos::RCP<std::vector<Element> > std_vec_;
303 mutable Teuchos::RCP<ConDualStdVector<Real> > dual_vec_;
307 ConStdVector(
const Teuchos::RCP<std::vector<Element> > & std_vec) : std_vec_(std_vec), dual_vec_(Teuchos::null) {}
311 Teuchos::RCP<const vector> xvalptr = ex.
getVector();
312 uint dimension = std_vec_->size();
313 for (
uint i=0; i<dimension; i++) {
314 (*std_vec_)[i] += (*xvalptr)[i];
319 uint dimension = std_vec_->size();
320 for (
uint i=0; i<dimension; i++) {
321 (*std_vec_)[i] *= alpha;
328 Teuchos::RCP<const vector> xvalptr = ex.
getVector();
330 uint dimension = std_vec_->size();
331 for (
uint i=0; i<dimension; i++) {
332 val += (*std_vec_)[i]*(*xvalptr)[i];
339 val = std::sqrt( dot(*
this) );
343 Teuchos::RCP<Vector<Real> >
clone()
const {
344 return Teuchos::rcp(
new ConStdVector( Teuchos::rcp(
new vector(std_vec_->size())) ) );
347 Teuchos::RCP<const std::vector<Element> >
getVector()
const {
355 Teuchos::RCP<Vector<Real> >
basis(
const int i )
const {
356 Teuchos::RCP<vector> e_rcp = Teuchos::rcp(
new vector(std_vec_->size(),0.0) );
357 Teuchos::RCP<ConStdVector> e = Teuchos::rcp(
new ConStdVector( e_rcp) );
362 int dimension()
const {
return static_cast<int>(std_vec_->size());}
365 dual_vec_ = Teuchos::rcp(
new ConDualStdVector<Real>( Teuchos::rcp(
new std::vector<Element>(*std_vec_) ) ) );
373 template <
class Real,
class Element>
377 typedef typename vector::size_type
uint;
381 Teuchos::RCP<std::vector<Element> > std_vec_;
382 mutable Teuchos::RCP<ConStdVector<Real> > dual_vec_;
386 ConDualStdVector(
const Teuchos::RCP<std::vector<Element> > & std_vec) : std_vec_(std_vec), dual_vec_(Teuchos::null) {}
390 Teuchos::RCP<const vector> xvalptr = ex.
getVector();
391 uint dimension = std_vec_->size();
392 for (
uint i=0; i<dimension; i++) {
393 (*std_vec_)[i] += (*xvalptr)[i];
398 uint dimension = std_vec_->size();
399 for (
uint i=0; i<dimension; i++) {
400 (*std_vec_)[i] *= alpha;
407 Teuchos::RCP<const vector> xvalptr = ex.
getVector();
408 uint dimension = std_vec_->size();
409 for (
uint i=0; i<dimension; i++) {
410 val += (*std_vec_)[i]*(*xvalptr)[i];
417 val = std::sqrt( dot(*
this) );
421 Teuchos::RCP<Vector<Real> >
clone()
const {
422 return Teuchos::rcp(
new ConDualStdVector( Teuchos::rcp(
new std::vector<Element>(std_vec_->size())) ) );
425 Teuchos::RCP<const std::vector<Element> >
getVector()
const {
433 Teuchos::RCP<Vector<Real> >
basis(
const int i )
const {
434 Teuchos::RCP<vector> e_rcp = Teuchos::rcp(
new vector(std_vec_->size(),0.0) );
435 Teuchos::RCP<ConDualStdVector> e = Teuchos::rcp(
new ConDualStdVector( e_rcp ) );
440 int dimension()
const {
return static_cast<int>(std_vec_->size());}
443 dual_vec_ = Teuchos::rcp(
new ConStdVector<Real>( Teuchos::rcp(
new std::vector<Element>(*std_vec_) ) ) );
454 template<
class Real,
class XPrim=StdVector<Real>,
class XDual=StdVector<Real> >
458 typedef typename vector::size_type
uint;
472 Teuchos::RCP<const std::vector<Real> >
Vp_;
474 Teuchos::RCP<FiniteDifference<Real> >
fd_;
481 using Teuchos::RCP;
using Teuchos::dyn_cast;
484 RCP<const vector> vp = dyn_cast<
const XPrim>(v).getVector();
487 RCP<vector> kvp = dyn_cast<XDual>(kv).getVector();
491 (*kvp)[0] = (2.0*(*vp)[0]-(*vp)[1])/dx2;
493 for(
uint i=1;i<nx_-1;++i) {
494 (*kvp)[i] = (2.0*(*vp)[i]-(*vp)[i-1]-(*vp)[i+1])/dx2;
497 (*kvp)[nx_-1] = (2.0*(*vp)[nx_-1]-(*vp)[nx_-2])/dx2;
504 Vp_((Teuchos::dyn_cast<const
StdVector<Real> >(
V)).getVector()), fd_(fd) {
507 dx_ = (1.0/(1.0+nx_));
517 using Teuchos::RCP;
using Teuchos::rcp;
using Teuchos::dyn_cast;
520 RCP<const vector> psip = dyn_cast<
const XPrim>(psi).getVector();
523 RCP<vector> kpsip = rcp(
new vector(nx_, 0.0) );
524 XDual kpsi(kpsip,fd_);
530 for(
uint i=0;i<nx_;++i) {
531 J += (*psip)[i]*(*kpsip)[i] + (*Vp_)[i]*pow((*psip)[i],2) + g_*pow((*psip)[i],4);
544 using Teuchos::RCP;
using Teuchos::rcp;
using Teuchos::dyn_cast;
547 RCP<const vector> psip = dyn_cast<
const XPrim>(psi).getVector();
550 RCP<vector> gp = dyn_cast<XDual>(g).getVector();
553 RCP<vector> kpsip = rcp(
new vector(nx_, 0.0) );
554 XDual kpsi(kpsip,fd_);
558 for(
uint i=0;i<nx_;++i) {
559 (*gp)[i] = ((*kpsip)[i] + (*Vp_)[i]*(*psip)[i] + 2.0*g_*pow((*psip)[i],3))*dx_;
570 using Teuchos::RCP;
using Teuchos::dyn_cast;
573 RCP<const vector> psip = dyn_cast<
const XPrim>(psi).getVector();
576 RCP<const vector> vp = dyn_cast<
const XPrim>(v).getVector();
579 RCP<vector> hvp = dyn_cast<XDual>(hv).getVector();
583 for(
uint i=0;i<nx_;++i) {
585 (*hvp)[i] += ( (*Vp_)[i] + 6.0*g_*pow((*psip)[i],2) )*(*vp)[i]*dx_;
594 template<
class Real,
class XPrim=StdVector<Real>,
class XDual=StdVector<Real>,
class CPrim=StdVector<Real>,
class CDual=StdVector<Real> >
598 typedef typename vector::size_type
uint;
603 Teuchos::RCP<FiniteDifference<Real> >
fd_;
608 nx_(n), dx_(dx), fd_(fd), exactsolve_(exactsolve) {}
617 using Teuchos::RCP;
using Teuchos::dyn_cast;
620 RCP<vector> cp = dyn_cast<CPrim>(c).getVector();
623 RCP<const vector> psip = dyn_cast<
const XPrim>(psi).getVector();
626 for(
uint i=0;i<nx_;++i) {
627 (*cp)[0] += dx_*pow((*psip)[i],2);
636 using Teuchos::RCP;
using Teuchos::dyn_cast;
639 RCP<vector> jvp = dyn_cast<CPrim>(jv).getVector();
642 RCP<const vector> vp = dyn_cast<
const XPrim>(v).getVector();
645 RCP<const vector> psip = dyn_cast<
const XPrim>(psi).getVector();
648 for(
uint i=0;i<nx_;++i) {
649 (*jvp)[0] += 2.0*dx_*(*psip)[i]*(*vp)[i];
658 using Teuchos::RCP;
using Teuchos::dyn_cast;
661 RCP<vector> ajvp = dyn_cast<XDual>(ajv).getVector();
664 RCP<const vector> vp = dyn_cast<
const CDual>(v).getVector();
667 RCP<const vector> psip = dyn_cast<
const XPrim>(psi).getVector();
669 for(
uint i=0;i<nx_;++i) {
670 (*ajvp)[i] = 2.0*dx_*(*psip)[i]*(*vp)[0];
682 using Teuchos::RCP;
using Teuchos::dyn_cast;
685 RCP<vector> ahuvp = dyn_cast<XDual>(ahuv).getVector();
688 RCP<const vector> up = dyn_cast<
const CDual>(u).getVector();
691 RCP<const vector> vp = dyn_cast<
const XPrim>(v).getVector();
694 RCP<const vector> psip = dyn_cast<
const XPrim>(psi).getVector();
696 for(
uint i=0;i<nx_;++i) {
697 (*ahuvp)[i] = 2.0*dx_*(*vp)[i]*(*up)[0];
709 using Teuchos::RCP;
using Teuchos::rcp;
using Teuchos::dyn_cast;
712 RCP<vector> v1p = dyn_cast<XPrim>(v1).getVector();
713 RCP<vector> v2p = dyn_cast<CDual>(v2).getVector();
714 RCP<const vector> b1p = dyn_cast<
const XDual>(b1).getVector();
715 RCP<const vector> b2p = dyn_cast<
const CPrim>(b2).getVector();
716 RCP<const vector> psip = dyn_cast<
const XPrim>(psi).getVector();
718 RCP<vector> jacp = rcp(
new vector(nx_, 0.0) );
719 RCP<vector> b1dp = rcp(
new vector(nx_, 0.0) );
721 for(
uint i=0;i<nx_;++i) {
722 (*jacp)[i] = (*psip)[i];
723 (*b1dp)[i] = (*b1p)[i];
734 Real d = 1.0/jac.dot(jac);
735 Real p = jac.dot(b1d);
737 (*v2p)[0] = d*(p-(*b2p)[0]);
740 v1.
scale(-(*v2p)[0]);
743 return std::vector<Real>(0);
std::vector< Element > vector
Provides the interface to evaluate objective functions.
Teuchos::RCP< FiniteDifference< Real > > fd_
Real value(const Vector< Real > &psi, Real &tol)
Evaluate .
virtual void scale(const Real alpha)=0
Compute where .
void scale(const Real alpha)
Compute where .
virtual void plus(const Vector &x)=0
Compute , where .
void applyAdjointJacobian(Vector< Real > &ajv, const Vector< Real > &v, const Vector< Real > &psi, Real &tol)
Evaluate .
std::vector< Element > vector
ConStdVector(const Teuchos::RCP< std::vector< Element > > &std_vec)
Objective_GrossPitaevskii(const Real &g, const Vector< Real > &V, Teuchos::RCP< FiniteDifference< Real > > fd)
std::vector< Real > vector
void value(Vector< Real > &c, const Vector< Real > &psi, Real &tol)
Evaluate .
Real norm() const
Returns where .
void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &psi, Real &tol)
Evaluate .
Teuchos::RCP< std::vector< Element > > getVector()
Teuchos::RCP< const std::vector< Element > > getVector() const
Teuchos::RCP< const vector > getVector() const
OptStdVector(const Teuchos::RCP< std::vector< Element > > &std_vec, Teuchos::RCP< FiniteDifference< Real > >fd)
Teuchos::RCP< Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
Teuchos::RCP< Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
Teuchos::RCP< FiniteDifference< Real > > fd_
std::vector< Element > vector
Defines the linear algebra or vector space interface.
std::vector< Real > vector
Teuchos::RCP< Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
void applyJacobian(Vector< Real > &jv, const Vector< Real > &v, const Vector< Real > &psi, Real &tol)
Evaluate .
Teuchos::RCP< Vector< Real > > basis(const int i) const
Return i-th basis vector.
Teuchos::RCP< Vector< Real > > basis(const int i) const
Return i-th basis vector.
Real dot(const Vector< Real > &x) const
Compute where .
std::vector< Element > vector
const Vector< Real > & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
Real dot(const Vector< Real > &x) const
Modify the dot product between primal variables to be .
const Vector< Real > & dual() const
Modify the dual of vector u to be .
std::vector< Real > solveAugmentedSystem(Vector< Real > &v1, Vector< Real > &v2, const Vector< Real > &b1, const Vector< Real > &b2, const Vector< Real > &psi, Real &tol)
int dimension() const
Return dimension of the vector space.
virtual std::vector< Real > solveAugmentedSystem(Vector< Real > &v1, Vector< Real > &v2, const Vector< Real > &b1, const Vector< Real > &b2, const Vector< Real > &x, Real &tol)
Approximately solves the augmented system where , , , , is an identity or Riesz operator...
Defines the equality constraint operator interface.
Teuchos::RCP< Vector< Real > > basis(const int i) const
Return i-th basis vector.
Provides the std::vector implementation of the ROL::Vector interface.
void scale(const Real alpha)
Compute where .
Teuchos::RCP< Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
void scale(const Real alpha)
Compute where .
Teuchos::RCP< vector > getVector()
Real norm() const
Returns where .
Teuchos::RCP< const std::vector< Element > > getVector() const
const Vector< Real > & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
OptDualStdVector(const Teuchos::RCP< std::vector< Element > > &std_vec, Teuchos::RCP< FiniteDifference< Real > >fd)
void applyK(const Vector< Real > &v, Vector< Real > &kv)
Apply finite difference operator.
Teuchos::RCP< const std::vector< Real > > Vp_
void plus(const Vector< Real > &x)
Compute , where .
Real dot(const Vector< Real > &x) const
Compute where .
int dimension() const
Return dimension of the vector space.
Real dot(const Vector< Real > &x) const
Compute where .
Teuchos::RCP< const std::vector< Element > > getVector() const
Real norm() const
Returns where .
void plus(const Vector< Real > &x)
Compute , where .
Teuchos::RCP< FiniteDifference< Real > > fd_
int dimension() const
Return dimension of the vector space.
ConDualStdVector(const Teuchos::RCP< std::vector< Element > > &std_vec)
Teuchos::RCP< Vector< Real > > basis(const int i) const
Return i-th basis vector.
void gradient(Vector< Real > &g, const Vector< Real > &psi, Real &tol)
Evaluate .
Teuchos::RCP< std::vector< Element > > getVector()
virtual void set(const Vector &x)
Set where .
Normalization_Constraint(int n, Real dx, Teuchos::RCP< FiniteDifference< Real > > fd, bool exactsolve)
void plus(const Vector< Real > &x)
Compute , where .
Real norm() const
Returns where .
void plus(const Vector< Real > &x)
Compute , where .
Teuchos::RCP< std::vector< Element > > getVector()
Teuchos::RCP< const std::vector< Element > > getVector() const
void applyAdjointHessian(Vector< Real > &ahuv, const Vector< Real > &u, const Vector< Real > &v, const Vector< Real > &psi, Real &tol)
Evaluate .
void scale(const Real alpha)
Compute where .
Teuchos::RCP< FiniteDifference< Real > > fd_
int dimension() const
Return dimension of the vector space.
const Vector< Real > & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...