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 "ROL_Stream.hpp"
85#include "Teuchos_GlobalMPISession.hpp"
86
87#include "ROL_ParameterList.hpp"
88#include "ROL_StdVector.hpp"
89#include "ROL_Objective.hpp"
90#include "ROL_Constraint.hpp"
91#include "ROL_Algorithm.hpp"
92#include "ROL_CompositeStep.hpp"
94
96
97
98using namespace ROL;
99
100template <class Real, class Element=Real>
101class OptStdVector; // Optimization space.
102
103template <class Real, class Element=Real>
104class OptDualStdVector; // Dual optimization space.
105
106template <class Real, class Element=Real>
107class ConStdVector; // Constraint space.
108
109template <class Real, class Element=Real>
110class ConDualStdVector; // Dual constraint space.
111
112// Vector space definitions:
113
114
115// Optimization space.
116template <class Real, class Element>
117class OptStdVector : public Vector<Real> {
118
119 typedef std::vector<Element> vector;
120 typedef typename vector::size_type uint;
121
122private:
123ROL::Ptr<std::vector<Element> > std_vec_;
124mutable ROL::Ptr<OptDualStdVector<Real> > dual_vec_;
125
126ROL::Ptr<FiniteDifference<Real> > fd_;
127
128
129public:
130
131OptStdVector(const ROL::Ptr<std::vector<Element> > & std_vec, ROL::Ptr<FiniteDifference<Real> >fd) :
132 std_vec_(std_vec), dual_vec_(ROL::nullPtr), fd_(fd) {}
133
134void plus( const Vector<Real> &x ) {
135 const OptStdVector &ex = dynamic_cast<const OptStdVector&>(x);
136 ROL::Ptr<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
143void 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
152Real dot( const Vector<Real> &x ) const {
153 Real val = 0;
154 const OptStdVector<Real, Element> & ex = dynamic_cast<const OptStdVector&>(x);
155 ROL::Ptr<const vector> xvalptr = ex.getVector();
156
157 ROL::Ptr<vector> kxvalptr = ROL::makePtr<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
168Real norm() const {
169 Real val = 0;
170 val = std::sqrt( dot(*this) );
171 return val;
172}
173
174ROL::Ptr<Vector<Real> > clone() const {
175 return ROL::makePtr<OptStdVector>( ROL::makePtr<vector>(std_vec_->size()), fd_ );
176}
177
178ROL::Ptr<const vector> getVector() const {
179 return std_vec_;
180}
181
182ROL::Ptr<vector> getVector() {
183 return std_vec_;
184}
185
186ROL::Ptr<Vector<Real> > basis( const int i ) const {
187 ROL::Ptr<vector> e_ptr = ROL::makePtr<vector>(std_vec_->size(),0.0);
188 ROL::Ptr<OptStdVector> e = ROL::makePtr<OptStdVector>( e_ptr, fd_ );
189 (*e_ptr)[i]= 1.0;
190 return e;
191}
192
193int dimension() const {return static_cast<int>(std_vec_->size());}
194
195
197const Vector<Real> & dual() const {
198 ROL::Ptr<vector> dual_vecp = ROL::makePtr<vector>(*std_vec_);
199 dual_vec_ = ROL::makePtr<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.
208template <class Real, class Element>
209class OptDualStdVector : public Vector<Real> {
210
211 typedef std::vector<Element> vector;
212 typedef typename vector::size_type uint;
213
214private:
215ROL::Ptr<std::vector<Element> > std_vec_;
216mutable ROL::Ptr<OptStdVector<Real> > dual_vec_;
217ROL::Ptr<FiniteDifference<Real> > fd_;
218
219public:
220
221OptDualStdVector(const ROL::Ptr<std::vector<Element> > & std_vec, ROL::Ptr<FiniteDifference<Real> >fd) :
222 std_vec_(std_vec), dual_vec_(ROL::nullPtr), fd_(fd) {}
223
224void plus( const Vector<Real> &x ) {
225 const OptDualStdVector &ex = dynamic_cast<const OptDualStdVector&>(x);
226 ROL::Ptr<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
233void 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
240Real dot( const Vector<Real> &x ) const {
241 Real val = 0;
242 const OptDualStdVector<Real, Element> & ex = dynamic_cast<const OptDualStdVector<Real, Element>&>(x);
243 ROL::Ptr<const vector> kxvalptr = ex.getVector();
244 ROL::Ptr<vector> xvalptr = ROL::makePtr<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
254Real norm() const {
255 Real val = 0;
256 val = std::sqrt( dot(*this) );
257 return val;
258}
259
260ROL::Ptr<Vector<Real> > clone() const {
261 return ROL::makePtr<OptDualStdVector>( ROL::makePtr<std::vector<Element>>(std_vec_->size()), fd_ );
262}
263
264ROL::Ptr<const std::vector<Element> > getVector() const {
265 return std_vec_;
266}
267
268ROL::Ptr<std::vector<Element> > getVector() {
269 return std_vec_;
270}
271
272ROL::Ptr<Vector<Real> > basis( const int i ) const {
273 ROL::Ptr<vector> e_ptr = ROL::makePtr<vector>(std_vec_->size(), 0.0 );
274 ROL::Ptr<OptDualStdVector> e = ROL::makePtr<OptDualStdVector>( e_ptr,fd_ );
275 (*e_ptr)[i] = 1.0;
276 return e;
277}
278
279int dimension() const {return static_cast<int>(std_vec_->size());}
280
281const Vector<Real> & dual() const {
282 ROL::Ptr<vector> dual_vecp = ROL::makePtr<vector>(*std_vec_);
283 dual_vec_ = ROL::makePtr<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.
295template <class Real, class Element>
296class ConStdVector : public Vector<Real> {
297
298 typedef std::vector<Element> vector;
299 typedef typename vector::size_type uint;
300
301private:
302ROL::Ptr<std::vector<Element> > std_vec_;
303mutable ROL::Ptr<ConDualStdVector<Real> > dual_vec_;
304
305public:
306
307ConStdVector(const ROL::Ptr<std::vector<Element> > & std_vec) : std_vec_(std_vec), dual_vec_(ROL::nullPtr) {}
308
309void plus( const Vector<Real> &x ) {
310 const ConStdVector &ex = dynamic_cast<const ConStdVector&>(x);
311 ROL::Ptr<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
318void 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
325Real dot( const Vector<Real> &x ) const {
326 Real val = 0;
327 const ConStdVector<Real, Element> & ex = dynamic_cast<const ConStdVector<Real, Element>&>(x);
328 ROL::Ptr<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
337Real norm() const {
338 Real val = 0;
339 val = std::sqrt( dot(*this) );
340 return val;
341}
342
343ROL::Ptr<Vector<Real> > clone() const {
344 return ROL::makePtr<ConStdVector>( ROL::makePtr<vector>(std_vec_->size()));
345}
346
347ROL::Ptr<const std::vector<Element> > getVector() const {
348 return std_vec_;
349}
350
351ROL::Ptr<std::vector<Element> > getVector() {
352 return std_vec_;
353}
354
355ROL::Ptr<Vector<Real> > basis( const int i ) const {
356 ROL::Ptr<vector> e_ptr = ROL::makePtr<vector>(std_vec_->size(),0.0);
357 ROL::Ptr<ConStdVector> e = ROL::makePtr<ConStdVector>( e_ptr);
358 (*e_ptr)[i] = 1.0;
359 return e;
360}
361
362int dimension() const {return static_cast<int>(std_vec_->size());}
363
364const Vector<Real> & dual() const {
365 dual_vec_ = ROL::makePtr<ConDualStdVector<Real>>( ROL::makePtr<std::vector<Element>>(*std_vec_) );
366 return *dual_vec_;
367}
368
369}; // class ConStdVector
370
371
372// Dual constraint space.
373template <class Real, class Element>
374class ConDualStdVector : public Vector<Real> {
375
376 typedef std::vector<Element> vector;
377 typedef typename vector::size_type uint;
378
379private:
380
381ROL::Ptr<std::vector<Element> > std_vec_;
382mutable ROL::Ptr<ConStdVector<Real> > dual_vec_;
383
384public:
385
386ConDualStdVector(const ROL::Ptr<std::vector<Element> > & std_vec) : std_vec_(std_vec), dual_vec_(ROL::nullPtr) {}
387
388void plus( const Vector<Real> &x ) {
389 const ConDualStdVector &ex = dynamic_cast<const ConDualStdVector&>(x);
390 ROL::Ptr<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
397void 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
404Real dot( const Vector<Real> &x ) const {
405 Real val = 0;
406 const ConDualStdVector<Real, Element> & ex = dynamic_cast<const ConDualStdVector<Real, Element>&>(x);
407 ROL::Ptr<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
415Real norm() const {
416 Real val = 0;
417 val = std::sqrt( dot(*this) );
418 return val;
419}
420
421ROL::Ptr<Vector<Real> > clone() const {
422 return ROL::makePtr<ConDualStdVector>( ROL::makePtr<std::vector<Element>>(std_vec_->size()));
423}
424
425ROL::Ptr<const std::vector<Element> > getVector() const {
426 return std_vec_;
427}
428
429ROL::Ptr<std::vector<Element> > getVector() {
430 return std_vec_;
431}
432
433ROL::Ptr<Vector<Real> > basis( const int i ) const {
434 ROL::Ptr<vector> e_ptr = ROL::makePtr<vector>(std_vec_->size(),0.0);
435 ROL::Ptr<ConDualStdVector> e = ROL::makePtr<ConDualStdVector>( e_ptr );
436 (*e_ptr)[i] = 1.0;
437 return e;
438}
439
440int dimension() const {return static_cast<int>(std_vec_->size());}
441
442const Vector<Real> & dual() const {
443 dual_vec_ = ROL::makePtr<ConStdVector<Real>>( ROL::makePtr<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
454template<class Real, class XPrim=StdVector<Real>, class XDual=StdVector<Real> >
455class 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 ROL::Ptr<const std::vector<Real> > Vp_;
473
474 ROL::Ptr<FiniteDifference<Real> > fd_;
475
477
479 void applyK(const Vector<Real> &v, Vector<Real> &kv) {
480
481
482
483 // Pointer to direction vector
484 ROL::Ptr<const vector> vp = dynamic_cast<const XPrim&>(v).getVector();
485
486 // Pointer to action of Hessian on direction vector
487 ROL::Ptr<vector> kvp = dynamic_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, ROL::Ptr<FiniteDifference<Real> > fd) : g_(g),
504 Vp_((dynamic_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
518
519 // Pointer to opt vector
520 ROL::Ptr<const vector> psip = dynamic_cast<const XPrim&>(psi).getVector();
521
522 // Pointer to K applied to opt vector
523 ROL::Ptr<vector> kpsip = ROL::makePtr<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
545
546 // Pointer to opt vector
547 ROL::Ptr<const vector> psip = dynamic_cast<const XPrim&>(psi).getVector();
548
549 // Pointer to gradient vector
550 ROL::Ptr<vector> gp = dynamic_cast<XDual&>(g).getVector();
551
552 // Pointer to K applied to opt vector
553 ROL::Ptr<vector> kpsip = ROL::makePtr<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
571
572 // Pointer to opt vector
573 ROL::Ptr<const vector> psip = dynamic_cast<const XPrim&>(psi).getVector();
574
575 // Pointer to direction vector
576 ROL::Ptr<const vector> vp = dynamic_cast<const XPrim&>(v).getVector();
577
578 // Pointer to action of Hessian on direction vector
579 ROL::Ptr<vector> hvp = dynamic_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
594template<class Real, class XPrim=StdVector<Real>, class XDual=StdVector<Real>, class CPrim=StdVector<Real>, class CDual=StdVector<Real> >
595class Normalization_Constraint : public Constraint<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 ROL::Ptr<FiniteDifference<Real> > fd_;
605
606 public:
607 Normalization_Constraint(int n, Real dx, ROL::Ptr<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
618
619 // Pointer to constraint vector (only one element)
620 ROL::Ptr<vector> cp = dynamic_cast<CPrim&>(c).getVector();
621
622 // Pointer to opt vector
623 ROL::Ptr<const vector> psip = dynamic_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
637
638 // Pointer to action of Jacobian of constraint on direction vector (yields scalar)
639 ROL::Ptr<vector> jvp = dynamic_cast<CPrim&>(jv).getVector();
640
641 // Pointer to direction vector
642 ROL::Ptr<const vector> vp = dynamic_cast<const XPrim&>(v).getVector();
643
644 // Pointer to opt vector
645 ROL::Ptr<const vector> psip = dynamic_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
659
660 // Pointer to action of adjoint of Jacobian of constraint on direction vector (yields vector)
661 ROL::Ptr<vector> ajvp = dynamic_cast<XDual&>(ajv).getVector();
662
663 // Pointer to direction vector
664 ROL::Ptr<const vector> vp = dynamic_cast<const CDual&>(v).getVector();
665
666 // Pointer to opt vector
667 ROL::Ptr<const vector> psip = dynamic_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
683
684 // The pointer to action of constraint Hessian in u,v inner product
685 ROL::Ptr<vector> ahuvp = dynamic_cast<XDual&>(ahuv).getVector();
686
687 // Pointer to direction vector u
688 ROL::Ptr<const vector> up = dynamic_cast<const CDual&>(u).getVector();
689
690 // Pointer to direction vector v
691 ROL::Ptr<const vector> vp = dynamic_cast<const XPrim&>(v).getVector();
692
693 // Pointer to opt vector
694 ROL::Ptr<const vector> psip = dynamic_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
710
711 if(exactsolve_) {
712 ROL::Ptr<vector> v1p = dynamic_cast<XPrim&>(v1).getVector();
713 ROL::Ptr<vector> v2p = dynamic_cast<CDual&>(v2).getVector();
714 ROL::Ptr<const vector> b1p = dynamic_cast<const XDual&>(b1).getVector();
715 ROL::Ptr<const vector> b2p = dynamic_cast<const CPrim&>(b2).getVector();
716 ROL::Ptr<const vector> psip = dynamic_cast<const XPrim&>(psi).getVector();
717
718 ROL::Ptr<vector> jacp = ROL::makePtr<vector>(nx_, 0.0);
719 ROL::Ptr<vector> b1dp = ROL::makePtr<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 Constraint<Real>::solveAugmentedSystem(v1,v2,b1,b2,psi,tol);
747 }
748 }
749};
750
Defines a no-output stream class ROL::NullStream and a function makeStreamPtr which either wraps a re...
ROL::Ptr< const std::vector< Element > > getVector() const
ROL::Ptr< ConStdVector< Real > > dual_vec_
Real dot(const Vector< Real > &x) const
Compute where .
ROL::Ptr< Vector< Real > > basis(const int i) const
Return i-th basis vector.
std::vector< Element > vector
int dimension() const
Return dimension of the vector space.
void plus(const Vector< Real > &x)
Compute , where .
Real norm() const
Returns where .
ConDualStdVector(const ROL::Ptr< std::vector< Element > > &std_vec)
void scale(const Real alpha)
Compute where .
const Vector< Real > & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis,...
ROL::Ptr< std::vector< Element > > getVector()
ROL::Ptr< std::vector< Element > > std_vec_
Real dot(const ROL::Vector< Real > &x) const
Compute where .
ROL::Ptr< Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
Real dot(const Vector< Real > &x) const
Compute where .
Real norm() const
Returns where .
ROL::Ptr< ConDualStdVector< Real > > dual_vec_
ROL::Ptr< std::vector< Element > > std_vec_
Real dot(const ROL::Vector< Real > &x) const
Compute where .
ROL::Ptr< Vector< Real > > basis(const int i) const
Return i-th basis vector.
void scale(const Real alpha)
Compute where .
std::vector< Element > vector
ConStdVector(const ROL::Ptr< std::vector< Element > > &std_vec)
ROL::Ptr< Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
ROL::Ptr< std::vector< Element > > getVector()
ROL::Ptr< const std::vector< Element > > getVector() const
void plus(const Vector< Real > &x)
Compute , where .
const Vector< Real > & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis,...
int dimension() const
Return dimension of the vector space.
ROL::Ptr< FiniteDifference< Real > > fd_
std::vector< Real > solveAugmentedSystem(Vector< Real > &v1, Vector< Real > &v2, const Vector< Real > &b1, const Vector< Real > &b2, const Vector< Real > &psi, Real &tol)
Normalization_Constraint(int n, Real dx, ROL::Ptr< FiniteDifference< Real > > fd, bool exactsolve)
void applyJacobian(Vector< Real > &jv, const Vector< Real > &v, const Vector< Real > &psi, Real &tol)
Evaluate .
void value(Vector< Real > &c, const Vector< Real > &psi, Real &tol)
Evaluate .
void applyAdjointHessian(Vector< Real > &ahuv, const Vector< Real > &u, const Vector< Real > &v, const Vector< Real > &psi, Real &tol)
Evaluate .
ROL::Ptr< const vector > getVector(const V &x)
void applyAdjointJacobian(Vector< Real > &ajv, const Vector< Real > &v, const Vector< Real > &psi, Real &tol)
Evaluate .
void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &psi, Real &tol)
Evaluate .
Objective_GrossPitaevskii(const Real &g, const Vector< Real > &V, ROL::Ptr< FiniteDifference< Real > > fd)
ROL::Ptr< const std::vector< Real > > Vp_
ROL::Ptr< FiniteDifference< Real > > fd_
ROL::Ptr< const vector > getVector(const V &x)
void applyK(const Vector< Real > &v, Vector< Real > &kv)
Apply finite difference operator.
void gradient(Vector< Real > &g, const Vector< Real > &psi, Real &tol)
Evaluate .
Real value(const Vector< Real > &psi, Real &tol)
Evaluate .
Real norm() const
Returns where .
int dimension() const
Return dimension of the vector space.
Real dot(const Vector< Real > &x) const
Compute where .
ROL::Ptr< Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
void scale(const Real alpha)
Compute where .
Real dot(const ROL::Vector< Real > &x) const
Compute where .
ROL::Ptr< std::vector< Element > > std_vec_
ROL::Ptr< const std::vector< Element > > getVector() const
ROL::Ptr< std::vector< Element > > getVector()
ROL::Ptr< Vector< Real > > basis(const int i) const
Return i-th basis vector.
const Vector< Real > & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis,...
void plus(const Vector< Real > &x)
Compute , where .
OptDualStdVector(const ROL::Ptr< std::vector< Element > > &std_vec, ROL::Ptr< FiniteDifference< Real > >fd)
std::vector< Element > vector
ROL::Ptr< OptStdVector< Real > > dual_vec_
ROL::Ptr< FiniteDifference< Real > > fd_
OptStdVector(const ROL::Ptr< std::vector< Element > > &std_vec, ROL::Ptr< FiniteDifference< Real > >fd)
void plus(const Vector< Real > &x)
Compute , where .
ROL::Ptr< vector > getVector()
ROL::Ptr< const std::vector< Element > > getVector() const
ROL::Ptr< Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
Real dot(const ROL::Vector< Real > &x) const
Compute where .
int dimension() const
Return dimension of the vector space.
ROL::Ptr< Vector< Real > > basis(const int i) const
Return i-th basis vector.
const Vector< Real > & dual() const
Modify the dual of vector u to be .
Real dot(const Vector< Real > &x) const
Modify the dot product between primal variables to be .
ROL::Ptr< FiniteDifference< Real > > fd_
Real norm() const
Returns where .
ROL::Ptr< const vector > getVector() const
ROL::Ptr< OptDualStdVector< Real > > dual_vec_
ROL::Ptr< std::vector< Element > > std_vec_
void scale(const Real alpha)
Compute where .
std::vector< Element > vector
Defines the general constraint operator interface.
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
Provides the interface to evaluate objective functions.
Provides the ROL::Vector interface for scalar values, to be used, for example, with scalar constraint...
Defines the linear algebra or vector space interface.
virtual void set(const Vector &x)
Set where .
virtual void scale(const Real alpha)=0
Compute where .
virtual void plus(const Vector &x)=0
Compute , where .