5#ifndef DUNE_ISTL_EIGENVALUE_ARPACKPP_HH
6#define DUNE_ISTL_EIGENVALUE_ARPACKPP_HH
8#if HAVE_ARPACKPP || defined DOXYGEN
15#include <dune/common/fvector.hh>
16#include <dune/common/exceptions.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;
134 typedef Dune::FieldVector<Real,dbvBlockSize> DomainBlockVectorBlock;
139 constexpr static int rbvBlockSize = mBlock;
140 typedef Dune::FieldVector<Real,rbvBlockSize> RangeBlockVectorBlock;
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)
274 title_(
" ArPackPlusPlus_Algorithms: "),
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;
318 const bool ivec =
true;
323 ARSymStdEig<Real,WrappedMatrix>
324 dprob(nrows, nev, &A, &WrappedMatrix::multMv, which, ncv, tol, maxit);
330 nconv = dprob.Eigenvalues(ev,ivec);
336 Real* x_raw = dprob.RawEigenvector(nev-1);
337 WrappedMatrix::arrayToDomainBlockVector(x_raw,x);
345 A.multMv(x_raw,Ax_raw);
346 WrappedMatrix::arrayToDomainBlockVector(Ax_raw,r);
348 const Real r_norm = r.two_norm();
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 ("
366 <<
"║residual║_2 = " << r_norm <<
"): "
367 <<
"λ = " << lambda << std::endl;
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;
420 const bool ivec =
true;
425 ARSymStdEig<Real,WrappedMatrix>
426 dprob(nrows, nev, &A, &WrappedMatrix::multMv, which, ncv, tol, maxit);
432 nconv = dprob.Eigenvalues(ev,ivec);
438 Real* x_raw = dprob.RawEigenvector(nev-1);
439 WrappedMatrix::arrayToDomainBlockVector(x_raw,x);
447 A.multMv(x_raw,Ax_raw);
448 WrappedMatrix::arrayToDomainBlockVector(Ax_raw,r);
450 const Real r_norm = r.two_norm();
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 ("
468 <<
"║residual║_2 = " << r_norm <<
"): "
469 <<
"λ = " << lambda << std::endl;
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;
521 const bool ivec =
true;
526 ARSymStdEig<Real,WrappedMatrix>
527 dprob(nrows, nev, &A, &WrappedMatrix::multMv, which, ncv, tol, maxit);
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);
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);
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 ("
584 <<
"║residual║_2 = {" << r_max_norm <<
","
585 << r_min_norm <<
"}, " <<
"λ = {"
586 << lambda_max <<
"," << lambda_min
587 <<
"}): cond_2 = " << cond_2 << std::endl;
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;
639 const bool ivec =
true;
644 ARSymStdEig<Real,WrappedMatrix>
645 dprob(ncols, nev, &A, &WrappedMatrix::multMtMv, which, ncv, tol, maxit);
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);
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);
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 ("
693 <<
"║residual║_2 = " << r_norm <<
"): "
694 <<
"σ = " << sigma << std::endl;
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;
751 const bool ivec =
true;
756 ARSymStdEig<Real,WrappedMatrix>
757 dprob(ncols, nev, &A, &WrappedMatrix::multMtMv, which, ncv, tol, maxit);
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);
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);
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 ("
805 <<
"║residual║_2 = " << r_norm <<
"): "
806 <<
"σ = " << sigma << std::endl;
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;
859 const bool ivec =
true;
864 ARSymStdEig<Real,WrappedMatrix>
865 dprob(ncols, nev, &A, &WrappedMatrix::multMtMv, which, ncv, tol, maxit);
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);
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;
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 ("
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;
Helper functions for determining the vector/matrix block level.
Some generic functions for pretty printing vectors and matrices.
This file implements a vector space as a tensor product of a given vector space. The number of compon...
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
Definition allocator.hh:11
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
Wrapper to use a range of ARPACK++ eigenvalue solvers.
Definition arpackpp.hh:245
const unsigned int verbosity_level_
Definition arpackpp.hh:958
unsigned int getIterationCount() const
Return the number of iterations in last application of an algorithm.
Definition arpackpp.hh:944
const std::string title_
Definition arpackpp.hh:966
BlockVector::field_type Real
Definition arpackpp.hh:247
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
const BCRSMatrix & m_
Definition arpackpp.hh:954
const unsigned int nIterationsMax_
Definition arpackpp.hh:955
const std::string blank_
Definition arpackpp.hh:967
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
unsigned int nIterations_
Definition arpackpp.hh:963
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
derive error class from the base class in common
Definition istlexception.hh:19