ROL
ROL_Lanczos.hpp
Go to the documentation of this file.
1// @HEADER
2// ************************************************************************
3//
4// Rapid Optimization Library (ROL) Package
5// Copyright (2014) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact lead developers:
38// Drew Kouri (dpkouri@sandia.gov) and
39// Denis Ridzal (dridzal@sandia.gov)
40//
41// ************************************************************************
42// @HEADER
43
44#ifndef ROL_LANCZOS_H
45#define ROL_LANCZOS_H
46
47#include "ROL_Krylov.hpp"
49#include "ROL_Vector.hpp"
50#include "ROL_Types.hpp"
51#include "ROL_LAPACK.hpp"
52
53namespace ROL {
54
63template<class Real>
64class Lanczos {
65
66 template <typename T> using ROL::Ptr = ROL::Ptr<T>;
67 template <typename T> using vector = std::vector<T>;
68
69 template typename vector<Real> size_type uint;
70
71 typedef Vector<Real> V;
73
74 typedef ROL::ParameterList PL;
75
76private:
77
78 ROL::LAPACK<int,Real> lapack_;
79
80 vector<ROL::Ptr<V> > Q_; // Orthogonal basis
81 vector<Real> alpha_; // Diagonal recursion coefficients
82 vector<Real> beta_; // Sub/super-diagonal recursion coefficients
83
84 // Temporary vectors for factorizations, linear solves, and eigenvalue calculations
89 vector<Real> y_; // Arnoldi expansion coefficients
90
91 vector<Real> work_; // Scratch space for eigenvalue decomposition
92 vector<int> ipiv_; // Pivots for LU
93
94 ROL::Ptr<V> u_; // An auxilliary vector
95
96 Real max_beta_; // maximum beta encountered
97 Real tol_beta_; // relative smallest beta allowed
98
99 Real tol_ortho_; // Maximum orthogonality loss tolerance
100
101 int maxit_; // Maximum number of vectors to store
102
103 int k_; // current iterate number
104
105
106 // Allocte memory for Arnoldi vectors and recurions coefficients
107 void allocate( void ) {
108
109 u_ = b.clone();
110
111 alpha_.reserve(maxit_);
112 beta_.reserve(maxit_);
113
114 dl_.reserve(maxit_);
115 d_.reserve(maxit_);
116 du_.reserve(maxit_);
117 du2_.reserve(maxit_);
118
119 work_.reserve(4*maxit_);
120
121 ipiv_.reserve(maxit_);
122
123 y_.reserve(maxit_);
124
125 alpha_.reserve(maxit_);
126 beta_.reserve(maxit_);
127
128 for( uint j=0; j<maxit_; ++j ) {
129 Q_.push_back(b.clone());
130 }
131 }
132
133public:
134
135 enum class FLAG_ITERATE : unsigned {
136 ITERATE_SUCCESS = 0,
137 ITERATE_SMALL_BETA, // Beta too small to continue
138 ITERATE_MAX_REACHED, // Reached maximum number of iterations
139 ITERATE_ORTHO_TOL, // Reached maximim orthogonality loss
141 };
142
143 enum class FLAG_SOLVE : unsigned {
144 SOLVE_SUCCESS = 0,
148 };
149
150
151 Lanczos( ROL::ParameterList &PL ) {
152 PL &krylovList = parlist.sublist("General").sublist("Krylov");
153 PL &lanczosList = krylovList.sublist("Lanczos");
154
155 Real tol_default = std::sqrt(ROL_EPSILON<Real>());
156
157 maxit_ = krylovList_.get("Iteration Limit",10);
158 tol_beta_ = lanczosList.get("Beta Relative Tolerance", tol_default);
159 tol_ortho_ = lanczosList.get("Orthogonality Tolerance", tol_default);
160
161 }
162
163 void initialize( const V& b ) {
164 allocate();
165 reset(b);
166
167 }
168
169 void initialize( const V &x0, const V &b, const LO &A, Real &tol ) {
170 allocate();
171 reset(x0,b,A,tol);
172
173 }
174
175
176 void reset( const V &b ) {
177 k_ = 0;
178 max_beta_ = 0;
179 Q_[0]->set(b);
180 beta_[0] = Q_[0]->norm();
181 max_beta_ = std::max(max_beta_,beta_[0]);
182 Q_[0]->scale(1.0/beta_[0]);
183 }
184
185 void reset( const V &x0, const V &b, const LO &A, Real &tol ) {
186 k_ = 0;
187 max_beta_ = 0;
188 Q_[0]->set(b);
189 A.apply(*u_,x0,tol);
190 Q_[0]->axpy(-1.0,*u_);
191 beta_[0] = Q_[0]->norm();
192 max_beta_ = std::max(max_beta_,beta_[0]);
193 Q_[0]->scale(1.0/beta_[0]);
194 }
195
196 FLAG_ITERATE iterate( const OP &A, Real &tol ) {
197
198 if( k_ == maxit_ ) {
199 return ITERATE_MAX_REACHED;
200 }
201
202 A.apply(*u_,*(Q_[k]),tol);
203 Real delta;
204
205 if( k_>0 ) {
206 u_->axpy(-beta_[k],V_[k_-1]);
207 }
208 alpha_[k] = u_->dot(*(V_[k]));
209 u_->axpy(alpha_[k],V_[k_]);
210
211 if( k_>0 ) {
212 delta = u_->dot(*(V_[k-1]));
213 u_->axpy(-delta,*(V_[k-1]));
214 }
215 delta = u_->dot(*(V_[k]));
216 alpha_[k] += delta;
217
218 if( k_ < maxit_-1 ) {
219 u_->axpy(-delta,*(V_[k_]));
220 beta_[k+1] = u_->norm();
221 max_beta_ = std::max(max_beta_,beta_[k+1]);
222 if(beta_[k+1] < tol_beta_*max_beta_) {
223 return ITERATE_SMALL_BETA;
224 }
225
226 V_[k+1]->set(*u_);
227 V_[k+1]->scale(1.0/beta_[k+1]);
228
229 // Check orthogonality
230 Real dotprod = V_[k+1]->dot(*(V_[0]));
231
232 if( std::sqrt(dotprod) > tol_ortho_ ) {
233 return ITERATE_ORTHO_TOL;
234 }
235 }
236
237 ++k_;
238 return ITERATE_SUCCESS;
239 }
240
241 // Compute the eigenvalues of the tridiagonal matrix T
242 void eigenvalues( std::vector<Real> &E ) {
243
244 std::vector<Real> Z(1,0); // Don't compute eigenvectors
245
246 int INFO;
247 int LDZ = 0;
248 const char COMPZ = 'N':
249 d_ = alpha_;
250 du_ = beta_;
251
252 lapack_->STEQR(COMPZ,k_,&d_[0],&du_[0],&Z[0],LDZ,&work_[0],&INFO);
253
254 if( INFO < 0 ) {
255 return SOLVE_ILLEGAL_VALUE;
256 }
257 else if( INFO > 0 ) {
258 return SOLVE_SINGULAR_U;
259 }
260
261 }
262
263 FLAG_SOLVE solve( V &x, Real tau=0 ) {
264
265 const char TRANS = 'N';
266 const int NRHS = 1;
267 int INFO;
268
269 // Fill arrays
270 for(uint j=0;j<k_;++j) {
271 d_[j] = alpha_[j]+tau;
272 }
273
274 dl_ = beta_;
275 du_ = beta_;
276 du2_.assign(maxit_,0);
277
278 // Do Tridiagonal LU factorization
279 lapack_->GTTRF(k_,&dl_[0],&d_[0],&du_[0],&du2_[0],&ipiv_[0],&INFO);
280
281 // Solve the factorized system for Arnoldi expansion coefficients
282 lapack_->GTTRS(TRANS,k_,1,&dl[0],&d_[0],&du_[0],&du2_[0],&ipiv_[0],&y_[0],k_,&INFO);
283
284 }
285
286
287
288}; // class LanczosFactorization
289} // namespace ROL
290
291#endif // ROL_LANCZOS_H
typename PV< Real >::size_type size_type
Contains definitions of custom data types in ROL.
Interface for computing the Lanczos vectors and approximate solutions to symmetric indefinite linear ...
vector< Real > alpha_
ROL::Ptr< V > u_
vector< Real > du_
std::vector< T > vector
vector< int > ipiv_
vector< Real > dl_
void allocate(void)
template vector< Real > size_type uint
FLAG_SOLVE solve(V &x, Real tau=0)
vector< Real > d_
vector< Real > du2_
vector< Real > beta_
vector< Real > y_
vector< ROL::Ptr< V > > Q_
void reset(const V &x0, const V &b, const LO &A, Real &tol)
void eigenvalues(std::vector< Real > &E)
ROL::LAPACK< int, Real > lapack_
void initialize(const V &b)
ROL::ParameterList PL
vector< Real > work_
void reset(const V &b)
void initialize(const V &x0, const V &b, const LO &A, Real &tol)
FLAG_ITERATE iterate(const OP &A, Real &tol)
Lanczos(ROL::ParameterList &PL)
LinearOperator< Real > OP
Vector< Real > V
Provides the interface to apply a linear operator.
virtual void apply(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const =0
Apply linear operator.
Defines the linear algebra or vector space interface.