DUNE-FEM (unstable)

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