3#ifndef DUNE_ISTL_EIGENVALUE_ARPACKPP_HH
4#define DUNE_ISTL_EIGENVALUE_ARPACKPP_HH
17#include <dune/istl/istlexception.hh>
48 template <
class BCRSMatrix>
49 class ArPackPlusPlus_BCRSMatrixWrapper
57 ArPackPlusPlus_BCRSMatrixWrapper (
const BCRSMatrix& A)
59 m_(A_.M() * mBlock), n_(A_.N() * nBlock)
64 "Only BCRSMatrices with blocklevel 2 are supported.");
68 domainBlockVector.resize(A_.N(),
false);
69 rangeBlockVector.resize(A_.M(),
false);
73 inline void multMv (Real* v, Real* w)
76 arrayToDomainBlockVector(v,domainBlockVector);
79 A_.mv(domainBlockVector,rangeBlockVector);
82 rangeBlockVectorToArray(rangeBlockVector,w);
86 inline void multMtMv (Real* v, Real* w)
89 arrayToDomainBlockVector(v,domainBlockVector);
92 A_.mv(domainBlockVector,rangeBlockVector);
93 A_.mtv(rangeBlockVector,domainBlockVector);
96 domainBlockVectorToArray(domainBlockVector,w);
100 inline void multMMtv (Real* v, Real* w)
103 arrayToRangeBlockVector(v,rangeBlockVector);
106 A_.mtv(rangeBlockVector,domainBlockVector);
107 A_.mv(domainBlockVector,rangeBlockVector);
110 rangeBlockVectorToArray(rangeBlockVector,w);
114 inline int nrows ()
const {
return m_; }
117 inline int ncols ()
const {
return n_; }
121 constexpr static int mBlock = BCRSMatrix::block_type::rows;
122 constexpr static int nBlock = BCRSMatrix::block_type::cols;
126 constexpr static int dbvBlockSize = nBlock;
132 constexpr static int rbvBlockSize = mBlock;
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;
145 domainBlockVectorToArray (
const DomainBlockVector& dbv, Real* v)
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];
155 rangeBlockVectorToArray (
const RangeBlockVector& rbv, Real* v)
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];
165 static inline void arrayToDomainBlockVector (
const Real* v,
166 DomainBlockVector& dbv)
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];
175 static inline void arrayToRangeBlockVector (
const Real* v,
176 RangeBlockVector& rbv)
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];
185 const BCRSMatrix& A_;
192 mutable DomainBlockVector domainBlockVector;
193 mutable RangeBlockVector rangeBlockVector;
234 template <
typename BCRSMatrix,
typename BlockVector>
235 class ArPackPlusPlus_Algorithms
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),
265 title_(
" ArPackPlusPlus_Algorithms: "),
266 blank_(title_.length(),
' ')
280 inline void computeSymMaxMagnitude (
const Real& epsilon,
281 BlockVector& x, Real& lambda)
const
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;
291 typedef ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
295 const int nrows = A.nrows();
296 const int ncols = A.ncols();
301 << nrows <<
"x" << ncols <<
").");
306 const Real tol = epsilon;
307 const int maxit = nIterationsMax_*nev;
308 Real* ev =
new Real[nev];
309 const bool ivec =
true;
314 ARSymStdEig<Real,WrappedMatrix>
315 dprob(nrows, nev, &A, &WrappedMatrix::multMv, which, ncv, tol, maxit);
318 if (verbosity_level_ > 3) dprob.Trace();
321 nconv = dprob.Eigenvalues(ev,ivec);
327 Real* x_raw = dprob.RawEigenvector(nev-1);
328 WrappedMatrix::arrayToDomainBlockVector(x_raw,x);
331 nIterations_ = dprob.GetIter();
335 Real* Ax_raw =
new Real[nrows];
336 A.multMv(x_raw,Ax_raw);
337 WrappedMatrix::arrayToDomainBlockVector(Ax_raw,r);
339 const Real r_norm = r.two_norm();
342 if (verbosity_level_ > 0)
344 if (verbosity_level_ > 1)
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;
355 std::cout << blank_ <<
"Result ("
356 <<
"#iterations = " << nIterations_ <<
", "
357 <<
"║residual║_2 = " << r_norm <<
"): "
358 <<
"λ = " << lambda << std::endl;
359 if (verbosity_level_ > 2)
382 inline void computeSymMinMagnitude (
const Real& epsilon,
383 BlockVector& x, Real& lambda)
const
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;
393 typedef ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
397 const int nrows = A.nrows();
398 const int ncols = A.ncols();
403 << nrows <<
"x" << ncols <<
").");
408 const Real tol = epsilon;
409 const int maxit = nIterationsMax_*nev;
410 Real* ev =
new Real[nev];
411 const bool ivec =
true;
416 ARSymStdEig<Real,WrappedMatrix>
417 dprob(nrows, nev, &A, &WrappedMatrix::multMv, which, ncv, tol, maxit);
420 if (verbosity_level_ > 3) dprob.Trace();
423 nconv = dprob.Eigenvalues(ev,ivec);
429 Real* x_raw = dprob.RawEigenvector(nev-1);
430 WrappedMatrix::arrayToDomainBlockVector(x_raw,x);
433 nIterations_ = dprob.GetIter();
437 Real* Ax_raw =
new Real[nrows];
438 A.multMv(x_raw,Ax_raw);
439 WrappedMatrix::arrayToDomainBlockVector(Ax_raw,r);
441 const Real r_norm = r.two_norm();
444 if (verbosity_level_ > 0)
446 if (verbosity_level_ > 1)
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;
457 std::cout << blank_ <<
"Result ("
458 <<
"#iterations = " << nIterations_ <<
", "
459 <<
"║residual║_2 = " << r_norm <<
"): "
460 <<
"λ = " << lambda << std::endl;
461 if (verbosity_level_ > 2)
484 inline void computeSymCond2 (
const Real& epsilon, Real& cond_2)
const
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;
494 typedef ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
498 const int nrows = A.nrows();
499 const int ncols = A.ncols();
504 << nrows <<
"x" << ncols <<
").");
509 const Real tol = epsilon;
510 const int maxit = nIterationsMax_*nev;
511 Real* ev =
new Real[nev];
512 const bool ivec =
true;
517 ARSymStdEig<Real,WrappedMatrix>
518 dprob(nrows, nev, &A, &WrappedMatrix::multMv, which, ncv, tol, maxit);
521 if (verbosity_level_ > 3) dprob.Trace();
524 nconv = dprob.Eigenvalues(ev,ivec);
527 const Real& lambda_max = ev[nev-1];
528 const Real& lambda_min = ev[0];
531 Real* x_max_raw = dprob.RawEigenvector(nev-1);
532 Real* x_min_raw = dprob.RawEigenvector(0);
535 cond_2 = std::abs(lambda_max / lambda_min);
538 nIterations_ = dprob.GetIter();
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)
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);
552 r_max_norm = std::sqrt(r_max_norm);
553 r_min_norm = std::sqrt(r_min_norm);
556 if (verbosity_level_ > 0)
558 if (verbosity_level_ > 1)
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;
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;
600 inline void computeNonSymMax (
const Real& epsilon,
601 BlockVector& x, Real& sigma)
const
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;
611 typedef ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
615 const int nrows = A.nrows();
616 const int ncols = A.ncols();
621 <<
"columns (" << nrows <<
"x" << ncols <<
")."
622 <<
" This case is not implemented, yet.");
627 const Real tol = epsilon;
628 const int maxit = nIterationsMax_*nev;
629 Real* ev =
new Real[nev];
630 const bool ivec =
true;
635 ARSymStdEig<Real,WrappedMatrix>
636 dprob(ncols, nev, &A, &WrappedMatrix::multMtMv, which, ncv, tol, maxit);
639 if (verbosity_level_ > 3) dprob.Trace();
642 nconv = dprob.Eigenvalues(ev,ivec);
645 const Real& lambda = ev[nev-1];
648 Real* x_raw = dprob.RawEigenvector(nev-1);
649 WrappedMatrix::arrayToDomainBlockVector(x_raw,x);
652 nIterations_ = dprob.GetIter();
656 Real* AtAx_raw =
new Real[ncols];
657 A.multMtMv(x_raw,AtAx_raw);
658 WrappedMatrix::arrayToDomainBlockVector(AtAx_raw,r);
660 const Real r_norm = r.two_norm();
664 sigma = std::sqrt(lambda);
667 if (verbosity_level_ > 0)
669 if (verbosity_level_ > 1)
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;
682 std::cout << blank_ <<
"Result ("
683 <<
"#iterations = " << nIterations_ <<
", "
684 <<
"║residual║_2 = " << r_norm <<
"): "
685 <<
"σ = " << sigma << std::endl;
686 if (verbosity_level_ > 2)
712 inline void computeNonSymMin (
const Real& epsilon,
713 BlockVector& x, Real& sigma)
const
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;
723 typedef ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
727 const int nrows = A.nrows();
728 const int ncols = A.ncols();
733 <<
"columns (" << nrows <<
"x" << ncols <<
")."
734 <<
" This case is not implemented, yet.");
739 const Real tol = epsilon;
740 const int maxit = nIterationsMax_*nev;
741 Real* ev =
new Real[nev];
742 const bool ivec =
true;
747 ARSymStdEig<Real,WrappedMatrix>
748 dprob(ncols, nev, &A, &WrappedMatrix::multMtMv, which, ncv, tol, maxit);
751 if (verbosity_level_ > 3) dprob.Trace();
754 nconv = dprob.Eigenvalues(ev,ivec);
757 const Real& lambda = ev[nev-1];
760 Real* x_raw = dprob.RawEigenvector(nev-1);
761 WrappedMatrix::arrayToDomainBlockVector(x_raw,x);
764 nIterations_ = dprob.GetIter();
768 Real* AtAx_raw =
new Real[ncols];
769 A.multMtMv(x_raw,AtAx_raw);
770 WrappedMatrix::arrayToDomainBlockVector(AtAx_raw,r);
772 const Real r_norm = r.two_norm();
776 sigma = std::sqrt(lambda);
779 if (verbosity_level_ > 0)
781 if (verbosity_level_ > 1)
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;
794 std::cout << blank_ <<
"Result ("
795 <<
"#iterations = " << nIterations_ <<
", "
796 <<
"║residual║_2 = " << r_norm <<
"): "
797 <<
"σ = " << sigma << std::endl;
798 if (verbosity_level_ > 2)
821 inline void computeNonSymCond2 (
const Real& epsilon, Real& cond_2)
const
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;
831 typedef ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
835 const int nrows = A.nrows();
836 const int ncols = A.ncols();
841 <<
"columns (" << nrows <<
"x" << ncols <<
")."
842 <<
" This case is not implemented, yet.");
847 const Real tol = epsilon;
848 const int maxit = nIterationsMax_*nev;
849 Real* ev =
new Real[nev];
850 const bool ivec =
true;
855 ARSymStdEig<Real,WrappedMatrix>
856 dprob(ncols, nev, &A, &WrappedMatrix::multMtMv, which, ncv, tol, maxit);
859 if (verbosity_level_ > 3) dprob.Trace();
862 nconv = dprob.Eigenvalues(ev,ivec);
865 const Real& lambda_max = ev[nev-1];
866 const Real& lambda_min = ev[0];
869 Real* x_max_raw = dprob.RawEigenvector(nev-1);
870 Real* x_min_raw = dprob.RawEigenvector(0);
873 nIterations_ = dprob.GetIter();
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)
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);
887 r_max_norm = std::sqrt(r_max_norm);
888 r_min_norm = std::sqrt(r_min_norm);
891 const Real sigma_max = std::sqrt(lambda_max);
892 const Real sigma_min = std::sqrt(lambda_min);
895 cond_2 = sigma_max / sigma_min;
898 if (verbosity_level_ > 0)
900 if (verbosity_level_ > 1)
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;
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;
926 delete[] AtAx_min_raw;
927 delete[] AtAx_max_raw;
935 inline unsigned int getIterationCount ()
const
937 if (nIterations_ == 0)
945 const BCRSMatrix& m_;
946 const unsigned int nIterationsMax_;
949 const unsigned int verbosity_level_;
954 mutable unsigned int nIterations_;
957 const std::string title_;
958 const std::string blank_;
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