Dune Core Modules (unstable)

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