3 #ifndef DUNE_ISTL_EIGENVALUE_ARPACKPP_HH
4 #define DUNE_ISTL_EIGENVALUE_ARPACKPP_HH
6 #if HAVE_ARPACKPP || defined DOXYGEN
13 #include <dune/common/fvector.hh>
14 #include <dune/common/exceptions.hh>
53 template <
class BCRSMatrix>
54 class ArPackPlusPlus_BCRSMatrixWrapper
62 ArPackPlusPlus_BCRSMatrixWrapper (
const BCRSMatrix& A)
64 m_(A_.M() * mBlock), n_(A_.N() * nBlock)
68 (blockLevel<BCRSMatrix>() == 2,
69 "Only BCRSMatrices with blocklevel 2 are supported.");
73 domainBlockVector.resize(A_.N());
74 rangeBlockVector.resize(A_.M());
78 inline void multMv (Real* v, Real* w)
81 arrayToDomainBlockVector(v,domainBlockVector);
84 A_.mv(domainBlockVector,rangeBlockVector);
87 rangeBlockVectorToArray(rangeBlockVector,w);
91 inline void multMtMv (Real* v, Real* w)
94 arrayToDomainBlockVector(v,domainBlockVector);
97 A_.mv(domainBlockVector,rangeBlockVector);
98 A_.mtv(rangeBlockVector,domainBlockVector);
101 domainBlockVectorToArray(domainBlockVector,w);
105 inline void multMMtv (Real* v, Real* w)
108 arrayToRangeBlockVector(v,rangeBlockVector);
111 A_.mtv(rangeBlockVector,domainBlockVector);
112 A_.mv(domainBlockVector,rangeBlockVector);
115 rangeBlockVectorToArray(rangeBlockVector,w);
119 inline int nrows ()
const {
return m_; }
122 inline int ncols ()
const {
return n_; }
126 constexpr
static int mBlock = BCRSMatrix::block_type::rows;
127 constexpr
static int nBlock = BCRSMatrix::block_type::cols;
131 constexpr
static int dbvBlockSize = nBlock;
132 typedef Dune::FieldVector<Real,dbvBlockSize> DomainBlockVectorBlock;
137 constexpr
static int rbvBlockSize = mBlock;
138 typedef Dune::FieldVector<Real,rbvBlockSize> RangeBlockVectorBlock;
144 typedef typename DomainBlockVectorBlock::size_type dbvb_size_type;
145 typedef typename RangeBlockVectorBlock::size_type rbvb_size_type;
150 domainBlockVectorToArray (
const DomainBlockVector& dbv, Real* v)
152 for (dbv_size_type block = 0; block < dbv.N(); ++block)
153 for (dbvb_size_type iBlock = 0; iBlock < dbvBlockSize; ++iBlock)
154 v[block*dbvBlockSize + iBlock] = dbv[block][iBlock];
160 rangeBlockVectorToArray (
const RangeBlockVector& rbv, Real* v)
162 for (rbv_size_type block = 0; block < rbv.N(); ++block)
163 for (rbvb_size_type iBlock = 0; iBlock < rbvBlockSize; ++iBlock)
164 v[block*rbvBlockSize + iBlock] = rbv[block][iBlock];
170 static inline void arrayToDomainBlockVector (
const Real* v,
171 DomainBlockVector& dbv)
173 for (dbv_size_type block = 0; block < dbv.N(); ++block)
174 for (dbvb_size_type iBlock = 0; iBlock < dbvBlockSize; ++iBlock)
175 dbv[block][iBlock] = v[block*dbvBlockSize + iBlock];
180 static inline void arrayToRangeBlockVector (
const Real* v,
181 RangeBlockVector& rbv)
183 for (rbv_size_type block = 0; block < rbv.N(); ++block)
184 for (rbvb_size_type iBlock = 0; iBlock < rbvBlockSize; ++iBlock)
185 rbv[block][iBlock] = v[block*rbvBlockSize + iBlock];
190 const BCRSMatrix& A_;
197 mutable DomainBlockVector domainBlockVector;
198 mutable RangeBlockVector rangeBlockVector;
241 template <
typename BCRSMatrix,
typename BlockVector>
267 const unsigned int nIterationsMax = 100000,
268 const unsigned int verbosity_level = 0)
272 title_(
" ArPackPlusPlus_Algorithms: "),
292 std::cout <<
title_ <<
"Computing an approximation of "
293 <<
"the dominant eigenvalue of a matrix which "
294 <<
"is assumed to be symmetric." << std::endl;
298 typedef Impl::ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
302 const int nrows = A.nrows();
303 const int ncols = A.ncols();
308 << nrows <<
"x" << ncols <<
").");
312 int ncv = std::min(20, nrows);
313 const Real tol = epsilon;
316 const bool ivec =
true;
321 ARSymStdEig<Real,WrappedMatrix>
322 dprob(nrows, nev, &A, &WrappedMatrix::multMv, which, ncv, tol, maxit);
328 nconv = dprob.Eigenvalues(ev,ivec);
334 Real* x_raw = dprob.RawEigenvector(nev-1);
335 WrappedMatrix::arrayToDomainBlockVector(x_raw,x);
343 A.multMv(x_raw,Ax_raw);
344 WrappedMatrix::arrayToDomainBlockVector(Ax_raw,r);
346 const Real r_norm = r.two_norm();
354 std::cout <<
blank_ <<
"Obtained eigenvalues of A by solving "
355 <<
"A*x = λ*x using the ARPACK++ class ARSym"
356 <<
"StdEig:" << std::endl;
357 std::cout <<
blank_ <<
" converged eigenvalues of A: "
358 << nconv <<
" / " << nev << std::endl;
359 std::cout <<
blank_ <<
" dominant eigenvalue of A: "
360 << lambda << std::endl;
362 std::cout <<
blank_ <<
"Result ("
364 <<
"║residual║_2 = " << r_norm <<
"): "
365 <<
"λ = " << lambda << std::endl;
394 std::cout <<
title_ <<
"Computing an approximation of the "
395 <<
"least dominant eigenvalue of a matrix which "
396 <<
"is assumed to be symmetric." << std::endl;
400 typedef Impl::ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
404 const int nrows = A.nrows();
405 const int ncols = A.ncols();
410 << nrows <<
"x" << ncols <<
").");
414 int ncv = std::min(20, nrows);
415 const Real tol = epsilon;
418 const bool ivec =
true;
423 ARSymStdEig<Real,WrappedMatrix>
424 dprob(nrows, nev, &A, &WrappedMatrix::multMv, which, ncv, tol, maxit);
430 nconv = dprob.Eigenvalues(ev,ivec);
436 Real* x_raw = dprob.RawEigenvector(nev-1);
437 WrappedMatrix::arrayToDomainBlockVector(x_raw,x);
445 A.multMv(x_raw,Ax_raw);
446 WrappedMatrix::arrayToDomainBlockVector(Ax_raw,r);
448 const Real r_norm = r.two_norm();
456 std::cout <<
blank_ <<
"Obtained eigenvalues of A by solving "
457 <<
"A*x = λ*x using the ARPACK++ class ARSym"
458 <<
"StdEig:" << std::endl;
459 std::cout <<
blank_ <<
" converged eigenvalues of A: "
460 << nconv <<
" / " << nev << std::endl;
461 std::cout <<
blank_ <<
" least dominant eigenvalue of A: "
462 << lambda << std::endl;
464 std::cout <<
blank_ <<
"Result ("
466 <<
"║residual║_2 = " << r_norm <<
"): "
467 <<
"λ = " << lambda << std::endl;
495 std::cout <<
title_ <<
"Computing an approximation of the "
496 <<
"spectral condition number of a matrix which "
497 <<
"is assumed to be symmetric." << std::endl;
501 typedef Impl::ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
505 const int nrows = A.nrows();
506 const int ncols = A.ncols();
511 << nrows <<
"x" << ncols <<
").");
515 int ncv = std::min(20, nrows);
516 const Real tol = epsilon;
519 const bool ivec =
true;
524 ARSymStdEig<Real,WrappedMatrix>
525 dprob(nrows, nev, &A, &WrappedMatrix::multMv, which, ncv, tol, maxit);
531 nconv = dprob.Eigenvalues(ev,ivec);
534 const Real& lambda_max = ev[nev-1];
535 const Real& lambda_min = ev[0];
538 Real* x_max_raw = dprob.RawEigenvector(nev-1);
539 Real* x_min_raw = dprob.RawEigenvector(0);
542 cond_2 = std::abs(lambda_max / lambda_min);
550 A.multMv(x_max_raw,Ax_max_raw);
551 A.multMv(x_min_raw,Ax_min_raw);
552 Real r_max_norm = 0.0;
553 Real r_min_norm = 0.0;
554 for (
int i = 0; i < nrows; ++i)
556 r_max_norm += std::pow(Ax_max_raw[i] - lambda_max * x_max_raw[i],2);
557 r_min_norm += std::pow(Ax_min_raw[i] - lambda_min * x_min_raw[i],2);
559 r_max_norm = std::sqrt(r_max_norm);
560 r_min_norm = std::sqrt(r_min_norm);
568 std::cout <<
blank_ <<
"Obtained eigenvalues of A by solving "
569 <<
"A*x = λ*x using the ARPACK++ class ARSym"
570 <<
"StdEig:" << std::endl;
571 std::cout <<
blank_ <<
" converged eigenvalues of A: "
572 << nconv <<
" / " << nev << std::endl;
573 std::cout <<
blank_ <<
" dominant eigenvalue of A: "
574 << lambda_max << std::endl;
575 std::cout <<
blank_ <<
" least dominant eigenvalue of A: "
576 << lambda_min << std::endl;
577 std::cout <<
blank_ <<
" spectral condition number of A: "
578 << cond_2 << std::endl;
580 std::cout <<
blank_ <<
"Result ("
582 <<
"║residual║_2 = {" << r_max_norm <<
","
583 << r_min_norm <<
"}, " <<
"λ = {"
584 << lambda_max <<
"," << lambda_min
585 <<
"}): cond_2 = " << cond_2 << std::endl;
612 std::cout <<
title_ <<
"Computing an approximation of the "
613 <<
"largest singular value of a matrix which "
614 <<
"is assumed to be nonsymmetric." << std::endl;
618 typedef Impl::ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
622 const int nrows = A.nrows();
623 const int ncols = A.ncols();
628 <<
"columns (" << nrows <<
"x" << ncols <<
")."
629 <<
" This case is not implemented, yet.");
633 int ncv = std::min(20, nrows);
634 const Real tol = epsilon;
637 const bool ivec =
true;
642 ARSymStdEig<Real,WrappedMatrix>
643 dprob(ncols, nev, &A, &WrappedMatrix::multMtMv, which, ncv, tol, maxit);
649 nconv = dprob.Eigenvalues(ev,ivec);
652 const Real& lambda = ev[nev-1];
655 Real* x_raw = dprob.RawEigenvector(nev-1);
656 WrappedMatrix::arrayToDomainBlockVector(x_raw,x);
664 A.multMtMv(x_raw,AtAx_raw);
665 WrappedMatrix::arrayToDomainBlockVector(AtAx_raw,r);
667 const Real r_norm = r.two_norm();
671 sigma = std::sqrt(lambda);
679 std::cout <<
blank_ <<
"Obtained singular values of A by sol"
680 <<
"ving (A^T*A)*x = σ²*x using the ARPACK++ "
681 <<
"class ARSymStdEig:" << std::endl;
682 std::cout <<
blank_ <<
" converged eigenvalues of A^T*A: "
683 << nconv <<
" / " << nev << std::endl;
684 std::cout <<
blank_ <<
" largest eigenvalue of A^T*A: "
685 << lambda << std::endl;
686 std::cout <<
blank_ <<
" => largest singular value of A: "
687 << sigma << std::endl;
689 std::cout <<
blank_ <<
"Result ("
691 <<
"║residual║_2 = " << r_norm <<
"): "
692 <<
"σ = " << sigma << std::endl;
724 std::cout <<
title_ <<
"Computing an approximation of the "
725 <<
"smallest singular value of a matrix which "
726 <<
"is assumed to be nonsymmetric." << std::endl;
730 typedef Impl::ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
734 const int nrows = A.nrows();
735 const int ncols = A.ncols();
740 <<
"columns (" << nrows <<
"x" << ncols <<
")."
741 <<
" This case is not implemented, yet.");
745 int ncv = std::min(20, nrows);
746 const Real tol = epsilon;
749 const bool ivec =
true;
754 ARSymStdEig<Real,WrappedMatrix>
755 dprob(ncols, nev, &A, &WrappedMatrix::multMtMv, which, ncv, tol, maxit);
761 nconv = dprob.Eigenvalues(ev,ivec);
764 const Real& lambda = ev[nev-1];
767 Real* x_raw = dprob.RawEigenvector(nev-1);
768 WrappedMatrix::arrayToDomainBlockVector(x_raw,x);
776 A.multMtMv(x_raw,AtAx_raw);
777 WrappedMatrix::arrayToDomainBlockVector(AtAx_raw,r);
779 const Real r_norm = r.two_norm();
783 sigma = std::sqrt(lambda);
791 std::cout <<
blank_ <<
"Obtained singular values of A by sol"
792 <<
"ving (A^T*A)*x = σ²*x using the ARPACK++ "
793 <<
"class ARSymStdEig:" << std::endl;
794 std::cout <<
blank_ <<
" converged eigenvalues of A^T*A: "
795 << nconv <<
" / " << nev << std::endl;
796 std::cout <<
blank_ <<
" smallest eigenvalue of A^T*A: "
797 << lambda << std::endl;
798 std::cout <<
blank_ <<
" => smallest singular value of A: "
799 << sigma << std::endl;
801 std::cout <<
blank_ <<
"Result ("
803 <<
"║residual║_2 = " << r_norm <<
"): "
804 <<
"σ = " << sigma << std::endl;
832 std::cout <<
title_ <<
"Computing an approximation of the "
833 <<
"spectral condition number of a matrix which "
834 <<
"is assumed to be nonsymmetric." << std::endl;
838 typedef Impl::ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
842 const int nrows = A.nrows();
843 const int ncols = A.ncols();
848 <<
"columns (" << nrows <<
"x" << ncols <<
")."
849 <<
" This case is not implemented, yet.");
853 int ncv = std::min(20, nrows);
854 const Real tol = epsilon;
857 const bool ivec =
true;
862 ARSymStdEig<Real,WrappedMatrix>
863 dprob(ncols, nev, &A, &WrappedMatrix::multMtMv, which, ncv, tol, maxit);
869 nconv = dprob.Eigenvalues(ev,ivec);
872 const Real& lambda_max = ev[nev-1];
873 const Real& lambda_min = ev[0];
876 Real* x_max_raw = dprob.RawEigenvector(nev-1);
877 Real* x_min_raw = dprob.RawEigenvector(0);
883 Real* AtAx_max_raw =
new Real[ncols];
884 Real* AtAx_min_raw =
new Real[ncols];
885 A.multMtMv(x_max_raw,AtAx_max_raw);
886 A.multMtMv(x_min_raw,AtAx_min_raw);
887 Real r_max_norm = 0.0;
888 Real r_min_norm = 0.0;
889 for (
int i = 0; i < ncols; ++i)
891 r_max_norm += std::pow(AtAx_max_raw[i] - lambda_max * x_max_raw[i],2);
892 r_min_norm += std::pow(AtAx_min_raw[i] - lambda_min * x_min_raw[i],2);
894 r_max_norm = std::sqrt(r_max_norm);
895 r_min_norm = std::sqrt(r_min_norm);
898 const Real sigma_max = std::sqrt(lambda_max);
899 const Real sigma_min = std::sqrt(lambda_min);
902 cond_2 = sigma_max / sigma_min;
910 std::cout <<
blank_ <<
"Obtained singular values of A by sol"
911 <<
"ving (A^T*A)*x = σ²*x using the ARPACK++ "
912 <<
"class ARSymStdEig:" << std::endl;
913 std::cout <<
blank_ <<
" converged eigenvalues of A^T*A: "
914 << nconv <<
" / " << nev << std::endl;
915 std::cout <<
blank_ <<
" largest eigenvalue of A^T*A: "
916 << lambda_max << std::endl;
917 std::cout <<
blank_ <<
" smallest eigenvalue of A^T*A: "
918 << lambda_min << std::endl;
919 std::cout <<
blank_ <<
" => largest singular value of A: "
920 << sigma_max << std::endl;
921 std::cout <<
blank_ <<
" => smallest singular value of A: "
922 << sigma_min << std::endl;
924 std::cout <<
blank_ <<
"Result ("
926 <<
"║residual║_2 = {" << r_max_norm <<
","
927 << r_min_norm <<
"}, " <<
"σ = {"
928 << sigma_max <<
"," << sigma_min
929 <<
"}): cond_2 = " << cond_2 << std::endl;
933 delete[] AtAx_min_raw;
934 delete[] AtAx_max_raw;
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...
Some generic functions for pretty printing vectors and matrices.
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:86
Definition: allocator.hh:9
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:464
typename Imp::BlockTraits< B >::field_type field_type
export the type representing the field
Definition: bcrsmatrix.hh:486
A vector of blocks with memory management.
Definition: bvector.hh:393
typename Imp::BlockTraits< B >::field_type field_type
export the type representing the field
Definition: bvector.hh:399
A::size_type size_type
The type for the index access.
Definition: bvector.hh:408
Wrapper to use a range of ARPACK++ eigenvalue solvers.
Definition: arpackpp.hh:243
const unsigned int verbosity_level_
Definition: arpackpp.hh:956
unsigned int getIterationCount() const
Return the number of iterations in last application of an algorithm.
Definition: arpackpp.hh:942
const std::string title_
Definition: arpackpp.hh:964
BlockVector::field_type Real
Definition: arpackpp.hh:245
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:287
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:389
const BCRSMatrix & m_
Definition: arpackpp.hh:952
const unsigned int nIterationsMax_
Definition: arpackpp.hh:953
const std::string blank_
Definition: arpackpp.hh:965
ArPackPlusPlus_Algorithms(const BCRSMatrix &m, const unsigned int nIterationsMax=100000, const unsigned int verbosity_level=0)
Construct from required parameters.
Definition: arpackpp.hh:266
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:491
unsigned int nIterations_
Definition: arpackpp.hh:961
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:719
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:828
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:607
derive error class from the base class in common
Definition: istlexception.hh:17