1 #include "ROL_StdVector.hpp"
2 #include "Teuchos_LAPACK.hpp"
3 
4 using namespace ROL;
5 
6 template<class Real>
7 class FiniteDifference {
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 
FiniteDifference(int n,double dx)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 
44         //! Given f, compute -u''=f
solve(ROL::Ptr<const std::vector<Real>> fp,ROL::Ptr<std::vector<Real>> up)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 
52         //! Same as above but with overwrite in place
solve(ROL::Ptr<std::vector<Real>> up)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 
57         //! Given u, compute f = -u''
apply(ROL::Ptr<const std::vector<Real>> up,ROL::Ptr<std::vector<Real>> fp)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 
67         //! Same as above but with overwrite in place
apply(ROL::Ptr<std::vector<Real>> fp)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 
86