dune-istl  2.5.1
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
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/bvector.hh> // provides Dune::BlockVector
17 #include <dune/istl/istlexception.hh> // provides Dune::ISTLError
18 #include <dune/istl/io.hh> // provides Dune::printvector(...)
19 
20 #ifdef Status
21 #undef Status // prevent preprocessor from damaging the ARPACK++
22  // code when "X11/Xlib.h" is included (the latter
23  // defines Status as "#define Status int" and
24  // ARPACK++ provides a class with a method called
25  // Status)
26 #endif
27 #include "arssym.h" // provides ARSymStdEig
28 
29 
30 namespace Dune
31 {
32 
48  template <class BCRSMatrix>
49  class ArPackPlusPlus_BCRSMatrixWrapper
50  {
51  public:
53  typedef typename BCRSMatrix::field_type Real;
54 
55  public:
57  ArPackPlusPlus_BCRSMatrixWrapper (const BCRSMatrix& A)
58  : A_(A),
59  m_(A_.M() * mBlock), n_(A_.N() * nBlock)
60  {
61  // assert that BCRSMatrix type has blocklevel 2
62  static_assert
64  "Only BCRSMatrices with blocklevel 2 are supported.");
65 
66  // allocate memory for auxiliary block vector objects
67  // which are compatible to matrix rows / columns
68  domainBlockVector.resize(A_.N(),false);
69  rangeBlockVector.resize(A_.M(),false);
70  }
71 
73  inline void multMv (Real* v, Real* w)
74  {
75  // get vector v as an object of appropriate type
76  arrayToDomainBlockVector(v,domainBlockVector);
77 
78  // perform matrix-vector product
79  A_.mv(domainBlockVector,rangeBlockVector);
80 
81  // get vector w from object of appropriate type
82  rangeBlockVectorToArray(rangeBlockVector,w);
83  };
84 
86  inline void multMtMv (Real* v, Real* w)
87  {
88  // get vector v as an object of appropriate type
89  arrayToDomainBlockVector(v,domainBlockVector);
90 
91  // perform matrix-vector product
92  A_.mv(domainBlockVector,rangeBlockVector);
93  A_.mtv(rangeBlockVector,domainBlockVector);
94 
95  // get vector w from object of appropriate type
96  domainBlockVectorToArray(domainBlockVector,w);
97  };
98 
100  inline void multMMtv (Real* v, Real* w)
101  {
102  // get vector v as an object of appropriate type
103  arrayToRangeBlockVector(v,rangeBlockVector);
104 
105  // perform matrix-vector product
106  A_.mtv(rangeBlockVector,domainBlockVector);
107  A_.mv(domainBlockVector,rangeBlockVector);
108 
109  // get vector w from object of appropriate type
110  rangeBlockVectorToArray(rangeBlockVector,w);
111  };
112 
114  inline int nrows () const { return m_; }
115 
117  inline int ncols () const { return n_; }
118 
119  protected:
120  // Number of rows and columns in each block of the matrix
121  constexpr static int mBlock = BCRSMatrix::block_type::rows;
122  constexpr static int nBlock = BCRSMatrix::block_type::cols;
123 
124  // Type of vectors in the domain of the linear map associated with
125  // the matrix, i.e. block vectors compatible to matrix rows
126  constexpr static int dbvBlockSize = nBlock;
127  typedef Dune::FieldVector<Real,dbvBlockSize> DomainBlockVectorBlock;
128  typedef Dune::BlockVector<DomainBlockVectorBlock> DomainBlockVector;
129 
130  // Type of vectors in the range of the linear map associated with
131  // the matrix, i.e. block vectors compatible to matrix columns
132  constexpr static int rbvBlockSize = mBlock;
133  typedef Dune::FieldVector<Real,rbvBlockSize> RangeBlockVectorBlock;
134  typedef Dune::BlockVector<RangeBlockVectorBlock> RangeBlockVector;
135 
136  // Types for vector index access
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;
141 
142  // Get vector v from a block vector object which is compatible to
143  // matrix rows
144  static inline void
145  domainBlockVectorToArray (const DomainBlockVector& dbv, Real* v)
146  {
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];
150  }
151 
152  // Get vector v from a block vector object which is compatible to
153  // matrix columns
154  static inline void
155  rangeBlockVectorToArray (const RangeBlockVector& rbv, Real* v)
156  {
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];
160  }
161 
162  public:
165  static inline void arrayToDomainBlockVector (const Real* v,
166  DomainBlockVector& dbv)
167  {
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];
171  }
172 
175  static inline void arrayToRangeBlockVector (const Real* v,
176  RangeBlockVector& rbv)
177  {
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];
181  }
182 
183  protected:
184  // The DUNE-ISTL BCRSMatrix
185  const BCRSMatrix& A_;
186 
187  // Number of rows and columns in the matrix
188  const int m_, n_;
189 
190  // Auxiliary block vector objects which are
191  // compatible to matrix rows / columns
192  mutable DomainBlockVector domainBlockVector;
193  mutable RangeBlockVector rangeBlockVector;
194  };
195 
196 
234  template <typename BCRSMatrix, typename BlockVector>
235  class ArPackPlusPlus_Algorithms
236  {
237  public:
238  typedef typename BlockVector::field_type Real;
239 
240  public:
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),
264  nIterations_(0),
265  title_(" ArPackPlusPlus_Algorithms: "),
266  blank_(title_.length(),' ')
267  {}
268 
280  inline void computeSymMaxMagnitude (const Real& epsilon,
281  BlockVector& x, Real& lambda) const
282  {
283  // print verbosity information
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;
288 
289  // use type ArPackPlusPlus_BCRSMatrixWrapper to store matrix information
290  // and to perform the product A*v (LU decomposition is not used)
291  typedef ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
292  WrappedMatrix A(m_);
293 
294  // get number of rows and columns in A
295  const int nrows = A.nrows();
296  const int ncols = A.ncols();
297 
298  // assert that A is square
299  if (nrows != ncols)
300  DUNE_THROW(Dune::ISTLError,"Matrix is not square ("
301  << nrows << "x" << ncols << ").");
302 
303  // allocate memory for variables, set parameters
304  const int nev = 1; // Number of eigenvalues to compute
305  const int ncv = 20; // Number of Arnoldi vectors generated at each iteration (0 == auto)
306  const Real tol = epsilon; // Stopping tolerance (relative accuracy of Ritz values) (0 == machine precision)
307  const int maxit = nIterationsMax_*nev; // Maximum number of Arnoldi update iterations allowed (0 == 100*nev)
308  Real* ev = new Real[nev]; // Computed eigenvalues of A
309  const bool ivec = true; // Flag deciding if eigenvectors shall be determined
310  int nconv; // Number of converged eigenvalues
311 
312  // define what we need: eigenvalues with largest magnitude
313  char which[] = "LM";
314  ARSymStdEig<Real,WrappedMatrix>
315  dprob(nrows, nev, &A, &WrappedMatrix::multMv, which, ncv, tol, maxit);
316 
317  // set ARPACK verbosity mode if requested
318  if (verbosity_level_ > 3) dprob.Trace();
319 
320  // find eigenvalues and eigenvectors of A, obtain the eigenvalues
321  nconv = dprob.Eigenvalues(ev,ivec);
322 
323  // obtain approximated dominant eigenvalue of A
324  lambda = ev[nev-1];
325 
326  // obtain associated approximated eigenvector of A
327  Real* x_raw = dprob.RawEigenvector(nev-1);
328  WrappedMatrix::arrayToDomainBlockVector(x_raw,x);
329 
330  // obtain number of Arnoldi update iterations actually taken
331  nIterations_ = dprob.GetIter();
332 
333  // compute residual norm
334  BlockVector r(x);
335  Real* Ax_raw = new Real[nrows];
336  A.multMv(x_raw,Ax_raw);
337  WrappedMatrix::arrayToDomainBlockVector(Ax_raw,r);
338  r.axpy(-lambda,x);
339  const Real r_norm = r.two_norm();
340 
341  // print verbosity information
342  if (verbosity_level_ > 0)
343  {
344  if (verbosity_level_ > 1)
345  {
346  // print some information about the problem
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;
354  }
355  std::cout << blank_ << "Result ("
356  << "#iterations = " << nIterations_ << ", "
357  << "║residual║_2 = " << r_norm << "): "
358  << "λ = " << lambda << std::endl;
359  if (verbosity_level_ > 2)
360  {
361  // print approximated eigenvector via DUNE-ISTL I/O methods
362  Dune::printvector(std::cout,x,blank_+"x",blank_+"row");
363  }
364  }
365 
366  // free dynamically allocated memory
367  delete[] Ax_raw;
368  delete[] ev;
369  }
370 
382  inline void computeSymMinMagnitude (const Real& epsilon,
383  BlockVector& x, Real& lambda) const
384  {
385  // print verbosity information
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;
390 
391  // use type ArPackPlusPlus_BCRSMatrixWrapper to store matrix information
392  // and to perform the product A*v (LU decomposition is not used)
393  typedef ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
394  WrappedMatrix A(m_);
395 
396  // get number of rows and columns in A
397  const int nrows = A.nrows();
398  const int ncols = A.ncols();
399 
400  // assert that A is square
401  if (nrows != ncols)
402  DUNE_THROW(Dune::ISTLError,"Matrix is not square ("
403  << nrows << "x" << ncols << ").");
404 
405  // allocate memory for variables, set parameters
406  const int nev = 1; // Number of eigenvalues to compute
407  const int ncv = 20; // Number of Arnoldi vectors generated at each iteration (0 == auto)
408  const Real tol = epsilon; // Stopping tolerance (relative accuracy of Ritz values) (0 == machine precision)
409  const int maxit = nIterationsMax_*nev; // Maximum number of Arnoldi update iterations allowed (0 == 100*nev)
410  Real* ev = new Real[nev]; // Computed eigenvalues of A
411  const bool ivec = true; // Flag deciding if eigenvectors shall be determined
412  int nconv; // Number of converged eigenvalues
413 
414  // define what we need: eigenvalues with smallest magnitude
415  char which[] = "SM";
416  ARSymStdEig<Real,WrappedMatrix>
417  dprob(nrows, nev, &A, &WrappedMatrix::multMv, which, ncv, tol, maxit);
418 
419  // set ARPACK verbosity mode if requested
420  if (verbosity_level_ > 3) dprob.Trace();
421 
422  // find eigenvalues and eigenvectors of A, obtain the eigenvalues
423  nconv = dprob.Eigenvalues(ev,ivec);
424 
425  // obtain approximated least dominant eigenvalue of A
426  lambda = ev[nev-1];
427 
428  // obtain associated approximated eigenvector of A
429  Real* x_raw = dprob.RawEigenvector(nev-1);
430  WrappedMatrix::arrayToDomainBlockVector(x_raw,x);
431 
432  // obtain number of Arnoldi update iterations actually taken
433  nIterations_ = dprob.GetIter();
434 
435  // compute residual norm
436  BlockVector r(x);
437  Real* Ax_raw = new Real[nrows];
438  A.multMv(x_raw,Ax_raw);
439  WrappedMatrix::arrayToDomainBlockVector(Ax_raw,r);
440  r.axpy(-lambda,x);
441  const Real r_norm = r.two_norm();
442 
443  // print verbosity information
444  if (verbosity_level_ > 0)
445  {
446  if (verbosity_level_ > 1)
447  {
448  // print some information about the problem
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;
456  }
457  std::cout << blank_ << "Result ("
458  << "#iterations = " << nIterations_ << ", "
459  << "║residual║_2 = " << r_norm << "): "
460  << "λ = " << lambda << std::endl;
461  if (verbosity_level_ > 2)
462  {
463  // print approximated eigenvector via DUNE-ISTL I/O methods
464  Dune::printvector(std::cout,x,blank_+"x",blank_+"row");
465  }
466  }
467 
468  // free dynamically allocated memory
469  delete[] Ax_raw;
470  delete[] ev;
471  }
472 
484  inline void computeSymCond2 (const Real& epsilon, Real& cond_2) const
485  {
486  // print verbosity information
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;
491 
492  // use type ArPackPlusPlus_BCRSMatrixWrapper to store matrix information
493  // and to perform the product A*v (LU decomposition is not used)
494  typedef ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
495  WrappedMatrix A(m_);
496 
497  // get number of rows and columns in A
498  const int nrows = A.nrows();
499  const int ncols = A.ncols();
500 
501  // assert that A is square
502  if (nrows != ncols)
503  DUNE_THROW(Dune::ISTLError,"Matrix is not square ("
504  << nrows << "x" << ncols << ").");
505 
506  // allocate memory for variables, set parameters
507  const int nev = 2; // Number of eigenvalues to compute
508  const int ncv = 20; // Number of Arnoldi vectors generated at each iteration (0 == auto)
509  const Real tol = epsilon; // Stopping tolerance (relative accuracy of Ritz values) (0 == machine precision)
510  const int maxit = nIterationsMax_*nev; // Maximum number of Arnoldi update iterations allowed (0 == 100*nev)
511  Real* ev = new Real[nev]; // Computed eigenvalues of A
512  const bool ivec = true; // Flag deciding if eigenvectors shall be determined
513  int nconv; // Number of converged eigenvalues
514 
515  // define what we need: eigenvalues from both ends of the spectrum
516  char which[] = "BE";
517  ARSymStdEig<Real,WrappedMatrix>
518  dprob(nrows, nev, &A, &WrappedMatrix::multMv, which, ncv, tol, maxit);
519 
520  // set ARPACK verbosity mode if requested
521  if (verbosity_level_ > 3) dprob.Trace();
522 
523  // find eigenvalues and eigenvectors of A, obtain the eigenvalues
524  nconv = dprob.Eigenvalues(ev,ivec);
525 
526  // obtain approximated dominant and least dominant eigenvalue of A
527  const Real& lambda_max = ev[nev-1];
528  const Real& lambda_min = ev[0];
529 
530  // obtain associated approximated eigenvectors of A
531  Real* x_max_raw = dprob.RawEigenvector(nev-1);
532  Real* x_min_raw = dprob.RawEigenvector(0);
533 
534  // obtain approximated spectral condition number of A
535  cond_2 = std::abs(lambda_max / lambda_min);
536 
537  // obtain number of Arnoldi update iterations actually taken
538  nIterations_ = dprob.GetIter();
539 
540  // compute each residual norm
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)
548  {
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);
551  }
552  r_max_norm = std::sqrt(r_max_norm);
553  r_min_norm = std::sqrt(r_min_norm);
554 
555  // print verbosity information
556  if (verbosity_level_ > 0)
557  {
558  if (verbosity_level_ > 1)
559  {
560  // print some information about the problem
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;
572  }
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;
579  }
580 
581  // free dynamically allocated memory
582  delete[] Ax_min_raw;
583  delete[] Ax_max_raw;
584  delete[] ev;
585  }
586 
600  inline void computeNonSymMax (const Real& epsilon,
601  BlockVector& x, Real& sigma) const
602  {
603  // print verbosity information
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;
608 
609  // use type ArPackPlusPlus_BCRSMatrixWrapper to store matrix information
610  // and to perform the product A^T*A*v (LU decomposition is not used)
611  typedef ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
612  WrappedMatrix A(m_);
613 
614  // get number of rows and columns in A
615  const int nrows = A.nrows();
616  const int ncols = A.ncols();
617 
618  // assert that A has more rows than columns (extend code later to the opposite case!)
619  if (nrows < ncols)
620  DUNE_THROW(Dune::ISTLError,"Matrix has less rows than "
621  << "columns (" << nrows << "x" << ncols << ")."
622  << " This case is not implemented, yet.");
623 
624  // allocate memory for variables, set parameters
625  const int nev = 1; // Number of eigenvalues to compute
626  const int ncv = 20; // Number of Arnoldi vectors generated at each iteration (0 == auto)
627  const Real tol = epsilon; // Stopping tolerance (relative accuracy of Ritz values) (0 == machine precision)
628  const int maxit = nIterationsMax_*nev; // Maximum number of Arnoldi update iterations allowed (0 == 100*nev)
629  Real* ev = new Real[nev]; // Computed eigenvalues of A^T*A
630  const bool ivec = true; // Flag deciding if eigenvectors shall be determined
631  int nconv; // Number of converged eigenvalues
632 
633  // define what we need: eigenvalues with largest algebraic value
634  char which[] = "LA";
635  ARSymStdEig<Real,WrappedMatrix>
636  dprob(ncols, nev, &A, &WrappedMatrix::multMtMv, which, ncv, tol, maxit);
637 
638  // set ARPACK verbosity mode if requested
639  if (verbosity_level_ > 3) dprob.Trace();
640 
641  // find eigenvalues and eigenvectors of A^T*A, obtain the eigenvalues
642  nconv = dprob.Eigenvalues(ev,ivec);
643 
644  // obtain approximated largest eigenvalue of A^T*A
645  const Real& lambda = ev[nev-1];
646 
647  // obtain associated approximated eigenvector of A^T*A
648  Real* x_raw = dprob.RawEigenvector(nev-1);
649  WrappedMatrix::arrayToDomainBlockVector(x_raw,x);
650 
651  // obtain number of Arnoldi update iterations actually taken
652  nIterations_ = dprob.GetIter();
653 
654  // compute residual norm
655  BlockVector r(x);
656  Real* AtAx_raw = new Real[ncols];
657  A.multMtMv(x_raw,AtAx_raw);
658  WrappedMatrix::arrayToDomainBlockVector(AtAx_raw,r);
659  r.axpy(-lambda,x);
660  const Real r_norm = r.two_norm();
661 
662  // calculate largest singular value of A (note that
663  // x is right-singular / left-singular vector of A)
664  sigma = std::sqrt(lambda);
665 
666  // print verbosity information
667  if (verbosity_level_ > 0)
668  {
669  if (verbosity_level_ > 1)
670  {
671  // print some information about the problem
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;
681  }
682  std::cout << blank_ << "Result ("
683  << "#iterations = " << nIterations_ << ", "
684  << "║residual║_2 = " << r_norm << "): "
685  << "σ = " << sigma << std::endl;
686  if (verbosity_level_ > 2)
687  {
688  // print approximated right-singular / left-singular vector
689  // via DUNE-ISTL I/O methods
690  Dune::printvector(std::cout,x,blank_+"x",blank_+"row");
691  }
692  }
693 
694  // free dynamically allocated memory
695  delete[] AtAx_raw;
696  delete[] ev;
697  }
698 
712  inline void computeNonSymMin (const Real& epsilon,
713  BlockVector& x, Real& sigma) const
714  {
715  // print verbosity information
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;
720 
721  // use type ArPackPlusPlus_BCRSMatrixWrapper to store matrix information
722  // and to perform the product A^T*A*v (LU decomposition is not used)
723  typedef ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
724  WrappedMatrix A(m_);
725 
726  // get number of rows and columns in A
727  const int nrows = A.nrows();
728  const int ncols = A.ncols();
729 
730  // assert that A has more rows than columns (extend code later to the opposite case!)
731  if (nrows < ncols)
732  DUNE_THROW(Dune::ISTLError,"Matrix has less rows than "
733  << "columns (" << nrows << "x" << ncols << ")."
734  << " This case is not implemented, yet.");
735 
736  // allocate memory for variables, set parameters
737  const int nev = 1; // Number of eigenvalues to compute
738  const int ncv = 20; // Number of Arnoldi vectors generated at each iteration (0 == auto)
739  const Real tol = epsilon; // Stopping tolerance (relative accuracy of Ritz values) (0 == machine precision)
740  const int maxit = nIterationsMax_*nev; // Maximum number of Arnoldi update iterations allowed (0 == 100*nev)
741  Real* ev = new Real[nev]; // Computed eigenvalues of A^T*A
742  const bool ivec = true; // Flag deciding if eigenvectors shall be determined
743  int nconv; // Number of converged eigenvalues
744 
745  // define what we need: eigenvalues with smallest algebraic value
746  char which[] = "SA";
747  ARSymStdEig<Real,WrappedMatrix>
748  dprob(ncols, nev, &A, &WrappedMatrix::multMtMv, which, ncv, tol, maxit);
749 
750  // set ARPACK verbosity mode if requested
751  if (verbosity_level_ > 3) dprob.Trace();
752 
753  // find eigenvalues and eigenvectors of A^T*A, obtain the eigenvalues
754  nconv = dprob.Eigenvalues(ev,ivec);
755 
756  // obtain approximated smallest eigenvalue of A^T*A
757  const Real& lambda = ev[nev-1];
758 
759  // obtain associated approximated eigenvector of A^T*A
760  Real* x_raw = dprob.RawEigenvector(nev-1);
761  WrappedMatrix::arrayToDomainBlockVector(x_raw,x);
762 
763  // obtain number of Arnoldi update iterations actually taken
764  nIterations_ = dprob.GetIter();
765 
766  // compute residual norm
767  BlockVector r(x);
768  Real* AtAx_raw = new Real[ncols];
769  A.multMtMv(x_raw,AtAx_raw);
770  WrappedMatrix::arrayToDomainBlockVector(AtAx_raw,r);
771  r.axpy(-lambda,x);
772  const Real r_norm = r.two_norm();
773 
774  // calculate smallest singular value of A (note that
775  // x is right-singular / left-singular vector of A)
776  sigma = std::sqrt(lambda);
777 
778  // print verbosity information
779  if (verbosity_level_ > 0)
780  {
781  if (verbosity_level_ > 1)
782  {
783  // print some information about the problem
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;
793  }
794  std::cout << blank_ << "Result ("
795  << "#iterations = " << nIterations_ << ", "
796  << "║residual║_2 = " << r_norm << "): "
797  << "σ = " << sigma << std::endl;
798  if (verbosity_level_ > 2)
799  {
800  // print approximated right-singular / left-singular vector
801  // via DUNE-ISTL I/O methods
802  Dune::printvector(std::cout,x,blank_+"x",blank_+"row");
803  }
804  }
805 
806  // free dynamically allocated memory
807  delete[] AtAx_raw;
808  delete[] ev;
809  }
810 
821  inline void computeNonSymCond2 (const Real& epsilon, Real& cond_2) const
822  {
823  // print verbosity information
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;
828 
829  // use type ArPackPlusPlus_BCRSMatrixWrapper to store matrix information
830  // and to perform the product A^T*A*v (LU decomposition is not used)
831  typedef ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
832  WrappedMatrix A(m_);
833 
834  // get number of rows and columns in A
835  const int nrows = A.nrows();
836  const int ncols = A.ncols();
837 
838  // assert that A has more rows than columns (extend code later to the opposite case!)
839  if (nrows < ncols)
840  DUNE_THROW(Dune::ISTLError,"Matrix has less rows than "
841  << "columns (" << nrows << "x" << ncols << ")."
842  << " This case is not implemented, yet.");
843 
844  // allocate memory for variables, set parameters
845  const int nev = 2; // Number of eigenvalues to compute
846  const int ncv = 20; // Number of Arnoldi vectors generated at each iteration (0 == auto)
847  const Real tol = epsilon; // Stopping tolerance (relative accuracy of Ritz values) (0 == machine precision)
848  const int maxit = nIterationsMax_*nev; // Maximum number of Arnoldi update iterations allowed (0 == 100*nev)
849  Real* ev = new Real[nev]; // Computed eigenvalues of A^T*A
850  const bool ivec = true; // Flag deciding if eigenvectors shall be determined
851  int nconv; // Number of converged eigenvalues
852 
853  // define what we need: eigenvalues from both ends of the spectrum
854  char which[] = "BE";
855  ARSymStdEig<Real,WrappedMatrix>
856  dprob(ncols, nev, &A, &WrappedMatrix::multMtMv, which, ncv, tol, maxit);
857 
858  // set ARPACK verbosity mode if requested
859  if (verbosity_level_ > 3) dprob.Trace();
860 
861  // find eigenvalues and eigenvectors of A^T*A, obtain the eigenvalues
862  nconv = dprob.Eigenvalues(ev,ivec);
863 
864  // obtain approximated largest and smallest eigenvalue of A^T*A
865  const Real& lambda_max = ev[nev-1];
866  const Real& lambda_min = ev[0];
867 
868  // obtain associated approximated eigenvectors of A^T*A
869  Real* x_max_raw = dprob.RawEigenvector(nev-1);
870  Real* x_min_raw = dprob.RawEigenvector(0);
871 
872  // obtain number of Arnoldi update iterations actually taken
873  nIterations_ = dprob.GetIter();
874 
875  // compute each residual norm
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)
883  {
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);
886  }
887  r_max_norm = std::sqrt(r_max_norm);
888  r_min_norm = std::sqrt(r_min_norm);
889 
890  // calculate largest and smallest singular value of A
891  const Real sigma_max = std::sqrt(lambda_max);
892  const Real sigma_min = std::sqrt(lambda_min);
893 
894  // obtain approximated spectral condition number of A
895  cond_2 = sigma_max / sigma_min;
896 
897  // print verbosity information
898  if (verbosity_level_ > 0)
899  {
900  if (verbosity_level_ > 1)
901  {
902  // print some information about the problem
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;
916  }
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;
923  }
924 
925  // free dynamically allocated memory
926  delete[] AtAx_min_raw;
927  delete[] AtAx_max_raw;
928  delete[] ev;
929  }
930 
935  inline unsigned int getIterationCount () const
936  {
937  if (nIterations_ == 0)
938  DUNE_THROW(Dune::ISTLError,"No algorithm applied, yet.");
939 
940  return nIterations_;
941  }
942 
943  protected:
944  // parameters related to iterative eigenvalue algorithms
945  const BCRSMatrix& m_;
946  const unsigned int nIterationsMax_;
947 
948  // verbosity setting
949  const unsigned int verbosity_level_;
950 
951  // memory for storing temporary variables (mutable as they shall
952  // just be effectless auxiliary variables of the const apply*(...)
953  // methods)
954  mutable unsigned int nIterations_;
955 
956  // constants for printing verbosity information
957  const std::string title_;
958  const std::string blank_;
959  };
960 
961 } // namespace Dune
962 
963 
964 #endif // HAVE_ARPACKPP
965 
966 #endif // DUNE_ISTL_EIGENVALUE_ARPACKPP_HH
A vector of blocks with memory management.
Definition: bvector.hh:312
B::field_type field_type
export the type representing the field
Definition: bcrsmatrix.hh:445
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:102
Definition: basearray.hh:19
This file implements a vector space as a tensor product of a given vector space. The number of compon...
The number of blocklevels the matrix contains.
Definition: bcrsmatrix.hh:465
B::field_type field_type
export the type representing the field
Definition: bvector.hh:319
derive error class from the base class in common
Definition: istlexception.hh:16