ROL
FiniteDifference.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Rapid Optimization Library (ROL) Package
4 //
5 // Copyright 2014 NTESS and the ROL contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #include "ROL_StdVector.hpp"
11 #include "Teuchos_LAPACK.hpp"
12 
13 using namespace ROL;
14 
15 template<class Real>
17 
18  private:
19  int n_;
20  double dx2_;
21  Teuchos::LAPACK<int,Real> lapack_;
22 
23  // subdiagonal is -1/dx^2
24  std::vector<Real> dl_;
25 
26  // diagonal is 2/dx^2
27  std::vector<Real> d_;
28 
29  // superdiagonal is -1/dx^2
30  std::vector<Real> du_;
31 
32  // second superdiagonal (du2 = 0)
33  std::vector<Real> du2_;
34 
35  // Pivot indices
36  std::vector<int> ipiv_;
37 
38  int info_;
39 
40 
41  public:
42 
43  FiniteDifference(int n, double dx) : n_(n),dx2_(dx*dx),
44  dl_(n_-1,-1.0/dx2_),
45  d_(n_,2.0/dx2_),
46  du_(n_-1,-1.0/dx2_),
47  du2_(n_-2,0.0),
48  ipiv_(n_,0) {
49 
50  lapack_.GTTRF(n_,&dl_[0],&d_[0],&du_[0],&du2_[0],&ipiv_[0],&info_);
51  }
52 
54  void solve(ROL::Ptr<const std::vector<Real> > fp, ROL::Ptr<std::vector<Real> > up) {
55  for(int i=0;i<n_;++i) {
56  (*up)[i] = (*fp)[i];
57  }
58  lapack_.GTTRS('n',n_,1,&dl_[0],&d_[0],&du_[0],&du2_[0],&ipiv_[0],&(*up)[0],n_,&info_);
59  }
60 
62  void solve(ROL::Ptr<std::vector<Real> > up) {
63  lapack_.GTTRS('n',n_,1,&dl_[0],&d_[0],&du_[0],&du2_[0],&ipiv_[0],&(*up)[0],n_,&info_);
64  }
65 
67  void apply(ROL::Ptr<const std::vector<Real> > up, ROL::Ptr<std::vector<Real> > fp) {
68  (*fp)[0] = (2.0*(*up)[0]-(*up)[1])/dx2_;
69 
70  for(int i=1;i<n_-1;++i) {
71  (*fp)[i] = (2.0*(*up)[i]-(*up)[i-1]-(*up)[i+1])/dx2_;
72  }
73  (*fp)[n_-1] = (2.0*(*up)[n_-1]-(*up)[n_-2])/dx2_;
74  }
75 
77  void apply(ROL::Ptr<std::vector<Real> > fp) {
78 
79  ROL::Ptr<std::vector<Real> > up = ROL::makePtr<std::vector<Real>>(n_, 0.0);
80  for(int i=0;i<n_;++i) {
81  (*up)[i] = (*fp)[i];
82  }
83 
84  (*fp)[0] = (2.0*(*up)[0]-(*up)[1])/dx2_;
85  for(int i=1;i<n_-1;++i) {
86  (*fp)[i] = (2.0*(*up)[i]-(*up)[i-1]-(*up)[i+1])/dx2_;
87  }
88  (*fp)[n_-1] = (2.0*(*up)[n_-1]-(*up)[n_-2])/dx2_;
89 
90  }
91 };
92 
93 
94 
std::vector< Real > dl_
std::vector< Real > du_
FiniteDifference(int n, double dx)
void solve(ROL::Ptr< std::vector< Real > > up)
Same as above but with overwrite in place.
std::vector< Real > d_
void solve(ROL::Ptr< const std::vector< Real > > fp, ROL::Ptr< std::vector< Real > > up)
Given f, compute -u&#39;&#39;=f.
Teuchos::LAPACK< int, Real > lapack_
std::vector< int > ipiv_
void apply(ROL::Ptr< std::vector< Real > > fp)
Same as above but with overwrite in place.
void apply(ROL::Ptr< const std::vector< Real > > up, ROL::Ptr< std::vector< Real > > fp)
Given u, compute f = -u&#39;&#39;.
std::vector< Real > du2_