ROL
ROL_StdLinearOperator.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 #ifndef ROL_STDLINEAROPERATOR_H
11 #define ROL_STDLINEAROPERATOR_H
12 
13 #include "ROL_LinearOperator.hpp"
14 #include "ROL_StdVector.hpp"
15 #include "ROL_LAPACK.hpp"
16 
17 
28 namespace ROL {
29 
30 template <class Real>
31 class StdLinearOperator : public LinearOperator<Real> {
32 
34 
35  typedef std::vector<Real> vector;
36 
37 private:
38 
39  ROL::Ptr<std::vector<Real> > A_;
40  int N_;
41  int INFO_;
42 
43  mutable vector PLU_;
44  mutable std::vector<int> ipiv_;
45 
46  ROL::LAPACK<int,Real> lapack_;
47 
48 public:
49 
51 
52  StdLinearOperator( ROL::Ptr<std::vector<Real> > &A ) : A_(A) {
53  int N2 = A_->size();
54  N_ = (std::round(std::sqrt(N2)));
55  bool isSquare = N_*N_ == N2;
56  ROL_TEST_FOR_EXCEPTION( !isSquare, std::invalid_argument,
57  "Error: vector representation of matrix must have a square "
58  "number of elements.");
59  ipiv_.resize(N_);
60  }
61 
62  virtual ~StdLinearOperator() {}
63 
65  void update( const Vector<Real> &x, bool flag = true, int iter = -1 ) {
66  ROL::Ptr<const vector> xp = dynamic_cast<const SV&>(x).getVector();
67  update(*xp,flag,iter);
68  }
69 
70  virtual void update( const std::vector<Real> &x, bool flag = true, int iter = -1 ) {}
71 
72  // Matrix multiplication
74  void apply( Vector<Real> &Hv, const Vector<Real> &v, Real &tol ) const {
75 
76  ROL::Ptr<vector> Hvp = dynamic_cast<SV&>(Hv).getVector();
77  ROL::Ptr<const vector> vp = dynamic_cast<const SV&>(v).getVector();
78  apply(*Hvp,*vp,tol);
79  }
80 
81  virtual void apply( std::vector<Real> &Hv, const std::vector<Real> &v, Real &tol ) const {
82  for( int i=0; i<N_; ++i ) {
83  Hv[i] = Real(0);
84  for( int j=0; j<N_; ++j ) {
85  Hv.at(i) += A_->at(N_*j+i)*v.at(j);
86  }
87  }
88  }
89 
90  // Matrix multiplication with transpose
92  void applyAdjoint( Vector<Real> &Hv, const Vector<Real> &v, Real &tol ) const {
93 
94  ROL::Ptr<vector> Hvp = dynamic_cast<SV&>(Hv).getVector();
95  ROL::Ptr<const vector> vp = dynamic_cast<const SV&>(v).getVector();
96  applyAdjoint(*Hvp,*vp,tol);
97  }
98 
99  virtual void applyAdjoint( std::vector<Real> &Hv, const std::vector<Real> &v, Real &tol ) const {
100  for( int i=0; i<N_; ++i ) {
101  Hv[i] = Real(0);
102  for( int j=0; j<N_; ++j ) {
103  Hv.at(i) += A_->at(N_*i+j)*v.at(j);
104  }
105  }
106  }
107 
108 
109  // Solve the system
110 
112  void applyInverse( Vector<Real> &Hv, const Vector<Real> &v, Real &tol ) const {
113 
114  ROL::Ptr<vector> Hvp = dynamic_cast<SV&>(Hv).getVector();
115  ROL::Ptr<const vector> vp = dynamic_cast<const SV&>(v).getVector();
116  applyInverse(*Hvp,*vp,tol);
117  }
118 
119  virtual void applyInverse( std::vector<Real> &Hv, const std::vector<Real> &v, Real &tol ) const {
120 
121  const int LDA = N_;
122  const int LDB = N_;
123  int INFO;
124  int NRHS = 1;
125 
126  Hv = v;
127  PLU_ = *A_;
128 
129  // Do LU factorization
130  lapack_.GETRF(N_,N_,&PLU_[0],LDA,&ipiv_[0],&INFO);
131 
132  ROL_TEST_FOR_EXCEPTION(INFO>0,std::logic_error,"Error in StdLinearOperator::applyInverse(): "
133  "Zero diagonal element encountered in matrix factor U(" << INFO << "," << INFO << ").");
134 
135  ROL_TEST_FOR_EXCEPTION(INFO<0,std::logic_error,"Error in StdLinearOperator::applyInverse(): "
136  "Illegal value encountered in element " << -INFO << " when performing LU factorization.");
137 
138  // Solve factored system
139  lapack_.GETRS('N',N_,NRHS,&PLU_[0],LDA,&ipiv_[0],&Hv[0],LDB,&INFO);
140 
141  ROL_TEST_FOR_EXCEPTION(INFO<0,std::logic_error,"Error in StdLinearOperator::applyInverse(): "
142  "Illegal value encountered in element " << -INFO << " when solving the factorized system. ");
143 
144  }
145 
146  // Solve the system with transposed matrix
147 
149  void applyAdjointInverse( Vector<Real> &Hv, const Vector<Real> &v, Real &tol ) const {
150 
151  ROL::Ptr<vector> Hvp = dynamic_cast<SV&>(Hv).getVector();
152  ROL::Ptr<const vector> vp = dynamic_cast<const SV&>(v).getVector();
153  applyAdjointInverse(*Hvp,*vp,tol);
154  }
155 
156  virtual void applyAdjointInverse( std::vector<Real> &Hv, const std::vector<Real> &v, Real &tol ) const {
157 
158  const int LDA = N_;
159  const int LDB = N_;
160  int INFO;
161  int NRHS = 1;
162 
163  Hv = v;
164  PLU_ = *A_;
165 
166  // Do LU factorization
167  lapack_.GETRF(N_,N_,&PLU_[0],LDA,&ipiv_[0],&INFO);
168 
169  ROL_TEST_FOR_EXCEPTION(INFO>0,std::logic_error,"Error in StdLinearOperator::applyAdjointInverse(): "
170  "Zero diagonal element encountered in matrix factor U(" << INFO << "," << INFO << ").");
171 
172  ROL_TEST_FOR_EXCEPTION(INFO<0,std::logic_error,"Error in StdLinearOperator::applyAdjointInverse(): "
173  "Illegal value encountered in element " << -INFO << " when performing LU factorization.");
174 
175  // Solve factored system
176  lapack_.GETRS('T',N_,NRHS,&PLU_[0],LDA,&ipiv_[0],&Hv[0],LDB,&INFO);
177 
178  ROL_TEST_FOR_EXCEPTION(INFO<0,std::logic_error,"Error in StdLinearOperator::applyAdjointInverse(): "
179  "Illegal value encountered in element " << -INFO << " when solving the factorized system. ");
180 
181  }
182 
183 }; // class LinearOperator
184 
185 } // namespace ROL
186 
187 #endif
void applyInverse(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const
Apply inverse of linear operator.
void apply(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const
Apply linear operator.
virtual void apply(std::vector< Real > &Hv, const std::vector< Real > &v, Real &tol) const
Provides the std::vector implementation to apply a linear operator, which is a std::vector representa...
virtual void applyAdjointInverse(std::vector< Real > &Hv, const std::vector< Real > &v, Real &tol) const
virtual void applyAdjoint(std::vector< Real > &Hv, const std::vector< Real > &v, Real &tol) const
ROL::LAPACK< int, Real > lapack_
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:46
virtual void applyInverse(std::vector< Real > &Hv, const std::vector< Real > &v, Real &tol) const
void update(const Vector< Real > &x, bool flag=true, int iter=-1)
Update linear operator.
Provides the interface to apply a linear operator.
void applyAdjointInverse(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const
Apply adjoint of the inverse linear operator.
virtual void update(const std::vector< Real > &x, bool flag=true, int iter=-1)
void applyAdjoint(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const
Apply adjoint of linear operator.
StdLinearOperator(ROL::Ptr< std::vector< Real > > &A)
ROL::Ptr< std::vector< Real > > A_