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
19template<class Real>
20class 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
74template<class Real>
75InnerProductMatrix<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
90template<class Real>
91InnerProductMatrix<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
106template<class Real>
109
110template<class Real>
111void 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
121template<class Real>
122void 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
131template<class Real>
132void 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
141template<class Real>
142void 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
159template<class Real>
160Real 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
176template<class Real>
177class 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
212template<class Real>
213InnerProductMatrixSolver<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
232template<class Real>
233InnerProductMatrixSolver<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
253template<class Real>
254void 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
267template<class Real>
268Real 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
V
Vector< Real > V
Definition ROL_BlockOperator2Diagonal.hpp:46
ROL_StdVector.hpp
InnerProductMatrixSolver
This class adds a solve method.
Definition InnerProductMatrix.hpp:177
InnerProductMatrixSolver::solve
void solve(ROL::Ptr< const std::vector< Real > > bp, ROL::Ptr< std::vector< Real > > xp)
solve for
Definition InnerProductMatrix.hpp:254
InnerProductMatrixSolver::lapack_
ROL::Ptr< Teuchos::LAPACK< int, Real > > lapack_
Definition InnerProductMatrix.hpp:180
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
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::applyaddtimes
void applyaddtimes(ROL::Ptr< const std::vector< Real > > xp, ROL::Ptr< std::vector< Real > > bp, Real factor)
Definition InnerProductMatrix.hpp:132
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
InnerProductMatrix::InnerProductMatrix
InnerProductMatrix(InnerProductMatrix< Real > *ipm)
InnerProductMatrix::U_
const std::vector< Real > U_
Definition InnerProductMatrix.hpp:63
InnerProductMatrix::w_
const std::vector< Real > w_
Definition InnerProductMatrix.hpp:65
InnerProductMatrix::apply
void apply(ROL::Ptr< const std::vector< Real > > xp, ROL::Ptr< std::vector< Real > > bp)
Definition InnerProductMatrix.hpp:111
InnerProductMatrix::solve
virtual void solve(ROL::Ptr< const std::vector< Real > > bp, ROL::Ptr< std::vector< Real > > xp)
Definition InnerProductMatrix.hpp:52
InnerProductMatrix::M_
std::vector< Real > M_
Definition InnerProductMatrix.hpp:67
InnerProductMatrix::update
void update(const std::vector< Real > &a)
Definition InnerProductMatrix.hpp:142
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
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
InnerProductMatrix::V_
const std::vector< Real > V_
Definition InnerProductMatrix.hpp:64
InnerProductMatrix::applyadd
void applyadd(ROL::Ptr< const std::vector< Real > > xp, ROL::Ptr< std::vector< Real > > bp)
Definition InnerProductMatrix.hpp:122
Generated by doxygen 1.9.8