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