ROL
InnerProductMatrix.hpp
Go to the documentation of this file.
1 #include<iostream>
2 #include<vector>
3 #include"ROL_StdVector.hpp"
4 #include"Teuchos_LAPACK.hpp"
5 
6 #ifndef __INNER_PRODUCT_MATRIX__
7 #define __INNER_PRODUCT_MATRIX__
8 
9 
10 template<class Real>
11 class InnerProductMatrix{
12  public:
13  InnerProductMatrix(const std::vector<Real> &U,
14  const std::vector<Real> &V,
15  const std::vector<Real> &w,
16  const int a=1);
17 
18  InnerProductMatrix(const std::vector<Real> &U,
19  const std::vector<Real> &V,
20  const std::vector<Real> &w,
21  const std::vector<Real> &a);
22 
23  InnerProductMatrix(InnerProductMatrix<Real> *ipm);
24 
25  void update(const std::vector<Real> &a);
26 
27  virtual ~InnerProductMatrix();
28 
29  void apply(ROL::Ptr<const std::vector<Real> > xp,
30  ROL::Ptr<std::vector<Real> > bp);
31 
32  void applyadd(ROL::Ptr<const std::vector<Real> > xp,
33  ROL::Ptr<std::vector<Real> > bp);
34 
35  void applyaddtimes(ROL::Ptr<const std::vector<Real> > xp,
36  ROL::Ptr<std::vector<Real> > bp, Real factor);
37 
38 
39  Real inner(ROL::Ptr<const std::vector<Real> > up,
40  ROL::Ptr<const std::vector<Real> > vp);
41 
42  // This method does nothing in the base class
43  virtual void solve(ROL::Ptr<const std::vector<Real> > bp,
44  ROL::Ptr<std::vector<Real> > xp){};
45 
46  // This method does nothing in the base class
47  virtual Real inv_inner(ROL::Ptr<const std::vector<Real> > up,
48  ROL::Ptr<const std::vector<Real> > vp){
49  return 0;}
50 
51  protected:
52  const int nq_;
53  const int ni_;
54  const std::vector<Real> U_;
55  const std::vector<Real> V_;
56  const std::vector<Real> w_;
57 
58  std::vector<Real> M_;
59 
60 
61 };
62 
63 
64 
65 template<class Real>
66 InnerProductMatrix<Real>::InnerProductMatrix( const std::vector<Real> &U,
67  const std::vector<Real> &V,
68  const std::vector<Real> &w,
69  const int a) :
70  nq_(w.size()), ni_(U.size()/nq_),
71  U_(U),V_(V),w_(w),M_(ni_*ni_,0) {
72  for(int i=0;i<ni_;++i) {
73  for(int j=0;j<ni_;++j) {
74  for(int k=0;k<nq_;++k) {
75  M_[i+j*ni_] += a*w_[k]*U_[k+i*nq_]*V_[k+j*nq_];
76  }
77  }
78  }
79 }
80 
81 template<class Real>
82 InnerProductMatrix<Real>::InnerProductMatrix( const std::vector<Real> &U,
83  const std::vector<Real> &V,
84  const std::vector<Real> &w,
85  const std::vector<Real> &a ) :
86  nq_(w.size()), ni_(U.size()/nq_),
87  U_(U),V_(V),w_(w),M_(ni_*ni_,0) {
88  for(int i=0;i<ni_;++i) {
89  for(int j=0;j<ni_;++j) {
90  for(int k=0;k<nq_;++k) {
91  M_[i+j*ni_] += a[k]*w_[k]*U_[k+i*nq_]*V_[k+j*nq_];
92  }
93  }
94  }
95 }
96 
97 template<class Real>
98 InnerProductMatrix<Real>::~InnerProductMatrix(){
99 }
100 
101 template<class Real>
102 void InnerProductMatrix<Real>::apply(ROL::Ptr<const std::vector<Real> > xp,
103  ROL::Ptr<std::vector<Real> > bp ) {
104  for(int i=0;i<ni_;++i) {
105  (*bp)[i] = 0;
106  for(int j=0;j<ni_;++j ) {
107  (*bp)[i] += M_[i+ni_*j]*(*xp)[j];
108  }
109  }
110 }
111 
112 template<class Real>
113 void InnerProductMatrix<Real>::applyadd(ROL::Ptr<const std::vector<Real> > xp,
114  ROL::Ptr<std::vector<Real> > bp ) {
115  for(int i=0;i<ni_;++i) {
116  for(int j=0;j<ni_;++j ) {
117  (*bp)[i] += M_[i+ni_*j]*(*xp)[j];
118  }
119  }
120 }
121 
122 template<class Real>
123 void InnerProductMatrix<Real>::applyaddtimes(ROL::Ptr<const std::vector<Real> > xp,
124  ROL::Ptr<std::vector<Real> > bp, Real factor ) {
125  for(int i=0;i<ni_;++i) {
126  for(int j=0;j<ni_;++j ) {
127  (*bp)[i] += factor*M_[i+ni_*j]*(*xp)[j];
128  }
129  }
130 }
131 
132 template<class Real>
133 void InnerProductMatrix<Real>::update( const std::vector<Real> &a ){
134 
135  std::fill(M_.begin(),M_.end(),0);
136  for(int i=0;i<ni_;++i) {
137  for(int j=0;j<ni_;++j) {
138  for(int k=0;k<nq_;++k) {
139  M_[i+j*ni_] += a[k]*w_[k]*U_[k+i*nq_]*V_[k+j*nq_];
140  }
141  }
142  }
143 }
144 
145 
146 
147 
148 
150 template<class Real>
151 Real InnerProductMatrix<Real>::inner( ROL::Ptr<const std::vector<Real> > up,
152  ROL::Ptr<const std::vector<Real> > vp ) {
153  Real J = 0;
154  ROL::Ptr<std::vector<Real> > Mvp = ROL::makePtr<std::vector<Real>>(ni_,0);
155  this->apply(vp,Mvp);
156  for(int i=0;i<ni_;++i) {
157  J += (*up)[i]*(*Mvp)[i];
158  }
159  return J;
160 }
161 
162 
163 
164 
165 
167 template<class Real>
168 class InnerProductMatrixSolver : public InnerProductMatrix<Real> {
169 
170  private:
171  ROL::Ptr<Teuchos::LAPACK<int,Real> > lapack_;
172  const int ni_;
173  const int nq_;
174  std::vector<Real> M_;
175  const char TRANS_;
176  std::vector<int> ipiv_;
177  std::vector<Real> PLU_;
178  const int nrhs_;
179  int info_;
180 
181  // Solve the system Ax=b for x
182  public:
183  InnerProductMatrixSolver(ROL::Ptr<Teuchos::LAPACK<int,Real> > lapack,
184  const std::vector<Real> &U=std::vector<Real>(),
185  const std::vector<Real> &V=std::vector<Real>(),
186  const std::vector<Real> &w=std::vector<Real>(),
187  const int a=1);
188 
189  InnerProductMatrixSolver(ROL::Ptr<Teuchos::LAPACK<int,Real> > lapack,
190  const std::vector<Real> &U=std::vector<Real>(),
191  const std::vector<Real> &V=std::vector<Real>(),
192  const std::vector<Real> &w=std::vector<Real>(),
193  const std::vector<Real> &a=std::vector<Real>());
194 
195  void solve(ROL::Ptr<const std::vector<Real> > bp,
196  ROL::Ptr<std::vector<Real> > xp);
197 
198  Real inv_inner(ROL::Ptr<const std::vector<Real> > up,
199  ROL::Ptr<const std::vector<Real> > vp);
200 };
201 
202 
203 template<class Real>
204 InnerProductMatrixSolver<Real>::InnerProductMatrixSolver(ROL::Ptr<Teuchos::LAPACK<int,Real> > lapack,
205  const std::vector<Real> &U,
206  const std::vector<Real> &V,
207  const std::vector<Real> &w,
208  const int a):
209  InnerProductMatrix<Real>(U,V,w,a),
210  lapack_(lapack),
211  ni_(InnerProductMatrix<Real>::ni_),
212  nq_(InnerProductMatrix<Real>::nq_),
213  M_(InnerProductMatrix<Real>::M_),
214  TRANS_('N'), ipiv_(ni_,0), PLU_(ni_*ni_,0),
215  nrhs_(1),info_(0){
216  PLU_ = M_;
217 
218  // Do matrix factorization
219  lapack->GETRF(ni_,ni_,&PLU_[0],ni_,&ipiv_[0],&info_);
220 
221 }
222 
223 template<class Real>
224 InnerProductMatrixSolver<Real>::InnerProductMatrixSolver(ROL::Ptr<Teuchos::LAPACK<int,Real> > lapack,
225  const std::vector<Real> &U,
226  const std::vector<Real> &V,
227  const std::vector<Real> &w,
228  const std::vector<Real> &a):
229  InnerProductMatrix<Real>(U,V,w,a),
230  lapack_(lapack),
231  ni_(InnerProductMatrix<Real>::ni_),
232  nq_(InnerProductMatrix<Real>::nq_),
233  M_(InnerProductMatrix<Real>::M_),
234  TRANS_('N'), ipiv_(ni_,0), PLU_(ni_*ni_,0),
235  nrhs_(1),info_(0){
236  PLU_ = M_;
237 
238  // Do matrix factorization
239  lapack->GETRF(ni_,ni_,&PLU_[0],ni_,&ipiv_[0],&info_);
240 
241 }
242 
244 template<class Real>
245 void InnerProductMatrixSolver<Real>::solve(ROL::Ptr<const std::vector<Real> > bp,
246  ROL::Ptr<std::vector<Real> > xp){
247 
248  int nrhs = bp->size()/ni_;
249 
250  *xp = *bp;
251 
252  // Solve LU-factored system
253  lapack_->GETRS(TRANS_,ni_,nrhs,&PLU_[0],ni_,&ipiv_[0],&(*xp)[0],ni_,&info_);
254 
255 }
256 
258 template<class Real>
259 Real InnerProductMatrixSolver<Real>::inv_inner( ROL::Ptr<const std::vector<Real> > up,
260  ROL::Ptr<const std::vector<Real> > vp ) {
261  Real J = 0;
262  ROL::Ptr<std::vector<Real> > Mivp = ROL::makePtr<std::vector<Real>>(ni_,0);
263  this->solve(vp,Mivp);
264  for(int i=0;i<ni_;++i) {
265  //std::cout << (*up)[i] << " " << (*vp)[i] << " " << (*Mivp)[i] << " \n";
266  J += (*up)[i]*(*Mivp)[i];
267  }
268  return J;
269 }
270 
271 #endif
272 
InnerProductMatrixSolver::lapack_
ROL::Ptr< Teuchos::LAPACK< int, Real > > lapack_
Definition: InnerProductMatrix.hpp:171
InnerProductMatrix::inv_inner
virtual Real inv_inner(ROL::Ptr< const std::vector< Real > > up, ROL::Ptr< const std::vector< Real > > vp)
Definition: InnerProductMatrix.hpp:47
InnerProductMatrixSolver::solve
void solve(ROL::Ptr< const std::vector< Real > > bp, ROL::Ptr< std::vector< Real > > xp)
solve for
Definition: InnerProductMatrix.hpp:245
InnerProductMatrix::M_
std::vector< Real > M_
Definition: InnerProductMatrix.hpp:58
InnerProductMatrix::U_
const std::vector< Real > U_
Definition: InnerProductMatrix.hpp:54
InnerProductMatrix::inner
Real inner(ROL::Ptr< const std::vector< Real > > up, ROL::Ptr< const std::vector< Real > > vp)
Compute the inner product .
Definition: InnerProductMatrix.hpp:151
InnerProductMatrixSolver::InnerProductMatrixSolver
InnerProductMatrixSolver(ROL::Ptr< Teuchos::LAPACK< int, Real > > lapack, const std::vector< Real > &U=std::vector< Real >(), const std::vector< Real > &V=std::vector< Real >(), const std::vector< Real > &w=std::vector< Real >(), const int a=1)
Definition: InnerProductMatrix.hpp:204
ROL::Vector
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:80
InnerProductMatrix::apply
void apply(ROL::Ptr< const std::vector< Real > > xp, ROL::Ptr< std::vector< Real > > bp)
Definition: InnerProductMatrix.hpp:102
InnerProductMatrix::w_
const std::vector< Real > w_
Definition: InnerProductMatrix.hpp:56
InnerProductMatrix::applyadd
void applyadd(ROL::Ptr< const std::vector< Real > > xp, ROL::Ptr< std::vector< Real > > bp)
Definition: InnerProductMatrix.hpp:113
InnerProductMatrix::InnerProductMatrix
InnerProductMatrix(const std::vector< Real > &U, const std::vector< Real > &V, const std::vector< Real > &w, const int a=1)
Definition: InnerProductMatrix.hpp:66
ROL_StdVector.hpp
InnerProductMatrixSolver::inv_inner
Real inv_inner(ROL::Ptr< const std::vector< Real > > up, ROL::Ptr< const std::vector< Real > > vp)
Compute the inner product .
Definition: InnerProductMatrix.hpp:259
InnerProductMatrix::V_
const std::vector< Real > V_
Definition: InnerProductMatrix.hpp:55
ROL::solve
virtual void solve(Vector< Real > &c, Vector< Real > &u, const Vector< Real > &z) override
Definition: ROL_Objective_TimeSimOpt.hpp:184
InnerProductMatrixSolver
This class adds a solve method.
Definition: InnerProductMatrix.hpp:168
InnerProductMatrix::applyaddtimes
void applyaddtimes(ROL::Ptr< const std::vector< Real > > xp, ROL::Ptr< std::vector< Real > > bp, Real factor)
Definition: InnerProductMatrix.hpp:123
ROL::apply
ROL::DiagonalOperator apply
InnerProductMatrix::update
void update(const std::vector< Real > &a)
Definition: InnerProductMatrix.hpp:133
InnerProductMatrix::solve
virtual void solve(ROL::Ptr< const std::vector< Real > > bp, ROL::Ptr< std::vector< Real > > xp)
Definition: InnerProductMatrix.hpp:43
Generated by   doxygen 1.8.5