5#ifndef DUNE_ISTL_EIGENVALUE_TEST_MATRIXINFO_HH
6#define DUNE_ISTL_EIGENVALUE_TEST_MATRIXINFO_HH
23#include "../arpackpp.hh"
24#include "../poweriteration.hh"
41template <
typename BCRSMatrix>
46 typedef typename BCRSMatrix::field_type
Real;
64 const bool verbose =
false,
65 const unsigned int arppp_a_verbosity_level = 0,
66 const unsigned int pia_verbosity_level = 0)
69 arppp_a_verbosity_level_(arppp_a_verbosity_level*verbose),
70 pia_verbosity_level_(pia_verbosity_level*verbose),
71 cond_2_(-1.0), symmetricity_assumed_(false)
75 (Dune::blockLevel<BCRSMatrix>() == 2,
76 "Only BCRSMatrices with blocklevel 2 are supported.");
80 (BCRSMatrix::block_type::rows == BCRSMatrix::block_type::cols,
81 "Only BCRSMatrices with square blocks are supported.");
84 const int nrows = m_.M() * BCRSMatrix::block_type::rows;
85 const int ncols = m_.N() * BCRSMatrix::block_type::cols;
88 << nrows <<
"x" << ncols <<
").");
94 if (cond_2_ == -1.0 || symmetricity_assumed_ != assume_symmetric)
97 std::cout <<
" MatrixInfo: Computing 2-norm condition number"
98 << (assume_symmetric ?
" (assuming that matrix is symmetric)." :
".")
101 if (assume_symmetric)
106 symmetricity_assumed_ = assume_symmetric;
127 Real lambda_max{}, lambda_min{}, cond_2{};
136 const ARPPP_A arppp_a(m_,100000,arppp_a_verbosity_level_);
141 PIA pia(m_,20000,pia_verbosity_level_);
142 static const bool avoidLinSolverCrime =
true;
148 const unsigned int piaLS_verbosity = 0;
149 PIALS piaLS(pia.getIterationMatrix(),piaLS_verbosity);
154 typename PIA::IterationOperator::domain_type,
155 typename PIA::IterationOperator::range_type> PIAPC;
156 PIAPC piaPC(pia.getIterationMatrix(),2,1.0);
157 const double piaLS_reduction = 1e-02;
158 const unsigned int piaLS_max_iter = 1000;
159 const unsigned int piaLS_verbosity = 0;
161 PIALS piaLS(pia.getIterationOperator(),piaPC,
162 piaLS_reduction,piaLS_max_iter,piaLS_verbosity);
169 const Real epsilon = 0.0;
171 arppp_a.computeSymMaxMagnitude(epsilon,x,lambda_max);
177 const Real epsilonPI = 1e-02;
178 const Real epsilonTLIME = 1e-08;
182 pia.applyPowerIteration(epsilonPI,x,lambda_max);
184 const Real gamma = m_.infinity_norm();
185 const Real eta = 0.0;
186 const Real delta = 1e-03;
188 pia.template applyTLIMEIteration<PIALS,avoidLinSolverCrime>
189 (gamma,eta,epsilonTLIME,piaLS,delta,2,external,x,lambda_max);
197 const Real epsilon = std::sqrt(std::numeric_limits<Real>::epsilon());
201 const Real gamma = 0.0;
202 const Real eta = 0.0;
203 const Real delta = 1e-03;
205 pia.template applyTLIMEIteration<PIALS,avoidLinSolverCrime>
206 (gamma,eta,epsilon,piaLS,delta,2,external,x,lambda_min);
213 if (std::abs(lambda_max) > m_.infinity_norm())
215 <<
"largest magnitude eigenvalue is greater than "
216 <<
"infinity norm of the matrix!");
220 std::cout <<
" Largest magnitude eigenvalue λ_max = "
221 << lambda_max << std::endl;
225 std::cout <<
" Smallest magnitude eigenvalue λ_min = "
226 << lambda_min << std::endl;
230 cond_2 = std::abs(lambda_max / lambda_min);
234 std::cout <<
" 2-norm condition number cond_2 = "
235 << cond_2 << std::endl;
247 Real sigma_max, sigma_min, cond_2;
256 const ARPPP_A arppp_a(m_,100000,arppp_a_verbosity_level_);
265 Real lambda_max{}, lambda_min{};
270 PIA pia(mtm,20000,pia_verbosity_level_);
271 static const bool avoidLinSolverCrime =
true;
277 const unsigned int piaLS_verbosity = 0;
278 PIALS piaLS(pia.getIterationMatrix(),piaLS_verbosity);
283 typename PIA::IterationOperator::domain_type,
284 typename PIA::IterationOperator::range_type> PIAPC;
285 PIAPC piaPC(pia.getIterationMatrix(),2,1.0);
286 const double piaLS_reduction = 1e-02;
287 const unsigned int piaLS_max_iter = 1000;
288 const unsigned int piaLS_verbosity = 0;
290 PIALS piaLS(pia.getIterationOperator(),piaPC,
291 piaLS_reduction,piaLS_max_iter,piaLS_verbosity);
297 const Real epsilon = 0.0;
299 arppp_a.computeNonSymMax(epsilon,x,sigma_max);
306 const Real epsilonPI = 1e-02;
307 const Real epsilonTLIME = 1e-08;
311 pia.applyPowerIteration(epsilonPI,x,lambda_max);
313 const Real gamma = mtm.infinity_norm();
314 const Real eta = 0.0;
315 const Real delta = 1e-03;
317 pia.template applyTLIMEIteration<PIALS,avoidLinSolverCrime>
318 (gamma,eta,epsilonTLIME,piaLS,delta,2,external,x,lambda_max);
321 sigma_max = std::sqrt(lambda_max);
328 const Real epsilon = std::sqrt(std::numeric_limits<Real>::epsilon());
332 const Real gamma = 0.0;
333 const Real eta = 0.0;
334 const Real delta = 1e-03;
336 pia.template applyTLIMEIteration<PIALS,avoidLinSolverCrime>
337 (gamma,eta,epsilon,piaLS,delta,2,external,x,lambda_min);
340 sigma_min = std::sqrt(lambda_min);
346 if (std::abs(lambda_max) > mtm.infinity_norm())
348 <<
"largest magnitude eigenvalue is greater than "
349 <<
"infinity norm of the matrix!");
353 std::cout <<
" Largest singular value σ_max = "
354 << sigma_max << std::endl;
358 std::cout <<
" Smallest singular value σ_min = "
359 << sigma_min << std::endl;
362 cond_2 = sigma_max / sigma_min;
366 std::cout <<
" 2-norm condition number cond_2 = "
367 << cond_2 << std::endl;
375 const BCRSMatrix& m_;
379 const unsigned int arppp_a_verbosity_level_;
380 const unsigned int pia_verbosity_level_;
384 mutable Real cond_2_;
385 mutable bool symmetricity_assumed_;
Helper functions for determining the vector/matrix block level.
This file implements a vector space as a tensor product of a given vector space. The number of compon...
Wrapper to use a range of ARPACK++ eigenvalue solvers.
Definition: arpackpp.hh:245
Bi-conjugate Gradient Stabilized (BiCG-STAB)
Definition: solvers.hh:419
A vector of blocks with memory management.
Definition: bvector.hh:392
Base class for Dune-Exceptions.
Definition: exceptions.hh:96
vector space out of a tensor product of fields.
Definition: fvector.hh:91
Iterative eigenvalue algorithms based on power iteration.
Definition: poweriteration.hh:176
Sequential SOR preconditioner.
Definition: preconditioners.hh:262
SuperLu Solver.
Definition: superlu.hh:271
Class template which yields information related to a square matrix like its spectral (i....
Definition: matrixinfo.hh:43
MatrixInfo(const BCRSMatrix &m, const bool verbose=false, const unsigned int arppp_a_verbosity_level=0, const unsigned int pia_verbosity_level=0)
Construct from required parameters.
Definition: matrixinfo.hh:63
Real computeSymCond2() const
Definition: matrixinfo.hh:123
Real getCond2(const bool assume_symmetric=true) const
return spectral (i.e. 2-norm) condition number of the matrix
Definition: matrixinfo.hh:92
BCRSMatrix::field_type Real
Type of the underlying field of the matrix.
Definition: matrixinfo.hh:46
static const int bvBlockSize
Definition: matrixinfo.hh:114
Real computeNonSymCond2() const
Definition: matrixinfo.hh:243
A few common exception classes.
Implements a vector constructed from a given type representing a field and a compile-time given size.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
void transposeMatMultMat(BCRSMatrix< FieldMatrix< T, n, m >, A > &res, const BCRSMatrix< FieldMatrix< T, k, n >, A1 > &mat, const BCRSMatrix< FieldMatrix< T, k, m >, A2 > &matt, bool tryHard=false)
Calculate product of a transposed sparse matrix with another sparse matrices ( ).
Definition: matrixmatrix.hh:590
provides functions for sparse matrix matrix multiplication.
Define general preconditioner interface.
Implementations of the inverse operator interface.
Classes for using SuperLU with ISTL matrices.