3#ifndef DUNE_ISTL_EIGENVALUE_TEST_MATRIXINFO_HH
4#define DUNE_ISTL_EIGENVALUE_TEST_MATRIXINFO_HH
20#include "../arpackpp.hh"
21#include "../poweriteration.hh"
38template <
typename BCRSMatrix>
43 typedef typename BCRSMatrix::field_type
Real;
61 const bool verbose =
false,
62 const unsigned int arppp_a_verbosity_level = 0,
63 const unsigned int pia_verbosity_level = 0)
66 arppp_a_verbosity_level_(arppp_a_verbosity_level*verbose),
67 pia_verbosity_level_(pia_verbosity_level*verbose),
68 cond_2_(-1.0), symmetricity_assumed_(false)
72 (BCRSMatrix::blocklevel == 2,
73 "Only BCRSMatrices with blocklevel 2 are supported.");
77 (BCRSMatrix::block_type::rows == BCRSMatrix::block_type::cols,
78 "Only BCRSMatrices with square blocks are supported.");
81 const int nrows = m_.M() * BCRSMatrix::block_type::rows;
82 const int ncols = m_.N() * BCRSMatrix::block_type::cols;
85 << nrows <<
"x" << ncols <<
").");
91 if (cond_2_ == -1.0 || symmetricity_assumed_ != assume_symmetric)
94 std::cout <<
" MatrixInfo: Computing 2-norm condition number"
95 << (assume_symmetric ?
" (assuming that matrix is symmetric)." :
".")
103 symmetricity_assumed_ = assume_symmetric;
124 Real lambda_max, lambda_min, cond_2;
132 typedef Dune::ArPackPlusPlus_Algorithms<BCRSMatrix,BlockVector> ARPPP_A;
133 const ARPPP_A arppp_a(m_,100000,arppp_a_verbosity_level_);
138 PIA pia(m_,20000,pia_verbosity_level_);
139 static const bool avoidLinSolverCrime =
true;
144 typedef Dune::SuperLU<BCRSMatrix> PIALS;
145 const unsigned int piaLS_verbosity = 0;
146 PIALS piaLS(pia.getIterationMatrix(),piaLS_verbosity);
151 typename PIA::IterationOperator::domain_type,
152 typename PIA::IterationOperator::range_type> PIAPC;
153 PIAPC piaPC(pia.getIterationMatrix(),2,1.0);
154 const double piaLS_reduction = 1e-02;
155 const unsigned int piaLS_max_iter = 1000;
156 const unsigned int piaLS_verbosity = 0;
158 PIALS piaLS(pia.getIterationOperator(),piaPC,
159 piaLS_reduction,piaLS_max_iter,piaLS_verbosity);
166 const Real epsilon = 0.0;
168 arppp_a.computeSymMaxMagnitude(epsilon,x,lambda_max);
174 const Real epsilonPI = 1e-02;
175 const Real epsilonTLIME = 1e-08;
179 pia.applyPowerIteration(epsilonPI,x,lambda_max);
181 const Real gamma = m_.infinity_norm();
182 const Real eta = 0.0;
183 const Real delta = 1e-03;
185 pia.template applyTLIMEIteration<PIALS,avoidLinSolverCrime>
186 (gamma,eta,epsilonTLIME,piaLS,delta,2,external,x,lambda_max);
194 const Real epsilon = 1e-11;
198 const Real gamma = 0.0;
199 const Real eta = 0.0;
200 const Real delta = 1e-03;
202 pia.template applyTLIMEIteration<PIALS,avoidLinSolverCrime>
203 (gamma,eta,epsilon,piaLS,delta,2,external,x,lambda_min);
210 if (std::abs(lambda_max) > m_.infinity_norm())
212 <<
"largest magnitude eigenvalue is greater than "
213 <<
"infinity norm of the matrix!");
217 std::cout <<
" Largest magnitude eigenvalue λ_max = "
218 << lambda_max << std::endl;
222 std::cout <<
" Smallest magnitude eigenvalue λ_min = "
223 << lambda_min << std::endl;
227 cond_2 = std::abs(lambda_max / lambda_min);
231 std::cout <<
" 2-norm condition number cond_2 = "
232 << cond_2 << std::endl;
244 Real sigma_max, sigma_min, cond_2;
252 typedef Dune::ArPackPlusPlus_Algorithms<BCRSMatrix,BlockVector> ARPPP_A;
253 const ARPPP_A arppp_a(m_,100000,arppp_a_verbosity_level_);
262 Real lambda_max, lambda_min;
267 PIA pia(mtm,20000,pia_verbosity_level_);
268 static const bool avoidLinSolverCrime =
true;
273 typedef Dune::SuperLU<BCRSMatrix> PIALS;
274 const unsigned int piaLS_verbosity = 0;
275 PIALS piaLS(pia.getIterationMatrix(),piaLS_verbosity);
280 typename PIA::IterationOperator::domain_type,
281 typename PIA::IterationOperator::range_type> PIAPC;
282 PIAPC piaPC(pia.getIterationMatrix(),2,1.0);
283 const double piaLS_reduction = 1e-02;
284 const unsigned int piaLS_max_iter = 1000;
285 const unsigned int piaLS_verbosity = 0;
287 PIALS piaLS(pia.getIterationOperator(),piaPC,
288 piaLS_reduction,piaLS_max_iter,piaLS_verbosity);
294 const Real epsilon = 0.0;
296 arppp_a.computeNonSymMax(epsilon,x,sigma_max);
303 const Real epsilonPI = 1e-02;
304 const Real epsilonTLIME = 1e-08;
308 pia.applyPowerIteration(epsilonPI,x,lambda_max);
311 const Real eta = 0.0;
312 const Real delta = 1e-03;
314 pia.template applyTLIMEIteration<PIALS,avoidLinSolverCrime>
315 (gamma,eta,epsilonTLIME,piaLS,delta,2,external,x,lambda_max);
318 sigma_max = std::sqrt(lambda_max);
325 const Real epsilon = 1e-11;
329 const Real gamma = 0.0;
330 const Real eta = 0.0;
331 const Real delta = 1e-03;
333 pia.template applyTLIMEIteration<PIALS,avoidLinSolverCrime>
334 (gamma,eta,epsilon,piaLS,delta,2,external,x,lambda_min);
337 sigma_min = std::sqrt(lambda_min);
343 if (std::abs(lambda_max) > mtm.infinity_norm())
345 <<
"largest magnitude eigenvalue is greater than "
346 <<
"infinity norm of the matrix!");
350 std::cout <<
" Largest singular value σ_max = "
351 << sigma_max << std::endl;
355 std::cout <<
" Smallest singular value σ_min = "
356 << sigma_min << std::endl;
359 cond_2 = sigma_max / sigma_min;
363 std::cout <<
" 2-norm condition number cond_2 = "
364 << cond_2 << std::endl;
372 const BCRSMatrix& m_;
376 const unsigned int arppp_a_verbosity_level_;
377 const unsigned int pia_verbosity_level_;
381 mutable Real cond_2_;
382 mutable bool symmetricity_assumed_;
This file implements a vector space as a tensor product of a given vector space. The number of compon...
Bi-conjugate Gradient Stabilized (BiCG-STAB)
Definition: solvers.hh:574
A vector of blocks with memory management.
Definition: bvector.hh:313
Base class for Dune-Exceptions.
Definition: exceptions.hh:94
vector space out of a tensor product of fields.
Definition: fvector.hh:93
A class template for performing some iterative eigenvalue algorithms based on power iteration.
Definition: poweriteration.hh:211
Sequential Gauss Seidel preconditioner.
Definition: preconditioners.hh:318
FieldTraits< ft >::real_type infinity_norm() const
infinity norm (maximum of absolute values of entries)
Definition: bvector.hh:211
Class template which yields information related to a square matrix like its spectral (i....
Definition: matrixinfo.hh:40
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:60
Real computeSymCond2() const
Definition: matrixinfo.hh:120
Real getCond2(const bool assume_symmetric=true) const
return spectral (i.e. 2-norm) condition number of the matrix
Definition: matrixinfo.hh:89
BCRSMatrix::field_type Real
Type of the underlying field of the matrix.
Definition: matrixinfo.hh:43
static const int bvBlockSize
Definition: matrixinfo.hh:111
Real computeNonSymCond2() const
Definition: matrixinfo.hh:240
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:216
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.