5#ifndef DUNE_ISTL_EIGENVALUE_ARPACKPP_HH
6#define DUNE_ISTL_EIGENVALUE_ARPACKPP_HH
8#if HAVE_ARPACKPP || defined DOXYGEN
20#include <dune/istl/istlexception.hh>
55 template <
class BCRSMatrix>
56 class ArPackPlusPlus_BCRSMatrixWrapper
64 ArPackPlusPlus_BCRSMatrixWrapper (
const BCRSMatrix& A)
66 m_(A_.M() * mBlock), n_(A_.N() * nBlock)
70 (blockLevel<BCRSMatrix>() == 2,
71 "Only BCRSMatrices with blocklevel 2 are supported.");
75 domainBlockVector.resize(A_.N());
76 rangeBlockVector.resize(A_.M());
80 inline void multMv (Real* v, Real* w)
83 arrayToDomainBlockVector(v,domainBlockVector);
86 A_.mv(domainBlockVector,rangeBlockVector);
89 rangeBlockVectorToArray(rangeBlockVector,w);
93 inline void multMtMv (Real* v, Real* w)
96 arrayToDomainBlockVector(v,domainBlockVector);
99 A_.mv(domainBlockVector,rangeBlockVector);
100 A_.mtv(rangeBlockVector,domainBlockVector);
103 domainBlockVectorToArray(domainBlockVector,w);
107 inline void multMMtv (Real* v, Real* w)
110 arrayToRangeBlockVector(v,rangeBlockVector);
113 A_.mtv(rangeBlockVector,domainBlockVector);
114 A_.mv(domainBlockVector,rangeBlockVector);
117 rangeBlockVectorToArray(rangeBlockVector,w);
121 inline int nrows ()
const {
return m_; }
124 inline int ncols ()
const {
return n_; }
128 constexpr static int mBlock = BCRSMatrix::block_type::rows;
129 constexpr static int nBlock = BCRSMatrix::block_type::cols;
133 constexpr static int dbvBlockSize = nBlock;
139 constexpr static int rbvBlockSize = mBlock;
146 typedef typename DomainBlockVectorBlock::size_type dbvb_size_type;
147 typedef typename RangeBlockVectorBlock::size_type rbvb_size_type;
152 domainBlockVectorToArray (
const DomainBlockVector& dbv, Real* v)
154 for (dbv_size_type block = 0; block < dbv.N(); ++block)
155 for (dbvb_size_type iBlock = 0; iBlock < dbvBlockSize; ++iBlock)
156 v[block*dbvBlockSize + iBlock] = dbv[block][iBlock];
162 rangeBlockVectorToArray (
const RangeBlockVector& rbv, Real* v)
164 for (rbv_size_type block = 0; block < rbv.N(); ++block)
165 for (rbvb_size_type iBlock = 0; iBlock < rbvBlockSize; ++iBlock)
166 v[block*rbvBlockSize + iBlock] = rbv[block][iBlock];
172 static inline void arrayToDomainBlockVector (
const Real* v,
173 DomainBlockVector& dbv)
175 for (dbv_size_type block = 0; block < dbv.N(); ++block)
176 for (dbvb_size_type iBlock = 0; iBlock < dbvBlockSize; ++iBlock)
177 dbv[block][iBlock] = v[block*dbvBlockSize + iBlock];
182 static inline void arrayToRangeBlockVector (
const Real* v,
183 RangeBlockVector& rbv)
185 for (rbv_size_type block = 0; block < rbv.N(); ++block)
186 for (rbvb_size_type iBlock = 0; iBlock < rbvBlockSize; ++iBlock)
187 rbv[block][iBlock] = v[block*rbvBlockSize + iBlock];
192 const BCRSMatrix& A_;
199 mutable DomainBlockVector domainBlockVector;
200 mutable RangeBlockVector rangeBlockVector;
243 template <
typename BCRSMatrix,
typename BlockVector>
269 const unsigned int nIterationsMax = 100000,
270 const unsigned int verbosity_level = 0)
271 : m_(m), nIterationsMax_(nIterationsMax),
272 verbosity_level_(verbosity_level),
274 title_(
" ArPackPlusPlus_Algorithms: "),
275 blank_(title_.length(),
' ')
293 if (verbosity_level_ > 0)
294 std::cout << title_ <<
"Computing an approximation of "
295 <<
"the dominant eigenvalue of a matrix which "
296 <<
"is assumed to be symmetric." << std::endl;
300 typedef Impl::ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
304 const int nrows = A.nrows();
305 const int ncols = A.ncols();
310 << nrows <<
"x" << ncols <<
").");
314 int ncv = std::min(20, nrows);
315 const Real tol = epsilon;
316 const int maxit = nIterationsMax_*nev;
317 Real* ev =
new Real[nev];
318 const bool ivec =
true;
323 ARSymStdEig<Real,WrappedMatrix>
324 dprob(nrows, nev, &A, &WrappedMatrix::multMv, which, ncv, tol, maxit);
327 if (verbosity_level_ > 3) dprob.Trace();
330 nconv = dprob.Eigenvalues(ev,ivec);
336 Real* x_raw = dprob.RawEigenvector(nev-1);
337 WrappedMatrix::arrayToDomainBlockVector(x_raw,x);
340 nIterations_ = dprob.GetIter();
344 Real* Ax_raw =
new Real[nrows];
345 A.multMv(x_raw,Ax_raw);
346 WrappedMatrix::arrayToDomainBlockVector(Ax_raw,r);
348 const Real r_norm = r.two_norm();
351 if (verbosity_level_ > 0)
353 if (verbosity_level_ > 1)
356 std::cout << blank_ <<
"Obtained eigenvalues of A by solving "
357 <<
"A*x = λ*x using the ARPACK++ class ARSym"
358 <<
"StdEig:" << std::endl;
359 std::cout << blank_ <<
" converged eigenvalues of A: "
360 << nconv <<
" / " << nev << std::endl;
361 std::cout << blank_ <<
" dominant eigenvalue of A: "
362 << lambda << std::endl;
364 std::cout << blank_ <<
"Result ("
365 <<
"#iterations = " << nIterations_ <<
", "
366 <<
"║residual║_2 = " << r_norm <<
"): "
367 <<
"λ = " << lambda << std::endl;
368 if (verbosity_level_ > 2)
395 if (verbosity_level_ > 0)
396 std::cout << title_ <<
"Computing an approximation of the "
397 <<
"least dominant eigenvalue of a matrix which "
398 <<
"is assumed to be symmetric." << std::endl;
402 typedef Impl::ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
406 const int nrows = A.nrows();
407 const int ncols = A.ncols();
412 << nrows <<
"x" << ncols <<
").");
416 int ncv = std::min(20, nrows);
417 const Real tol = epsilon;
418 const int maxit = nIterationsMax_*nev;
419 Real* ev =
new Real[nev];
420 const bool ivec =
true;
425 ARSymStdEig<Real,WrappedMatrix>
426 dprob(nrows, nev, &A, &WrappedMatrix::multMv, which, ncv, tol, maxit);
429 if (verbosity_level_ > 3) dprob.Trace();
432 nconv = dprob.Eigenvalues(ev,ivec);
438 Real* x_raw = dprob.RawEigenvector(nev-1);
439 WrappedMatrix::arrayToDomainBlockVector(x_raw,x);
442 nIterations_ = dprob.GetIter();
446 Real* Ax_raw =
new Real[nrows];
447 A.multMv(x_raw,Ax_raw);
448 WrappedMatrix::arrayToDomainBlockVector(Ax_raw,r);
450 const Real r_norm = r.two_norm();
453 if (verbosity_level_ > 0)
455 if (verbosity_level_ > 1)
458 std::cout << blank_ <<
"Obtained eigenvalues of A by solving "
459 <<
"A*x = λ*x using the ARPACK++ class ARSym"
460 <<
"StdEig:" << std::endl;
461 std::cout << blank_ <<
" converged eigenvalues of A: "
462 << nconv <<
" / " << nev << std::endl;
463 std::cout << blank_ <<
" least dominant eigenvalue of A: "
464 << lambda << std::endl;
466 std::cout << blank_ <<
"Result ("
467 <<
"#iterations = " << nIterations_ <<
", "
468 <<
"║residual║_2 = " << r_norm <<
"): "
469 <<
"λ = " << lambda << std::endl;
470 if (verbosity_level_ > 2)
496 if (verbosity_level_ > 0)
497 std::cout << title_ <<
"Computing an approximation of the "
498 <<
"spectral condition number of a matrix which "
499 <<
"is assumed to be symmetric." << std::endl;
503 typedef Impl::ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
507 const int nrows = A.nrows();
508 const int ncols = A.ncols();
513 << nrows <<
"x" << ncols <<
").");
517 int ncv = std::min(20, nrows);
518 const Real tol = epsilon;
519 const int maxit = nIterationsMax_*nev;
520 Real* ev =
new Real[nev];
521 const bool ivec =
true;
526 ARSymStdEig<Real,WrappedMatrix>
527 dprob(nrows, nev, &A, &WrappedMatrix::multMv, which, ncv, tol, maxit);
530 if (verbosity_level_ > 3) dprob.Trace();
533 nconv = dprob.Eigenvalues(ev,ivec);
536 const Real& lambda_max = ev[nev-1];
537 const Real& lambda_min = ev[0];
540 Real* x_max_raw = dprob.RawEigenvector(nev-1);
541 Real* x_min_raw = dprob.RawEigenvector(0);
544 cond_2 = std::abs(lambda_max / lambda_min);
547 nIterations_ = dprob.GetIter();
550 Real* Ax_max_raw =
new Real[nrows];
551 Real* Ax_min_raw =
new Real[nrows];
552 A.multMv(x_max_raw,Ax_max_raw);
553 A.multMv(x_min_raw,Ax_min_raw);
554 Real r_max_norm = 0.0;
555 Real r_min_norm = 0.0;
556 for (
int i = 0; i < nrows; ++i)
558 r_max_norm += std::pow(Ax_max_raw[i] - lambda_max * x_max_raw[i],2);
559 r_min_norm += std::pow(Ax_min_raw[i] - lambda_min * x_min_raw[i],2);
561 r_max_norm = std::sqrt(r_max_norm);
562 r_min_norm = std::sqrt(r_min_norm);
565 if (verbosity_level_ > 0)
567 if (verbosity_level_ > 1)
570 std::cout << blank_ <<
"Obtained eigenvalues of A by solving "
571 <<
"A*x = λ*x using the ARPACK++ class ARSym"
572 <<
"StdEig:" << std::endl;
573 std::cout << blank_ <<
" converged eigenvalues of A: "
574 << nconv <<
" / " << nev << std::endl;
575 std::cout << blank_ <<
" dominant eigenvalue of A: "
576 << lambda_max << std::endl;
577 std::cout << blank_ <<
" least dominant eigenvalue of A: "
578 << lambda_min << std::endl;
579 std::cout << blank_ <<
" spectral condition number of A: "
580 << cond_2 << std::endl;
582 std::cout << blank_ <<
"Result ("
583 <<
"#iterations = " << nIterations_ <<
", "
584 <<
"║residual║_2 = {" << r_max_norm <<
","
585 << r_min_norm <<
"}, " <<
"λ = {"
586 << lambda_max <<
"," << lambda_min
587 <<
"}): cond_2 = " << cond_2 << std::endl;
613 if (verbosity_level_ > 0)
614 std::cout << title_ <<
"Computing an approximation of the "
615 <<
"largest singular value of a matrix which "
616 <<
"is assumed to be nonsymmetric." << std::endl;
620 typedef Impl::ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
624 const int nrows = A.nrows();
625 const int ncols = A.ncols();
630 <<
"columns (" << nrows <<
"x" << ncols <<
")."
631 <<
" This case is not implemented, yet.");
635 int ncv = std::min(20, nrows);
636 const Real tol = epsilon;
637 const int maxit = nIterationsMax_*nev;
638 Real* ev =
new Real[nev];
639 const bool ivec =
true;
644 ARSymStdEig<Real,WrappedMatrix>
645 dprob(ncols, nev, &A, &WrappedMatrix::multMtMv, which, ncv, tol, maxit);
648 if (verbosity_level_ > 3) dprob.Trace();
651 nconv = dprob.Eigenvalues(ev,ivec);
654 const Real& lambda = ev[nev-1];
657 Real* x_raw = dprob.RawEigenvector(nev-1);
658 WrappedMatrix::arrayToDomainBlockVector(x_raw,x);
661 nIterations_ = dprob.GetIter();
665 Real* AtAx_raw =
new Real[ncols];
666 A.multMtMv(x_raw,AtAx_raw);
667 WrappedMatrix::arrayToDomainBlockVector(AtAx_raw,r);
669 const Real r_norm = r.two_norm();
673 sigma = std::sqrt(lambda);
676 if (verbosity_level_ > 0)
678 if (verbosity_level_ > 1)
681 std::cout << blank_ <<
"Obtained singular values of A by sol"
682 <<
"ving (A^T*A)*x = σ²*x using the ARPACK++ "
683 <<
"class ARSymStdEig:" << std::endl;
684 std::cout << blank_ <<
" converged eigenvalues of A^T*A: "
685 << nconv <<
" / " << nev << std::endl;
686 std::cout << blank_ <<
" largest eigenvalue of A^T*A: "
687 << lambda << std::endl;
688 std::cout << blank_ <<
" => largest singular value of A: "
689 << sigma << std::endl;
691 std::cout << blank_ <<
"Result ("
692 <<
"#iterations = " << nIterations_ <<
", "
693 <<
"║residual║_2 = " << r_norm <<
"): "
694 <<
"σ = " << sigma << std::endl;
695 if (verbosity_level_ > 2)
725 if (verbosity_level_ > 0)
726 std::cout << title_ <<
"Computing an approximation of the "
727 <<
"smallest singular value of a matrix which "
728 <<
"is assumed to be nonsymmetric." << std::endl;
732 typedef Impl::ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
736 const int nrows = A.nrows();
737 const int ncols = A.ncols();
742 <<
"columns (" << nrows <<
"x" << ncols <<
")."
743 <<
" This case is not implemented, yet.");
747 int ncv = std::min(20, nrows);
748 const Real tol = epsilon;
749 const int maxit = nIterationsMax_*nev;
750 Real* ev =
new Real[nev];
751 const bool ivec =
true;
756 ARSymStdEig<Real,WrappedMatrix>
757 dprob(ncols, nev, &A, &WrappedMatrix::multMtMv, which, ncv, tol, maxit);
760 if (verbosity_level_ > 3) dprob.Trace();
763 nconv = dprob.Eigenvalues(ev,ivec);
766 const Real& lambda = ev[nev-1];
769 Real* x_raw = dprob.RawEigenvector(nev-1);
770 WrappedMatrix::arrayToDomainBlockVector(x_raw,x);
773 nIterations_ = dprob.GetIter();
777 Real* AtAx_raw =
new Real[ncols];
778 A.multMtMv(x_raw,AtAx_raw);
779 WrappedMatrix::arrayToDomainBlockVector(AtAx_raw,r);
781 const Real r_norm = r.two_norm();
785 sigma = std::sqrt(lambda);
788 if (verbosity_level_ > 0)
790 if (verbosity_level_ > 1)
793 std::cout << blank_ <<
"Obtained singular values of A by sol"
794 <<
"ving (A^T*A)*x = σ²*x using the ARPACK++ "
795 <<
"class ARSymStdEig:" << std::endl;
796 std::cout << blank_ <<
" converged eigenvalues of A^T*A: "
797 << nconv <<
" / " << nev << std::endl;
798 std::cout << blank_ <<
" smallest eigenvalue of A^T*A: "
799 << lambda << std::endl;
800 std::cout << blank_ <<
" => smallest singular value of A: "
801 << sigma << std::endl;
803 std::cout << blank_ <<
"Result ("
804 <<
"#iterations = " << nIterations_ <<
", "
805 <<
"║residual║_2 = " << r_norm <<
"): "
806 <<
"σ = " << sigma << std::endl;
807 if (verbosity_level_ > 2)
833 if (verbosity_level_ > 0)
834 std::cout << title_ <<
"Computing an approximation of the "
835 <<
"spectral condition number of a matrix which "
836 <<
"is assumed to be nonsymmetric." << std::endl;
840 typedef Impl::ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
844 const int nrows = A.nrows();
845 const int ncols = A.ncols();
850 <<
"columns (" << nrows <<
"x" << ncols <<
")."
851 <<
" This case is not implemented, yet.");
855 int ncv = std::min(20, nrows);
856 const Real tol = epsilon;
857 const int maxit = nIterationsMax_*nev;
858 Real* ev =
new Real[nev];
859 const bool ivec =
true;
864 ARSymStdEig<Real,WrappedMatrix>
865 dprob(ncols, nev, &A, &WrappedMatrix::multMtMv, which, ncv, tol, maxit);
868 if (verbosity_level_ > 3) dprob.Trace();
871 nconv = dprob.Eigenvalues(ev,ivec);
874 const Real& lambda_max = ev[nev-1];
875 const Real& lambda_min = ev[0];
878 Real* x_max_raw = dprob.RawEigenvector(nev-1);
879 Real* x_min_raw = dprob.RawEigenvector(0);
882 nIterations_ = dprob.GetIter();
885 Real* AtAx_max_raw =
new Real[ncols];
886 Real* AtAx_min_raw =
new Real[ncols];
887 A.multMtMv(x_max_raw,AtAx_max_raw);
888 A.multMtMv(x_min_raw,AtAx_min_raw);
889 Real r_max_norm = 0.0;
890 Real r_min_norm = 0.0;
891 for (
int i = 0; i < ncols; ++i)
893 r_max_norm += std::pow(AtAx_max_raw[i] - lambda_max * x_max_raw[i],2);
894 r_min_norm += std::pow(AtAx_min_raw[i] - lambda_min * x_min_raw[i],2);
896 r_max_norm = std::sqrt(r_max_norm);
897 r_min_norm = std::sqrt(r_min_norm);
900 const Real sigma_max = std::sqrt(lambda_max);
901 const Real sigma_min = std::sqrt(lambda_min);
904 cond_2 = sigma_max / sigma_min;
907 if (verbosity_level_ > 0)
909 if (verbosity_level_ > 1)
912 std::cout << blank_ <<
"Obtained singular values of A by sol"
913 <<
"ving (A^T*A)*x = σ²*x using the ARPACK++ "
914 <<
"class ARSymStdEig:" << std::endl;
915 std::cout << blank_ <<
" converged eigenvalues of A^T*A: "
916 << nconv <<
" / " << nev << std::endl;
917 std::cout << blank_ <<
" largest eigenvalue of A^T*A: "
918 << lambda_max << std::endl;
919 std::cout << blank_ <<
" smallest eigenvalue of A^T*A: "
920 << lambda_min << std::endl;
921 std::cout << blank_ <<
" => largest singular value of A: "
922 << sigma_max << std::endl;
923 std::cout << blank_ <<
" => smallest singular value of A: "
924 << sigma_min << std::endl;
926 std::cout << blank_ <<
"Result ("
927 <<
"#iterations = " << nIterations_ <<
", "
928 <<
"║residual║_2 = {" << r_max_norm <<
","
929 << r_min_norm <<
"}, " <<
"σ = {"
930 << sigma_max <<
"," << sigma_min
931 <<
"}): cond_2 = " << cond_2 << std::endl;
935 delete[] AtAx_min_raw;
936 delete[] AtAx_max_raw;
946 if (nIterations_ == 0)
955 const unsigned int nIterationsMax_;
958 const unsigned int verbosity_level_;
963 mutable unsigned int nIterations_;
966 const std::string title_;
967 const std::string blank_;
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
unsigned int getIterationCount() const
Return the number of iterations in last application of an algorithm.
Definition: arpackpp.hh:944
void computeSymMaxMagnitude(const Real &epsilon, BlockVector &x, Real &lambda) const
Assume the matrix to be square, symmetric and perform IRLM to compute an approximation lambda of its ...
Definition: arpackpp.hh:289
void computeSymMinMagnitude(const Real &epsilon, BlockVector &x, Real &lambda) const
Assume the matrix to be square, symmetric and perform IRLM to compute an approximation lambda of its ...
Definition: arpackpp.hh:391
ArPackPlusPlus_Algorithms(const BCRSMatrix &m, const unsigned int nIterationsMax=100000, const unsigned int verbosity_level=0)
Construct from required parameters.
Definition: arpackpp.hh:268
void computeSymCond2(const Real &epsilon, Real &cond_2) const
Assume the matrix to be square, symmetric and perform IRLM to compute an approximation of its spectra...
Definition: arpackpp.hh:493
void computeNonSymMin(const Real &epsilon, BlockVector &x, Real &sigma) const
Assume the matrix to be nonsymmetric and perform IRLM to compute an approximation sigma of its smalle...
Definition: arpackpp.hh:721
void computeNonSymCond2(const Real &epsilon, Real &cond_2) const
Assume the matrix to be nonsymmetric and perform IRLM to compute an approximation of its spectral con...
Definition: arpackpp.hh:830
void computeNonSymMax(const Real &epsilon, BlockVector &x, Real &sigma) const
Assume the matrix to be nonsymmetric and perform IRLM to compute an approximation sigma of its larges...
Definition: arpackpp.hh:609
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:466
typename Imp::BlockTraits< B >::field_type field_type
export the type representing the field
Definition: bcrsmatrix.hh:488
A vector of blocks with memory management.
Definition: bvector.hh:392
typename Imp::BlockTraits< B >::field_type field_type
export the type representing the field
Definition: bvector.hh:398
A::size_type size_type
The type for the index access.
Definition: bvector.hh:407
vector space out of a tensor product of fields.
Definition: fvector.hh:91
derive error class from the base class in common
Definition: istlexception.hh:19
A few common exception classes.
Some generic functions for pretty printing vectors and matrices.
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 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:89
Dune namespace.
Definition: alignedallocator.hh:13