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