Dune Core Modules (2.5.0)

arpackpp.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_ARPACKPP_HH
4#define DUNE_ISTL_EIGENVALUE_ARPACKPP_HH
5
6#if HAVE_ARPACKPP
7
8#include <cmath> // provides std::abs, std::pow, std::sqrt
9
10#include <iostream> // provides std::cout, std::endl
11#include <string> // provides std::string
12
13#include <dune/common/fvector.hh> // provides Dune::FieldVector
14#include <dune/common/exceptions.hh> // provides DUNE_THROW(...)
15
16#include <dune/istl/bvector.hh> // provides Dune::BlockVector
17#include <dune/istl/istlexception.hh> // provides Dune::ISTLError
18#include <dune/istl/io.hh> // provides Dune::printvector(...)
19
20#ifdef Status
21#undef Status // prevent preprocessor from damaging the ARPACK++
22 // code when "X11/Xlib.h" is included (the latter
23 // defines Status as "#define Status int" and
24 // ARPACK++ provides a class with a method called
25 // Status)
26#endif
27#include "arssym.h" // provides ARSymStdEig
28
29
30namespace Dune
31{
32
48 template <class BCRSMatrix>
49 class ArPackPlusPlus_BCRSMatrixWrapper
50 {
51 public:
53 typedef typename BCRSMatrix::field_type Real;
54
55 public:
57 ArPackPlusPlus_BCRSMatrixWrapper (const BCRSMatrix& A)
58 : A_(A),
59 m_(A_.M() * mBlock), n_(A_.N() * nBlock)
60 {
61 // assert that BCRSMatrix type has blocklevel 2
62 static_assert
64 "Only BCRSMatrices with blocklevel 2 are supported.");
65
66 // allocate memory for auxiliary block vector objects
67 // which are compatible to matrix rows / columns
68 domainBlockVector.resize(A_.N(),false);
69 rangeBlockVector.resize(A_.M(),false);
70 }
71
73 inline void multMv (Real* v, Real* w)
74 {
75 // get vector v as an object of appropriate type
76 arrayToDomainBlockVector(v,domainBlockVector);
77
78 // perform matrix-vector product
79 A_.mv(domainBlockVector,rangeBlockVector);
80
81 // get vector w from object of appropriate type
82 rangeBlockVectorToArray(rangeBlockVector,w);
83 };
84
86 inline void multMtMv (Real* v, Real* w)
87 {
88 // get vector v as an object of appropriate type
89 arrayToDomainBlockVector(v,domainBlockVector);
90
91 // perform matrix-vector product
92 A_.mv(domainBlockVector,rangeBlockVector);
93 A_.mtv(rangeBlockVector,domainBlockVector);
94
95 // get vector w from object of appropriate type
96 domainBlockVectorToArray(domainBlockVector,w);
97 };
98
100 inline void multMMtv (Real* v, Real* w)
101 {
102 // get vector v as an object of appropriate type
103 arrayToRangeBlockVector(v,rangeBlockVector);
104
105 // perform matrix-vector product
106 A_.mtv(rangeBlockVector,domainBlockVector);
107 A_.mv(domainBlockVector,rangeBlockVector);
108
109 // get vector w from object of appropriate type
110 rangeBlockVectorToArray(rangeBlockVector,w);
111 };
112
114 inline int nrows () const { return m_; }
115
117 inline int ncols () const { return n_; }
118
119 protected:
120 // Number of rows and columns in each block of the matrix
121 constexpr static int mBlock = BCRSMatrix::block_type::rows;
122 constexpr static int nBlock = BCRSMatrix::block_type::cols;
123
124 // Type of vectors in the domain of the linear map associated with
125 // the matrix, i.e. block vectors compatible to matrix rows
126 constexpr static int dbvBlockSize = nBlock;
127 typedef Dune::FieldVector<Real,dbvBlockSize> DomainBlockVectorBlock;
128 typedef Dune::BlockVector<DomainBlockVectorBlock> DomainBlockVector;
129
130 // Type of vectors in the range of the linear map associated with
131 // the matrix, i.e. block vectors compatible to matrix columns
132 constexpr static int rbvBlockSize = mBlock;
133 typedef Dune::FieldVector<Real,rbvBlockSize> RangeBlockVectorBlock;
134 typedef Dune::BlockVector<RangeBlockVectorBlock> RangeBlockVector;
135
136 // Types for vector index access
137 typedef typename DomainBlockVector::size_type dbv_size_type;
138 typedef typename RangeBlockVector::size_type rbv_size_type;
139 typedef typename DomainBlockVectorBlock::size_type dbvb_size_type;
140 typedef typename RangeBlockVectorBlock::size_type rbvb_size_type;
141
142 // Get vector v from a block vector object which is compatible to
143 // matrix rows
144 static inline void
145 domainBlockVectorToArray (const DomainBlockVector& dbv, Real* v)
146 {
147 for (dbv_size_type block = 0; block < dbv.N(); ++block)
148 for (dbvb_size_type iBlock = 0; iBlock < dbvBlockSize; ++iBlock)
149 v[block*dbvBlockSize + iBlock] = dbv[block][iBlock];
150 }
151
152 // Get vector v from a block vector object which is compatible to
153 // matrix columns
154 static inline void
155 rangeBlockVectorToArray (const RangeBlockVector& rbv, Real* v)
156 {
157 for (rbv_size_type block = 0; block < rbv.N(); ++block)
158 for (rbvb_size_type iBlock = 0; iBlock < rbvBlockSize; ++iBlock)
159 v[block*rbvBlockSize + iBlock] = rbv[block][iBlock];
160 }
161
162 public:
165 static inline void arrayToDomainBlockVector (const Real* v,
166 DomainBlockVector& dbv)
167 {
168 for (dbv_size_type block = 0; block < dbv.N(); ++block)
169 for (dbvb_size_type iBlock = 0; iBlock < dbvBlockSize; ++iBlock)
170 dbv[block][iBlock] = v[block*dbvBlockSize + iBlock];
171 }
172
175 static inline void arrayToRangeBlockVector (const Real* v,
176 RangeBlockVector& rbv)
177 {
178 for (rbv_size_type block = 0; block < rbv.N(); ++block)
179 for (rbvb_size_type iBlock = 0; iBlock < rbvBlockSize; ++iBlock)
180 rbv[block][iBlock] = v[block*rbvBlockSize + iBlock];
181 }
182
183 protected:
184 // The DUNE-ISTL BCRSMatrix
185 const BCRSMatrix& A_;
186
187 // Number of rows and columns in the matrix
188 const int m_, n_;
189
190 // Auxiliary block vector objects which are
191 // compatible to matrix rows / columns
192 mutable DomainBlockVector domainBlockVector;
193 mutable RangeBlockVector rangeBlockVector;
194 };
195
196
234 template <typename BCRSMatrix, typename BlockVector>
235 class ArPackPlusPlus_Algorithms
236 {
237 public:
238 typedef typename BlockVector::field_type Real;
239
240 public:
259 ArPackPlusPlus_Algorithms (const BCRSMatrix& m,
260 const unsigned int nIterationsMax = 100000,
261 const unsigned int verbosity_level = 0)
262 : m_(m), nIterationsMax_(nIterationsMax),
263 verbosity_level_(verbosity_level),
264 nIterations_(0),
265 title_(" ArPackPlusPlus_Algorithms: "),
266 blank_(title_.length(),' ')
267 {}
268
280 inline void computeSymMaxMagnitude (const Real& epsilon,
281 BlockVector& x, Real& lambda) const
282 {
283 // print verbosity information
284 if (verbosity_level_ > 0)
285 std::cout << title_ << "Computing an approximation of "
286 << "the dominant eigenvalue of a matrix which "
287 << "is assumed to be symmetric." << std::endl;
288
289 // use type ArPackPlusPlus_BCRSMatrixWrapper to store matrix information
290 // and to perform the product A*v (LU decomposition is not used)
291 typedef ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
292 WrappedMatrix A(m_);
293
294 // get number of rows and columns in A
295 const int nrows = A.nrows();
296 const int ncols = A.ncols();
297
298 // assert that A is square
299 if (nrows != ncols)
300 DUNE_THROW(Dune::ISTLError,"Matrix is not square ("
301 << nrows << "x" << ncols << ").");
302
303 // allocate memory for variables, set parameters
304 const int nev = 1; // Number of eigenvalues to compute
305 const int ncv = 20; // Number of Arnoldi vectors generated at each iteration (0 == auto)
306 const Real tol = epsilon; // Stopping tolerance (relative accuracy of Ritz values) (0 == machine precision)
307 const int maxit = nIterationsMax_*nev; // Maximum number of Arnoldi update iterations allowed (0 == 100*nev)
308 Real* ev = new Real[nev]; // Computed eigenvalues of A
309 const bool ivec = true; // Flag deciding if eigenvectors shall be determined
310 int nconv; // Number of converged eigenvalues
311
312 // define what we need: eigenvalues with largest magnitude
313 char which[] = "LM";
314 ARSymStdEig<Real,WrappedMatrix>
315 dprob(nrows, nev, &A, &WrappedMatrix::multMv, which, ncv, tol, maxit);
316
317 // set ARPACK verbosity mode if requested
318 if (verbosity_level_ > 3) dprob.Trace();
319
320 // find eigenvalues and eigenvectors of A, obtain the eigenvalues
321 nconv = dprob.Eigenvalues(ev,ivec);
322
323 // obtain approximated dominant eigenvalue of A
324 lambda = ev[nev-1];
325
326 // obtain associated approximated eigenvector of A
327 Real* x_raw = dprob.RawEigenvector(nev-1);
328 WrappedMatrix::arrayToDomainBlockVector(x_raw,x);
329
330 // obtain number of Arnoldi update iterations actually taken
331 nIterations_ = dprob.GetIter();
332
333 // compute residual norm
334 BlockVector r(x);
335 Real* Ax_raw = new Real[nrows];
336 A.multMv(x_raw,Ax_raw);
337 WrappedMatrix::arrayToDomainBlockVector(Ax_raw,r);
338 r.axpy(-lambda,x);
339 const Real r_norm = r.two_norm();
340
341 // print verbosity information
342 if (verbosity_level_ > 0)
343 {
344 if (verbosity_level_ > 1)
345 {
346 // print some information about the problem
347 std::cout << blank_ << "Obtained eigenvalues of A by solving "
348 << "A*x = λ*x using the ARPACK++ class ARSym"
349 << "StdEig:" << std::endl;
350 std::cout << blank_ << " converged eigenvalues of A: "
351 << nconv << " / " << nev << std::endl;
352 std::cout << blank_ << " dominant eigenvalue of A: "
353 << lambda << std::endl;
354 }
355 std::cout << blank_ << "Result ("
356 << "#iterations = " << nIterations_ << ", "
357 << "║residual║_2 = " << r_norm << "): "
358 << "λ = " << lambda << std::endl;
359 if (verbosity_level_ > 2)
360 {
361 // print approximated eigenvector via DUNE-ISTL I/O methods
362 Dune::printvector(std::cout,x,blank_+"x",blank_+"row");
363 }
364 }
365
366 // free dynamically allocated memory
367 delete[] Ax_raw;
368 delete[] ev;
369 }
370
382 inline void computeSymMinMagnitude (const Real& epsilon,
383 BlockVector& x, Real& lambda) const
384 {
385 // print verbosity information
386 if (verbosity_level_ > 0)
387 std::cout << title_ << "Computing an approximation of the "
388 << "least dominant eigenvalue of a matrix which "
389 << "is assumed to be symmetric." << std::endl;
390
391 // use type ArPackPlusPlus_BCRSMatrixWrapper to store matrix information
392 // and to perform the product A*v (LU decomposition is not used)
393 typedef ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
394 WrappedMatrix A(m_);
395
396 // get number of rows and columns in A
397 const int nrows = A.nrows();
398 const int ncols = A.ncols();
399
400 // assert that A is square
401 if (nrows != ncols)
402 DUNE_THROW(Dune::ISTLError,"Matrix is not square ("
403 << nrows << "x" << ncols << ").");
404
405 // allocate memory for variables, set parameters
406 const int nev = 1; // Number of eigenvalues to compute
407 const int ncv = 20; // Number of Arnoldi vectors generated at each iteration (0 == auto)
408 const Real tol = epsilon; // Stopping tolerance (relative accuracy of Ritz values) (0 == machine precision)
409 const int maxit = nIterationsMax_*nev; // Maximum number of Arnoldi update iterations allowed (0 == 100*nev)
410 Real* ev = new Real[nev]; // Computed eigenvalues of A
411 const bool ivec = true; // Flag deciding if eigenvectors shall be determined
412 int nconv; // Number of converged eigenvalues
413
414 // define what we need: eigenvalues with smallest magnitude
415 char which[] = "SM";
416 ARSymStdEig<Real,WrappedMatrix>
417 dprob(nrows, nev, &A, &WrappedMatrix::multMv, which, ncv, tol, maxit);
418
419 // set ARPACK verbosity mode if requested
420 if (verbosity_level_ > 3) dprob.Trace();
421
422 // find eigenvalues and eigenvectors of A, obtain the eigenvalues
423 nconv = dprob.Eigenvalues(ev,ivec);
424
425 // obtain approximated least dominant eigenvalue of A
426 lambda = ev[nev-1];
427
428 // obtain associated approximated eigenvector of A
429 Real* x_raw = dprob.RawEigenvector(nev-1);
430 WrappedMatrix::arrayToDomainBlockVector(x_raw,x);
431
432 // obtain number of Arnoldi update iterations actually taken
433 nIterations_ = dprob.GetIter();
434
435 // compute residual norm
436 BlockVector r(x);
437 Real* Ax_raw = new Real[nrows];
438 A.multMv(x_raw,Ax_raw);
439 WrappedMatrix::arrayToDomainBlockVector(Ax_raw,r);
440 r.axpy(-lambda,x);
441 const Real r_norm = r.two_norm();
442
443 // print verbosity information
444 if (verbosity_level_ > 0)
445 {
446 if (verbosity_level_ > 1)
447 {
448 // print some information about the problem
449 std::cout << blank_ << "Obtained eigenvalues of A by solving "
450 << "A*x = λ*x using the ARPACK++ class ARSym"
451 << "StdEig:" << std::endl;
452 std::cout << blank_ << " converged eigenvalues of A: "
453 << nconv << " / " << nev << std::endl;
454 std::cout << blank_ << " least dominant eigenvalue of A: "
455 << lambda << std::endl;
456 }
457 std::cout << blank_ << "Result ("
458 << "#iterations = " << nIterations_ << ", "
459 << "║residual║_2 = " << r_norm << "): "
460 << "λ = " << lambda << std::endl;
461 if (verbosity_level_ > 2)
462 {
463 // print approximated eigenvector via DUNE-ISTL I/O methods
464 Dune::printvector(std::cout,x,blank_+"x",blank_+"row");
465 }
466 }
467
468 // free dynamically allocated memory
469 delete[] Ax_raw;
470 delete[] ev;
471 }
472
484 inline void computeSymCond2 (const Real& epsilon, Real& cond_2) const
485 {
486 // print verbosity information
487 if (verbosity_level_ > 0)
488 std::cout << title_ << "Computing an approximation of the "
489 << "spectral condition number of a matrix which "
490 << "is assumed to be symmetric." << std::endl;
491
492 // use type ArPackPlusPlus_BCRSMatrixWrapper to store matrix information
493 // and to perform the product A*v (LU decomposition is not used)
494 typedef ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
495 WrappedMatrix A(m_);
496
497 // get number of rows and columns in A
498 const int nrows = A.nrows();
499 const int ncols = A.ncols();
500
501 // assert that A is square
502 if (nrows != ncols)
503 DUNE_THROW(Dune::ISTLError,"Matrix is not square ("
504 << nrows << "x" << ncols << ").");
505
506 // allocate memory for variables, set parameters
507 const int nev = 2; // Number of eigenvalues to compute
508 const int ncv = 20; // Number of Arnoldi vectors generated at each iteration (0 == auto)
509 const Real tol = epsilon; // Stopping tolerance (relative accuracy of Ritz values) (0 == machine precision)
510 const int maxit = nIterationsMax_*nev; // Maximum number of Arnoldi update iterations allowed (0 == 100*nev)
511 Real* ev = new Real[nev]; // Computed eigenvalues of A
512 const bool ivec = true; // Flag deciding if eigenvectors shall be determined
513 int nconv; // Number of converged eigenvalues
514
515 // define what we need: eigenvalues from both ends of the spectrum
516 char which[] = "BE";
517 ARSymStdEig<Real,WrappedMatrix>
518 dprob(nrows, nev, &A, &WrappedMatrix::multMv, which, ncv, tol, maxit);
519
520 // set ARPACK verbosity mode if requested
521 if (verbosity_level_ > 3) dprob.Trace();
522
523 // find eigenvalues and eigenvectors of A, obtain the eigenvalues
524 nconv = dprob.Eigenvalues(ev,ivec);
525
526 // obtain approximated dominant and least dominant eigenvalue of A
527 const Real& lambda_max = ev[nev-1];
528 const Real& lambda_min = ev[0];
529
530 // obtain associated approximated eigenvectors of A
531 Real* x_max_raw = dprob.RawEigenvector(nev-1);
532 Real* x_min_raw = dprob.RawEigenvector(0);
533
534 // obtain approximated spectral condition number of A
535 cond_2 = std::abs(lambda_max / lambda_min);
536
537 // obtain number of Arnoldi update iterations actually taken
538 nIterations_ = dprob.GetIter();
539
540 // compute each residual norm
541 Real* Ax_max_raw = new Real[nrows];
542 Real* Ax_min_raw = new Real[nrows];
543 A.multMv(x_max_raw,Ax_max_raw);
544 A.multMv(x_min_raw,Ax_min_raw);
545 Real r_max_norm = 0.0;
546 Real r_min_norm = 0.0;
547 for (int i = 0; i < nrows; ++i)
548 {
549 r_max_norm += std::pow(Ax_max_raw[i] - lambda_max * x_max_raw[i],2);
550 r_min_norm += std::pow(Ax_min_raw[i] - lambda_min * x_min_raw[i],2);
551 }
552 r_max_norm = std::sqrt(r_max_norm);
553 r_min_norm = std::sqrt(r_min_norm);
554
555 // print verbosity information
556 if (verbosity_level_ > 0)
557 {
558 if (verbosity_level_ > 1)
559 {
560 // print some information about the problem
561 std::cout << blank_ << "Obtained eigenvalues of A by solving "
562 << "A*x = λ*x using the ARPACK++ class ARSym"
563 << "StdEig:" << std::endl;
564 std::cout << blank_ << " converged eigenvalues of A: "
565 << nconv << " / " << nev << std::endl;
566 std::cout << blank_ << " dominant eigenvalue of A: "
567 << lambda_max << std::endl;
568 std::cout << blank_ << " least dominant eigenvalue of A: "
569 << lambda_min << std::endl;
570 std::cout << blank_ << " spectral condition number of A: "
571 << cond_2 << std::endl;
572 }
573 std::cout << blank_ << "Result ("
574 << "#iterations = " << nIterations_ << ", "
575 << "║residual║_2 = {" << r_max_norm << ","
576 << r_min_norm << "}, " << "λ = {"
577 << lambda_max << "," << lambda_min
578 << "}): cond_2 = " << cond_2 << std::endl;
579 }
580
581 // free dynamically allocated memory
582 delete[] Ax_min_raw;
583 delete[] Ax_max_raw;
584 delete[] ev;
585 }
586
600 inline void computeNonSymMax (const Real& epsilon,
601 BlockVector& x, Real& sigma) const
602 {
603 // print verbosity information
604 if (verbosity_level_ > 0)
605 std::cout << title_ << "Computing an approximation of the "
606 << "largest singular value of a matrix which "
607 << "is assumed to be nonsymmetric." << std::endl;
608
609 // use type ArPackPlusPlus_BCRSMatrixWrapper to store matrix information
610 // and to perform the product A^T*A*v (LU decomposition is not used)
611 typedef ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
612 WrappedMatrix A(m_);
613
614 // get number of rows and columns in A
615 const int nrows = A.nrows();
616 const int ncols = A.ncols();
617
618 // assert that A has more rows than columns (extend code later to the opposite case!)
619 if (nrows < ncols)
620 DUNE_THROW(Dune::ISTLError,"Matrix has less rows than "
621 << "columns (" << nrows << "x" << ncols << ")."
622 << " This case is not implemented, yet.");
623
624 // allocate memory for variables, set parameters
625 const int nev = 1; // Number of eigenvalues to compute
626 const int ncv = 20; // Number of Arnoldi vectors generated at each iteration (0 == auto)
627 const Real tol = epsilon; // Stopping tolerance (relative accuracy of Ritz values) (0 == machine precision)
628 const int maxit = nIterationsMax_*nev; // Maximum number of Arnoldi update iterations allowed (0 == 100*nev)
629 Real* ev = new Real[nev]; // Computed eigenvalues of A^T*A
630 const bool ivec = true; // Flag deciding if eigenvectors shall be determined
631 int nconv; // Number of converged eigenvalues
632
633 // define what we need: eigenvalues with largest algebraic value
634 char which[] = "LA";
635 ARSymStdEig<Real,WrappedMatrix>
636 dprob(ncols, nev, &A, &WrappedMatrix::multMtMv, which, ncv, tol, maxit);
637
638 // set ARPACK verbosity mode if requested
639 if (verbosity_level_ > 3) dprob.Trace();
640
641 // find eigenvalues and eigenvectors of A^T*A, obtain the eigenvalues
642 nconv = dprob.Eigenvalues(ev,ivec);
643
644 // obtain approximated largest eigenvalue of A^T*A
645 const Real& lambda = ev[nev-1];
646
647 // obtain associated approximated eigenvector of A^T*A
648 Real* x_raw = dprob.RawEigenvector(nev-1);
649 WrappedMatrix::arrayToDomainBlockVector(x_raw,x);
650
651 // obtain number of Arnoldi update iterations actually taken
652 nIterations_ = dprob.GetIter();
653
654 // compute residual norm
655 BlockVector r(x);
656 Real* AtAx_raw = new Real[ncols];
657 A.multMtMv(x_raw,AtAx_raw);
658 WrappedMatrix::arrayToDomainBlockVector(AtAx_raw,r);
659 r.axpy(-lambda,x);
660 const Real r_norm = r.two_norm();
661
662 // calculate largest singular value of A (note that
663 // x is right-singular / left-singular vector of A)
664 sigma = std::sqrt(lambda);
665
666 // print verbosity information
667 if (verbosity_level_ > 0)
668 {
669 if (verbosity_level_ > 1)
670 {
671 // print some information about the problem
672 std::cout << blank_ << "Obtained singular values of A by sol"
673 << "ving (A^T*A)*x = σ²*x using the ARPACK++ "
674 << "class ARSymStdEig:" << std::endl;
675 std::cout << blank_ << " converged eigenvalues of A^T*A: "
676 << nconv << " / " << nev << std::endl;
677 std::cout << blank_ << " largest eigenvalue of A^T*A: "
678 << lambda << std::endl;
679 std::cout << blank_ << " => largest singular value of A: "
680 << sigma << std::endl;
681 }
682 std::cout << blank_ << "Result ("
683 << "#iterations = " << nIterations_ << ", "
684 << "║residual║_2 = " << r_norm << "): "
685 << "σ = " << sigma << std::endl;
686 if (verbosity_level_ > 2)
687 {
688 // print approximated right-singular / left-singular vector
689 // via DUNE-ISTL I/O methods
690 Dune::printvector(std::cout,x,blank_+"x",blank_+"row");
691 }
692 }
693
694 // free dynamically allocated memory
695 delete[] AtAx_raw;
696 delete[] ev;
697 }
698
712 inline void computeNonSymMin (const Real& epsilon,
713 BlockVector& x, Real& sigma) const
714 {
715 // print verbosity information
716 if (verbosity_level_ > 0)
717 std::cout << title_ << "Computing an approximation of the "
718 << "smallest singular value of a matrix which "
719 << "is assumed to be nonsymmetric." << std::endl;
720
721 // use type ArPackPlusPlus_BCRSMatrixWrapper to store matrix information
722 // and to perform the product A^T*A*v (LU decomposition is not used)
723 typedef ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
724 WrappedMatrix A(m_);
725
726 // get number of rows and columns in A
727 const int nrows = A.nrows();
728 const int ncols = A.ncols();
729
730 // assert that A has more rows than columns (extend code later to the opposite case!)
731 if (nrows < ncols)
732 DUNE_THROW(Dune::ISTLError,"Matrix has less rows than "
733 << "columns (" << nrows << "x" << ncols << ")."
734 << " This case is not implemented, yet.");
735
736 // allocate memory for variables, set parameters
737 const int nev = 1; // Number of eigenvalues to compute
738 const int ncv = 20; // Number of Arnoldi vectors generated at each iteration (0 == auto)
739 const Real tol = epsilon; // Stopping tolerance (relative accuracy of Ritz values) (0 == machine precision)
740 const int maxit = nIterationsMax_*nev; // Maximum number of Arnoldi update iterations allowed (0 == 100*nev)
741 Real* ev = new Real[nev]; // Computed eigenvalues of A^T*A
742 const bool ivec = true; // Flag deciding if eigenvectors shall be determined
743 int nconv; // Number of converged eigenvalues
744
745 // define what we need: eigenvalues with smallest algebraic value
746 char which[] = "SA";
747 ARSymStdEig<Real,WrappedMatrix>
748 dprob(ncols, nev, &A, &WrappedMatrix::multMtMv, which, ncv, tol, maxit);
749
750 // set ARPACK verbosity mode if requested
751 if (verbosity_level_ > 3) dprob.Trace();
752
753 // find eigenvalues and eigenvectors of A^T*A, obtain the eigenvalues
754 nconv = dprob.Eigenvalues(ev,ivec);
755
756 // obtain approximated smallest eigenvalue of A^T*A
757 const Real& lambda = ev[nev-1];
758
759 // obtain associated approximated eigenvector of A^T*A
760 Real* x_raw = dprob.RawEigenvector(nev-1);
761 WrappedMatrix::arrayToDomainBlockVector(x_raw,x);
762
763 // obtain number of Arnoldi update iterations actually taken
764 nIterations_ = dprob.GetIter();
765
766 // compute residual norm
767 BlockVector r(x);
768 Real* AtAx_raw = new Real[ncols];
769 A.multMtMv(x_raw,AtAx_raw);
770 WrappedMatrix::arrayToDomainBlockVector(AtAx_raw,r);
771 r.axpy(-lambda,x);
772 const Real r_norm = r.two_norm();
773
774 // calculate smallest singular value of A (note that
775 // x is right-singular / left-singular vector of A)
776 sigma = std::sqrt(lambda);
777
778 // print verbosity information
779 if (verbosity_level_ > 0)
780 {
781 if (verbosity_level_ > 1)
782 {
783 // print some information about the problem
784 std::cout << blank_ << "Obtained singular values of A by sol"
785 << "ving (A^T*A)*x = σ²*x using the ARPACK++ "
786 << "class ARSymStdEig:" << std::endl;
787 std::cout << blank_ << " converged eigenvalues of A^T*A: "
788 << nconv << " / " << nev << std::endl;
789 std::cout << blank_ << " smallest eigenvalue of A^T*A: "
790 << lambda << std::endl;
791 std::cout << blank_ << " => smallest singular value of A: "
792 << sigma << std::endl;
793 }
794 std::cout << blank_ << "Result ("
795 << "#iterations = " << nIterations_ << ", "
796 << "║residual║_2 = " << r_norm << "): "
797 << "σ = " << sigma << std::endl;
798 if (verbosity_level_ > 2)
799 {
800 // print approximated right-singular / left-singular vector
801 // via DUNE-ISTL I/O methods
802 Dune::printvector(std::cout,x,blank_+"x",blank_+"row");
803 }
804 }
805
806 // free dynamically allocated memory
807 delete[] AtAx_raw;
808 delete[] ev;
809 }
810
821 inline void computeNonSymCond2 (const Real& epsilon, Real& cond_2) const
822 {
823 // print verbosity information
824 if (verbosity_level_ > 0)
825 std::cout << title_ << "Computing an approximation of the "
826 << "spectral condition number of a matrix which "
827 << "is assumed to be nonsymmetric." << std::endl;
828
829 // use type ArPackPlusPlus_BCRSMatrixWrapper to store matrix information
830 // and to perform the product A^T*A*v (LU decomposition is not used)
831 typedef ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
832 WrappedMatrix A(m_);
833
834 // get number of rows and columns in A
835 const int nrows = A.nrows();
836 const int ncols = A.ncols();
837
838 // assert that A has more rows than columns (extend code later to the opposite case!)
839 if (nrows < ncols)
840 DUNE_THROW(Dune::ISTLError,"Matrix has less rows than "
841 << "columns (" << nrows << "x" << ncols << ")."
842 << " This case is not implemented, yet.");
843
844 // allocate memory for variables, set parameters
845 const int nev = 2; // Number of eigenvalues to compute
846 const int ncv = 20; // Number of Arnoldi vectors generated at each iteration (0 == auto)
847 const Real tol = epsilon; // Stopping tolerance (relative accuracy of Ritz values) (0 == machine precision)
848 const int maxit = nIterationsMax_*nev; // Maximum number of Arnoldi update iterations allowed (0 == 100*nev)
849 Real* ev = new Real[nev]; // Computed eigenvalues of A^T*A
850 const bool ivec = true; // Flag deciding if eigenvectors shall be determined
851 int nconv; // Number of converged eigenvalues
852
853 // define what we need: eigenvalues from both ends of the spectrum
854 char which[] = "BE";
855 ARSymStdEig<Real,WrappedMatrix>
856 dprob(ncols, nev, &A, &WrappedMatrix::multMtMv, which, ncv, tol, maxit);
857
858 // set ARPACK verbosity mode if requested
859 if (verbosity_level_ > 3) dprob.Trace();
860
861 // find eigenvalues and eigenvectors of A^T*A, obtain the eigenvalues
862 nconv = dprob.Eigenvalues(ev,ivec);
863
864 // obtain approximated largest and smallest eigenvalue of A^T*A
865 const Real& lambda_max = ev[nev-1];
866 const Real& lambda_min = ev[0];
867
868 // obtain associated approximated eigenvectors of A^T*A
869 Real* x_max_raw = dprob.RawEigenvector(nev-1);
870 Real* x_min_raw = dprob.RawEigenvector(0);
871
872 // obtain number of Arnoldi update iterations actually taken
873 nIterations_ = dprob.GetIter();
874
875 // compute each residual norm
876 Real* AtAx_max_raw = new Real[ncols];
877 Real* AtAx_min_raw = new Real[ncols];
878 A.multMtMv(x_max_raw,AtAx_max_raw);
879 A.multMtMv(x_min_raw,AtAx_min_raw);
880 Real r_max_norm = 0.0;
881 Real r_min_norm = 0.0;
882 for (int i = 0; i < ncols; ++i)
883 {
884 r_max_norm += std::pow(AtAx_max_raw[i] - lambda_max * x_max_raw[i],2);
885 r_min_norm += std::pow(AtAx_min_raw[i] - lambda_min * x_min_raw[i],2);
886 }
887 r_max_norm = std::sqrt(r_max_norm);
888 r_min_norm = std::sqrt(r_min_norm);
889
890 // calculate largest and smallest singular value of A
891 const Real sigma_max = std::sqrt(lambda_max);
892 const Real sigma_min = std::sqrt(lambda_min);
893
894 // obtain approximated spectral condition number of A
895 cond_2 = sigma_max / sigma_min;
896
897 // print verbosity information
898 if (verbosity_level_ > 0)
899 {
900 if (verbosity_level_ > 1)
901 {
902 // print some information about the problem
903 std::cout << blank_ << "Obtained singular values of A by sol"
904 << "ving (A^T*A)*x = σ²*x using the ARPACK++ "
905 << "class ARSymStdEig:" << std::endl;
906 std::cout << blank_ << " converged eigenvalues of A^T*A: "
907 << nconv << " / " << nev << std::endl;
908 std::cout << blank_ << " largest eigenvalue of A^T*A: "
909 << lambda_max << std::endl;
910 std::cout << blank_ << " smallest eigenvalue of A^T*A: "
911 << lambda_min << std::endl;
912 std::cout << blank_ << " => largest singular value of A: "
913 << sigma_max << std::endl;
914 std::cout << blank_ << " => smallest singular value of A: "
915 << sigma_min << std::endl;
916 }
917 std::cout << blank_ << "Result ("
918 << "#iterations = " << nIterations_ << ", "
919 << "║residual║_2 = {" << r_max_norm << ","
920 << r_min_norm << "}, " << "σ = {"
921 << sigma_max << "," << sigma_min
922 << "}): cond_2 = " << cond_2 << std::endl;
923 }
924
925 // free dynamically allocated memory
926 delete[] AtAx_min_raw;
927 delete[] AtAx_max_raw;
928 delete[] ev;
929 }
930
935 inline unsigned int getIterationCount () const
936 {
937 if (nIterations_ == 0)
938 DUNE_THROW(Dune::ISTLError,"No algorithm applied, yet.");
939
940 return nIterations_;
941 }
942
943 protected:
944 // parameters related to iterative eigenvalue algorithms
945 const BCRSMatrix& m_;
946 const unsigned int nIterationsMax_;
947
948 // verbosity setting
949 const unsigned int verbosity_level_;
950
951 // memory for storing temporary variables (mutable as they shall
952 // just be effectless auxiliary variables of the const apply*(...)
953 // methods)
954 mutable unsigned int nIterations_;
955
956 // constants for printing verbosity information
957 const std::string title_;
958 const std::string blank_;
959 };
960
961} // namespace Dune
962
963
964#endif // HAVE_ARPACKPP
965
966#endif // DUNE_ISTL_EIGENVALUE_ARPACKPP_HH
This file implements a vector space as a tensor product of a given vector space. The number of compon...
B::field_type field_type
export the type representing the field
Definition: bcrsmatrix.hh:445
@ blocklevel
The number of blocklevels the matrix contains.
Definition: bcrsmatrix.hh:465
A vector of blocks with memory management.
Definition: bvector.hh:313
B::field_type field_type
export the type representing the field
Definition: bvector.hh:319
vector space out of a tensor product of fields.
Definition: fvector.hh:93
derive error class from the base class in common
Definition: istlexception.hh:16
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 printvector(std::ostream &s, const V &v, std::string title, std::string rowtext, int columns=1, int width=10, int precision=2)
Print an ISTL vector.
Definition: io.hh:102
Some generic functions for pretty printing vectors and matrices.
Dune namespace.
Definition: alignment.hh:11
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 13, 23:29, 2024)