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>
21 #undef Status // prevent preprocessor from damaging the ARPACK++
52 template <
class BCRSMatrix>
53 class ArPackPlusPlus_BCRSMatrixWrapper
61 ArPackPlusPlus_BCRSMatrixWrapper (
const BCRSMatrix& A)
63 m_(A_.M() * mBlock), n_(A_.N() * nBlock)
68 "Only BCRSMatrices with blocklevel 2 are supported.");
72 domainBlockVector.resize(A_.N());
73 rangeBlockVector.resize(A_.M());
77 inline void multMv (Real* v, Real* w)
80 arrayToDomainBlockVector(v,domainBlockVector);
83 A_.mv(domainBlockVector,rangeBlockVector);
86 rangeBlockVectorToArray(rangeBlockVector,w);
90 inline void multMtMv (Real* v, Real* w)
93 arrayToDomainBlockVector(v,domainBlockVector);
96 A_.mv(domainBlockVector,rangeBlockVector);
97 A_.mtv(rangeBlockVector,domainBlockVector);
100 domainBlockVectorToArray(domainBlockVector,w);
104 inline void multMMtv (Real* v, Real* w)
107 arrayToRangeBlockVector(v,rangeBlockVector);
110 A_.mtv(rangeBlockVector,domainBlockVector);
111 A_.mv(domainBlockVector,rangeBlockVector);
114 rangeBlockVectorToArray(rangeBlockVector,w);
118 inline int nrows ()
const {
return m_; }
121 inline int ncols ()
const {
return n_; }
125 constexpr
static int mBlock = BCRSMatrix::block_type::rows;
126 constexpr
static int nBlock = BCRSMatrix::block_type::cols;
130 constexpr
static int dbvBlockSize = nBlock;
131 typedef Dune::FieldVector<Real,dbvBlockSize> DomainBlockVectorBlock;
136 constexpr
static int rbvBlockSize = mBlock;
137 typedef Dune::FieldVector<Real,rbvBlockSize> RangeBlockVectorBlock;
143 typedef typename DomainBlockVectorBlock::size_type dbvb_size_type;
144 typedef typename RangeBlockVectorBlock::size_type rbvb_size_type;
149 domainBlockVectorToArray (
const DomainBlockVector& dbv, Real* v)
151 for (dbv_size_type block = 0; block < dbv.N(); ++block)
152 for (dbvb_size_type iBlock = 0; iBlock < dbvBlockSize; ++iBlock)
153 v[block*dbvBlockSize + iBlock] = dbv[block][iBlock];
159 rangeBlockVectorToArray (
const RangeBlockVector& rbv, Real* v)
161 for (rbv_size_type block = 0; block < rbv.N(); ++block)
162 for (rbvb_size_type iBlock = 0; iBlock < rbvBlockSize; ++iBlock)
163 v[block*rbvBlockSize + iBlock] = rbv[block][iBlock];
169 static inline void arrayToDomainBlockVector (
const Real* v,
170 DomainBlockVector& dbv)
172 for (dbv_size_type block = 0; block < dbv.N(); ++block)
173 for (dbvb_size_type iBlock = 0; iBlock < dbvBlockSize; ++iBlock)
174 dbv[block][iBlock] = v[block*dbvBlockSize + iBlock];
179 static inline void arrayToRangeBlockVector (
const Real* v,
180 RangeBlockVector& rbv)
182 for (rbv_size_type block = 0; block < rbv.N(); ++block)
183 for (rbvb_size_type iBlock = 0; iBlock < rbvBlockSize; ++iBlock)
184 rbv[block][iBlock] = v[block*rbvBlockSize + iBlock];
189 const BCRSMatrix& A_;
196 mutable DomainBlockVector domainBlockVector;
197 mutable RangeBlockVector rangeBlockVector;
240 template <
typename BCRSMatrix,
typename BlockVector>
266 const unsigned int nIterationsMax = 100000,
267 const unsigned int verbosity_level = 0)
271 title_(
" ArPackPlusPlus_Algorithms: "),
291 std::cout <<
title_ <<
"Computing an approximation of "
292 <<
"the dominant eigenvalue of a matrix which "
293 <<
"is assumed to be symmetric." << std::endl;
297 typedef Impl::ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
301 const int nrows = A.nrows();
302 const int ncols = A.ncols();
307 << nrows <<
"x" << ncols <<
").");
311 int ncv = std::min(20, nrows);
312 const Real tol = epsilon;
315 const bool ivec =
true;
320 ARSymStdEig<Real,WrappedMatrix>
321 dprob(nrows, nev, &A, &WrappedMatrix::multMv, which, ncv, tol, maxit);
327 nconv = dprob.Eigenvalues(ev,ivec);
333 Real* x_raw = dprob.RawEigenvector(nev-1);
334 WrappedMatrix::arrayToDomainBlockVector(x_raw,x);
342 A.multMv(x_raw,Ax_raw);
343 WrappedMatrix::arrayToDomainBlockVector(Ax_raw,r);
345 const Real r_norm = r.two_norm();
353 std::cout <<
blank_ <<
"Obtained eigenvalues of A by solving "
354 <<
"A*x = λ*x using the ARPACK++ class ARSym"
355 <<
"StdEig:" << std::endl;
356 std::cout <<
blank_ <<
" converged eigenvalues of A: "
357 << nconv <<
" / " << nev << std::endl;
358 std::cout <<
blank_ <<
" dominant eigenvalue of A: "
359 << lambda << std::endl;
361 std::cout <<
blank_ <<
"Result ("
363 <<
"║residual║_2 = " << r_norm <<
"): "
364 <<
"λ = " << lambda << std::endl;
393 std::cout <<
title_ <<
"Computing an approximation of the "
394 <<
"least dominant eigenvalue of a matrix which "
395 <<
"is assumed to be symmetric." << std::endl;
399 typedef Impl::ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
403 const int nrows = A.nrows();
404 const int ncols = A.ncols();
409 << nrows <<
"x" << ncols <<
").");
413 int ncv = std::min(20, nrows);
414 const Real tol = epsilon;
417 const bool ivec =
true;
422 ARSymStdEig<Real,WrappedMatrix>
423 dprob(nrows, nev, &A, &WrappedMatrix::multMv, which, ncv, tol, maxit);
429 nconv = dprob.Eigenvalues(ev,ivec);
435 Real* x_raw = dprob.RawEigenvector(nev-1);
436 WrappedMatrix::arrayToDomainBlockVector(x_raw,x);
444 A.multMv(x_raw,Ax_raw);
445 WrappedMatrix::arrayToDomainBlockVector(Ax_raw,r);
447 const Real r_norm = r.two_norm();
455 std::cout <<
blank_ <<
"Obtained eigenvalues of A by solving "
456 <<
"A*x = λ*x using the ARPACK++ class ARSym"
457 <<
"StdEig:" << std::endl;
458 std::cout <<
blank_ <<
" converged eigenvalues of A: "
459 << nconv <<
" / " << nev << std::endl;
460 std::cout <<
blank_ <<
" least dominant eigenvalue of A: "
461 << lambda << std::endl;
463 std::cout <<
blank_ <<
"Result ("
465 <<
"║residual║_2 = " << r_norm <<
"): "
466 <<
"λ = " << lambda << std::endl;
494 std::cout <<
title_ <<
"Computing an approximation of the "
495 <<
"spectral condition number of a matrix which "
496 <<
"is assumed to be symmetric." << std::endl;
500 typedef Impl::ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
504 const int nrows = A.nrows();
505 const int ncols = A.ncols();
510 << nrows <<
"x" << ncols <<
").");
514 int ncv = std::min(20, nrows);
515 const Real tol = epsilon;
518 const bool ivec =
true;
523 ARSymStdEig<Real,WrappedMatrix>
524 dprob(nrows, nev, &A, &WrappedMatrix::multMv, which, ncv, tol, maxit);
530 nconv = dprob.Eigenvalues(ev,ivec);
533 const Real& lambda_max = ev[nev-1];
534 const Real& lambda_min = ev[0];
537 Real* x_max_raw = dprob.RawEigenvector(nev-1);
538 Real* x_min_raw = dprob.RawEigenvector(0);
541 cond_2 = std::abs(lambda_max / lambda_min);
549 A.multMv(x_max_raw,Ax_max_raw);
550 A.multMv(x_min_raw,Ax_min_raw);
551 Real r_max_norm = 0.0;
552 Real r_min_norm = 0.0;
553 for (
int i = 0; i < nrows; ++i)
555 r_max_norm += std::pow(Ax_max_raw[i] - lambda_max * x_max_raw[i],2);
556 r_min_norm += std::pow(Ax_min_raw[i] - lambda_min * x_min_raw[i],2);
558 r_max_norm = std::sqrt(r_max_norm);
559 r_min_norm = std::sqrt(r_min_norm);
567 std::cout <<
blank_ <<
"Obtained eigenvalues of A by solving "
568 <<
"A*x = λ*x using the ARPACK++ class ARSym"
569 <<
"StdEig:" << std::endl;
570 std::cout <<
blank_ <<
" converged eigenvalues of A: "
571 << nconv <<
" / " << nev << std::endl;
572 std::cout <<
blank_ <<
" dominant eigenvalue of A: "
573 << lambda_max << std::endl;
574 std::cout <<
blank_ <<
" least dominant eigenvalue of A: "
575 << lambda_min << std::endl;
576 std::cout <<
blank_ <<
" spectral condition number of A: "
577 << cond_2 << std::endl;
579 std::cout <<
blank_ <<
"Result ("
581 <<
"║residual║_2 = {" << r_max_norm <<
","
582 << r_min_norm <<
"}, " <<
"λ = {"
583 << lambda_max <<
"," << lambda_min
584 <<
"}): cond_2 = " << cond_2 << std::endl;
611 std::cout <<
title_ <<
"Computing an approximation of the "
612 <<
"largest singular value of a matrix which "
613 <<
"is assumed to be nonsymmetric." << std::endl;
617 typedef Impl::ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
621 const int nrows = A.nrows();
622 const int ncols = A.ncols();
627 <<
"columns (" << nrows <<
"x" << ncols <<
")."
628 <<
" This case is not implemented, yet.");
632 int ncv = std::min(20, nrows);
633 const Real tol = epsilon;
636 const bool ivec =
true;
641 ARSymStdEig<Real,WrappedMatrix>
642 dprob(ncols, nev, &A, &WrappedMatrix::multMtMv, which, ncv, tol, maxit);
648 nconv = dprob.Eigenvalues(ev,ivec);
651 const Real& lambda = ev[nev-1];
654 Real* x_raw = dprob.RawEigenvector(nev-1);
655 WrappedMatrix::arrayToDomainBlockVector(x_raw,x);
663 A.multMtMv(x_raw,AtAx_raw);
664 WrappedMatrix::arrayToDomainBlockVector(AtAx_raw,r);
666 const Real r_norm = r.two_norm();
670 sigma = std::sqrt(lambda);
678 std::cout <<
blank_ <<
"Obtained singular values of A by sol"
679 <<
"ving (A^T*A)*x = σ²*x using the ARPACK++ "
680 <<
"class ARSymStdEig:" << std::endl;
681 std::cout <<
blank_ <<
" converged eigenvalues of A^T*A: "
682 << nconv <<
" / " << nev << std::endl;
683 std::cout <<
blank_ <<
" largest eigenvalue of A^T*A: "
684 << lambda << std::endl;
685 std::cout <<
blank_ <<
" => largest singular value of A: "
686 << sigma << std::endl;
688 std::cout <<
blank_ <<
"Result ("
690 <<
"║residual║_2 = " << r_norm <<
"): "
691 <<
"σ = " << sigma << std::endl;
723 std::cout <<
title_ <<
"Computing an approximation of the "
724 <<
"smallest singular value of a matrix which "
725 <<
"is assumed to be nonsymmetric." << std::endl;
729 typedef Impl::ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
733 const int nrows = A.nrows();
734 const int ncols = A.ncols();
739 <<
"columns (" << nrows <<
"x" << ncols <<
")."
740 <<
" This case is not implemented, yet.");
744 int ncv = std::min(20, nrows);
745 const Real tol = epsilon;
748 const bool ivec =
true;
753 ARSymStdEig<Real,WrappedMatrix>
754 dprob(ncols, nev, &A, &WrappedMatrix::multMtMv, which, ncv, tol, maxit);
760 nconv = dprob.Eigenvalues(ev,ivec);
763 const Real& lambda = ev[nev-1];
766 Real* x_raw = dprob.RawEigenvector(nev-1);
767 WrappedMatrix::arrayToDomainBlockVector(x_raw,x);
775 A.multMtMv(x_raw,AtAx_raw);
776 WrappedMatrix::arrayToDomainBlockVector(AtAx_raw,r);
778 const Real r_norm = r.two_norm();
782 sigma = std::sqrt(lambda);
790 std::cout <<
blank_ <<
"Obtained singular values of A by sol"
791 <<
"ving (A^T*A)*x = σ²*x using the ARPACK++ "
792 <<
"class ARSymStdEig:" << std::endl;
793 std::cout <<
blank_ <<
" converged eigenvalues of A^T*A: "
794 << nconv <<
" / " << nev << std::endl;
795 std::cout <<
blank_ <<
" smallest eigenvalue of A^T*A: "
796 << lambda << std::endl;
797 std::cout <<
blank_ <<
" => smallest singular value of A: "
798 << sigma << std::endl;
800 std::cout <<
blank_ <<
"Result ("
802 <<
"║residual║_2 = " << r_norm <<
"): "
803 <<
"σ = " << sigma << std::endl;
831 std::cout <<
title_ <<
"Computing an approximation of the "
832 <<
"spectral condition number of a matrix which "
833 <<
"is assumed to be nonsymmetric." << std::endl;
837 typedef Impl::ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
841 const int nrows = A.nrows();
842 const int ncols = A.ncols();
847 <<
"columns (" << nrows <<
"x" << ncols <<
")."
848 <<
" This case is not implemented, yet.");
852 int ncv = std::min(20, nrows);
853 const Real tol = epsilon;
856 const bool ivec =
true;
861 ARSymStdEig<Real,WrappedMatrix>
862 dprob(ncols, nev, &A, &WrappedMatrix::multMtMv, which, ncv, tol, maxit);
868 nconv = dprob.Eigenvalues(ev,ivec);
871 const Real& lambda_max = ev[nev-1];
872 const Real& lambda_min = ev[0];
875 Real* x_max_raw = dprob.RawEigenvector(nev-1);
876 Real* x_min_raw = dprob.RawEigenvector(0);
882 Real* AtAx_max_raw =
new Real[ncols];
883 Real* AtAx_min_raw =
new Real[ncols];
884 A.multMtMv(x_max_raw,AtAx_max_raw);
885 A.multMtMv(x_min_raw,AtAx_min_raw);
886 Real r_max_norm = 0.0;
887 Real r_min_norm = 0.0;
888 for (
int i = 0; i < ncols; ++i)
890 r_max_norm += std::pow(AtAx_max_raw[i] - lambda_max * x_max_raw[i],2);
891 r_min_norm += std::pow(AtAx_min_raw[i] - lambda_min * x_min_raw[i],2);
893 r_max_norm = std::sqrt(r_max_norm);
894 r_min_norm = std::sqrt(r_min_norm);
897 const Real sigma_max = std::sqrt(lambda_max);
898 const Real sigma_min = std::sqrt(lambda_min);
901 cond_2 = sigma_max / sigma_min;
909 std::cout <<
blank_ <<
"Obtained singular values of A by sol"
910 <<
"ving (A^T*A)*x = σ²*x using the ARPACK++ "
911 <<
"class ARSymStdEig:" << std::endl;
912 std::cout <<
blank_ <<
" converged eigenvalues of A^T*A: "
913 << nconv <<
" / " << nev << std::endl;
914 std::cout <<
blank_ <<
" largest eigenvalue of A^T*A: "
915 << lambda_max << std::endl;
916 std::cout <<
blank_ <<
" smallest eigenvalue of A^T*A: "
917 << lambda_min << std::endl;
918 std::cout <<
blank_ <<
" => largest singular value of A: "
919 << sigma_max << std::endl;
920 std::cout <<
blank_ <<
" => smallest singular value of A: "
921 << sigma_min << std::endl;
923 std::cout <<
blank_ <<
"Result ("
925 <<
"║residual║_2 = {" << r_max_norm <<
","
926 << r_min_norm <<
"}, " <<
"σ = {"
927 << sigma_max <<
"," << sigma_min
928 <<
"}): cond_2 = " << cond_2 << std::endl;
932 delete[] AtAx_min_raw;
933 delete[] AtAx_max_raw;
971 #endif // HAVE_ARPACKPP
973 #endif // DUNE_ISTL_EIGENVALUE_ARPACKPP_HH