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
10template<class Real>
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
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
65template<class Real>
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
81template<class Real>
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
97template<class Real>
99}
100
101template<class Real>
102void 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
112template<class Real>
113void 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
122template<class Real>
123void 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
132template<class Real>
133void 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
150template<class Real>
151Real 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
167template<class 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
203template<class Real>
204InnerProductMatrixSolver<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
223template<class Real>
224InnerProductMatrixSolver<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
244template<class Real>
245void 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
258template<class Real>
259Real 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
Vector< Real > V
This class adds a solve method.
void solve(ROL::Ptr< const std::vector< Real > > bp, ROL::Ptr< std::vector< Real > > xp)
solve for
ROL::Ptr< Teuchos::LAPACK< int, Real > > lapack_
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)
Real inv_inner(ROL::Ptr< const std::vector< Real > > up, ROL::Ptr< const std::vector< Real > > vp)
Compute the inner product .
void applyaddtimes(ROL::Ptr< const std::vector< Real > > xp, ROL::Ptr< std::vector< Real > > bp, Real factor)
virtual Real inv_inner(ROL::Ptr< const std::vector< Real > > up, ROL::Ptr< const std::vector< Real > > vp)
InnerProductMatrix(InnerProductMatrix< Real > *ipm)
const std::vector< Real > U_
const std::vector< Real > w_
void apply(ROL::Ptr< const std::vector< Real > > xp, ROL::Ptr< std::vector< Real > > bp)
virtual void solve(ROL::Ptr< const std::vector< Real > > bp, ROL::Ptr< std::vector< Real > > xp)
std::vector< Real > M_
void update(const std::vector< Real > &a)
Real inner(ROL::Ptr< const std::vector< Real > > up, ROL::Ptr< const std::vector< Real > > vp)
Compute the inner product .
InnerProductMatrix(const std::vector< Real > &U, const std::vector< Real > &V, const std::vector< Real > &w, const int a=1)
const std::vector< Real > V_
void applyadd(ROL::Ptr< const std::vector< Real > > xp, ROL::Ptr< std::vector< Real > > bp)