Extension: Matrix functions
#include "dg/matrix/matrix.h"
lanczos.h
Go to the documentation of this file.
1#pragma once
2#include <cusp/dia_matrix.h>
3#include <cusp/coo_matrix.h>
4
5#include "dg/algorithm.h"
6#include "tridiaginv.h"
7#include "matrixfunction.h"
8
9
14namespace dg{
15namespace mat{
16
66template< class ContainerType >
68{
69 public:
71 using HCooMatrix = cusp::coo_matrix<int, value_type, cusp::host_memory>;
72 using HDiaMatrix = cusp::dia_matrix<int, value_type, cusp::host_memory>;
73 using HVec = dg::HVec;
82 UniversalLanczos( const ContainerType& copyable, unsigned max_iterations)
83 {
84 m_v = m_vp = m_vm = copyable;
85 m_max_iter = max_iterations;
86 m_iter = max_iterations;
87 //sub matrix and vector
88 set_iter( max_iterations);
89 }
91 template<class ...Params>
92 void construct( Params&& ...ps)
93 {
94 //construct and swap
95 *this = UniversalLanczos( std::forward<Params>( ps)...);
96 }
97
100 void set_max( unsigned new_max) {
101 m_max_iter = new_max;
102 set_iter( new_max);
103 }
104
107 unsigned get_max() const {return m_max_iter;}
108
111 void set_verbose( bool verbose){ m_verbose = verbose;}
112
115 value_type get_bnorm() const{return m_bnorm;}
116
144 template < class MatrixType, class ContainerType0, class ContainerType1,
145 class ContainerType2, class FuncTe1>
146 unsigned solve(ContainerType0& x, FuncTe1 f,
147 MatrixType&& A, const ContainerType1& b,
148 const ContainerType2& weights, value_type eps,
149 value_type nrmb_correction = 1.,
150 std::string error_norm = "universal",
151 value_type res_fac = 1.,
152 unsigned q = 1 )
153 {
154 tridiag( f, std::forward<MatrixType>(A), b, weights, eps,
155 nrmb_correction, error_norm, res_fac, q);
156 if( "residual" == error_norm)
157 m_yH = f( m_TH);
158 //Compute x = |b|_M V f(T) e1
159 normMbVy(std::forward<MatrixType>(A), m_TH, m_yH, x, b,
160 m_bnorm);
161 return m_iter;
162 }
183 template< class MatrixType, class ContainerType0, class ContainerType1>
184 const HDiaMatrix& tridiag( MatrixType&& A, const ContainerType0& b,
185 const ContainerType1& weights, value_type eps = 1e-12,
186 value_type nrmb_correction = 1.,
187 std::string error_norm = "universal", value_type res_fac = 1.,
188 unsigned q = 1
189 )
190 {
191 auto op = make_Linear_Te1( -1);
192 tridiag( op, std::forward<MatrixType>(A), b, weights, eps,
193 nrmb_correction, error_norm, res_fac, q);
194 return m_TH;
195 }
196
200 unsigned get_iter() const {return m_iter;}
201 private:
202
215 template< class MatrixType, class DiaMatrixType, class ContainerType0,
216 class ContainerType1,class ContainerType2>
217 void normMbVy( MatrixType&& A,
218 const DiaMatrixType& T,
219 const ContainerType0& y,
220 ContainerType1& x,
221 const ContainerType2& b, value_type bnorm)
222 {
223 dg::blas1::copy(0., x);
224 if( 0 == bnorm )
225 {
226 return;
227 }
228 dg::blas1::axpby(1./bnorm, b, 0.0, m_v); //m_v[1] = b/||b||
229 dg::blas1::copy(0., m_vm);
230 // check if (potentially) all higher elements in y are zero
231 unsigned less_iter = 0;
232 for( unsigned i=0; i<y.size(); i++)
233 if( y[i] != 0)
234 less_iter = i+1;
235 dg::blas1::axpby( y[0]*bnorm, m_v, 1., x); //Compute b= |b| V y
236 for ( unsigned i=0; i<less_iter-1; i++)
237 {
238 dg::blas2::symv( std::forward<MatrixType>(A), m_v, m_vp);
240 -T.values(i,0)/T.values(i,2), m_vm,
241 -T.values(i,1)/T.values(i,2), m_v,
242 1.0/T.values(i,2), m_vp);
243 dg::blas1::axpby( y[i+1]*bnorm, m_vp, 1., x); //Compute b= |b| V y
244 m_vm.swap( m_v);
245 m_v.swap( m_vp);
246
247 }
248 }
249 template < class MatrixType, class ContainerType1,
250 class ContainerType2, class UnaryOp>
251 void tridiag(UnaryOp f,
252 MatrixType&& A, const ContainerType1& b,
253 const ContainerType2& weights, value_type eps,
254 value_type nrmb_correction,
255 std::string error_norm = "residual",
256 value_type res_fac = 1.,
257 unsigned q = 1 )
258 {
259#ifdef MPI_VERSION
260 int rank;
261 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
262#endif //MPI
263 m_bnorm = sqrt(dg::blas2::dot(b, weights, b));
264 if( m_verbose)
265 {
266 DG_RANK0 std::cout << "# Norm of b "<<m_bnorm <<"\n";
267 DG_RANK0 std::cout << "# Res factor "<<res_fac <<"\n";
268 DG_RANK0 std::cout << "# Residual errors: \n";
269 }
270 if( m_bnorm == 0)
271 {
272 set_iter(1);
273 return;
274 }
275 value_type residual;
276 dg::blas1::axpby(1./m_bnorm, b, 0.0, m_v); //m_v[1] = x/||x||
277 value_type betaip = 0.;
278 value_type alphai = 0.;
279 for( unsigned i=0; i<m_max_iter; i++)
280 {
281 m_TH.values(i,0) = betaip; // -1 diagonal
282 dg::blas2::symv(std::forward<MatrixType>(A), m_v, m_vp);
283 dg::blas1::axpby(-betaip, m_vm, 1.0, m_vp); // only - if i>0, therefore no if (i>0)
284 alphai = dg::blas2::dot(m_vp, weights, m_v);
285 m_TH.values(i,1) = alphai;
286 dg::blas1::axpby(-alphai, m_v, 1.0, m_vp);
287 betaip = sqrt(dg::blas2::dot(m_vp, weights, m_vp));
288 if (betaip == 0)
289 {
290 if( m_verbose)
291 DG_RANK0 std::cout << "beta["<<i+1 <<"]=0 encountered\n";
292 set_iter(i+1);
293 break;
294 }
295 m_TH.values(i,2) = betaip; // +1 diagonal
296
297 value_type xnorm = 0.;
298 if( "residual" == error_norm)
299 {
300 residual = compute_residual_error( m_TH, i)*m_bnorm;
301 xnorm = m_bnorm;
302 }
303 else
304 {
305 if( i>=q &&( (i<=10) || (i>10 && i%10 == 0) ))
306 {
307 residual = compute_universal_error( m_TH, i, q, f,
308 m_yH)*m_bnorm;
309 xnorm = dg::fast_l2norm( m_yH)*m_bnorm;
310 }
311 else
312 {
313 residual = 1e10;
314 xnorm = m_bnorm;
315 }
316 }
317 if( m_verbose)
318 DG_RANK0 std::cout << "# ||r||_W = " << residual << "\tat i = " << i << "\n";
319 if (res_fac*residual< eps*(xnorm + nrmb_correction) )
320 {
321 set_iter(i+1);
322 break;
323 }
324 dg::blas1::scal(m_vp, 1./betaip);
325 m_vm.swap(m_v);
326 m_v.swap( m_vp);
327 set_iter( m_max_iter);
328 }
329 }
330 value_type compute_residual_error( const HDiaMatrix& TH, unsigned iter)
331 {
332 value_type T1 = compute_Tinv_m1( TH, iter+1);
333 return TH.values(iter,2)*fabs(T1); //Tinv_i1
334 }
335 template<class UnaryOp>
336 value_type compute_universal_error( const HDiaMatrix& TH, unsigned iter,
337 unsigned q, UnaryOp f, HVec& yH)
338 {
339 unsigned new_iter = iter + 1 + q;
340 set_iter( iter+1);
341 HDiaMatrix THtilde( new_iter, new_iter, 3*new_iter-2, 3);
342 THtilde.diagonal_offsets[0] = -1;
343 THtilde.diagonal_offsets[1] = 0;
344 THtilde.diagonal_offsets[2] = 1;
345 for( unsigned u=0; u<iter+1; u++)
346 {
347 THtilde.values(u,0) = TH.values(u,0);
348 THtilde.values(u,1) = TH.values(u,1);
349 THtilde.values(u,2) = TH.values(u,2);
350 }
351 for( unsigned u=1; u<=q; u++)
352 {
353 THtilde.values( iter+u, 0) = u==1 ? TH.values(iter,2) :
354 TH.values( iter+1-u, 1);
355 THtilde.values( iter+u, 1) = TH.values( iter-u, 1);
356 THtilde.values( iter+u, 2) = TH.values( iter-u, 0);
357 }
358 yH = f( TH);
359 HVec yHtilde = f( THtilde);
360 for( unsigned u=0; u<yH.size(); u++)
361 yHtilde[u] -= yH[u];
362 value_type norm = dg::fast_l2norm( yHtilde);
363 return norm;
364 }
365
368 void set_iter( unsigned new_iter) {
369 // The alignment (which is the pitch of the underlying values)
370 // of m_max_iter preserves the existing elements
371 m_TH.resize(new_iter, new_iter, 3*new_iter-2, 3, m_max_iter);
372 m_TH.diagonal_offsets[0] = -1;
373 m_TH.diagonal_offsets[1] = 0;
374 m_TH.diagonal_offsets[2] = 1;
375 m_iter = new_iter;
376 }
377 ContainerType m_v, m_vp, m_vm;
378 HDiaMatrix m_TH;
379 HVec m_yH;
380 unsigned m_iter, m_max_iter;
381 bool m_verbose = false;
382 value_type m_bnorm = 0.;
383};
384
385
386} //namespace mat
387} //namespace dg
388
Tridiagonalize and approximate via Lanczos algorithm. A is self-adjoint in the weights .
Definition: lanczos.h:68
cusp::dia_matrix< int, value_type, cusp::host_memory > HDiaMatrix
Definition: lanczos.h:72
value_type get_bnorm() const
Norm of b from last call to operator()
Definition: lanczos.h:115
get_value_type< ContainerType > value_type
value type of the ContainerType class
Definition: lanczos.h:70
cusp::coo_matrix< int, value_type, cusp::host_memory > HCooMatrix
Definition: lanczos.h:71
UniversalLanczos(const ContainerType &copyable, unsigned max_iterations)
Allocate memory for the method.
Definition: lanczos.h:82
void set_verbose(bool verbose)
Set or unset debugging output during iterations.
Definition: lanczos.h:111
dg::HVec HVec
Definition: lanczos.h:73
unsigned get_max() const
Get the current maximum number of iterations.
Definition: lanczos.h:107
unsigned get_iter() const
Get the number of iterations in the last call to tridiag or solve (same as T.num_rows)
Definition: lanczos.h:200
UniversalLanczos()
Allocate nothing, Call construct method before usage.
Definition: lanczos.h:75
unsigned solve(ContainerType0 &x, FuncTe1 f, MatrixType &&A, const ContainerType1 &b, const ContainerType2 &weights, value_type eps, value_type nrmb_correction=1., std::string error_norm="universal", value_type res_fac=1., unsigned q=1)
via Lanczos and matrix function computation. A is self-adjoint in the weights .
Definition: lanczos.h:146
void set_max(unsigned new_max)
Set the maximum number of iterations.
Definition: lanczos.h:100
const HDiaMatrix & tridiag(MatrixType &&A, const ContainerType0 &b, const ContainerType1 &weights, value_type eps=1e-12, value_type nrmb_correction=1., std::string error_norm="universal", value_type res_fac=1., unsigned q=1)
Tridiagonalization of A using Lanczos method.
Definition: lanczos.h:184
void construct(Params &&...ps)
Perfect forward parameters to one of the constructors.
Definition: lanczos.h:92
void copy(const ContainerTypeIn &source, ContainerTypeOut &target)
void axpbypgz(get_value_type< ContainerType > alpha, const ContainerType1 &x, get_value_type< ContainerType > beta, const ContainerType2 &y, get_value_type< ContainerType > gamma, ContainerType &z)
void axpby(get_value_type< ContainerType > alpha, const ContainerType1 &x, get_value_type< ContainerType > beta, ContainerType &y)
void scal(ContainerType &x, get_value_type< ContainerType > alpha)
void symv(MatrixType &&M, const ContainerType1 &x, ContainerType2 &y)
get_value_type< MatrixType > dot(const ContainerType1 &x, const MatrixType &m, const ContainerType2 &y)
MPI_Vector< thrust::host_vector< real_type > > weights(const aRealMPITopology2d< real_type > &g)
value_type compute_Tinv_m1(const cusp::dia_matrix< int, value_type, cusp::host_memory > &T, unsigned size)
Computes the value of via a Thomas algorithm.
Definition: tridiaginv.h:32
static auto fast_l2norm
typename TensorTraits< std::decay_t< Vector > >::value_type get_value_type
static auto make_Linear_Te1(int exp)
Create a functor that computes directly.
Definition: matrixfunction.h:190
thrust::host_vector< double > HVec
Classes for Krylov space approximations of a Matrix-Vector product.
#define DG_RANK0