VirtualEigenvalueSolver.hxx
1 #ifndef SELDON_FILE_VIRTUAL_EIGENVALUE_SOLVER_HXX
2 
3 #ifdef SELDON_WITH_SLEPC
4 #include <petscconf.h>
5 
6 #if defined(PETSC_USE_COMPLEX)
7 #define Petsc_Scalar complex<double>
8 #else
9 #define Petsc_Scalar double
10 #endif
11 
12 #define Petsc_Error_Code int
13 #define Petsc_Int int
14 
15 #endif
16 
17 namespace Seldon
18 {
19 
21  class AnasaziParam
22  {
23  public:
25 
30  enum {SOLVER_LOBPCG, SOLVER_BKS, SOLVER_BD};
31 
33  enum {ORTHO_DGKS, ORTHO_SVQB};
34 
35  protected:
37  int ortho_manager;
38 
40  int nb_blocks;
41 
43  int restart_number;
44 
46  int type_solver;
47 
48  public:
49  AnasaziParam();
50 
51  int GetNbBlocks() const;
52  void SetNbBlocks(int);
53  int GetNbMaximumRestarts() const;
54  void SetNbMaximumRestarts(int);
55  int GetOrthoManager() const;
56 
57  int GetEigensolverType() const;
58  void SetEigensolverType(int type);
59 
60  };
61 
62 
64  class SlepcParam
65  {
66  protected:
68  int type_solver;
69  int block_size;
70  int nstep;
71  int type_extraction;
72  int quadrature_rule;
73  int nstep_inner, nstep_outer;
74  int npoints;
75  int moment_size, block_max_size, npart;
76  double delta_rank, delta_spur;
77  int borth;
78  int double_exp;
79  int init_size;
80  int krylov_restart, restart_number, restart_add;
81  int nb_conv_vector, nb_conv_vector_proj;
82  bool non_locking_variant;
83  double restart_ratio;
84  string method;
85  int shift_type;
86 
87  public:
88  enum {SHIFT_CONSTANT, SHIFT_RAYLEIGH, SHIFT_WILKINSON};
89 
90  enum {QUADRULE_TRAPEZE, QUADRULE_CHEBY};
91 
92  enum {EXTRACT_RITZ, EXTRACT_HANKEL};
93 
94  enum {POWER, SUBSPACE, ARNOLDI, LANCZOS, KRYLOVSCHUR, GD, JD,
95  RQCG, LOBPCG, CISS, LAPACK, ARPACK, BLZPACK, TRLAN, BLOPEX, PRIMME, FEAST};
96 
97  SlepcParam();
98 
99  const char* GetEigensolverChar() const;
100  int GetEigensolverType() const;
101  void SetEigensolverType(int type);
102 
103  void SetBlockSize(int n);
104  int GetBlockSize() const;
105  void SetMaximumBlockSize(int n);
106  int GetMaximumBlockSize() const;
107 
108  void SetNumberOfSteps(int n);
109  int GetNumberOfSteps() const;
110 
111  void SetExtractionType(int);
112  int GetExtractionType() const;
113 
114  void SetQuadratureRuleType(int);
115  int GetQuadratureRuleType() const;
116 
117  void SetInnerSteps(int);
118  void SetOuterSteps(int);
119  int GetInnerSteps() const;
120  int GetOuterSteps() const;
121 
122  int GetNumberIntegrationPoints() const;
123  void SetNumberIntegrationPoints(int);
124  int GetMomentSize() const;
125  void SetMomentSize(int);
126  int GetNumberPartitions() const;
127  void SetNumberPartitions(int);
128 
129  void SetThresholdRank(double);
130  double GetThresholdRank() const;
131  void SetThresholdSpurious(double);
132  double GetThresholdSpurious() const;
133 
134  int GetBorthogonalization() const;
135  void SetBorthogonalization(int);
136  int GetDoubleExpansion() const;
137  void SetDoubleExpansion(int);
138  int GetInitialSize() const;
139  void SetInitialSize(int);
140  int GetKrylovRestart() const;
141  void SetKrylovRestart(int);
142  int GetRestartNumber() const;
143  void SetRestartNumber(int);
144  int GetRestartNumberAdd() const;
145  void SetRestartNumberAdd(int);
146 
147  int GetNumberConvergedVectors() const;
148  void SetNumberConvergedVectors(int);
149  int GetNumberConvergedVectorsProjected() const;
150  void SetNumberConvergedVectorsProjected(int);
151 
152  bool UseNonLockingVariant() const;
153  void SetNonLockingVariant(bool n);
154  double GetRestartRatio() const;
155  void SetRestartRatio(double d);
156  string GetMethod() const;
157  void SetMethod(string s);
158  int GetShiftType() const;
159  void SetShiftType(int n);
160 
161  };
162 
163 
165  class FeastParam
166  {
167  protected:
169  bool evaluate_number_eigenval;
170 
172  int nb_points_quadrature;
173 
175  int type_integration;
176 
178  double emin_interval, emax_interval;
179 
181  complex<double> center_spectrum; double radius_spectrum;
182 
184  double ratio_ellipse, angle_ellipse;
185 
186  public:
187  FeastParam();
188 
189  void EnableEstimateNumberEigenval(bool p = true);
190  bool EstimateNumberEigenval() const;
191 
192  int GetNumOfQuadraturePoints() const;
193  void SetNumOfQuadraturePoints(int);
194 
195  int GetTypeIntegration() const;
196  void SetTypeIntegration(int);
197 
198  double GetCircleRadiusSpectrum() const;
199  complex<double> GetCircleCenterSpectrum() const;
200  double GetRatioEllipseSpectrum() const;
201  double GetAngleEllipseSpectrum() const;
202 
203  double GetLowerBoundInterval() const;
204  double GetUpperBoundInterval() const;
205 
206  void SetIntervalSpectrum(double, double);
207  void SetCircleSpectrum(const complex<double>& z, double r);
208  void SetEllipseSpectrum(const complex<double>& z, double r,
209  double ratio, double teta);
210 
211  };
212 
213  template<class T0, class T1>
214  void to_complex_eigen(const T0& x, T1& y);
215 
216  template<class T0, class T1>
217  void to_complex_eigen(const complex<T0>& x, T1& y);
218 
219  template<class T0, class T1>
220  void to_complex_eigen(const T0& x, complex<T1>& y);
221 
222  template<class T0, class T1>
223  void to_complex_eigen(const complex<T0>& x, complex<T1>& y);
224 
225 
227  template<class T>
228  class EigenvalueComparisonClass
229  {
230  public:
231  virtual int CompareEigenvalue(const T& Lr, const T& Li, const T& Lr2, const T& Li2) = 0;
232 
233  };
234 
236  class GeneralEigenProblem_Base
237  {
238  public:
240 
245  enum {SMALL_EIGENVALUES, LARGE_EIGENVALUES, CENTERED_EIGENVALUES};
246 
248  enum {SORTED_REAL, SORTED_IMAG, SORTED_MODULUS, SORTED_USER};
249 
250  protected:
252  int nb_eigenvalues_wanted;
253 
255  int nb_arnoldi_vectors;
256 
258  bool automatic_selection_arnoldi_vectors;
259 
261  int type_spectrum_wanted;
262 
264  int type_sort_eigenvalues;
265 
267  double stopping_criterion;
268 
270  int nb_maximum_iterations;
271 
273  int nb_linear_solves; int display_every;
274  int print_level;
275 
277  int nb_prod;
278 
280  int n_, nglob;
281 
282 #ifdef SELDON_WITH_MPI
283  MPI_Comm comm;
286 
288  MPI_Comm rci_comm;
289 
291  MPI_Comm comm_global;
292  bool global_comm_set;
293 #endif
294 
295  public:
296  GeneralEigenProblem_Base();
297  virtual ~GeneralEigenProblem_Base();
298 
299 #ifdef SELDON_WITH_MPI
300  void SetCommunicator(const MPI_Comm& comm_);
301  MPI_Comm& GetCommunicator();
302  MPI_Comm& GetRciCommunicator();
303 
304  void SetGlobalCommunicator(const MPI_Comm& comm_);
305  MPI_Comm& GetGlobalCommunicator();
306 #endif
307 
308  int GetRankCommunicator() const;
309  int GetGlobalRankCommunicator() const;
310 
311  int GetNbAskedEigenvalues() const;
312  void SetNbAskedEigenvalues(int n);
313 
314  int GetTypeSpectrum() const;
315  int GetTypeSorting() const;
316 
317  void SetStoppingCriterion(double eps);
318  double GetStoppingCriterion() const;
319  void SetNbMaximumIterations(int n);
320  int GetNbMaximumIterations() const;
321 
322  int GetNbMatrixVectorProducts() const;
323 
324  int GetNbArnoldiVectors() const;
325  void SetNbArnoldiVectors(int n);
326 
327  int GetM() const;
328  int GetGlobalM() const;
329  int GetN() const;
330 
331  int GetPrintLevel() const;
332  void SetPrintLevel(int lvl);
333 
334  void IncrementProdMatVect();
335 
336  int GetNbLinearSolves() const;
337  void IncrementLinearSolves();
338 
339  void Init(int n);
340 
341  virtual bool IsSymmetricProblem() const;
342  virtual bool IsHermitianProblem() const;
343 
344  };
345 
346 
348  template<class T>
349  class GeneralEigenProblem : public GeneralEigenProblem_Base
350  {
351  protected:
353  T shift, shift_imag;
354 
356  EigenvalueComparisonClass<T>* compar_eigenval;
357 
358  public:
359  typedef typename ClassComplexType<T>::Tcplx Tcplx;
360  typedef typename ClassComplexType<T>::Treal Treal;
361  GeneralEigenProblem();
362 
363  T GetShiftValue() const;
364  T GetImagShiftValue() const;
365  void SetShiftValue(const T&);
366  void SetImagShiftValue(const T&);
367 
368  void GetComplexShift(const Treal&, const Treal&, Tcplx&) const;
369  void GetComplexShift(const Tcplx&, const Tcplx&, Tcplx&) const;
370 
371  void SetTypeSpectrum(int type, const T& val,
372  int type_sort = SORTED_MODULUS);
373 
374  void SetTypeSpectrum(int type, const complex<T>& val,
375  int type_sort = SORTED_MODULUS);
376 
377  void SetUserComparisonClass(EigenvalueComparisonClass<T>* ev);
378 
379 #ifdef SELDON_WITH_SLEPC
380  template<class T0>
381  static Petsc_Error_Code
382  GetComparisonEigenvalueSlepcGen(Petsc_Scalar, Petsc_Scalar, Petsc_Scalar, Petsc_Scalar, T0,
383  Petsc_Int*, void*);
384 
385  static Petsc_Error_Code
386  GetComparisonEigenvalueSlepcGen(Petsc_Scalar, Petsc_Scalar, Petsc_Scalar, Petsc_Scalar, Petsc_Scalar,
387  Petsc_Int*, void*);
388 
389  static Petsc_Error_Code
390  GetComparisonEigenvalueSlepc(Petsc_Scalar, Petsc_Scalar, Petsc_Scalar, Petsc_Scalar,
391  Petsc_Int*, void*);
392 #endif
393 
394  void FillComplexEigenvectors(int m, const Tcplx& Emid, Treal eps,
395  const Vector<Tcplx>& lambda_cplx,
396  const Matrix<Tcplx, General, ColMajor>& Ecplx,
397  Vector<T>& Lr, Vector<T>& Li,
398  Matrix<T, General, ColMajor>& E);
399 
400  virtual void DistributeEigenvectors(Matrix<T, General, ColMajor>& eigen_vec);
401 
402  protected:
403 #ifdef SELDON_WITH_MPI
404  void AssembleEigenvectors(Matrix<T, General, ColMajor>&, const Vector<int>&,
405  Vector<int>*, Vector<Vector<int> >*, int, int, int);
406 #endif
407 
408  };
409 
411 
420  template<class T>
421  class EigenProblem_Base : public GeneralEigenProblem<T>
422  {
423  public :
425 
441  enum {REGULAR_MODE, SHIFTED_MODE, IMAG_SHIFTED_MODE, INVERT_MODE,
442  BUCKLING_MODE, CAYLEY_MODE};
443 
445  enum { REAL_PART, IMAG_PART, COMPLEX_PART};
446 
447  protected :
448  typedef typename ClassComplexType<T>::Tcplx Tcplx;
449  typedef typename ClassComplexType<T>::Treal Treal;
450 
452  int eigenvalue_computation_mode;
453 
455 
458  int nb_add_eigenvalues;
459 
461 
466  bool use_cholesky;
467 
469 
473  bool diagonal_mass;
474 
477  bool complex_system; int selected_part;
478 
480  AnasaziParam anasazi_param;
481 
483  SlepcParam slepc_param;
484 
486  FeastParam feast_param;
487 
488  public :
489  EigenProblem_Base();
490 
491  // basic functions
492  int GetComputationalMode() const;
493  void SetComputationalMode(int mode);
494 
495  int GetNbAdditionalEigenvalues() const;
496  void SetNbAdditionalEigenvalues(int n);
497 
498  AnasaziParam& GetAnasaziParameters();
499  SlepcParam& GetSlepcParameters();
500  FeastParam& GetFeastParameters();
501 
502  void SetCholeskyFactoForMass(bool chol = true);
503  bool UseCholeskyFactoForMass() const;
504 
505  void SetDiagonalMass(bool diag = true);
506  bool DiagonalMass() const;
507 
508  void PrintErrorInit() const;
509 
510  // mass matrix stuff
511  virtual void ComputeDiagonalMass() = 0;
512  virtual void FactorizeDiagonalMass() = 0;
513  virtual void GetSqrtDiagonal(Vector<T>&) = 0;
514 
515  virtual void MltInvSqrtDiagonalMass(Vector<Treal>& X) = 0;
516  virtual void MltSqrtDiagonalMass(Vector<Treal>& X) = 0;
517 
518  virtual void MltInvSqrtDiagonalMass(Vector<Tcplx>& X) = 0;
519  virtual void MltSqrtDiagonalMass(Vector<Tcplx>& X) = 0;
520 
521  virtual void ComputeMassForCholesky();
522  virtual void ComputeMassMatrix();
523 
524  virtual void MltMass(const Vector<Treal>& X, Vector<Treal>& Y) = 0;
525  virtual void MltMass(const Vector<Tcplx>& X, Vector<Tcplx>& Y) = 0;
526  virtual void MltMass(const SeldonTranspose&, const Vector<Treal>& X, Vector<Treal>& Y) = 0;
527  virtual void MltMass(const SeldonTranspose&, const Vector<Tcplx>& X, Vector<Tcplx>& Y) = 0;
528 
529  // stiffness matrix stuff
530  virtual void ComputeStiffnessMatrix();
531  virtual void ComputeStiffnessMatrix(const T& a, const T& b);
532 
533  virtual void MltStiffness(const Vector<Treal>& X, Vector<Treal>& Y) = 0;
534  virtual void MltStiffness(const Vector<Tcplx>& X, Vector<Tcplx>& Y) = 0;
535 
536  virtual void MltStiffness(const T& a, const T& b,
537  const Vector<Treal>& X, Vector<Treal>& Y) = 0;
538 
539  virtual void MltStiffness(const T& a, const T& b,
540  const Vector<Tcplx>& X, Vector<Tcplx>& Y) = 0;
541 
542  virtual void MltStiffness(const SeldonTranspose&, const Vector<Treal>& X, Vector<Treal>& Y) = 0;
543  virtual void MltStiffness(const SeldonTranspose&, const Vector<Tcplx>& X, Vector<Tcplx>& Y) = 0;
544 
545  // functions to overload (factorisation of mass and/or stiffness matrix)
546  virtual void ComputeAndFactorizeStiffnessMatrix(const Treal& a, const Treal& b,
547  int real_p = COMPLEX_PART);
548 
549  virtual void ComputeAndFactorizeStiffnessMatrix(const Tcplx& a, const Tcplx& b,
550  int real_p = COMPLEX_PART);
551 
552  virtual void ComputeSolution(const Vector<Treal>& X, Vector<Treal>& Y);
553  virtual void ComputeSolution(const Vector<Tcplx>& X, Vector<Tcplx>& Y);
554 
555  virtual void ComputeSolution(const SeldonTranspose&,
556  const Vector<Treal>& X, Vector<Treal>& Y);
557 
558  virtual void ComputeSolution(const SeldonTranspose&,
559  const Vector<Tcplx>& X, Vector<Tcplx>& Y);
560 
561  virtual void FactorizeCholeskyMass();
562 
563  virtual void MltCholeskyMass(const SeldonTranspose& transA, Vector<Treal>& X);
564  virtual void MltCholeskyMass(const SeldonTranspose& transA, Vector<Tcplx>& X);
565 
566  virtual void SolveCholeskyMass(const SeldonTranspose& transA, Vector<Treal>& X);
567  virtual void SolveCholeskyMass(const SeldonTranspose& transA, Vector<Tcplx>& X);
568 
569  virtual void Clear();
570 
571  // modification of eigenvectors to take into account
572  // the use of matrix D^-1/2 K D^-1/2 or L^-1 K L^-T instead of K
573  // => eigenvectors are recovered by multiplying them by matrix D^1/2 or by L^T
574  template<class T0, class Prop, class Storage>
575  friend void ApplyScalingEigenvec(EigenProblem_Base<T0>& var,
576  Vector<T0>& eigen_values,
577  Vector<T0>& lambda_imag,
578  Matrix<T0, Prop, Storage>& eigen_vectors,
579  const T0& shiftr, const T0& shifti);
580 
581  template<class T0, class Prop, class Storage>
582  friend void ApplyScalingEigenvec(EigenProblem_Base<complex<T0> >& var,
583  Vector<complex<T0> >& eigen_values,
584  Vector<complex<T0> >& lambda_imag,
585  Matrix<complex<T0>, Prop, Storage>& eigen_vectors,
586  const complex<T0>& shiftr,
587  const complex<T0>& shifti);
588 
589  };
590 
591 
593  template<class T, class StiffValue = T,
594  class MassValue = typename ClassComplexType<T>::Treal>
595  class VirtualEigenProblem : public EigenProblem_Base<T>
596  {
597  protected:
599  VirtualMatrix<MassValue>* Mh;
600 
602  VirtualMatrix<StiffValue>* Kh;
603 
605  Vector<MassValue> sqrt_diagonal_mass;
606 
607  typedef typename ClassComplexType<T>::Tcplx Tcplx;
608  typedef typename ClassComplexType<T>::Treal Treal;
609 
610  public:
611  VirtualEigenProblem();
612 
613  void InitMatrix(VirtualMatrix<StiffValue>& K, int n = -1);
614  void InitMatrix(VirtualMatrix<StiffValue>& K, VirtualMatrix<MassValue>& M, int n = -1);
615 
616  void SetMatrix(VirtualMatrix<StiffValue>& K, VirtualMatrix<MassValue>& M);
617 
618  void SetTypeSpectrum(int type, const T& val,
619  int type_sort = EigenProblem_Base<T>::SORTED_MODULUS);
620 
621  void SetTypeSpectrum(int type, const complex<T>& val,
622  int type_sort = EigenProblem_Base<T>::SORTED_MODULUS);
623 
624  bool IsSymmetricProblem() const;
625  bool IsHermitianProblem() const;
626 
627  void ComputeDiagonalMass();
628  void FactorizeDiagonalMass();
629  void GetSqrtDiagonal(Vector<T>&);
630 
631  void MltInvSqrtDiagonalMass(Vector<Treal>& X);
632  void MltSqrtDiagonalMass(Vector<Treal>& X);
633 
634  void MltInvSqrtDiagonalMass(Vector<Tcplx>& X);
635  void MltSqrtDiagonalMass(Vector<Tcplx>& X);
636 
637  void MltMass(const Vector<Treal>& X, Vector<Treal>& Y);
638  void MltMass(const Vector<Tcplx>& X, Vector<Tcplx>& Y);
639  void MltMass(const SeldonTranspose&, const Vector<Treal>& X, Vector<Treal>& Y);
640  void MltMass(const SeldonTranspose&, const Vector<Tcplx>& X, Vector<Tcplx>& Y);
641 
642  void MltStiffness(const Vector<Treal>& X, Vector<Treal>& Y);
643  void MltStiffness(const Vector<Tcplx>& X, Vector<Tcplx>& Y);
644 
645  void MltStiffness(const T& coef_mass, const T& coef_stiff,
646  const Vector<Treal>& X, Vector<Treal>& Y);
647 
648  void MltStiffness(const T& coef_mass, const T& coef_stiff,
649  const Vector<Tcplx>& X, Vector<Tcplx>& Y);
650 
651  void MltStiffness(const SeldonTranspose&, const Vector<Treal>& X, Vector<Treal>& Y);
652  void MltStiffness(const SeldonTranspose&, const Vector<Tcplx>& X, Vector<Tcplx>& Y);
653 
654  void Clear();
655 
656  };
657 
658 
659  // sorting eigenvalues
660  template<class T, class Storage1, class Storage2>
661  void SortEigenvalues(Vector<T>& lambda_r, Vector<T>& lambda_i,
662  Matrix<T, General, Storage1>& eigen_old,
663  Matrix<T, General, Storage2>& eigen_new,
664  int type_spectrum, int type_sort,
665  const T& shift_r, const T& shift_i);
666 
667  // sorting eigenvalues
668  template<class T, class Storage1, class Storage2>
669  void SortEigenvalues(Vector<complex<T> >& lambda_r,
670  Vector<complex<T> >& lambda_i,
671  Matrix<complex<T>, General, Storage1>& eigen_old,
672  Matrix<complex<T>, General, Storage2>& eigen_new,
673  int type_spectrum, int type_sort,
674  const complex<T>& shift_r, const complex<T>& shift_i);
675 
676 
678  template<class T, class Tstiff, class Prop, class Storage,
679  class Tmass = typename ClassComplexType<T>::Treal,
680  class PropM = Symmetric, class StorageM = RowSymPacked>
681  class DenseEigenProblem : public VirtualEigenProblem<T, Tstiff, Tmass>
682  {
683  protected :
684  typedef typename ClassComplexType<T>::Tcplx Tcplx;
685  typedef typename ClassComplexType<T>::Treal Treal;
686 
688  Matrix<Treal, Prop, Storage> mat_lu_real;
689 
691  Matrix<Tcplx, Prop, Storage> mat_lu_cplx;
692 
694  Matrix<Tmass, PropM, StorageM>* Mh;
695 
697  Matrix<Tstiff, Prop, Storage>* Kh;
698 
700  Vector<int> pivot;
701 
703  Matrix<Treal, PropM, StorageM> mat_chol;
704 
706  Vector<Treal> Xchol_real, Xchol_imag;
707 
708  public :
709 
710  DenseEigenProblem();
711 
712  void InitMatrix(Matrix<Tstiff, Prop, Storage>&);
713  void InitMatrix(Matrix<Tstiff, Prop, Storage>&, Matrix<Tmass, PropM, StorageM>& );
714 
715  void ComputeDiagonalMass();
716  void FactorizeCholeskyMass();
717 
718  void MltCholeskyMass(const SeldonTranspose& transA, Vector<Treal>& X);
719  void MltCholeskyMass(const SeldonTranspose& transA, Vector<Tcplx>& X);
720 
721  void SolveCholeskyMass(const SeldonTranspose& transA, Vector<Treal>& X);
722  void SolveCholeskyMass(const SeldonTranspose& transA, Vector<Tcplx>& X);
723 
724  void MltMass(const Vector<Treal>& X, Vector<Treal>& Y);
725  void MltMass(const Vector<Tcplx>& X, Vector<Tcplx>& Y);
726  void MltMass(const SeldonTranspose&, const Vector<Treal>& X, Vector<Treal>& Y);
727  void MltMass(const SeldonTranspose&, const Vector<Tcplx>& X, Vector<Tcplx>& Y);
728 
729  void MltStiffness(const Vector<Treal>& X, Vector<Treal>& Y);
730  void MltStiffness(const Vector<Tcplx>& X, Vector<Tcplx>& Y);
731  void MltStiffness(const T& coef_mass, const T& coef_stiff,
732  const Vector<Treal>& X, Vector<Treal>& Y);
733 
734  void MltStiffness(const T& coef_mass, const T& coef_stiff,
735  const Vector<Tcplx>& X, Vector<Tcplx>& Y);
736 
737  void MltStiffness(const SeldonTranspose&, const Vector<Treal>& X, Vector<Treal>& Y);
738  void MltStiffness(const SeldonTranspose&, const Vector<Tcplx>& X, Vector<Tcplx>& Y);
739 
740  void ComputeAndFactoRealMatrix(const Treal&, const Treal& a,
741  const Treal& b, int which);
742 
743  void ComputeAndFactoRealMatrix(const Tcplx&, const Treal& a,
744  const Treal& b, int which);
745 
746  void ComputeAndFactorizeStiffnessMatrix(const Treal& a, const Treal& b,
747  int which_part =
748  EigenProblem_Base<T>::COMPLEX_PART);
749 
750  void ComputeAndFactorizeStiffnessMatrix(const Tcplx& a, const Tcplx& b,
751  int which_part =
752  EigenProblem_Base<T>::COMPLEX_PART);
753 
754  void ComputeSolution(const Vector<Treal>& X, Vector<Treal>& Y);
755  void ComputeSolution(const Vector<Tcplx>& X, Vector<Tcplx>& Y);
756 
757  void ComputeSolution(const SeldonTranspose& transA,
758  const Vector<Treal>& X, Vector<Treal>& Y);
759 
760  void ComputeSolution(const SeldonTranspose& transA,
761  const Vector<Tcplx>& X, Vector<Tcplx>& Y);
762 
763  void Clear();
764 
765  };
766 
767 
769  template<class T, class MatStiff,
770  class MatMass = Matrix<typename ClassComplexType<T>::Treal,
771  Symmetric, ArrayRowSymSparse> >
772  class SparseEigenProblem
773  : public VirtualEigenProblem<T, typename MatStiff::entry_type,
774  typename MatMass::entry_type>
775  {
776  protected :
777  typedef typename ClassComplexType<T>::Tcplx Tcplx;
778  typedef typename ClassComplexType<T>::Treal Treal;
779 
780 #ifdef SELDON_WITH_MPI
781  SparseDistributedSolver<Treal> mat_lu_real;
783 
785  SparseDistributedSolver<Tcplx> mat_lu_cplx;
786 
788  DistributedCholeskySolver<Treal> chol_facto_mass_matrix;
789 #else
790  SparseDirectSolver<Treal> mat_lu_real;
792 
794  SparseDirectSolver<Tcplx> mat_lu_cplx;
795 
797  SparseCholeskySolver<Treal> chol_facto_mass_matrix;
798 #endif
799 
801  Vector<Treal> Xchol_real, Xchol_imag;
802 
803  MatMass* Mh;
804  MatStiff* Kh;
805 
806  bool distributed; int nloc;
807  Vector<int> local_col_numbers;
808  Vector<int>* ProcSharingRows;
809  Vector<Vector<int> >* SharingRowNumbers;
810  int nodl_scalar_, nb_unknowns_scal_;
811 
812  int RetrieveLocalNumbers(MatStiff& K);
813 
814  public :
815 
816  SparseEigenProblem();
817 
818  void SelectCholeskySolver(int type);
819 
820  void InitMatrix(MatStiff&);
821  void InitMatrix(MatStiff&, MatMass&);
822 
823  void ComputeDiagonalMass();
824  void FactorizeCholeskyMass();
825 
826  template<class Storage>
827  void FactorizeCholeskyMass(Matrix<Treal, Symmetric, Storage>& M);
828 
829  template<class T0, class Prop, class Storage>
830  void FactorizeCholeskyMass(Matrix<T0, Prop, Storage>& M);
831 
832 #ifdef SELDON_WITH_MPI
833  template<class Storage>
834  void FactorizeCholeskyMass(DistributedMatrix<Treal, Symmetric, Storage>& M);
835 
836  template<class T0, class Prop, class Storage>
837  void FactorizeCholeskyMass(DistributedMatrix<T0, Prop, Storage>& M);
838 #endif
839 
840  void MltCholeskyMass(const SeldonTranspose& transA, Vector<Treal>& X);
841  void MltCholeskyMass(const SeldonTranspose& transA, Vector<Tcplx>& X);
842 
843  void SolveCholeskyMass(const SeldonTranspose& transA, Vector<Treal>& X);
844  void SolveCholeskyMass(const SeldonTranspose& transA, Vector<Tcplx>& X);
845 
846  void MltMass(const Vector<Treal>& X, Vector<Treal>& Y);
847  void MltMass(const Vector<Tcplx>& X, Vector<Tcplx>& Y);
848  void MltMass(const SeldonTranspose&, const Vector<Treal>& X, Vector<Treal>& Y);
849  void MltMass(const SeldonTranspose&, const Vector<Tcplx>& X, Vector<Tcplx>& Y);
850 
851  void MltStiffness(const Vector<Treal>& X, Vector<Treal>& Y);
852  void MltStiffness(const Vector<Tcplx>& X, Vector<Tcplx>& Y);
853  void MltStiffness(const T& coef_mass, const T& coef_stiff,
854  const Vector<Treal>& X, Vector<Treal>& Y);
855 
856  void MltStiffness(const T& coef_mass, const T& coef_stiff,
857  const Vector<Tcplx>& X, Vector<Tcplx>& Y);
858 
859  void MltStiffness(const SeldonTranspose&, const Vector<Treal>& X, Vector<Treal>& Y);
860  void MltStiffness(const SeldonTranspose&, const Vector<Tcplx>& X, Vector<Tcplx>& Y);
861 
862  void ComputeAndFactoRealMatrix(const Treal&, const Treal& a,
863  const Treal& b, int which);
864 
865  void ComputeAndFactoRealMatrix(const Tcplx&, const Treal& a,
866  const Treal& b, int which);
867 
868  void ComputeAndFactorizeStiffnessMatrix(const Treal& a, const Treal& b,
869  int which =
870  EigenProblem_Base<T>::COMPLEX_PART);
871 
872  void ComputeAndFactorizeStiffnessMatrix(const Tcplx& a, const Tcplx& b,
873  int which =
874  EigenProblem_Base<T>::COMPLEX_PART);
875 
876  void ComputeSolution(const Vector<Treal>& X, Vector<Treal>& Y);
877  void ComputeSolution(const Vector<Tcplx>& X, Vector<Tcplx>& Y);
878 
879  void ComputeSolution(const SeldonTranspose& transA,
880  const Vector<Treal>& X, Vector<Treal>& Y);
881 
882  void ComputeSolution(const SeldonTranspose& transA,
883  const Vector<Tcplx>& X, Vector<Tcplx>& Y);
884 
885  void DistributeEigenvectors(Matrix<T, General, ColMajor>& eigen_vec);
886 
887  void Clear();
888 
889  };
890 
891 
893  class TypeEigenvalueSolver
894  {
895  public :
896  enum {DEFAULT, ARPACK, ANASAZI, FEAST, SLEPC};
897 
898  static int GetDefaultSolver();
899  static int GetDefaultPolynomialSolver();
900  static int GetDefaultNonLinearSolver();
901 
902  };
903 
904  template<class T, class Prop, class Storage>
905  void GetEigenvaluesEigenvectors(EigenProblem_Base<T>& var_eig,
906  Vector<T>& lambda, Vector<T>& lambda_imag,
907  Matrix<T, Prop, Storage>& eigen_vec,
908  int type_solver = TypeEigenvalueSolver::DEFAULT);
909 
910 }
911 
912 #define SELDON_FILE_VIRTUAL_EIGENVALUE_SOLVER_HXX
913 #endif
914 
Seldon::GeneralEigenProblem_Base::GetM
int GetM() const
returns number of rows
Definition: VirtualEigenvalueSolver.cxx:359
Seldon::GeneralEigenProblem_Base::GetNbArnoldiVectors
int GetNbArnoldiVectors() const
returns the number of Arnoldi vectors to use
Definition: VirtualEigenvalueSolver.cxx:344
Seldon::EigenProblem_Base::PrintErrorInit
void PrintErrorInit() const
prints error of initialization and aborts program
Definition: EigenvalueSolver.cxx:538
Seldon::VirtualEigenProblem::SetTypeSpectrum
void SetTypeSpectrum(int type, const T &val, int type_sort=EigenProblem_Base< T >::SORTED_MODULUS)
sets which eigenvalues are searched
Definition: VirtualEigenvalueSolver.cxx:1119
Seldon::AnasaziParam::nb_blocks
int nb_blocks
number of blocks for blocked solvers
Definition: VirtualEigenvalueSolver.hxx:40
Seldon::DenseEigenProblem
computation of a few eigenvalues for dense matrices
Definition: EigenvalueSolver.hxx:330
Seldon::DenseEigenProblem::mat_lu_real
Matrix< Treal, Prop, Storage > mat_lu_real
LU factorisation of a real matrix.
Definition: VirtualEigenvalueSolver.hxx:688
Seldon::AnasaziParam::SetEigensolverType
void SetEigensolverType(int type)
sets the solver used in Anasazi
Definition: VirtualEigenvalueSolver.cxx:70
Seldon::SparseEigenProblem::mat_lu_real
SparseDirectSolver< Treal > mat_lu_real
LU factorisation of sparse matrix.
Definition: VirtualEigenvalueSolver.hxx:791
Seldon::VirtualEigenProblem::MltInvSqrtDiagonalMass
void MltInvSqrtDiagonalMass(Vector< Treal > &X)
multiplication of X by D^-1/2
Definition: VirtualEigenvalueSolver.cxx:1244
Seldon::SparseEigenProblem::MltStiffness
void MltStiffness(const Vector< Treal > &X, Vector< Treal > &Y)
multiplication by stiffness matrix
Definition: VirtualEigenvalueSolver.cxx:2899
Seldon::DenseEigenProblem::MltCholeskyMass
void MltCholeskyMass(const TransStatus &transA, Vector< T0 > &X)
computation of L X or L^T X
Definition: EigenvalueSolver.cxx:1247
Seldon::SeldonTranspose
Definition: MatrixFlag.hxx:32
Seldon::VirtualEigenProblem::sqrt_diagonal_mass
Vector< MassValue > sqrt_diagonal_mass
diagonal D^1/2 if the mass matrix is diagonal positive
Definition: VirtualEigenvalueSolver.hxx:605
Seldon::GeneralEigenProblem_Base::GeneralEigenProblem_Base
GeneralEigenProblem_Base()
default constructor
Definition: VirtualEigenvalueSolver.cxx:172
Seldon::TypeEigenvalueSolver
list of availables eigenvalue solvers
Definition: EigenvalueSolver.hxx:512
Seldon::SparseEigenProblem::chol_facto_mass_matrix
SparseCholeskySolver< Treal > chol_facto_mass_matrix
Cholesky factorisation of mass matrix if required.
Definition: VirtualEigenvalueSolver.hxx:797
Seldon::SparseEigenProblem::chol_facto_mass_matrix
SparseCholeskySolver< double > chol_facto_mass_matrix
Cholesky factorisation of mass matrix if required.
Definition: EigenvalueSolver.hxx:412
Seldon::AnasaziParam::type_solver
int type_solver
which solver ?
Definition: VirtualEigenvalueSolver.hxx:46
Seldon::GeneralEigenProblem::shift
T shift
shift sigma (if type_spectrum = centered_eigenvalues)
Definition: VirtualEigenvalueSolver.hxx:353
Seldon::SparseEigenProblem
computation of a few eigenvalues for sparse matrices
Definition: EigenvalueSolver.hxx:396
Seldon::FeastParam::ratio_ellipse
double ratio_ellipse
parameters for an ellipse
Definition: VirtualEigenvalueSolver.hxx:184
Seldon::EigenProblem_Base::EigenProblem_Base
EigenProblem_Base()
default constructor
Definition: EigenvalueSolver.cxx:15
Seldon::GeneralEigenProblem::SetUserComparisonClass
void SetUserComparisonClass(EigenvalueComparisonClass< T > *ev)
sets the class where two eigenvalues can be compared
Definition: VirtualEigenvalueSolver.cxx:573
Seldon::SparseEigenProblem::mat_lu_cplx
SparseDirectSolver< Complexe > mat_lu_cplx
factorisation of complex system
Definition: EigenvalueSolver.hxx:409
Seldon::VirtualEigenProblem::MltStiffness
void MltStiffness(const Vector< Treal > &X, Vector< Treal > &Y)
matrix vector product with stifness matrix Y = K X
Definition: VirtualEigenvalueSolver.cxx:1345
Seldon::GeneralEigenProblem::GetComplexShift
void GetComplexShift(const Treal &, const Treal &, Tcplx &) const
forms the complex shift from real and imaginary part
Definition: VirtualEigenvalueSolver.cxx:522
Seldon::EigenProblem_Base::eigenvalue_computation_mode
int eigenvalue_computation_mode
mode used to find eigenvalues (regular, shifted, Cayley, etc)
Definition: EigenvalueSolver.hxx:68
Seldon::VirtualEigenProblem::IsHermitianProblem
bool IsHermitianProblem() const
returns true if the matrix is hermitian
Definition: VirtualEigenvalueSolver.cxx:1172
Seldon::AnasaziParam::restart_number
int restart_number
restart parameter for blocked solvers
Definition: VirtualEigenvalueSolver.hxx:43
Seldon::DenseEigenProblem::MltMass
void MltMass(const Vector< Treal > &X, Vector< Treal > &Y)
matrix-vector product Y = M X where M is the mass matrix
Definition: VirtualEigenvalueSolver.cxx:1993
Seldon::SparseEigenProblem::FactorizeCholeskyMass
void FactorizeCholeskyMass()
computes Cholesky factorisation of M from matrix M
Definition: EigenvalueSolver.cxx:1481
Seldon::GeneralEigenProblem_Base::type_sort_eigenvalues
int type_sort_eigenvalues
large eigenvalues because of their real part, imaginary part or magnitude ?
Definition: VirtualEigenvalueSolver.hxx:264
Seldon::DistributedCholeskySolver
Definition: DistributedCholeskySolver.hxx:27
Seldon::DenseEigenProblem::ComputeSolution
void ComputeSolution(const Vector< T0 > &X, Vector< T0 > &Y)
solution of (a M + b K) Y = X
Definition: EigenvalueSolver.cxx:1410
Seldon::SlepcParam
Parameters for Slepc package.
Definition: VirtualEigenvalueSolver.hxx:64
Seldon::Vector
Definition: SeldonHeader.hxx:207
Seldon::FeastParam::SetCircleSpectrum
void SetCircleSpectrum(const complex< double > &z, double r)
sets a circle where eigenvalues are searched
Definition: VirtualEigenvalueSolver.cxx:150
Seldon::GeneralEigenProblem_Base::GetTypeSpectrum
int GetTypeSpectrum() const
returns the spectrum desired (large, small eigenvalues, etc)
Definition: VirtualEigenvalueSolver.cxx:294
Seldon::DenseEigenProblem::Xchol_real
Vector< Treal > Xchol_real
temporary vectors for Cholesky
Definition: VirtualEigenvalueSolver.hxx:706
Seldon::EigenProblem_Base::ComputeAndFactorizeStiffnessMatrix
void ComputeAndFactorizeStiffnessMatrix(const T &a, const T &b)
computation of matrix a M + b K and factorisation of this matrix
Definition: EigenvalueSolver.cxx:758
Seldon::SparseDistributedSolver
general class for direct solver
Definition: DistributedSolver.hxx:15
Seldon::GeneralEigenProblem_Base::GetNbLinearSolves
int GetNbLinearSolves() const
returns the number of linear solves
Definition: VirtualEigenvalueSolver.cxx:408
Seldon::DenseEigenProblem::mat_chol
Matrix< Treal, PropM, StorageM > mat_chol
Cholesky factorisation of mass matrix.
Definition: VirtualEigenvalueSolver.hxx:703
Seldon::SparseEigenProblem::Clear
void Clear()
clears memory used by the object
Definition: EigenvalueSolver.cxx:1913
Seldon::EigenProblem_Base::GetComputationalMode
int GetComputationalMode() const
returns the spectral transformation used for evaluation of eigenvalues
Definition: EigenvalueSolver.cxx:133
Seldon::GeneralEigenProblem::compar_eigenval
EigenvalueComparisonClass< T > * compar_eigenval
class for comparing eigenvalues
Definition: VirtualEigenvalueSolver.hxx:356
Seldon::SortEigenvalues
void SortEigenvalues(Vector< T > &lambda_r, Vector< T > &lambda_i, Matrix< T, General, Storage1, Alloc1 > &eigen_old, Matrix< T, General, Storage2, Alloc2 > &eigen_new, int type_spectrum, int type_sort, const T &shift_r, const T &shift_i)
sorting eigenvalues
Definition: EigenvalueSolver.cxx:1040
Seldon::DenseEigenProblem::Clear
void Clear()
clearing variables used for eigenvalue resolution
Definition: EigenvalueSolver.cxx:1443
Seldon::VirtualEigenProblem::Mh
VirtualMatrix< MassValue > * Mh
mass matrix
Definition: VirtualEigenvalueSolver.hxx:599
Seldon::EigenProblem_Base::MltStiffness
void MltStiffness(const Vector< T > &X, Vector< T > &Y)
matrix vector product with stifness matrix Y = K X
Definition: EigenvalueSolver.cxx:709
Seldon::EigenProblem_Base::ComputeStiffnessMatrix
void ComputeStiffnessMatrix()
computation of stiffness matrix K
Definition: EigenvalueSolver.cxx:691
Seldon::Matrix
Definition: SeldonHeader.hxx:226
Seldon::FeastParam::nb_points_quadrature
int nb_points_quadrature
number of points in the contour
Definition: VirtualEigenvalueSolver.hxx:172
Seldon::AnasaziParam::GetNbMaximumRestarts
int GetNbMaximumRestarts() const
returns the restart parameter used in blocked solvers
Definition: VirtualEigenvalueSolver.cxx:42
Seldon::EigenProblem_Base::ComputeDiagonalMass
void ComputeDiagonalMass()
computation of diagonal of mass matrix
Definition: EigenvalueSolver.cxx:634
Seldon::EigenProblem_Base::ComputeMassMatrix
void ComputeMassMatrix()
computation of mass matrix M
Definition: EigenvalueSolver.cxx:663
Seldon::GeneralEigenProblem_Base::nb_prod
int nb_prod
number of matrix-vector products
Definition: VirtualEigenvalueSolver.hxx:277
Seldon::EigenProblem_Base::diagonal_mass
bool diagonal_mass
if true, the generalized problem is reduced to a standard problem
Definition: EigenvalueSolver.hxx:98
Seldon::SparseEigenProblem::SparseEigenProblem
SparseEigenProblem()
default constructor
Definition: EigenvalueSolver.cxx:1461
Seldon::GeneralEigenProblem_Base::nb_linear_solves
int nb_linear_solves
number of linear solves
Definition: VirtualEigenvalueSolver.hxx:273
Seldon::ArrayRowSymSparse
Definition: Storage.hxx:132
Seldon::GeneralEigenProblem_Base::GetNbMaximumIterations
int GetNbMaximumIterations() const
returns the maximal number of iterations allowed for the iterative algorithm
Definition: VirtualEigenvalueSolver.cxx:329
Seldon::SparseEigenProblem::MltCholeskyMass
void MltCholeskyMass(const TransStatus &transA, Vector< T0 > &X)
computes L X or L^T x if M = L L^T
Definition: EigenvalueSolver.cxx:1526
Seldon::DenseEigenProblem::Mh
Matrix< Tmass, PropM, StorageM > * Mh
mass matrix
Definition: VirtualEigenvalueSolver.hxx:694
Seldon::GeneralEigenProblem::SetImagShiftValue
void SetImagShiftValue(const T &)
Sets the imaginary part of shift value.
Definition: VirtualEigenvalueSolver.cxx:513
Seldon::SparseDirectSolver< Treal >
Seldon::GeneralEigenProblem_Base::SetPrintLevel
void SetPrintLevel(int lvl)
sets the level of verbosity
Definition: VirtualEigenvalueSolver.cxx:387
Seldon::SparseEigenProblem::ComputeSolution
void ComputeSolution(const Vector< T0 > &X, Vector< T0 > &Y)
solves (a M + b K) Y = X with stored factorization
Definition: EigenvalueSolver.cxx:1845
Seldon::DenseEigenProblem::Kh
Matrix< Tstiff, Prop, Storage > * Kh
stiffness matrix
Definition: VirtualEigenvalueSolver.hxx:697
Seldon::EigenProblem_Base::GetSlepcParameters
SlepcParam & GetSlepcParameters()
returns parameters specific to Slepc
Definition: VirtualEigenvalueSolver.cxx:830
Seldon::FeastParam::evaluate_number_eigenval
bool evaluate_number_eigenval
if true, we want to estimate the number of eigenvalues
Definition: VirtualEigenvalueSolver.hxx:169
Seldon::GeneralEigenProblem::GetShiftValue
T GetShiftValue() const
returns the shift value used
Definition: VirtualEigenvalueSolver.cxx:483
Seldon::GeneralEigenProblem::SetTypeSpectrum
void SetTypeSpectrum(int type, const T &val, int type_sort=SORTED_MODULUS)
sets which eigenvalues are searched
Definition: VirtualEigenvalueSolver.cxx:544
Seldon::FeastParam::GetLowerBoundInterval
double GetLowerBoundInterval() const
returns lower bound of the interval where eigenvalues are searched
Definition: VirtualEigenvalueSolver.cxx:128
Seldon::VirtualEigenProblem::GetSqrtDiagonal
void GetSqrtDiagonal(Vector< T > &)
fills D^1/2
Definition: VirtualEigenvalueSolver.cxx:1233
Seldon::TypeEigenvalueSolver::GetDefaultSolver
static int GetDefaultSolver()
returns the default eigenvalue solver (linear eigenproblem)
Definition: EigenvalueSolver.cxx:1936
Seldon::RowSymPacked
Definition: Storage.hxx:157
Seldon::EigenProblem_Base::MltInvSqrtDiagonalMass
void MltInvSqrtDiagonalMass(Vector< T0 > &X)
multiplication of X by D^-1/2
Definition: EigenvalueSolver.cxx:614
Seldon::SparseEigenProblem::Xchol_real
Vector< Treal > Xchol_real
temporary vectors for Cholesky
Definition: VirtualEigenvalueSolver.hxx:801
Seldon::EigenProblem_Base::ComputeSolution
void ComputeSolution(const Vector< T0 > &X, Vector< T0 > &Y)
solving the linear system (a M + b K) Y = X
Definition: EigenvalueSolver.cxx:781
Seldon::EigenProblem_Base::SetNbAdditionalEigenvalues
void SetNbAdditionalEigenvalues(int n)
sets the number of additional eigenvalues
Definition: EigenvalueSolver.cxx:245
Seldon::AnasaziParam::GetNbBlocks
int GetNbBlocks() const
returns the number of blocks used in blocked solvers
Definition: VirtualEigenvalueSolver.cxx:28
Seldon::SparseEigenProblem::mat_lu_cplx
SparseDirectSolver< Tcplx > mat_lu_cplx
factorisation of complex system
Definition: VirtualEigenvalueSolver.hxx:794
Seldon::EigenProblem_Base::SolveCholeskyMass
void SolveCholeskyMass(const TransStatus &transA, Vector< T > &X)
computation of L^-1 X or L^-T x if M = L L^T
Definition: EigenvalueSolver.cxx:819
Seldon::EigenProblem_Base::feast_param
FeastParam feast_param
additional parameters for Feast
Definition: VirtualEigenvalueSolver.hxx:486
Seldon::AnasaziParam::GetEigensolverType
int GetEigensolverType() const
returns the solver used in Anasazi
Definition: VirtualEigenvalueSolver.cxx:63
Seldon::GeneralEigenProblem_Base::GetRankCommunicator
int GetRankCommunicator() const
returns rank for the solver communicator
Definition: VirtualEigenvalueSolver.cxx:254
Seldon::GeneralEigenProblem_Base
Base class to solve a general eigenvalue problem (linear or not)
Definition: VirtualEigenvalueSolver.hxx:236
Seldon::EigenProblem_Base::Clear
void Clear()
memory release
Definition: EigenvalueSolver.cxx:827
Seldon::EigenProblem_Base::use_cholesky
bool use_cholesky
if true, the generalized problem is reduced to a standard problem
Definition: EigenvalueSolver.hxx:91
Seldon::DenseEigenProblem::mat_lu_cplx
Matrix< Tcplx, Prop, Storage > mat_lu_cplx
LU factorisation of a complex matrix.
Definition: VirtualEigenvalueSolver.hxx:691
Seldon::EigenProblem_Base::FactorizeDiagonalMass
void FactorizeDiagonalMass()
computation of D^1/2 from D
Definition: EigenvalueSolver.cxx:604
Seldon::FeastParam::emin_interval
double emin_interval
interval where eigenvalues are searched (real symmetric or hermitian)
Definition: VirtualEigenvalueSolver.hxx:178
Seldon::SparseEigenProblem::ComputeAndFactorizeStiffnessMatrix
void ComputeAndFactorizeStiffnessMatrix(const T &a, const T &b)
computes and factorizes a M + b K where M is the mass matrix and K the stiffness matrix
Definition: EigenvalueSolver.cxx:1635
Seldon::GeneralEigenProblem_Base::IncrementProdMatVect
void IncrementProdMatVect()
increment of the number of matrix vector products
Definition: VirtualEigenvalueSolver.cxx:401
Seldon::EigenProblem_Base::GetAnasaziParameters
AnasaziParam & GetAnasaziParameters()
returns parameters specific to Anasazi
Definition: VirtualEigenvalueSolver.cxx:822
Seldon::SparseEigenProblem::MltMass
void MltMass(const Vector< Treal > &X, Vector< Treal > &Y)
multiplication by mass matrix
Definition: VirtualEigenvalueSolver.cxx:2772
Seldon::FeastParam::center_spectrum
complex< double > center_spectrum
disk where eigenvalues are searched (non-symmetric)
Definition: VirtualEigenvalueSolver.hxx:181
Seldon::VirtualEigenProblem::Clear
void Clear()
memory release
Definition: VirtualEigenvalueSolver.cxx:1449
Seldon::EigenvalueComparisonClass
base class for setting comparison between eigenvalues
Definition: VirtualEigenvalueSolver.hxx:228
Seldon::General
Definition: Properties.hxx:26
Seldon::SparseCholeskySolver< Treal >
Seldon::AnasaziParam::GetOrthoManager
int GetOrthoManager() const
returns orthogonalization manager set in Anasazi
Definition: VirtualEigenvalueSolver.cxx:56
Seldon::EigenProblem_Base::MltMass
void MltMass(const Vector< T > &X, Vector< T > &Y)
matrix vector product with mass matrix Y = M X
Definition: EigenvalueSolver.cxx:672
Seldon::AnasaziParam
Parameters for Anasazi package.
Definition: VirtualEigenvalueSolver.hxx:21
Seldon::SlepcParam::type_solver
int type_solver
which solver ?
Definition: VirtualEigenvalueSolver.hxx:68
Seldon::EigenProblem_Base::SetCholeskyFactoForMass
void SetCholeskyFactoForMass(bool chol=true)
indicates the use of Cholesky factorisation in order to solve a standard eigenvalue problem instead o...
Definition: EigenvalueSolver.cxx:387
Seldon::DenseEigenProblem::ComputeDiagonalMass
void ComputeDiagonalMass()
computation of diagonal of mass matrix
Definition: VirtualEigenvalueSolver.cxx:1890
Seldon::GeneralEigenProblem_Base::IncrementLinearSolves
void IncrementLinearSolves()
increments the number of linear solves
Definition: VirtualEigenvalueSolver.cxx:415
Seldon::FeastParam::FeastParam
FeastParam()
default constructor
Definition: VirtualEigenvalueSolver.cxx:114
Seldon::EigenProblem_Base::DiagonalMass
bool DiagonalMass() const
returns true if the mass matrix is diagonal
Definition: EigenvalueSolver.cxx:412
Seldon::VirtualEigenProblem::IsSymmetricProblem
bool IsSymmetricProblem() const
returns true if the matrix is symmetric
Definition: VirtualEigenvalueSolver.cxx:1151
Seldon::EigenProblem_Base::anasazi_param
AnasaziParam anasazi_param
additional parameters for Anasazi
Definition: VirtualEigenvalueSolver.hxx:480
Seldon::TypeEigenvalueSolver::GetDefaultPolynomialSolver
static int GetDefaultPolynomialSolver()
returns the default eigenvalue solver (polynomial eigenproblem)
Definition: VirtualEigenvalueSolver.cxx:3649
Seldon::GeneralEigenProblem_Base::~GeneralEigenProblem_Base
virtual ~GeneralEigenProblem_Base()
destructor
Definition: VirtualEigenvalueSolver.cxx:202
Seldon::DenseEigenProblem::InitMatrix
void InitMatrix(Matrix< Tstiff, Prop, Storage > &)
Sets stiffness matrix, eigenproblem is K x = lambda x.
Definition: VirtualEigenvalueSolver.cxx:1865
Seldon::GeneralEigenProblem_Base::GetPrintLevel
int GetPrintLevel() const
returns level of verbosity
Definition: VirtualEigenvalueSolver.cxx:380
Seldon::GeneralEigenProblem_Base::GetNbAskedEigenvalues
int GetNbAskedEigenvalues() const
returns the number of eigenvalues asked by the user
Definition: VirtualEigenvalueSolver.cxx:280
Seldon::VirtualEigenProblem::MltMass
void MltMass(const Vector< Treal > &X, Vector< Treal > &Y)
matrix vector product with mass matrix Y = M X
Definition: VirtualEigenvalueSolver.cxx:1285
Seldon::FeastParam::SetIntervalSpectrum
void SetIntervalSpectrum(double, double)
sets the interval where eigenvalues are searched
Definition: VirtualEigenvalueSolver.cxx:142
Seldon::VirtualEigenProblem::VirtualEigenProblem
VirtualEigenProblem()
default constructor
Definition: VirtualEigenvalueSolver.cxx:1044
Seldon::DenseEigenProblem::DenseEigenProblem
DenseEigenProblem()
default constructor
Definition: EigenvalueSolver.cxx:1215
Seldon::EigenProblem_Base
Base class to solve an eigenvalue problem.
Definition: EigenvalueSolver.hxx:18
Seldon::SparseEigenProblem::ComputeAndFactoRealMatrix
void ComputeAndFactoRealMatrix(const Treal &, const Treal &a, const Treal &b, int which)
computes and factorizes a M + b K where M is the mass matrix and K the stiffness matrix
Definition: VirtualEigenvalueSolver.cxx:3132
Seldon::EigenProblem_Base::nb_add_eigenvalues
int nb_add_eigenvalues
additional number of eigenvalues
Definition: EigenvalueSolver.hxx:77
Seldon::VirtualEigenProblem::Kh
VirtualMatrix< StiffValue > * Kh
stiffness matrix
Definition: VirtualEigenvalueSolver.hxx:602
Seldon::FeastParam::type_integration
int type_integration
Method of integration.
Definition: VirtualEigenvalueSolver.hxx:175
Seldon::VirtualEigenProblem::MltSqrtDiagonalMass
void MltSqrtDiagonalMass(Vector< Treal > &X)
multiplication of X by D^1/2
Definition: VirtualEigenvalueSolver.cxx:1254
Seldon::DenseEigenProblem::FactorizeCholeskyMass
void FactorizeCholeskyMass()
Cholesky factorisation of mass matrix.
Definition: EigenvalueSolver.cxx:1225
Seldon::GeneralEigenProblem::GeneralEigenProblem
GeneralEigenProblem()
default constructor
Definition: VirtualEigenvalueSolver.cxx:469
Seldon::TypeEigenvalueSolver::GetDefaultNonLinearSolver
static int GetDefaultNonLinearSolver()
returns the default eigenvalue solver (non-linear eigenproblem)
Definition: VirtualEigenvalueSolver.cxx:3664
Seldon::DistributedMatrix
matrix distributed over all the processors
Definition: DistributedMatrix.hxx:506
Seldon::GeneralEigenProblem_Base::n_
int n_
size of the problem
Definition: VirtualEigenvalueSolver.hxx:280
Seldon::GeneralEigenProblem_Base::GetN
int GetN() const
returns number of columns
Definition: VirtualEigenvalueSolver.cxx:373
Seldon::AnasaziParam::SetNbBlocks
void SetNbBlocks(int)
returns the number of blocks used in blocked solvers
Definition: VirtualEigenvalueSolver.cxx:35
Seldon::EigenProblem_Base::GetNbAdditionalEigenvalues
int GetNbAdditionalEigenvalues() const
returns the additional number of eigenvalues
Definition: EigenvalueSolver.cxx:157
Seldon::EigenProblem_Base::UseCholeskyFactoForMass
bool UseCholeskyFactoForMass() const
returns true if Cholesky factorisation has to be used for mass matrix
Definition: EigenvalueSolver.cxx:396
Seldon::EigenProblem_Base::slepc_param
SlepcParam slepc_param
additional parameters for Slepc
Definition: VirtualEigenvalueSolver.hxx:483
Seldon::DenseEigenProblem::pivot
Vector< int > pivot
pivot used by the LU factorisation
Definition: EigenvalueSolver.hxx:339
Seldon::FeastParam::SetEllipseSpectrum
void SetEllipseSpectrum(const complex< double > &z, double r, double ratio, double teta)
sets an ellipse where eigenvalues are searched
Definition: VirtualEigenvalueSolver.cxx:158
Seldon::VirtualEigenProblem::FactorizeDiagonalMass
void FactorizeDiagonalMass()
computation of D^1/2 from D
Definition: VirtualEigenvalueSolver.cxx:1222
Seldon::SparseEigenProblem::DistributeEigenvectors
void DistributeEigenvectors(Matrix< T, General, ColMajor > &eigen_vec)
changes final eigenvectors if needed
Definition: VirtualEigenvalueSolver.cxx:3602
Seldon::GeneralEigenProblem_Base::automatic_selection_arnoldi_vectors
bool automatic_selection_arnoldi_vectors
if true nb_arnoldi_vectors is automatically computed
Definition: VirtualEigenvalueSolver.hxx:258
Seldon::DenseEigenProblem::MltStiffness
void MltStiffness(const Vector< Treal > &X, Vector< Treal > &Y)
matrix-vector product Y = K X where K is the stiffness matrix
Definition: VirtualEigenvalueSolver.cxx:2057
Seldon::EigenProblem_Base::SetDiagonalMass
void SetDiagonalMass(bool diag=true)
indicates that the mass matrix is diagonal
Definition: EigenvalueSolver.cxx:404
Seldon::GeneralEigenProblem::SetShiftValue
void SetShiftValue(const T &)
Sets the real part of shift value.
Definition: VirtualEigenvalueSolver.cxx:505
Seldon::GeneralEigenProblem_Base::SetNbMaximumIterations
void SetNbMaximumIterations(int n)
sets the maximal number of iterations allowed for the iterative algorithm
Definition: VirtualEigenvalueSolver.cxx:322
Seldon::GeneralEigenProblem_Base::nb_maximum_iterations
int nb_maximum_iterations
Maximal number of iterations.
Definition: VirtualEigenvalueSolver.hxx:270
Seldon::GeneralEigenProblem_Base::GetStoppingCriterion
double GetStoppingCriterion() const
returns the stopping criterion
Definition: VirtualEigenvalueSolver.cxx:315
Seldon::Symmetric
Definition: Properties.hxx:30
Seldon::VirtualEigenProblem::ComputeDiagonalMass
void ComputeDiagonalMass()
computation of diagonal of mass matrix
Definition: VirtualEigenvalueSolver.cxx:1203
Seldon::VirtualEigenProblem
base class for eigenvalue problems
Definition: VirtualEigenvalueSolver.hxx:595
Seldon::EigenProblem_Base::ApplyScalingEigenvec
friend void ApplyScalingEigenvec(EigenPb &var, Vector1 &eigen_values, Vector1 &lambda_imag, Matrix1 &eigen_vectors, const T0 &shiftr, const T0 &shifti)
modification of eigenvectors to take into account scaling by mass matrix
Definition: EigenvalueSolver.cxx:846
Seldon::SparseEigenProblem::SolveCholeskyMass
void SolveCholeskyMass(const TransStatus &transA, Vector< T0 > &X)
computes L^-1 X or L^-T x if M = L L^T
Definition: EigenvalueSolver.cxx:1536
Seldon::DenseEigenProblem::ComputeAndFactoRealMatrix
void ComputeAndFactoRealMatrix(const Treal &, const Treal &a, const Treal &b, int which)
computation and factorisation of a M + b K
Definition: VirtualEigenvalueSolver.cxx:2190
Seldon::DenseEigenProblem::SolveCholeskyMass
void SolveCholeskyMass(const TransStatus &transA, Vector< T0 > &X)
computation of L^-1 X or L^-T X
Definition: EigenvalueSolver.cxx:1258
Seldon::VirtualEigenProblem::SetMatrix
void SetMatrix(VirtualMatrix< StiffValue > &K, VirtualMatrix< MassValue > &M)
sets pointers to the stiffness and mass matrix
Definition: VirtualEigenvalueSolver.cxx:1105
Seldon::AnasaziParam::AnasaziParam
AnasaziParam()
default constructor
Definition: VirtualEigenvalueSolver.cxx:18
Seldon
Seldon namespace.
Definition: Array.cxx:24
Seldon::GeneralEigenProblem_Base::nb_eigenvalues_wanted
int nb_eigenvalues_wanted
number of eigenvalues to be computed
Definition: VirtualEigenvalueSolver.hxx:252
Seldon::SlepcParam::SlepcParam
SlepcParam()
Default constructor.
Definition: VirtualEigenvalueSolver.cxx:82
Seldon::EigenProblem_Base::MltCholeskyMass
void MltCholeskyMass(const TransStatus &transA, Vector< T > &X)
computation of L X or L^T x if M = L L^T
Definition: EigenvalueSolver.cxx:809
Seldon::EigenProblem_Base::GetFeastParameters
FeastParam & GetFeastParameters()
returns parameters specific to Feast
Definition: VirtualEigenvalueSolver.cxx:838
Seldon::EigenProblem_Base::complex_system
bool complex_system
if true consider Real( (a M + bK)^-1) or Imag( (a M + b K)^-1 ) or the whole system,...
Definition: EigenvalueSolver.hxx:128
Seldon::GeneralEigenProblem_Base::Init
void Init(int n)
initialisation of the size of the eigenvalue problem
Definition: VirtualEigenvalueSolver.cxx:434
Seldon::GeneralEigenProblem_Base::GetGlobalM
int GetGlobalM() const
returns global number of rows
Definition: VirtualEigenvalueSolver.cxx:366
Seldon::GeneralEigenProblem::GetImagShiftValue
T GetImagShiftValue() const
returns the imaginary part of shift value used
Definition: VirtualEigenvalueSolver.cxx:497
Seldon::EigenProblem_Base::ComputeMassForCholesky
void ComputeMassForCholesky()
computation of mass matrix
Definition: EigenvalueSolver.cxx:654
Seldon::DenseEigenProblem::ComputeAndFactorizeStiffnessMatrix
void ComputeAndFactorizeStiffnessMatrix(const T &a, const T &b)
computation and factorisation of a M + b K
Definition: EigenvalueSolver.cxx:1310
Seldon::GeneralEigenProblem_Base::GetNbMatrixVectorProducts
int GetNbMatrixVectorProducts() const
returns the number of matrix-vector products performed since last call to Init
Definition: VirtualEigenvalueSolver.cxx:337
Seldon::AnasaziParam::ortho_manager
int ortho_manager
orthogonalization manager
Definition: VirtualEigenvalueSolver.hxx:37
Seldon::GeneralEigenProblem_Base::nb_arnoldi_vectors
int nb_arnoldi_vectors
number of Arnoldi vectors
Definition: VirtualEigenvalueSolver.hxx:255
Seldon::GeneralEigenProblem_Base::GetGlobalRankCommunicator
int GetGlobalRankCommunicator() const
returns rank for the global communicator
Definition: VirtualEigenvalueSolver.cxx:267
Seldon::GeneralEigenProblem
class for solving a general eigenproblem with parameter T (double or complex)
Definition: VirtualEigenvalueSolver.hxx:349
Seldon::FeastParam
parameters for Feast package
Definition: VirtualEigenvalueSolver.hxx:165
Seldon::GeneralEigenProblem_Base::type_spectrum_wanted
int type_spectrum_wanted
which spectrum ? Near from Zero ? Near from Infinity ? or near from a value ?
Definition: VirtualEigenvalueSolver.hxx:261
Seldon::GeneralEigenProblem::DistributeEigenvectors
virtual void DistributeEigenvectors(Matrix< T, General, ColMajor > &eigen_vec)
changes final eigenvectors if needed
Definition: VirtualEigenvalueSolver.cxx:737
Seldon::FeastParam::GetUpperBoundInterval
double GetUpperBoundInterval() const
returns upper bound of the interval where eigenvalues are searched
Definition: VirtualEigenvalueSolver.cxx:135
Seldon::EigenProblem_Base::SetComputationalMode
void SetComputationalMode(int mode)
sets the spectral transformation used for evaluation of eigenvalues
Definition: EigenvalueSolver.cxx:141
Seldon::GeneralEigenProblem_Base::SetNbAskedEigenvalues
void SetNbAskedEigenvalues(int n)
sets the number of eigenvalues to compute
Definition: VirtualEigenvalueSolver.cxx:287
Seldon::AnasaziParam::SetNbMaximumRestarts
void SetNbMaximumRestarts(int)
sets the restart parameter used in blocked solvers
Definition: VirtualEigenvalueSolver.cxx:49
Seldon::SparseEigenProblem::SelectCholeskySolver
void SelectCholeskySolver(int type)
sets Cholesky solver to use
Definition: EigenvalueSolver.cxx:1472
Seldon::GeneralEigenProblem_Base::SetStoppingCriterion
void SetStoppingCriterion(double eps)
modifies the stopping critertion
Definition: VirtualEigenvalueSolver.cxx:308
Seldon::VirtualMatrix
Abstract base class for all matrices.
Definition: Matrix_Base.hxx:42
Seldon::GeneralEigenProblem_Base::SetNbArnoldiVectors
void SetNbArnoldiVectors(int n)
sets the number of Arnoldi vectors to use
Definition: VirtualEigenvalueSolver.cxx:351
Seldon::EigenProblem_Base::MltSqrtDiagonalMass
void MltSqrtDiagonalMass(Vector< T0 > &X)
multiplication of X by D^1/2
Definition: EigenvalueSolver.cxx:624
Seldon::SparseEigenProblem::ComputeDiagonalMass
void ComputeDiagonalMass()
computation of diagonal of mass matrix
Definition: VirtualEigenvalueSolver.cxx:2506
Seldon::GeneralEigenProblem_Base::GetTypeSorting
int GetTypeSorting() const
returns how eigenvalues are sorted (real, imaginary part or modulus)
Definition: VirtualEigenvalueSolver.cxx:301
Seldon::VirtualEigenProblem::InitMatrix
void InitMatrix(VirtualMatrix< StiffValue > &K, int n=-1)
initialization of a standard eigenvalue problem
Definition: VirtualEigenvalueSolver.cxx:1059
Seldon::EigenProblem_Base::FactorizeCholeskyMass
void FactorizeCholeskyMass()
computation of Cholesky factorisation of M from matrix M
Definition: EigenvalueSolver.cxx:799
Seldon::GeneralEigenProblem_Base::stopping_criterion
double stopping_criterion
threshold for Arpack's iterative process
Definition: VirtualEigenvalueSolver.hxx:267