ROL
gross-pitaevskii/example_02.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ************************************************************************
3 //
4 // Rapid Optimization Library (ROL) Package
5 // Copyright (2014) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact lead developers:
38 // Drew Kouri (dpkouri@sandia.gov) and
39 // Denis Ridzal (dridzal@sandia.gov)
40 //
41 // ************************************************************************
42 // @HEADER
43 
82 #include <iostream>
83 
84 #include "Teuchos_oblackholestream.hpp"
85 #include "Teuchos_GlobalMPISession.hpp"
86 #include "Teuchos_XMLParameterListHelpers.hpp"
87 
88 #include "ROL_StdVector.hpp"
89 #include "ROL_Objective.hpp"
91 #include "ROL_Algorithm.hpp"
92 #include "ROL_CompositeStep.hpp"
94 
96 
97 
98 using namespace ROL;
99 
100 template <class Real, class Element=Real>
101 class OptStdVector; // Optimization space.
102 
103 template <class Real, class Element=Real>
104 class OptDualStdVector; // Dual optimization space.
105 
106 template <class Real, class Element=Real>
107 class ConStdVector; // Constraint space.
108 
109 template <class Real, class Element=Real>
110 class ConDualStdVector; // Dual constraint space.
111 
112 // Vector space definitions:
113 
114 
115 // Optimization space.
116 template <class Real, class Element>
117 class OptStdVector : public Vector<Real> {
118 
119  typedef std::vector<Element> vector;
120  typedef typename vector::size_type uint;
121 
122 private:
123 Teuchos::RCP<std::vector<Element> > std_vec_;
124 mutable Teuchos::RCP<OptDualStdVector<Real> > dual_vec_;
125 
126 Teuchos::RCP<FiniteDifference<Real> > fd_;
127 
128 
129 public:
130 
131 OptStdVector(const Teuchos::RCP<std::vector<Element> > & std_vec, Teuchos::RCP<FiniteDifference<Real> >fd) :
132  std_vec_(std_vec), dual_vec_(Teuchos::null), fd_(fd) {}
133 
134 void plus( const Vector<Real> &x ) {
135  const OptStdVector &ex = Teuchos::dyn_cast<const OptStdVector>(x);
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];
140  }
141 }
142 
143 void scale( const Real alpha ) {
144  uint dimension = std_vec_->size();
145  for (uint i=0; i<dimension; i++) {
146  (*std_vec_)[i] *= alpha;
147  }
148 }
149 
150 
152 Real dot( const Vector<Real> &x ) const {
153  Real val = 0;
154  const OptStdVector<Real, Element> & ex = Teuchos::dyn_cast<const OptStdVector>(x);
155  Teuchos::RCP<const vector> xvalptr = ex.getVector();
156 
157  Teuchos::RCP<vector> kxvalptr = Teuchos::rcp( new vector(std_vec_->size(), 0.0) );
158 
159  fd_->apply(xvalptr,kxvalptr);
160 
161  uint dimension = std_vec_->size();
162  for (uint i=0; i<dimension; i++) {
163  val += (*std_vec_)[i]*(*kxvalptr)[i];
164  }
165  return val;
166 }
167 
168 Real norm() const {
169  Real val = 0;
170  val = std::sqrt( dot(*this) );
171  return val;
172 }
173 
174 Teuchos::RCP<Vector<Real> > clone() const {
175  return Teuchos::rcp( new OptStdVector( Teuchos::rcp( new vector(std_vec_->size()) ),fd_ ) );
176 }
177 
178 Teuchos::RCP<const vector> getVector() const {
179  return std_vec_;
180 }
181 
182 Teuchos::RCP<vector> getVector() {
183  return std_vec_;
184 }
185 
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_ ) );
189  (*e_rcp)[i]= 1.0;
190  return e;
191 }
192 
193 int dimension() const {return static_cast<int>(std_vec_->size());}
194 
195 
197 const Vector<Real> & dual() const {
198  Teuchos::RCP<vector> dual_vecp = Teuchos::rcp(new vector(*std_vec_));
199  dual_vec_ = Teuchos::rcp( new OptDualStdVector<Real>( dual_vecp, fd_ ) );
200  fd_->apply(dual_vecp);
201  return *dual_vec_;
202 }
203 
204 }; // class OptStdVector
205 
206 
207 // Dual optimization space.
208 template <class Real, class Element>
209 class OptDualStdVector : public Vector<Real> {
210 
211  typedef std::vector<Element> vector;
212  typedef typename vector::size_type uint;
213 
214 private:
215 Teuchos::RCP<std::vector<Element> > std_vec_;
216 mutable Teuchos::RCP<OptStdVector<Real> > dual_vec_;
217 Teuchos::RCP<FiniteDifference<Real> > fd_;
218 
219 public:
220 
221 OptDualStdVector(const Teuchos::RCP<std::vector<Element> > & std_vec, Teuchos::RCP<FiniteDifference<Real> >fd) :
222  std_vec_(std_vec), dual_vec_(Teuchos::null), fd_(fd) {}
223 
224 void plus( const Vector<Real> &x ) {
225  const OptDualStdVector &ex = Teuchos::dyn_cast<const OptDualStdVector>(x);
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];
230  }
231 }
232 
233 void scale( const Real alpha ) {
234  uint dimension = std_vec_->size();
235  for (uint i=0; i<dimension; i++) {
236  (*std_vec_)[i] *= alpha;
237  }
238 }
239 
240 Real dot( const Vector<Real> &x ) const {
241  Real val = 0;
242  const OptDualStdVector<Real, Element> & ex = Teuchos::dyn_cast<const OptDualStdVector<Real, Element> >(x);
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);
246 
247  uint dimension = std_vec_->size();
248  for (unsigned i=0; i<dimension; i++) {
249  val += (*std_vec_)[i]*(*xvalptr)[i];
250  }
251  return val;
252 }
253 
254 Real norm() const {
255  Real val = 0;
256  val = std::sqrt( dot(*this) );
257  return val;
258 }
259 
260 Teuchos::RCP<Vector<Real> > clone() const {
261  return Teuchos::rcp( new OptDualStdVector( Teuchos::rcp( new std::vector<Element>(std_vec_->size()) ), fd_ ) );
262 }
263 
264 Teuchos::RCP<const std::vector<Element> > getVector() const {
265  return std_vec_;
266 }
267 
268 Teuchos::RCP<std::vector<Element> > getVector() {
269  return std_vec_;
270 }
271 
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_ ) );
275  (*e_rcp)[i] = 1.0;
276  return e;
277 }
278 
279 int dimension() const {return static_cast<int>(std_vec_->size());}
280 
281 const Vector<Real> & dual() const {
282  Teuchos::RCP<vector> dual_vecp = Teuchos::rcp(new vector(*std_vec_));
283  dual_vec_ = Teuchos::rcp( new OptStdVector<Real>( dual_vecp, fd_ ) );
284 
285  fd_->solve(dual_vecp);
286  return *dual_vec_;
287 }
288 
289 }; // class OptDualStdVector
290 
291 
292 
293 
294 // Constraint space.
295 template <class Real, class Element>
296 class ConStdVector : public Vector<Real> {
297 
298  typedef std::vector<Element> vector;
299  typedef typename vector::size_type uint;
300 
301 private:
302 Teuchos::RCP<std::vector<Element> > std_vec_;
303 mutable Teuchos::RCP<ConDualStdVector<Real> > dual_vec_;
304 
305 public:
306 
307 ConStdVector(const Teuchos::RCP<std::vector<Element> > & std_vec) : std_vec_(std_vec), dual_vec_(Teuchos::null) {}
308 
309 void plus( const Vector<Real> &x ) {
310  const ConStdVector &ex = Teuchos::dyn_cast<const ConStdVector>(x);
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];
315  }
316 }
317 
318 void scale( const Real alpha ) {
319  uint dimension = std_vec_->size();
320  for (uint i=0; i<dimension; i++) {
321  (*std_vec_)[i] *= alpha;
322  }
323 }
324 
325 Real dot( const Vector<Real> &x ) const {
326  Real val = 0;
327  const ConStdVector<Real, Element> & ex = Teuchos::dyn_cast<const ConStdVector<Real, Element> >(x);
328  Teuchos::RCP<const vector> xvalptr = ex.getVector();
329 
330  uint dimension = std_vec_->size();
331  for (uint i=0; i<dimension; i++) {
332  val += (*std_vec_)[i]*(*xvalptr)[i];
333  }
334  return val;
335 }
336 
337 Real norm() const {
338  Real val = 0;
339  val = std::sqrt( dot(*this) );
340  return val;
341 }
342 
343 Teuchos::RCP<Vector<Real> > clone() const {
344  return Teuchos::rcp( new ConStdVector( Teuchos::rcp(new vector(std_vec_->size())) ) );
345 }
346 
347 Teuchos::RCP<const std::vector<Element> > getVector() const {
348  return std_vec_;
349 }
350 
351 Teuchos::RCP<std::vector<Element> > getVector() {
352  return std_vec_;
353 }
354 
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) );
358  (*e_rcp)[i] = 1.0;
359  return e;
360 }
361 
362 int dimension() const {return static_cast<int>(std_vec_->size());}
363 
364 const Vector<Real> & dual() const {
365  dual_vec_ = Teuchos::rcp( new ConDualStdVector<Real>( Teuchos::rcp( new std::vector<Element>(*std_vec_) ) ) );
366  return *dual_vec_;
367 }
368 
369 }; // class ConStdVector
370 
371 
372 // Dual constraint space.
373 template <class Real, class Element>
374 class ConDualStdVector : public Vector<Real> {
375 
376  typedef std::vector<Element> vector;
377  typedef typename vector::size_type uint;
378 
379 private:
380 
381 Teuchos::RCP<std::vector<Element> > std_vec_;
382 mutable Teuchos::RCP<ConStdVector<Real> > dual_vec_;
383 
384 public:
385 
386 ConDualStdVector(const Teuchos::RCP<std::vector<Element> > & std_vec) : std_vec_(std_vec), dual_vec_(Teuchos::null) {}
387 
388 void plus( const Vector<Real> &x ) {
389  const ConDualStdVector &ex = Teuchos::dyn_cast<const ConDualStdVector>(x);
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];
394  }
395 }
396 
397 void scale( const Real alpha ) {
398  uint dimension = std_vec_->size();
399  for (uint i=0; i<dimension; i++) {
400  (*std_vec_)[i] *= alpha;
401  }
402 }
403 
404 Real dot( const Vector<Real> &x ) const {
405  Real val = 0;
406  const ConDualStdVector<Real, Element> & ex = Teuchos::dyn_cast<const ConDualStdVector<Real, Element> >(x);
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];
411  }
412  return val;
413 }
414 
415 Real norm() const {
416  Real val = 0;
417  val = std::sqrt( dot(*this) );
418  return val;
419 }
420 
421 Teuchos::RCP<Vector<Real> > clone() const {
422  return Teuchos::rcp( new ConDualStdVector( Teuchos::rcp(new std::vector<Element>(std_vec_->size())) ) );
423 }
424 
425 Teuchos::RCP<const std::vector<Element> > getVector() const {
426  return std_vec_;
427 }
428 
429 Teuchos::RCP<std::vector<Element> > getVector() {
430  return std_vec_;
431 }
432 
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 ) );
436  (*e_rcp)[i] = 1.0;
437  return e;
438 }
439 
440 int dimension() const {return static_cast<int>(std_vec_->size());}
441 
442 const Vector<Real> & dual() const {
443  dual_vec_ = Teuchos::rcp( new ConStdVector<Real>( Teuchos::rcp( new std::vector<Element>(*std_vec_) ) ) );
444  return *dual_vec_;
445 }
446 
447 }; // class ConDualStdVector
448 
449 /*** End of declaration of four vector spaces. ***/
450 
451 
452 
454 template<class Real, class XPrim=StdVector<Real>, class XDual=StdVector<Real> >
455 class Objective_GrossPitaevskii : public Objective<Real> {
456 
457  typedef std::vector<Real> vector;
458  typedef typename vector::size_type uint;
459 
460  private:
461 
463  Real g_;
464 
466  uint nx_;
467 
469  Real dx_;
470 
472  Teuchos::RCP<const std::vector<Real> > Vp_;
473 
474  Teuchos::RCP<FiniteDifference<Real> > fd_;
475 
477 
479  void applyK(const Vector<Real> &v, Vector<Real> &kv) {
480 
481  using Teuchos::RCP; using Teuchos::dyn_cast;
482 
483  // Pointer to direction vector
484  RCP<const vector> vp = dyn_cast<const XPrim>(v).getVector();
485 
486  // Pointer to action of Hessian on direction vector
487  RCP<vector> kvp = dyn_cast<XDual>(kv).getVector();
488 
489  Real dx2 = dx_*dx_;
490 
491  (*kvp)[0] = (2.0*(*vp)[0]-(*vp)[1])/dx2;
492 
493  for(uint i=1;i<nx_-1;++i) {
494  (*kvp)[i] = (2.0*(*vp)[i]-(*vp)[i-1]-(*vp)[i+1])/dx2;
495  }
496 
497  (*kvp)[nx_-1] = (2.0*(*vp)[nx_-1]-(*vp)[nx_-2])/dx2;
498 
499  }
500 
501  public:
502 
503  Objective_GrossPitaevskii(const Real &g, const Vector<Real> &V, Teuchos::RCP<FiniteDifference<Real> > fd) : g_(g),
504  Vp_((Teuchos::dyn_cast<const StdVector<Real> >(V)).getVector()), fd_(fd) {
505 
506  nx_ = Vp_->size();
507  dx_ = (1.0/(1.0+nx_));
508  }
509 
511 
515  Real value( const Vector<Real> &psi, Real &tol ) {
516 
517  using Teuchos::RCP; using Teuchos::rcp; using Teuchos::dyn_cast;
518 
519  // Pointer to opt vector
520  RCP<const vector> psip = dyn_cast<const XPrim>(psi).getVector();
521 
522  // Pointer to K applied to opt vector
523  RCP<vector> kpsip = rcp( new vector(nx_, 0.0) );
524  XDual kpsi(kpsip,fd_);
525 
526  Real J = 0;
527 
528  applyK(psi,kpsi);
529 
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);
532  }
533 
534  // Rescale for numerical integration by trapezoidal rule
535  J *= 0.5*dx_;
536 
537  return J;
538  }
539 
541 
542  void gradient( Vector<Real> &g, const Vector<Real> &psi, Real &tol ) {
543 
544  using Teuchos::RCP; using Teuchos::rcp; using Teuchos::dyn_cast;
545 
546  // Pointer to opt vector
547  RCP<const vector> psip = dyn_cast<const XPrim>(psi).getVector();
548 
549  // Pointer to gradient vector
550  RCP<vector> gp = dyn_cast<XDual>(g).getVector();
551 
552  // Pointer to K applied to opt vector
553  RCP<vector> kpsip = rcp( new vector(nx_, 0.0) );
554  XDual kpsi(kpsip,fd_);
555 
556  applyK(psi,kpsi);
557 
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_;
560  }
561 
562  }
563 
564 
565 
567 
568  void hessVec( Vector<Real> &hv, const Vector<Real> &v, const Vector<Real> &psi, Real &tol ) {
569 
570  using Teuchos::RCP; using Teuchos::dyn_cast;
571 
572  // Pointer to opt vector
573  RCP<const vector> psip = dyn_cast<const XPrim>(psi).getVector();
574 
575  // Pointer to direction vector
576  RCP<const vector> vp = dyn_cast<const XPrim>(v).getVector();
577 
578  // Pointer to action of Hessian on direction vector
579  RCP<vector> hvp = dyn_cast<XDual>(hv).getVector();
580 
581  applyK(v,hv);
582 
583  for(uint i=0;i<nx_;++i) {
584  (*hvp)[i] *= dx_;
585  (*hvp)[i] += ( (*Vp_)[i] + 6.0*g_*pow((*psip)[i],2) )*(*vp)[i]*dx_;
586  }
587 
588  }
589 
590 };
591 
592 
594 template<class Real, class XPrim=StdVector<Real>, class XDual=StdVector<Real>, class CPrim=StdVector<Real>, class CDual=StdVector<Real> >
595 class Normalization_Constraint : public EqualityConstraint<Real> {
596 
597  typedef std::vector<Real> vector;
598  typedef typename vector::size_type uint;
599 
600  private:
601  uint nx_;
602  Real dx_;
603  Teuchos::RCP<FiniteDifference<Real> > fd_;
604  bool exactsolve_;
605 
606  public:
607  Normalization_Constraint(int n, Real dx, Teuchos::RCP<FiniteDifference<Real> > fd, bool exactsolve) :
608  nx_(n), dx_(dx), fd_(fd), exactsolve_(exactsolve) {}
609 
611 
615  void value(Vector<Real> &c, const Vector<Real> &psi, Real &tol){
616 
617  using Teuchos::RCP; using Teuchos::dyn_cast;
618 
619  // Pointer to constraint vector (only one element)
620  RCP<vector> cp = dyn_cast<CPrim>(c).getVector();
621 
622  // Pointer to opt vector
623  RCP<const vector> psip = dyn_cast<const XPrim>(psi).getVector();
624 
625  (*cp)[0] = -1.0;
626  for(uint i=0;i<nx_;++i) {
627  (*cp)[0] += dx_*pow((*psip)[i],2);
628  }
629  }
630 
632 
634  void applyJacobian(Vector<Real> &jv, const Vector<Real> &v, const Vector<Real> &psi, Real &tol){
635 
636  using Teuchos::RCP; using Teuchos::dyn_cast;
637 
638  // Pointer to action of Jacobian of constraint on direction vector (yields scalar)
639  RCP<vector> jvp = dyn_cast<CPrim>(jv).getVector();
640 
641  // Pointer to direction vector
642  RCP<const vector> vp = dyn_cast<const XPrim>(v).getVector();
643 
644  // Pointer to opt vector
645  RCP<const vector> psip = dyn_cast<const XPrim>(psi).getVector();
646 
647  (*jvp)[0] = 0;
648  for(uint i=0;i<nx_;++i) {
649  (*jvp)[0] += 2.0*dx_*(*psip)[i]*(*vp)[i];
650  }
651  }
652 
654 
656  void applyAdjointJacobian(Vector<Real> &ajv, const Vector<Real> &v, const Vector<Real> &psi, Real &tol){
657 
658  using Teuchos::RCP; using Teuchos::dyn_cast;
659 
660  // Pointer to action of adjoint of Jacobian of constraint on direction vector (yields vector)
661  RCP<vector> ajvp = dyn_cast<XDual>(ajv).getVector();
662 
663  // Pointer to direction vector
664  RCP<const vector> vp = dyn_cast<const CDual>(v).getVector();
665 
666  // Pointer to opt vector
667  RCP<const vector> psip = dyn_cast<const XPrim>(psi).getVector();
668 
669  for(uint i=0;i<nx_;++i) {
670  (*ajvp)[i] = 2.0*dx_*(*psip)[i]*(*vp)[0];
671  }
672  }
673 
675 
679  const Vector<Real> &psi, Real &tol){
680 
681 
682  using Teuchos::RCP; using Teuchos::dyn_cast;
683 
684  // The pointer to action of constraint Hessian in u,v inner product
685  RCP<vector> ahuvp = dyn_cast<XDual>(ahuv).getVector();
686 
687  // Pointer to direction vector u
688  RCP<const vector> up = dyn_cast<const CDual>(u).getVector();
689 
690  // Pointer to direction vector v
691  RCP<const vector> vp = dyn_cast<const XPrim>(v).getVector();
692 
693  // Pointer to opt vector
694  RCP<const vector> psip = dyn_cast<const XPrim>(psi).getVector();
695 
696  for(uint i=0;i<nx_;++i) {
697  (*ahuvp)[i] = 2.0*dx_*(*vp)[i]*(*up)[0];
698  }
699  }
700 
706  std::vector<Real> solveAugmentedSystem(Vector<Real> &v1, Vector<Real> &v2, const Vector<Real> &b1,
707  const Vector<Real> &b2, const Vector<Real> &psi, Real &tol) {
708 
709  using Teuchos::RCP; using Teuchos::rcp; using Teuchos::dyn_cast;
710 
711  if(exactsolve_) {
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();
717 
718  RCP<vector> jacp = rcp( new vector(nx_, 0.0) );
719  RCP<vector> b1dp = rcp( new vector(nx_, 0.0) );
720 
721  for(uint i=0;i<nx_;++i) {
722  (*jacp)[i] = (*psip)[i];
723  (*b1dp)[i] = (*b1p)[i];
724  }
725 
726  // The Jacobian of the constraint is \f$c'(\psi)=2dx\psi\f$
727  XDual jac(jacp,fd_);
728  jac.scale(2.0*dx_);
729 
730  // A Dual-in-name-only version of b1, so we can compute the desired inner products involving inv(K)
731  XDual b1d(b1dp,fd_);
732 
733  // \f$ (c'K^{-1}*c'^\ast)^{-1} \f$
734  Real d = 1.0/jac.dot(jac);
735  Real p = jac.dot(b1d);
736 
737  (*v2p)[0] = d*(p-(*b2p)[0]);
738 
739  v1.set(jac.dual());
740  v1.scale(-(*v2p)[0]);
741  v1.plus(b1d.dual());
742 
743  return std::vector<Real>(0);
744  }
745  else{
746  return EqualityConstraint<Real>::solveAugmentedSystem(v1,v2,b1,b2,psi,tol);
747  }
748  }
749 };
750 
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)
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.
Definition: ROL_Vector.hpp:74
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 .
Definition: ROL_Vector.hpp:196
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...