dune-istl  2.8.0
arpackpp.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 #ifndef DUNE_ISTL_EIGENVALUE_ARPACKPP_HH
4 #define DUNE_ISTL_EIGENVALUE_ARPACKPP_HH
5 
6 #if HAVE_ARPACKPP || defined DOXYGEN
7 
8 #include <cmath> // provides std::abs, std::pow, std::sqrt
9 
10 #include <iostream> // provides std::cout, std::endl
11 #include <string> // provides std::string
12 
13 #include <dune/common/fvector.hh> // provides Dune::FieldVector
14 #include <dune/common/exceptions.hh> // provides DUNE_THROW(...)
15 
16 #include <dune/istl/blocklevel.hh> // provides Dune::blockLevel
17 #include <dune/istl/bvector.hh> // provides Dune::BlockVector
18 #include <dune/istl/istlexception.hh> // provides Dune::ISTLError
19 #include <dune/istl/io.hh> // provides Dune::printvector(...)
20 
21 #ifdef Status
22 #undef Status // prevent preprocessor from damaging the ARPACK++
23  // code when "X11/Xlib.h" is included (the latter
24  // defines Status as "#define Status int" and
25  // ARPACK++ provides a class with a method called
26  // Status)
27 #endif
28 #include "arssym.h" // provides ARSymStdEig
29 
30 namespace Dune
31 {
32 
37  namespace Impl {
53  template <class BCRSMatrix>
54  class ArPackPlusPlus_BCRSMatrixWrapper
55  {
56  public:
58  typedef typename BCRSMatrix::field_type Real;
59 
60  public:
62  ArPackPlusPlus_BCRSMatrixWrapper (const BCRSMatrix& A)
63  : A_(A),
64  m_(A_.M() * mBlock), n_(A_.N() * nBlock)
65  {
66  // assert that BCRSMatrix type has blocklevel 2
67  static_assert
68  (blockLevel<BCRSMatrix>() == 2,
69  "Only BCRSMatrices with blocklevel 2 are supported.");
70 
71  // allocate memory for auxiliary block vector objects
72  // which are compatible to matrix rows / columns
73  domainBlockVector.resize(A_.N());
74  rangeBlockVector.resize(A_.M());
75  }
76 
78  inline void multMv (Real* v, Real* w)
79  {
80  // get vector v as an object of appropriate type
81  arrayToDomainBlockVector(v,domainBlockVector);
82 
83  // perform matrix-vector product
84  A_.mv(domainBlockVector,rangeBlockVector);
85 
86  // get vector w from object of appropriate type
87  rangeBlockVectorToArray(rangeBlockVector,w);
88  };
89 
91  inline void multMtMv (Real* v, Real* w)
92  {
93  // get vector v as an object of appropriate type
94  arrayToDomainBlockVector(v,domainBlockVector);
95 
96  // perform matrix-vector product
97  A_.mv(domainBlockVector,rangeBlockVector);
98  A_.mtv(rangeBlockVector,domainBlockVector);
99 
100  // get vector w from object of appropriate type
101  domainBlockVectorToArray(domainBlockVector,w);
102  };
103 
105  inline void multMMtv (Real* v, Real* w)
106  {
107  // get vector v as an object of appropriate type
108  arrayToRangeBlockVector(v,rangeBlockVector);
109 
110  // perform matrix-vector product
111  A_.mtv(rangeBlockVector,domainBlockVector);
112  A_.mv(domainBlockVector,rangeBlockVector);
113 
114  // get vector w from object of appropriate type
115  rangeBlockVectorToArray(rangeBlockVector,w);
116  };
117 
119  inline int nrows () const { return m_; }
120 
122  inline int ncols () const { return n_; }
123 
124  protected:
125  // Number of rows and columns in each block of the matrix
126  constexpr static int mBlock = BCRSMatrix::block_type::rows;
127  constexpr static int nBlock = BCRSMatrix::block_type::cols;
128 
129  // Type of vectors in the domain of the linear map associated with
130  // the matrix, i.e. block vectors compatible to matrix rows
131  constexpr static int dbvBlockSize = nBlock;
132  typedef Dune::FieldVector<Real,dbvBlockSize> DomainBlockVectorBlock;
133  typedef Dune::BlockVector<DomainBlockVectorBlock> DomainBlockVector;
134 
135  // Type of vectors in the range of the linear map associated with
136  // the matrix, i.e. block vectors compatible to matrix columns
137  constexpr static int rbvBlockSize = mBlock;
138  typedef Dune::FieldVector<Real,rbvBlockSize> RangeBlockVectorBlock;
139  typedef Dune::BlockVector<RangeBlockVectorBlock> RangeBlockVector;
140 
141  // Types for vector index access
142  typedef typename DomainBlockVector::size_type dbv_size_type;
143  typedef typename RangeBlockVector::size_type rbv_size_type;
144  typedef typename DomainBlockVectorBlock::size_type dbvb_size_type;
145  typedef typename RangeBlockVectorBlock::size_type rbvb_size_type;
146 
147  // Get vector v from a block vector object which is compatible to
148  // matrix rows
149  static inline void
150  domainBlockVectorToArray (const DomainBlockVector& dbv, Real* v)
151  {
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];
155  }
156 
157  // Get vector v from a block vector object which is compatible to
158  // matrix columns
159  static inline void
160  rangeBlockVectorToArray (const RangeBlockVector& rbv, Real* v)
161  {
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];
165  }
166 
167  public:
170  static inline void arrayToDomainBlockVector (const Real* v,
171  DomainBlockVector& dbv)
172  {
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];
176  }
177 
180  static inline void arrayToRangeBlockVector (const Real* v,
181  RangeBlockVector& rbv)
182  {
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];
186  }
187 
188  protected:
189  // The DUNE-ISTL BCRSMatrix
190  const BCRSMatrix& A_;
191 
192  // Number of rows and columns in the matrix
193  const int m_, n_;
194 
195  // Auxiliary block vector objects which are
196  // compatible to matrix rows / columns
197  mutable DomainBlockVector domainBlockVector;
198  mutable RangeBlockVector rangeBlockVector;
199  };
200  } // end namespace Impl
201 
241  template <typename BCRSMatrix, typename BlockVector>
243  {
244  public:
245  typedef typename BlockVector::field_type Real;
246 
247  public:
267  const unsigned int nIterationsMax = 100000,
268  const unsigned int verbosity_level = 0)
269  : m_(m), nIterationsMax_(nIterationsMax),
270  verbosity_level_(verbosity_level),
271  nIterations_(0),
272  title_(" ArPackPlusPlus_Algorithms: "),
273  blank_(title_.length(),' ')
274  {}
275 
287  inline void computeSymMaxMagnitude (const Real& epsilon,
288  BlockVector& x, Real& lambda) const
289  {
290  // print verbosity information
291  if (verbosity_level_ > 0)
292  std::cout << title_ << "Computing an approximation of "
293  << "the dominant eigenvalue of a matrix which "
294  << "is assumed to be symmetric." << std::endl;
295 
296  // use type ArPackPlusPlus_BCRSMatrixWrapper to store matrix information
297  // and to perform the product A*v (LU decomposition is not used)
298  typedef Impl::ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
299  WrappedMatrix A(m_);
300 
301  // get number of rows and columns in A
302  const int nrows = A.nrows();
303  const int ncols = A.ncols();
304 
305  // assert that A is square
306  if (nrows != ncols)
307  DUNE_THROW(Dune::ISTLError,"Matrix is not square ("
308  << nrows << "x" << ncols << ").");
309 
310  // allocate memory for variables, set parameters
311  const int nev = 1; // Number of eigenvalues to compute
312  int ncv = std::min(20, nrows); // Number of Arnoldi vectors generated at each iteration (0 == auto)
313  const Real tol = epsilon; // Stopping tolerance (relative accuracy of Ritz values) (0 == machine precision)
314  const int maxit = nIterationsMax_*nev; // Maximum number of Arnoldi update iterations allowed (0 == 100*nev)
315  Real* ev = new Real[nev]; // Computed eigenvalues of A
316  const bool ivec = true; // Flag deciding if eigenvectors shall be determined
317  int nconv; // Number of converged eigenvalues
318 
319  // define what we need: eigenvalues with largest magnitude
320  char which[] = "LM";
321  ARSymStdEig<Real,WrappedMatrix>
322  dprob(nrows, nev, &A, &WrappedMatrix::multMv, which, ncv, tol, maxit);
323 
324  // set ARPACK verbosity mode if requested
325  if (verbosity_level_ > 3) dprob.Trace();
326 
327  // find eigenvalues and eigenvectors of A, obtain the eigenvalues
328  nconv = dprob.Eigenvalues(ev,ivec);
329 
330  // obtain approximated dominant eigenvalue of A
331  lambda = ev[nev-1];
332 
333  // obtain associated approximated eigenvector of A
334  Real* x_raw = dprob.RawEigenvector(nev-1);
335  WrappedMatrix::arrayToDomainBlockVector(x_raw,x);
336 
337  // obtain number of Arnoldi update iterations actually taken
338  nIterations_ = dprob.GetIter();
339 
340  // compute residual norm
341  BlockVector r(x);
342  Real* Ax_raw = new Real[nrows];
343  A.multMv(x_raw,Ax_raw);
344  WrappedMatrix::arrayToDomainBlockVector(Ax_raw,r);
345  r.axpy(-lambda,x);
346  const Real r_norm = r.two_norm();
347 
348  // print verbosity information
349  if (verbosity_level_ > 0)
350  {
351  if (verbosity_level_ > 1)
352  {
353  // print some information about the problem
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;
361  }
362  std::cout << blank_ << "Result ("
363  << "#iterations = " << nIterations_ << ", "
364  << "║residual║_2 = " << r_norm << "): "
365  << "λ = " << lambda << std::endl;
366  if (verbosity_level_ > 2)
367  {
368  // print approximated eigenvector via DUNE-ISTL I/O methods
369  Dune::printvector(std::cout,x,blank_+"x",blank_+"row");
370  }
371  }
372 
373  // free dynamically allocated memory
374  delete[] Ax_raw;
375  delete[] ev;
376  }
377 
389  inline void computeSymMinMagnitude (const Real& epsilon,
390  BlockVector& x, Real& lambda) const
391  {
392  // print verbosity information
393  if (verbosity_level_ > 0)
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;
397 
398  // use type ArPackPlusPlus_BCRSMatrixWrapper to store matrix information
399  // and to perform the product A*v (LU decomposition is not used)
400  typedef Impl::ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
401  WrappedMatrix A(m_);
402 
403  // get number of rows and columns in A
404  const int nrows = A.nrows();
405  const int ncols = A.ncols();
406 
407  // assert that A is square
408  if (nrows != ncols)
409  DUNE_THROW(Dune::ISTLError,"Matrix is not square ("
410  << nrows << "x" << ncols << ").");
411 
412  // allocate memory for variables, set parameters
413  const int nev = 1; // Number of eigenvalues to compute
414  int ncv = std::min(20, nrows); // Number of Arnoldi vectors generated at each iteration (0 == auto)
415  const Real tol = epsilon; // Stopping tolerance (relative accuracy of Ritz values) (0 == machine precision)
416  const int maxit = nIterationsMax_*nev; // Maximum number of Arnoldi update iterations allowed (0 == 100*nev)
417  Real* ev = new Real[nev]; // Computed eigenvalues of A
418  const bool ivec = true; // Flag deciding if eigenvectors shall be determined
419  int nconv; // Number of converged eigenvalues
420 
421  // define what we need: eigenvalues with smallest magnitude
422  char which[] = "SM";
423  ARSymStdEig<Real,WrappedMatrix>
424  dprob(nrows, nev, &A, &WrappedMatrix::multMv, which, ncv, tol, maxit);
425 
426  // set ARPACK verbosity mode if requested
427  if (verbosity_level_ > 3) dprob.Trace();
428 
429  // find eigenvalues and eigenvectors of A, obtain the eigenvalues
430  nconv = dprob.Eigenvalues(ev,ivec);
431 
432  // obtain approximated least dominant eigenvalue of A
433  lambda = ev[nev-1];
434 
435  // obtain associated approximated eigenvector of A
436  Real* x_raw = dprob.RawEigenvector(nev-1);
437  WrappedMatrix::arrayToDomainBlockVector(x_raw,x);
438 
439  // obtain number of Arnoldi update iterations actually taken
440  nIterations_ = dprob.GetIter();
441 
442  // compute residual norm
443  BlockVector r(x);
444  Real* Ax_raw = new Real[nrows];
445  A.multMv(x_raw,Ax_raw);
446  WrappedMatrix::arrayToDomainBlockVector(Ax_raw,r);
447  r.axpy(-lambda,x);
448  const Real r_norm = r.two_norm();
449 
450  // print verbosity information
451  if (verbosity_level_ > 0)
452  {
453  if (verbosity_level_ > 1)
454  {
455  // print some information about the problem
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;
463  }
464  std::cout << blank_ << "Result ("
465  << "#iterations = " << nIterations_ << ", "
466  << "║residual║_2 = " << r_norm << "): "
467  << "λ = " << lambda << std::endl;
468  if (verbosity_level_ > 2)
469  {
470  // print approximated eigenvector via DUNE-ISTL I/O methods
471  Dune::printvector(std::cout,x,blank_+"x",blank_+"row");
472  }
473  }
474 
475  // free dynamically allocated memory
476  delete[] Ax_raw;
477  delete[] ev;
478  }
479 
491  inline void computeSymCond2 (const Real& epsilon, Real& cond_2) const
492  {
493  // print verbosity information
494  if (verbosity_level_ > 0)
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;
498 
499  // use type ArPackPlusPlus_BCRSMatrixWrapper to store matrix information
500  // and to perform the product A*v (LU decomposition is not used)
501  typedef Impl::ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
502  WrappedMatrix A(m_);
503 
504  // get number of rows and columns in A
505  const int nrows = A.nrows();
506  const int ncols = A.ncols();
507 
508  // assert that A is square
509  if (nrows != ncols)
510  DUNE_THROW(Dune::ISTLError,"Matrix is not square ("
511  << nrows << "x" << ncols << ").");
512 
513  // allocate memory for variables, set parameters
514  const int nev = 2; // Number of eigenvalues to compute
515  int ncv = std::min(20, nrows); // Number of Arnoldi vectors generated at each iteration (0 == auto)
516  const Real tol = epsilon; // Stopping tolerance (relative accuracy of Ritz values) (0 == machine precision)
517  const int maxit = nIterationsMax_*nev; // Maximum number of Arnoldi update iterations allowed (0 == 100*nev)
518  Real* ev = new Real[nev]; // Computed eigenvalues of A
519  const bool ivec = true; // Flag deciding if eigenvectors shall be determined
520  int nconv; // Number of converged eigenvalues
521 
522  // define what we need: eigenvalues from both ends of the spectrum
523  char which[] = "BE";
524  ARSymStdEig<Real,WrappedMatrix>
525  dprob(nrows, nev, &A, &WrappedMatrix::multMv, which, ncv, tol, maxit);
526 
527  // set ARPACK verbosity mode if requested
528  if (verbosity_level_ > 3) dprob.Trace();
529 
530  // find eigenvalues and eigenvectors of A, obtain the eigenvalues
531  nconv = dprob.Eigenvalues(ev,ivec);
532 
533  // obtain approximated dominant and least dominant eigenvalue of A
534  const Real& lambda_max = ev[nev-1];
535  const Real& lambda_min = ev[0];
536 
537  // obtain associated approximated eigenvectors of A
538  Real* x_max_raw = dprob.RawEigenvector(nev-1);
539  Real* x_min_raw = dprob.RawEigenvector(0);
540 
541  // obtain approximated spectral condition number of A
542  cond_2 = std::abs(lambda_max / lambda_min);
543 
544  // obtain number of Arnoldi update iterations actually taken
545  nIterations_ = dprob.GetIter();
546 
547  // compute each residual norm
548  Real* Ax_max_raw = new Real[nrows];
549  Real* Ax_min_raw = new Real[nrows];
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)
555  {
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);
558  }
559  r_max_norm = std::sqrt(r_max_norm);
560  r_min_norm = std::sqrt(r_min_norm);
561 
562  // print verbosity information
563  if (verbosity_level_ > 0)
564  {
565  if (verbosity_level_ > 1)
566  {
567  // print some information about the problem
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;
579  }
580  std::cout << blank_ << "Result ("
581  << "#iterations = " << nIterations_ << ", "
582  << "║residual║_2 = {" << r_max_norm << ","
583  << r_min_norm << "}, " << "λ = {"
584  << lambda_max << "," << lambda_min
585  << "}): cond_2 = " << cond_2 << std::endl;
586  }
587 
588  // free dynamically allocated memory
589  delete[] Ax_min_raw;
590  delete[] Ax_max_raw;
591  delete[] ev;
592  }
593 
607  inline void computeNonSymMax (const Real& epsilon,
608  BlockVector& x, Real& sigma) const
609  {
610  // print verbosity information
611  if (verbosity_level_ > 0)
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;
615 
616  // use type ArPackPlusPlus_BCRSMatrixWrapper to store matrix information
617  // and to perform the product A^T*A*v (LU decomposition is not used)
618  typedef Impl::ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
619  WrappedMatrix A(m_);
620 
621  // get number of rows and columns in A
622  const int nrows = A.nrows();
623  const int ncols = A.ncols();
624 
625  // assert that A has more rows than columns (extend code later to the opposite case!)
626  if (nrows < ncols)
627  DUNE_THROW(Dune::ISTLError,"Matrix has less rows than "
628  << "columns (" << nrows << "x" << ncols << ")."
629  << " This case is not implemented, yet.");
630 
631  // allocate memory for variables, set parameters
632  const int nev = 1; // Number of eigenvalues to compute
633  int ncv = std::min(20, nrows); // Number of Arnoldi vectors generated at each iteration (0 == auto)
634  const Real tol = epsilon; // Stopping tolerance (relative accuracy of Ritz values) (0 == machine precision)
635  const int maxit = nIterationsMax_*nev; // Maximum number of Arnoldi update iterations allowed (0 == 100*nev)
636  Real* ev = new Real[nev]; // Computed eigenvalues of A^T*A
637  const bool ivec = true; // Flag deciding if eigenvectors shall be determined
638  int nconv; // Number of converged eigenvalues
639 
640  // define what we need: eigenvalues with largest algebraic value
641  char which[] = "LA";
642  ARSymStdEig<Real,WrappedMatrix>
643  dprob(ncols, nev, &A, &WrappedMatrix::multMtMv, which, ncv, tol, maxit);
644 
645  // set ARPACK verbosity mode if requested
646  if (verbosity_level_ > 3) dprob.Trace();
647 
648  // find eigenvalues and eigenvectors of A^T*A, obtain the eigenvalues
649  nconv = dprob.Eigenvalues(ev,ivec);
650 
651  // obtain approximated largest eigenvalue of A^T*A
652  const Real& lambda = ev[nev-1];
653 
654  // obtain associated approximated eigenvector of A^T*A
655  Real* x_raw = dprob.RawEigenvector(nev-1);
656  WrappedMatrix::arrayToDomainBlockVector(x_raw,x);
657 
658  // obtain number of Arnoldi update iterations actually taken
659  nIterations_ = dprob.GetIter();
660 
661  // compute residual norm
662  BlockVector r(x);
663  Real* AtAx_raw = new Real[ncols];
664  A.multMtMv(x_raw,AtAx_raw);
665  WrappedMatrix::arrayToDomainBlockVector(AtAx_raw,r);
666  r.axpy(-lambda,x);
667  const Real r_norm = r.two_norm();
668 
669  // calculate largest singular value of A (note that
670  // x is right-singular / left-singular vector of A)
671  sigma = std::sqrt(lambda);
672 
673  // print verbosity information
674  if (verbosity_level_ > 0)
675  {
676  if (verbosity_level_ > 1)
677  {
678  // print some information about the problem
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;
688  }
689  std::cout << blank_ << "Result ("
690  << "#iterations = " << nIterations_ << ", "
691  << "║residual║_2 = " << r_norm << "): "
692  << "σ = " << sigma << std::endl;
693  if (verbosity_level_ > 2)
694  {
695  // print approximated right-singular / left-singular vector
696  // via DUNE-ISTL I/O methods
697  Dune::printvector(std::cout,x,blank_+"x",blank_+"row");
698  }
699  }
700 
701  // free dynamically allocated memory
702  delete[] AtAx_raw;
703  delete[] ev;
704  }
705 
719  inline void computeNonSymMin (const Real& epsilon,
720  BlockVector& x, Real& sigma) const
721  {
722  // print verbosity information
723  if (verbosity_level_ > 0)
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;
727 
728  // use type ArPackPlusPlus_BCRSMatrixWrapper to store matrix information
729  // and to perform the product A^T*A*v (LU decomposition is not used)
730  typedef Impl::ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
731  WrappedMatrix A(m_);
732 
733  // get number of rows and columns in A
734  const int nrows = A.nrows();
735  const int ncols = A.ncols();
736 
737  // assert that A has more rows than columns (extend code later to the opposite case!)
738  if (nrows < ncols)
739  DUNE_THROW(Dune::ISTLError,"Matrix has less rows than "
740  << "columns (" << nrows << "x" << ncols << ")."
741  << " This case is not implemented, yet.");
742 
743  // allocate memory for variables, set parameters
744  const int nev = 1; // Number of eigenvalues to compute
745  int ncv = std::min(20, nrows); // Number of Arnoldi vectors generated at each iteration (0 == auto)
746  const Real tol = epsilon; // Stopping tolerance (relative accuracy of Ritz values) (0 == machine precision)
747  const int maxit = nIterationsMax_*nev; // Maximum number of Arnoldi update iterations allowed (0 == 100*nev)
748  Real* ev = new Real[nev]; // Computed eigenvalues of A^T*A
749  const bool ivec = true; // Flag deciding if eigenvectors shall be determined
750  int nconv; // Number of converged eigenvalues
751 
752  // define what we need: eigenvalues with smallest algebraic value
753  char which[] = "SA";
754  ARSymStdEig<Real,WrappedMatrix>
755  dprob(ncols, nev, &A, &WrappedMatrix::multMtMv, which, ncv, tol, maxit);
756 
757  // set ARPACK verbosity mode if requested
758  if (verbosity_level_ > 3) dprob.Trace();
759 
760  // find eigenvalues and eigenvectors of A^T*A, obtain the eigenvalues
761  nconv = dprob.Eigenvalues(ev,ivec);
762 
763  // obtain approximated smallest eigenvalue of A^T*A
764  const Real& lambda = ev[nev-1];
765 
766  // obtain associated approximated eigenvector of A^T*A
767  Real* x_raw = dprob.RawEigenvector(nev-1);
768  WrappedMatrix::arrayToDomainBlockVector(x_raw,x);
769 
770  // obtain number of Arnoldi update iterations actually taken
771  nIterations_ = dprob.GetIter();
772 
773  // compute residual norm
774  BlockVector r(x);
775  Real* AtAx_raw = new Real[ncols];
776  A.multMtMv(x_raw,AtAx_raw);
777  WrappedMatrix::arrayToDomainBlockVector(AtAx_raw,r);
778  r.axpy(-lambda,x);
779  const Real r_norm = r.two_norm();
780 
781  // calculate smallest singular value of A (note that
782  // x is right-singular / left-singular vector of A)
783  sigma = std::sqrt(lambda);
784 
785  // print verbosity information
786  if (verbosity_level_ > 0)
787  {
788  if (verbosity_level_ > 1)
789  {
790  // print some information about the problem
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;
800  }
801  std::cout << blank_ << "Result ("
802  << "#iterations = " << nIterations_ << ", "
803  << "║residual║_2 = " << r_norm << "): "
804  << "σ = " << sigma << std::endl;
805  if (verbosity_level_ > 2)
806  {
807  // print approximated right-singular / left-singular vector
808  // via DUNE-ISTL I/O methods
809  Dune::printvector(std::cout,x,blank_+"x",blank_+"row");
810  }
811  }
812 
813  // free dynamically allocated memory
814  delete[] AtAx_raw;
815  delete[] ev;
816  }
817 
828  inline void computeNonSymCond2 (const Real& epsilon, Real& cond_2) const
829  {
830  // print verbosity information
831  if (verbosity_level_ > 0)
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;
835 
836  // use type ArPackPlusPlus_BCRSMatrixWrapper to store matrix information
837  // and to perform the product A^T*A*v (LU decomposition is not used)
838  typedef Impl::ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
839  WrappedMatrix A(m_);
840 
841  // get number of rows and columns in A
842  const int nrows = A.nrows();
843  const int ncols = A.ncols();
844 
845  // assert that A has more rows than columns (extend code later to the opposite case!)
846  if (nrows < ncols)
847  DUNE_THROW(Dune::ISTLError,"Matrix has less rows than "
848  << "columns (" << nrows << "x" << ncols << ")."
849  << " This case is not implemented, yet.");
850 
851  // allocate memory for variables, set parameters
852  const int nev = 2; // Number of eigenvalues to compute
853  int ncv = std::min(20, nrows); // Number of Arnoldi vectors generated at each iteration (0 == auto)
854  const Real tol = epsilon; // Stopping tolerance (relative accuracy of Ritz values) (0 == machine precision)
855  const int maxit = nIterationsMax_*nev; // Maximum number of Arnoldi update iterations allowed (0 == 100*nev)
856  Real* ev = new Real[nev]; // Computed eigenvalues of A^T*A
857  const bool ivec = true; // Flag deciding if eigenvectors shall be determined
858  int nconv; // Number of converged eigenvalues
859 
860  // define what we need: eigenvalues from both ends of the spectrum
861  char which[] = "BE";
862  ARSymStdEig<Real,WrappedMatrix>
863  dprob(ncols, nev, &A, &WrappedMatrix::multMtMv, which, ncv, tol, maxit);
864 
865  // set ARPACK verbosity mode if requested
866  if (verbosity_level_ > 3) dprob.Trace();
867 
868  // find eigenvalues and eigenvectors of A^T*A, obtain the eigenvalues
869  nconv = dprob.Eigenvalues(ev,ivec);
870 
871  // obtain approximated largest and smallest eigenvalue of A^T*A
872  const Real& lambda_max = ev[nev-1];
873  const Real& lambda_min = ev[0];
874 
875  // obtain associated approximated eigenvectors of A^T*A
876  Real* x_max_raw = dprob.RawEigenvector(nev-1);
877  Real* x_min_raw = dprob.RawEigenvector(0);
878 
879  // obtain number of Arnoldi update iterations actually taken
880  nIterations_ = dprob.GetIter();
881 
882  // compute each residual norm
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)
890  {
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);
893  }
894  r_max_norm = std::sqrt(r_max_norm);
895  r_min_norm = std::sqrt(r_min_norm);
896 
897  // calculate largest and smallest singular value of A
898  const Real sigma_max = std::sqrt(lambda_max);
899  const Real sigma_min = std::sqrt(lambda_min);
900 
901  // obtain approximated spectral condition number of A
902  cond_2 = sigma_max / sigma_min;
903 
904  // print verbosity information
905  if (verbosity_level_ > 0)
906  {
907  if (verbosity_level_ > 1)
908  {
909  // print some information about the problem
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;
923  }
924  std::cout << blank_ << "Result ("
925  << "#iterations = " << nIterations_ << ", "
926  << "║residual║_2 = {" << r_max_norm << ","
927  << r_min_norm << "}, " << "σ = {"
928  << sigma_max << "," << sigma_min
929  << "}): cond_2 = " << cond_2 << std::endl;
930  }
931 
932  // free dynamically allocated memory
933  delete[] AtAx_min_raw;
934  delete[] AtAx_max_raw;
935  delete[] ev;
936  }
937 
942  inline unsigned int getIterationCount () const
943  {
944  if (nIterations_ == 0)
945  DUNE_THROW(Dune::ISTLError,"No algorithm applied, yet.");
946 
947  return nIterations_;
948  }
949 
950  protected:
951  // parameters related to iterative eigenvalue algorithms
952  const BCRSMatrix& m_;
953  const unsigned int nIterationsMax_;
954 
955  // verbosity setting
956  const unsigned int verbosity_level_;
957 
958  // memory for storing temporary variables (mutable as they shall
959  // just be effectless auxiliary variables of the const apply*(...)
960  // methods)
961  mutable unsigned int nIterations_;
962 
963  // constants for printing verbosity information
964  const std::string title_;
965  const std::string blank_;
966  };
967 
970 } // namespace Dune
971 
972 #endif // HAVE_ARPACKPP
973 
974 #endif // DUNE_ISTL_EIGENVALUE_ARPACKPP_HH
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