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