Dune Core Modules (2.5.0)

matrixinfo.hh
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3#ifndef DUNE_ISTL_EIGENVALUE_TEST_MATRIXINFO_HH
4#define DUNE_ISTL_EIGENVALUE_TEST_MATRIXINFO_HH
5
6#include <cmath> // provides std::abs and std::sqrt
7#include <cassert> // provides assert
8
9#include <iostream> // provides std::cout, std::endl
10
11#include <dune/common/exceptions.hh> // provides DUNE_THROW(...), Dune::Exception
12#include <dune/common/fvector.hh> // provides Dune::FieldVector
13
14#include <dune/istl/bvector.hh> // provides Dune::BlockVector
15#include <dune/istl/superlu.hh> // provides Dune::SuperLU
16#include <dune/istl/preconditioners.hh> // provides Dune::SeqGS
17#include <dune/istl/solvers.hh> // provides Dune::BiCGSTABSolver
18#include <dune/istl/matrixmatrix.hh> // provides Dune::transposeMatMultMat(...)
19
20#include "../arpackpp.hh" // provides Dune::ArPackPlusPlus_Algorithms
21#include "../poweriteration.hh" // provides Dune::PowerIteration_Algorithms
22
23
38template <typename BCRSMatrix>
40{
41public:
43 typedef typename BCRSMatrix::field_type Real;
44
45public:
60 MatrixInfo (const BCRSMatrix& m,
61 const bool verbose = false,
62 const unsigned int arppp_a_verbosity_level = 0,
63 const unsigned int pia_verbosity_level = 0)
64 : m_(m),
65 verbose_(verbose),
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)
69 {
70 // assert that BCRSMatrix type has blocklevel 2
71 static_assert
72 (BCRSMatrix::blocklevel == 2,
73 "Only BCRSMatrices with blocklevel 2 are supported.");
74
75 // assert that BCRSMatrix type has square blocks
76 static_assert
77 (BCRSMatrix::block_type::rows == BCRSMatrix::block_type::cols,
78 "Only BCRSMatrices with square blocks are supported.");
79
80 // assert that m_ is square
81 const int nrows = m_.M() * BCRSMatrix::block_type::rows;
82 const int ncols = m_.N() * BCRSMatrix::block_type::cols;
83 if (nrows != ncols)
84 DUNE_THROW(Dune::Exception,"Matrix is not square ("
85 << nrows << "x" << ncols << ").");
86 }
87
89 inline Real getCond2 (const bool assume_symmetric = true) const
90 {
91 if (cond_2_ == -1.0 || symmetricity_assumed_ != assume_symmetric)
92 {
93 if (verbose_)
94 std::cout << " MatrixInfo: Computing 2-norm condition number"
95 << (assume_symmetric ? " (assuming that matrix is symmetric)." : ".")
96 << std::endl;
97
98 if (assume_symmetric)
99 cond_2_ = computeSymCond2(); // assume that m_ is symmetric
100 else
101 cond_2_ = computeNonSymCond2(); // don't assume that m_ is symmetric
102
103 symmetricity_assumed_ = assume_symmetric;
104 }
105 return cond_2_;
106 }
107
108protected:
111 static const int bvBlockSize = BCRSMatrix::block_type::rows;
114
115protected:
120 inline Real computeSymCond2 () const
121 {
122 // 1) allocate memory for largest and smallest magnitude eigenvalue
123 // as well as the spectral (i.e. 2-norm) condition number
124 Real lambda_max, lambda_min, cond_2;
125
126 // 2) allocate memory for starting vectors and approximated
127 // eigenvectors
128 BlockVector x(m_.M());
129
130#if HAVE_ARPACKPP
131 // 3) setup ARPACK++ eigenvalue algorithms
132 typedef Dune::ArPackPlusPlus_Algorithms<BCRSMatrix,BlockVector> ARPPP_A;
133 const ARPPP_A arppp_a(m_,100000,arppp_a_verbosity_level_);
134#endif // HAVE_ARPACKPP
135
136 // 4) setup power iteration based iterative eigenvalue algorithms
138 PIA pia(m_,20000,pia_verbosity_level_);
139 static const bool avoidLinSolverCrime = true;
140
141#if HAVE_SUPERLU
142 // 5) select a linear solver for power iteration based iterative
143 // eigenvalue algorithms
144 typedef Dune::SuperLU<BCRSMatrix> PIALS;
145 const unsigned int piaLS_verbosity = 0;
146 PIALS piaLS(pia.getIterationMatrix(),piaLS_verbosity);
147#else
148 // 5) select a linear solver for power iteration based iterative
149 // eigenvalue algorithms
150 typedef Dune::SeqGS<BCRSMatrix,
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);
160#endif // HAVE_SUPERLU
161
162#if HAVE_ARPACKPP
163 // 6) get largest magnitude eigenvalue via ARPACK++
164 // (assume that m_ is symmetric)
165 {
166 const Real epsilon = 0.0;
167 // x = 1.0; (not supported yet)
168 arppp_a.computeSymMaxMagnitude(epsilon,x,lambda_max);
169 }
170#else
171 // 6) get largest magnitude eigenvalue via a combination of
172 // power and TLIME iteration (assume that m_ is symmetric)
173 {
174 const Real epsilonPI = 1e-02;
175 const Real epsilonTLIME = 1e-08;
176 x = 1.0;
177 // 6.1) perform power iteration for largest magnitude
178 // eigenvalue (predictor)
179 pia.applyPowerIteration(epsilonPI,x,lambda_max);
180 // 6.2) perform TLIME iteration to improve result (corrector)
181 const Real gamma = m_.infinity_norm();
182 const Real eta = 0.0;
183 const Real delta = 1e-03;
184 bool external;
185 pia.template applyTLIMEIteration<PIALS,avoidLinSolverCrime>
186 (gamma,eta,epsilonTLIME,piaLS,delta,2,external,x,lambda_max);
187 assert(external);
188 }
189#endif // HAVE_ARPACKPP
190
191 // 7) get smallest magnitude eigenvalue via TLIME iteration
192 // (assume that m_ is symmetric)
193 {
194 const Real epsilon = 1e-11;
195 x = 1.0;
196 // 7.1) perform TLIME iteration for smallest magnitude
197 // eigenvalue
198 const Real gamma = 0.0;
199 const Real eta = 0.0;
200 const Real delta = 1e-03;
201 bool external;
202 pia.template applyTLIMEIteration<PIALS,avoidLinSolverCrime>
203 (gamma,eta,epsilon,piaLS,delta,2,external,x,lambda_min);
204 assert(external);
205 }
206
207 // 8) check largest magnitude eigenvalue (we have
208 // ||m|| >= |lambda| for each eigenvalue lambda
209 // of a matrix m and each matrix norm ||.||)
210 if (std::abs(lambda_max) > m_.infinity_norm())
211 DUNE_THROW(Dune::Exception,"Absolute value of approximated "
212 << "largest magnitude eigenvalue is greater than "
213 << "infinity norm of the matrix!");
214
215 // 9) output largest magnitude eigenvalue
216 if (verbose_)
217 std::cout << " Largest magnitude eigenvalue λ_max = "
218 << lambda_max << std::endl;
219
220 // 10) output smallest magnitude eigenvalue
221 if (verbose_)
222 std::cout << " Smallest magnitude eigenvalue λ_min = "
223 << lambda_min << std::endl;
224
225 // 11) compute spectral (i.e. 2-norm) condition number
226 // (assume that m_ is symmetric)
227 cond_2 = std::abs(lambda_max / lambda_min);
228
229 // 12) output spectral (i.e. 2-norm) condition number
230 if (verbose_)
231 std::cout << " 2-norm condition number cond_2 = "
232 << cond_2 << std::endl;
233
234 // 13) return spectral (i.e. 2-norm) condition number
235 return cond_2;
236 }
237
240 inline Real computeNonSymCond2 () const
241 {
242 // 1) allocate memory for largest and smallest singular value
243 // as well as the spectral (i.e. 2-norm) condition number
244 Real sigma_max, sigma_min, cond_2;
245
246 // 2) allocate memory for starting vectors and approximated
247 // eigenvectors respectively singular vectors
248 BlockVector x(m_.M());
249
250#if HAVE_ARPACKPP
251 // 3) setup ARPACK++ eigenvalue algorithms
252 typedef Dune::ArPackPlusPlus_Algorithms<BCRSMatrix,BlockVector> ARPPP_A;
253 const ARPPP_A arppp_a(m_,100000,arppp_a_verbosity_level_);
254#endif // HAVE_ARPACKPP
255
256 // 4) compute m^t*m
257 BCRSMatrix mtm;
258 Dune::transposeMatMultMat(mtm,m_,m_);
259
260 // 5) allocate memory for largest and smallest magnitude
261 // eigenvalue of m^t*m
262 Real lambda_max, lambda_min;
263
264 // 6) setup power iteration based iterative eigenvalue algorithms
265 // for m^t*m
267 PIA pia(mtm,20000,pia_verbosity_level_);
268 static const bool avoidLinSolverCrime = true;
269
270#if HAVE_SUPERLU
271 // 7) select a linear solver for power iteration based iterative
272 // eigenvalue algorithms for m^t*m
273 typedef Dune::SuperLU<BCRSMatrix> PIALS;
274 const unsigned int piaLS_verbosity = 0;
275 PIALS piaLS(pia.getIterationMatrix(),piaLS_verbosity);
276#else
277 // 7) select a linear solver for power iteration based iterative
278 // eigenvalue algorithms for m^t*m
279 typedef Dune::SeqGS<BCRSMatrix,
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);
289#endif // HAVE_SUPERLU
290
291#if HAVE_ARPACKPP
292 // 8) get largest singular value via ARPACK++
293 {
294 const Real epsilon = 0.0;
295 // x = 1.0; (not supported yet)
296 arppp_a.computeNonSymMax(epsilon,x,sigma_max);
297 }
298#else
299 // 8) get largest singular value as square root of the largest
300 // magnitude eigenvalue of m^t*m via a combination of power
301 // and TLIME iteration
302 {
303 const Real epsilonPI = 1e-02;
304 const Real epsilonTLIME = 1e-08;
305 x = 1.0;
306 // 8.1) perform power iteration for largest magnitude
307 // eigenvalue of m^t*m (predictor)
308 pia.applyPowerIteration(epsilonPI,x,lambda_max);
309 // 8.2) perform TLIME iteration to improve result (corrector)
310 const Real gamma = mtm.infinity_norm();
311 const Real eta = 0.0;
312 const Real delta = 1e-03;
313 bool external;
314 pia.template applyTLIMEIteration<PIALS,avoidLinSolverCrime>
315 (gamma,eta,epsilonTLIME,piaLS,delta,2,external,x,lambda_max);
316 assert(external);
317 // 8.3) get largest singular value
318 sigma_max = std::sqrt(lambda_max);
319 }
320#endif // HAVE_ARPACKPP
321
322 // 9) get smallest singular value as square root of the smallest
323 // magnitude eigenvalue of m^t*m via TLIME iteration
324 {
325 const Real epsilon = 1e-11;
326 x = 1.0;
327 // 9.1) perform TLIME iteration for smallest magnitude
328 // eigenvalue of m^t*m
329 const Real gamma = 0.0;
330 const Real eta = 0.0;
331 const Real delta = 1e-03;
332 bool external;
333 pia.template applyTLIMEIteration<PIALS,avoidLinSolverCrime>
334 (gamma,eta,epsilon,piaLS,delta,2,external,x,lambda_min);
335 assert(external);
336 // 9.2) get smallest singular value
337 sigma_min = std::sqrt(lambda_min);
338 }
339
340 // 10) check largest magnitude eigenvalue (we have
341 // ||m|| >= |lambda| for each eigenvalue lambda
342 // of a matrix m and each matrix norm ||.||)
343 if (std::abs(lambda_max) > mtm.infinity_norm())
344 DUNE_THROW(Dune::Exception,"Absolute value of approximated "
345 << "largest magnitude eigenvalue is greater than "
346 << "infinity norm of the matrix!");
347
348 // 11) output largest singular value
349 if (verbose_)
350 std::cout << " Largest singular value σ_max = "
351 << sigma_max << std::endl;
352
353 // 12) output smallest singular value
354 if (verbose_)
355 std::cout << " Smallest singular value σ_min = "
356 << sigma_min << std::endl;
357
358 // 13) compute spectral (i.e. 2-norm) condition number
359 cond_2 = sigma_max / sigma_min;
360
361 // 14) output spectral (i.e. 2-norm) condition number
362 if (verbose_)
363 std::cout << " 2-norm condition number cond_2 = "
364 << cond_2 << std::endl;
365
366 // 15) return spectral (i.e. 2-norm) condition number
367 return cond_2;
368 }
369
370protected:
371 // parameters related to computation of matrix information
372 const BCRSMatrix& m_;
373
374 // verbosity setting
375 const bool verbose_;
376 const unsigned int arppp_a_verbosity_level_;
377 const unsigned int pia_verbosity_level_;
378
379 // memory for storing matrix information
380 // (mutable as matrix information is computed on demand)
381 mutable Real cond_2_;
382 mutable bool symmetricity_assumed_;
383};
384
385
386#endif // DUNE_ISTL_EIGENVALUE_TEST_MATRIXINFO_HH
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.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 23, 23:29, 2024)