You are here: Concepts>Concepts Web>Class documentation

Class documentation of Concepts

Loading...
Searching...
No Matches
arpackpp.hh
Go to the documentation of this file.
1
6#ifndef ARPACKPP_HH_
7#define ARPACKPP_HH_
8
9/* The embedded Arpack includes have virtual overloaded functions
10 For this file, we shut this warning
11 */
12#include "basics/warnings/push.h"
13#include "basics/warnings/ignore_warning_overloaded_virtual.h"
14
15#include "ARPACK.hh"
16//classes for complex (general) eigenvalue problems
17#include <eigensolver/arpackpp/arscomp.h>
18#include <eigensolver/arpackpp/argcomp.h>
19//classes for real symmetric (general) eigenvalue problems
20#include <eigensolver/arpackpp/argsym.h>
21
22// Restore warning settings
23#include "basics/warnings/pop.h"
24
25#include "operator.hh"
26
27namespace eigensolver {
28
29 using concepts::Real;
30 using concepts::Cmplx;
31
35 template<class T>
36 class ArpackStdOperatorWrapper {
37 public:
38
46
52 template<class F>
53 void multOPx(F* v, F* w) {
54 concepts::Vector<F> vec1(OP_.dimX(), v);
55 concepts::Vector<F> vec2(OP_.dimY());
56 OP_(vec1, vec2);
57 for (uint i = 0; i < OP_.dimY(); ++i)
58 w[i] = vec2[i];
59 }
60
61 private:
62 concepts::Operator<T>& OP_;
63 };
64
68 template<class T, class U, class V>
69 class ArpackOperatorWrapper {
70 public:
71
82
83 //FIXME replace for-loop with memcpy
84 //TODO what happens if dimensions don't agree? Exception??
90 template<class S>
91 void multAxp(S* v, S* w) {
92 concepts::Vector<S> vec1(A_.dimX(), v);
93 concepts::Vector<S> vec2(A_.dimY());
94 A_(vec1, vec2);
95 for (uint i = 0; i < A_.dimY(); ++i)
96 w[i] = vec2[i];
97 }
98
104 template<class S>
105 void multBxp(S* v, S* w) {
106 concepts::Vector<S> vec1(B_.dimX(), v);
107 concepts::Vector<S> vec2(B_.dimY());
108 B_(vec1, vec2);
109 for (uint i = 0; i < B_.dimY(); ++i)
110 w[i] = vec2[i];
111 }
112
118 template<class S>
119 void multOPx(S* v, S* w) {
120 concepts::Vector<S> vec1(OP_.dimX(), v);
121 concepts::Vector<S> vec2(OP_.dimY());
122 OP_(vec1, vec2);
123 for (uint i = 0; i < OP_.dimY(); ++i)
124 w[i] = vec2[i];
125 }
126
133 template<class S>
134 void multOPxRegular(S* v, S* w) {
135
136 concepts::Vector<S> vec1(OP_.dimX(), v);
137 concepts::Vector<S> vec2(A_.dimY());
138 //vec2 = A*v
139 A_(vec1, vec2);
140 //vec1 = OP * vec2 = OP * A * v
141 OP_(vec2, vec1);
142 for (uint i = 0; i < OP_.dimY(); ++i)
143 w[i] = vec1[i];
144 }
145
146 private:
148 concepts::Operator<T>& OP_;
150 concepts::Operator<U>& A_;
152 concepts::Operator<V>& B_;
153 };
154
160 template<class T>
161 class ArPackppStd: public EigenSolver<concepts::Cmplx> {
162 public:
163
174 ArPackppStd(concepts::Operator<T>& OP,
175 int kmax = 1,
176 concepts::Cmplx sigma = 0.0,
177 Real tol = 0.0,
178 int maxiter = 300) :
179 whichEV_((char *) "LM"),
180 modus_(ArPack<Real>::SHIFTINV),
181 aow_(OP),
182 arpackSolver_(
183 OP.dimX(), kmax, &aow_, &ArpackStdOperatorWrapper<T>::multOPx, sigma,
184 whichEV_, 0, tol, maxiter),
185 computed_(false) {
186 if (kmax < 1 || kmax > ((int)OP.dimX() - 2))
187 throw conceptsException(concepts::ExceptionBase("number of eigenpairs to be computed has to be larger than or equal to one and smaller than or equal to the dimension of the matrix minus 2"));
188 ;
189 }
190 ;
191
206 ArPackppStd(concepts::Operator<T>& OP,
207 int kmax = 1,
208 char* which = (char*) "LM",
209 Real tol = 0.0,
210 int maxiter = 300) :
211 whichEV_(which),
212 modus_(ArPack<Real>::SHIFTINV),
213 aow_(OP),
214 arpackSolver_(
215 OP.dimX(), kmax, &aow_, &ArpackStdOperatorWrapper<T>::multOPx,
216 whichEV_, 0, tol, maxiter),
217 computed_(false) {
218 if (kmax < 1 || kmax > ((int)OP.dimX() - 2))
219 throw conceptsException(concepts::ExceptionBase("number of eigenpairs to be computed has to be larger than or equal to one and smaller than or equal to the dimension of the matrix minus 2"));
220 ;
221 }
222 ;
223
225 virtual ~ArPackppStd() {
226 for (uint i = 0; i < this->arpackSolver_.EigenvectorsFound(); ++i)
227 delete eigenVectors_[i];
228 }
229
234 virtual const concepts::Array<Cmplx>&
235 getEV() {
236 if (!computed_)
237 compute_();
238 return eigenValues_;
239 }
240
245 virtual concepts::Array<concepts::Vector<Cmplx>*>&
246 getEF() {
247 if (!computed_)
248 compute_();
249 return eigenVectors_;
250 }
251
253 virtual uint iterations() const {
254 return iter_;
255 }
256
258 virtual uint converged() const {
259 return convEigenvalues_;
260 }
261
263 inline concepts::Array<Cmplx> getRESID() {
264 concepts::Array<Cmplx> res(this->arpackSolver_.RawResidualVector(),
265 this->arpackSolver_.GetN());
266 return res;
267 }
268
269 protected:
270 virtual std::ostream&
271 info(std::ostream& os) const {
272 os << concepts::typeOf(*this) << "\n";
273 switch (modus_) {
274 case ArPack<concepts::Real>::REGINV:
275 os << "regular inverse mode \n";
276 break;
277 case ArPack<concepts::Real>::SHIFTINV:
278 os << "shift-invert mode \n";
279 break;
280 case ArPack<concepts::Real>::NORMAL:
281 os << "normal mode \n";
282 break;
283 case ArPack<concepts::Real>::SHIFTINVIMAG:
284 os << "shift-invert mode \n";
285 break;
286 }
287 return os;
288 }
289
290 private:
292 void compute_() {
293 int numOfVecs = this->arpackSolver_.FindEigenvectors();
294 eigenValues_.resize(this->arpackSolver_.ConvergedEigenvalues());
295
296 for (int i = 0; i < this->arpackSolver_.ConvergedEigenvalues(); ++i) {
297 eigenValues_[i] = this->arpackSolver_.Eigenvalue(i);
298 }
299 eigenVectors_.resize(this->arpackSolver_.GetN());
300 for (int i = 0; i < numOfVecs; ++i)
301 eigenVectors_[i] = new concepts::Vector<concepts::Cmplx>(
302 this->arpackSolver_.GetN(), this->arpackSolver_.RawEigenvectors()
303 + i * this->arpackSolver_.GetN());
304
305 computed_ = true;
306 iter_ = arpackSolver_.GetIter();
307 convEigenvalues_ = arpackSolver_.ConvergedEigenvalues();
308 }
309 ;
310
312 char* whichEV_;
313
315 ArPack<Real>::modus modus_;
316
317 // /// ???
318 //uint n_;
319
321 ArpackStdOperatorWrapper<T> aow_;
322
324 ARCompStdEig<Real, ArpackStdOperatorWrapper<T> > arpackSolver_;
325
327 uint iter_;
328
330 uint convEigenvalues_;
331
333 bool computed_;
334
336 concepts::Array<Cmplx> eigenValues_;
337
339 concepts::Array<concepts::Vector<Cmplx>*> eigenVectors_;
340
341 };
342
349 class ArPackppSymGen: public EigenSolver<Real> {
350
351 public:
352
369 ArPackppSymGen(concepts::Operator<Real>& OP,
370 concepts::Operator<Real>& A,
371 concepts::Operator<Real>& B,
372 int kmax = 1,
373 Real sigma = 0.1,
374 Real tol = 0.0,
375 int maxiter = 300);
376
398 ArPackppSymGen(char invertType,
399 concepts::Operator<Real>& OP,
400 concepts::Operator<Real>& A,
401 int kmax = 1,
402 Real sigma = 0.1,
403 Real tol = 0.0,
404 int maxiter = 300);
405
406// /** Regular mode constructor
407// @param OP multiplication operator, OP * x = inv(B) * x
408// @param A matrix of the left hand side
409// @param B matrix of the right hand side
410// @param kmax number of eigenpairs to be computed, default 1
411// @param which defines which sort of eigenvalues should be computed, i.e.
412// eigenvalues of largest magnitude ("LM"), smallest magnitude ("SM"),
413// largest real part ("LR"), smallest real part ("SR"),
414// largest imaginary part ("LI") or smallest imaginary part ("SI"), default
415// "LM"
416// @param tol convergence tolerance for the eigenpairs, default 0.0 (is
417// replaced by LAPACK's \c DLAMCH('EPS'))
418// @param maxiter maximum number of Arnoldi iterations, default 300
419// @warning \c A has to be real and symmetric.
420// @warning \c B has to be real, symmetric and positive definite.
421// @warning \c kmax has to be larger than or equal to 1 and smaller than
422// or equal to the dimension of the matrix minus 2.
423// */
424// ArPackppSymGen(concepts::Operator<Real>& OP,
425// concepts::Operator<Real>& A,
426// concepts::Operator<Real>& B,
427// int kmax = 1,
428// char* which = (char*) "LM",
429// Real tol = 0.0,
430// int maxiter = 300);
431
433 virtual ~ArPackppSymGen() {
434 for (uint i = 0; i < this->arpackSolver_.EigenvectorsFound(); ++i)
435 delete eigenVectors_[i];
436 }
437
441 virtual const concepts::Array<Real>&
442 getEV() {
443 if (!computed_)
444 compute_();
445 return eigenValues_;
446 }
447
451 virtual concepts::Array<concepts::Vector<Real>*>&
452 getEF() {
453 if (!computed_)
454 compute_();
455 return eigenVectors_;
456 }
457
459 virtual uint iterations() const {
460 return iter_;
461 }
462
464 virtual uint converged() const {
465 return convEigenvalues_;
466 }
467
469 inline concepts::Array<Real> getRESID() {
470 concepts::Array<Real> res(this->arpackSolver_.RawResidualVector(),
471 this->arpackSolver_.GetN());
472 return res;
473 }
474
475 protected:
476 virtual std::ostream& info(std::ostream& os) const {
477 os << concepts::typeOf(*this) << std::endl ;
478 switch (modus_) {
479 case ArPack<concepts::Real>::REGINV:
480 os << "regular inverse mode" << std::endl;
481 break;
482 case ArPack<concepts::Real>::SHIFTINV:
483 os << "shift-invert mode" << std::endl;
484 break;
485 case ArPack<concepts::Real>::NORMAL:
486 os << "normal mode" << std::endl;
487 break;
488 case ArPack<concepts::Real>::SHIFTINVIMAG:
489 os << "shift-invert mode" << std::endl;
490 break;
491 }
492 return os;
493 }
494
495 private:
496
498 void compute_();
499
501 char* whichEV_;
502
504 ArPack<Real>::modus modus_;
505
507 ArpackOperatorWrapper<concepts::Real, concepts::Real, concepts::Real> aow_;
508
510 ARSymGenEig<Real, ArpackOperatorWrapper<concepts::Real, concepts::Real,
511 concepts::Real> , ArpackOperatorWrapper<concepts::Real, concepts::Real,
512 concepts::Real> > arpackSolver_;
513
514 // /// ???
515 //uint n_;
516
518 uint iter_;
519
521 uint convEigenvalues_;
522
524 bool computed_;
525
527 concepts::Array<Real> eigenValues_;
528
530 concepts::Array<concepts::Vector<Real>*> eigenVectors_;
531
532 };
533
534
540 template<class F, class G, class H = concepts::Real>
541 class ArPackppGen: public EigenSolver<concepts::Cmplx> {
542 public:
543
557 ArPackppGen(concepts::Operator<F>& OP,
558 concepts::Operator<G>& A,
559 concepts::Operator<H>& B,
560 int kmax = 1,
561 concepts::Cmplx sigma = 0.0,
562 Real tol = 0.0,
563 int maxiter = 300) :
564 whichEV_((char*) "LM"),
565 modus_(ArPack<Real>::SHIFTINV),
566 aow_(OP, A, B),
567 arpackSolver_(
568 A.dimX(),
569 kmax,
570 &aow_,
571 &ArpackOperatorWrapper<F, G, H>::multOPx,
572 &aow_,
573 &ArpackOperatorWrapper<F, G, H>::multBxp,
574 sigma,
575 whichEV_,
576 0,
577 tol,
578 maxiter),
579 computed_(false) {
580 if (kmax < 1 || kmax > ((int)OP.dimX() - 2))
581 throw conceptsException(concepts::ExceptionBase("number of eigenpairs to be computed has to be larger than or equal to one and smaller than or equal to the dimension of the matrix minus 2"));
582 ;
583 }
584 ;
585
603 ArPackppGen(concepts::Operator<F>& OP,
604 concepts::Operator<G>& A,
605 concepts::Operator<H>& B,
606 int kmax = 1,
607 char* which = (char*) "LM",
608 Real tol = 0.0,
609 int maxiter = 300) :
610 whichEV_(which),
611 modus_(ArPack<Real>::SHIFTINV),
612 aow_(OP, A, B),
613 arpackSolver_(
614 A.dimX(),
615 kmax,
616 &aow_,
617 &ArpackOperatorWrapper<F, G, H>::multOPxRegular,
618 &aow_,
619 &ArpackOperatorWrapper<F, G, H>::multBxp,
620 whichEV_,
621 0,
622 tol,
623 maxiter),
624 computed_(false) {
625 if (kmax < 1 || kmax > ((int)OP.dimX() - 2))
626 throw conceptsException(concepts::ExceptionBase("number of eigenpairs to be computed has to be larger than or equal to one and smaller than or equal to the dimension of the matrix minus 2"));
627 ;
628 }
629 ;
630
632 virtual ~ArPackppGen() {
633 for (uint i = 0; i < this->arpackSolver_.EigenvectorsFound(); ++i)
634 delete eigenVectors_[i];
635 }
636
641 virtual const concepts::Array<Cmplx>&
642 getEV() {
643 if (!computed_)
644 compute_();
645 return eigenValues_;
646 }
647
652 virtual concepts::Array<concepts::Vector<Cmplx>*>&
653 getEF() {
654 if (!computed_)
655 compute_();
656 return eigenVectors_;
657 }
658
660 virtual uint iterations() const {
661 return iter_;
662 }
663
665 virtual uint converged() const {
666 return convEigenvalues_;
667 }
668
670 inline concepts::Array<Cmplx> getRESID() {
671 concepts::Array<Cmplx> res(this->arpackSolver_.RawResidualVector(),
672 this->arpackSolver_.GetN());
673 return res;
674 }
675 protected:
676 virtual std::ostream&
677 info(std::ostream& os) const {
678 os << concepts::typeOf(*this) << std::endl;
679 switch (modus_) {
680 case ArPack<concepts::Real>::REGINV:
681 os << "regular inverse mode" << std::endl;
682 break;
683 case ArPack<concepts::Real>::SHIFTINV:
684 os << "shift-invert mode" << std::endl;
685 break;
686 case ArPack<concepts::Real>::NORMAL:
687 os << "normal mode" << std::endl;
688 break;
689 case ArPack<concepts::Real>::SHIFTINVIMAG:
690 os << "shift-invert mode" << std::endl;
691 break;
692 }
693 return os;
694 }
695 private:
696
698 void compute_() {
699 int numOfVecs = this->arpackSolver_.FindEigenvectors();
700 eigenValues_.resize(this->arpackSolver_.ConvergedEigenvalues());
701
702 for (int i = 0; i < this->arpackSolver_.ConvergedEigenvalues(); ++i) {
703 eigenValues_[i] = this->arpackSolver_.Eigenvalue(i);
704 }
705 eigenVectors_.resize(this->arpackSolver_.GetN());
706 for (int i = 0; i < numOfVecs; ++i)
707 eigenVectors_[i] = new concepts::Vector<concepts::Cmplx>(
708 this->arpackSolver_.GetN(), this->arpackSolver_.RawEigenvectors()
709 + i * this->arpackSolver_.GetN());
710
711 computed_ = true;
712 iter_ = arpackSolver_.GetIter();
713 convEigenvalues_ = arpackSolver_.ConvergedEigenvalues();
714 }
715 ;
716
718 char* whichEV_;
719
721 ArPack<Real>::modus modus_;
722
723 // /// ???
724 //uint n_;
725
727 ArpackOperatorWrapper<F, G, H> aow_;
728
730 ARCompGenEig<Real, ArpackOperatorWrapper<F, G, H> , ArpackOperatorWrapper<
731 F, G, H> > arpackSolver_;
732
734 uint iter_;
735
737 uint convEigenvalues_;
738
740 bool computed_;
741
743 concepts::Array<Cmplx> eigenValues_;
744
746 concepts::Array<concepts::Vector<Cmplx>*> eigenVectors_;
747
748 };
749
750} // namespace eigensolver
751
752#endif /* ARPACKPP_HH_ */
conceptsException
#define conceptsException(exc)
Definition exceptions.hh:344
concepts::Array::resize
void resize(const uint sz)
Definition array.hh:281
concepts::Operator::dimY
virtual const uint dimY() const
Definition compositions.hh:98
concepts::Operator::dimX
virtual const uint dimX() const
Definition compositions.hh:93
eigensolver::ArPackppGen::getRESID
concepts::Array< Cmplx > getRESID()
Returns the RESID vector.
Definition arpackpp.hh:670
eigensolver::ArPackppGen::~ArPackppGen
virtual ~ArPackppGen()
Deconstructor.
Definition arpackpp.hh:632
eigensolver::ArPackppGen::iterations
virtual uint iterations() const
Returns the actual number of Arnoldi iterations.
Definition arpackpp.hh:660
eigensolver::ArPackppGen::converged
virtual uint converged() const
Returns the number of converged eigenpairs.
Definition arpackpp.hh:665
eigensolver::ArPackppGen::ArPackppGen
ArPackppGen(concepts::Operator< F > &OP, concepts::Operator< G > &A, concepts::Operator< H > &B, int kmax=1, concepts::Cmplx sigma=0.0, Real tol=0.0, int maxiter=300)
Definition arpackpp.hh:557
eigensolver::ArPackppGen::info
virtual std::ostream & info(std::ostream &os) const
Returns information in an output stream.
Definition arpackpp.hh:677
eigensolver::ArPackppGen::ArPackppGen
ArPackppGen(concepts::Operator< F > &OP, concepts::Operator< G > &A, concepts::Operator< H > &B, int kmax=1, char *which=(char *) "LM", Real tol=0.0, int maxiter=300)
Definition arpackpp.hh:603
eigensolver::ArPackppGen::getEV
virtual const concepts::Array< Cmplx > & getEV()
Definition arpackpp.hh:642
eigensolver::ArPackppGen::getEF
virtual concepts::Array< concepts::Vector< Cmplx > * > & getEF()
Definition arpackpp.hh:653
eigensolver::ArPackppStd::converged
virtual uint converged() const
Returns the number of converged eigenpairs.
Definition arpackpp.hh:258
eigensolver::ArPackppStd::~ArPackppStd
virtual ~ArPackppStd()
Deconstructor.
Definition arpackpp.hh:225
eigensolver::ArPackppStd::ArPackppStd
ArPackppStd(concepts::Operator< T > &OP, int kmax=1, concepts::Cmplx sigma=0.0, Real tol=0.0, int maxiter=300)
Definition arpackpp.hh:174
eigensolver::ArPackppStd::ArPackppStd
ArPackppStd(concepts::Operator< T > &OP, int kmax=1, char *which=(char *) "LM", Real tol=0.0, int maxiter=300)
Definition arpackpp.hh:206
eigensolver::ArPackppStd::getEF
virtual concepts::Array< concepts::Vector< Cmplx > * > & getEF()
Definition arpackpp.hh:246
eigensolver::ArPackppStd::getEV
virtual const concepts::Array< Cmplx > & getEV()
Definition arpackpp.hh:235
eigensolver::ArPackppStd::getRESID
concepts::Array< Cmplx > getRESID()
Returns the RESID vector.
Definition arpackpp.hh:263
eigensolver::ArPackppStd::iterations
virtual uint iterations() const
Returns the actual number of Arnoldi iterations.
Definition arpackpp.hh:253
eigensolver::ArPackppStd::info
virtual std::ostream & info(std::ostream &os) const
Returns information in an output stream.
Definition arpackpp.hh:271
eigensolver::ArPackppSymGen::converged
virtual uint converged() const
Returns the number of converged eigenpairs.
Definition arpackpp.hh:464
eigensolver::ArPackppSymGen::info
virtual std::ostream & info(std::ostream &os) const
Returns information in an output stream.
Definition arpackpp.hh:476
eigensolver::ArPackppSymGen::iterations
virtual uint iterations() const
Returns the actual number of Arnoldi iterations.
Definition arpackpp.hh:459
eigensolver::ArPackppSymGen::~ArPackppSymGen
virtual ~ArPackppSymGen()
Deconstructor.
Definition arpackpp.hh:433
eigensolver::ArPackppSymGen::getRESID
concepts::Array< Real > getRESID()
Returns the RESID vector.
Definition arpackpp.hh:469
eigensolver::ArPackppSymGen::getEV
virtual const concepts::Array< Real > & getEV()
Definition arpackpp.hh:442
eigensolver::ArPackppSymGen::ArPackppSymGen
ArPackppSymGen(char invertType, concepts::Operator< Real > &OP, concepts::Operator< Real > &A, int kmax=1, Real sigma=0.1, Real tol=0.0, int maxiter=300)
eigensolver::ArPackppSymGen::getEF
virtual concepts::Array< concepts::Vector< Real > * > & getEF()
Definition arpackpp.hh:452
eigensolver::ArPackppSymGen::ArPackppSymGen
ArPackppSymGen(concepts::Operator< Real > &OP, concepts::Operator< Real > &A, concepts::Operator< Real > &B, int kmax=1, Real sigma=0.1, Real tol=0.0, int maxiter=300)
eigensolver::ArpackOperatorWrapper::ArpackOperatorWrapper
ArpackOperatorWrapper(concepts::Operator< T > &OP, concepts::Operator< U > &A, concepts::Operator< V > &B)
Definition arpackpp.hh:78
eigensolver::ArpackStdOperatorWrapper::ArpackStdOperatorWrapper
ArpackStdOperatorWrapper(concepts::Operator< T > &OP)
Definition arpackpp.hh:43
concepts::typeOf
std::string typeOf(const T &t)
Definition output.hh:43
concepts::Real
double Real
Definition typedefs.hh:39
concepts::Cmplx
std::complex< Real > Cmplx
Type for a complex number. It also depends on the setting of Real.
Definition typedefs.hh:42
Generated on Wed Sep 13 2023 21:06:18 for Concepts by doxygen 1.9.8