Reference documentation for deal.II version Git 193422c69f 2020-07-08 17:07:46 +0200
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
arpack_solver.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2010 - 2019 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_arpack_solver_h
17 #define dealii_arpack_solver_h
18 
19 #include <deal.II/base/config.h>
20 
22 
24 
25 #include <cstring>
26 
27 
28 #ifdef DEAL_II_WITH_ARPACK
29 
31 
32 
33 extern "C" void
34 dnaupd_(int * ido,
35  char * bmat,
36  unsigned int *n,
37  char * which,
38  unsigned int *nev,
39  const double *tol,
40  double * resid,
41  int * ncv,
42  double * v,
43  int * ldv,
44  int * iparam,
45  int * ipntr,
46  double * workd,
47  double * workl,
48  int * lworkl,
49  int * info);
50 
51 extern "C" void
52 dsaupd_(int * ido,
53  char * bmat,
54  unsigned int *n,
55  char * which,
56  unsigned int *nev,
57  double * tol,
58  double * resid,
59  int * ncv,
60  double * v,
61  int * ldv,
62  int * iparam,
63  int * ipntr,
64  double * workd,
65  double * workl,
66  int * lworkl,
67  int * info);
68 
69 extern "C" void
70 dneupd_(int * rvec,
71  char * howmany,
72  int * select,
73  double * d,
74  double * di,
75  double * z,
76  int * ldz,
77  double * sigmar,
78  double * sigmai,
79  double * workev,
80  char * bmat,
81  unsigned int *n,
82  char * which,
83  unsigned int *nev,
84  double * tol,
85  double * resid,
86  int * ncv,
87  double * v,
88  int * ldv,
89  int * iparam,
90  int * ipntr,
91  double * workd,
92  double * workl,
93  int * lworkl,
94  int * info);
95 
96 extern "C" void
97 dseupd_(int * rvec,
98  char * howmany,
99  int * select,
100  double * d,
101  double * z,
102  int * ldz,
103  double * sigmar,
104  char * bmat,
105  unsigned int *n,
106  char * which,
107  unsigned int *nev,
108  double * tol,
109  double * resid,
110  int * ncv,
111  double * v,
112  int * ldv,
113  int * iparam,
114  int * ipntr,
115  double * workd,
116  double * workl,
117  int * lworkl,
118  int * info);
119 
167 class ArpackSolver : public Subscriptor
168 {
169 public:
174 
175 
181  {
221  };
222 
227  {
233  explicit AdditionalData(
234  const unsigned int number_of_arnoldi_vectors = 15,
236  const bool symmetric = false);
237 
243  const unsigned int number_of_arnoldi_vectors;
244 
249 
253  const bool symmetric;
254  };
255 
259  SolverControl &
260  control() const;
261 
266  const AdditionalData &data = AdditionalData());
267 
271  template <typename VectorType>
272  void
273  set_initial_vector(const VectorType &vec);
274 
280  void
281  set_shift(const std::complex<double> sigma);
282 
332  template <typename VectorType,
333  typename MatrixType1,
334  typename MatrixType2,
335  typename INVERSE>
336  void
337  solve(const MatrixType1 & A,
338  const MatrixType2 & B,
339  const INVERSE & inverse,
340  std::vector<std::complex<double>> &eigenvalues,
341  std::vector<VectorType> & eigenvectors,
342  const unsigned int n_eigenvalues = 0);
343 
344 protected:
350 
355 
360  std::vector<double> resid;
361 
365  double sigmar;
366 
370  double sigmai;
371 
372 
373 private:
378  int,
379  int,
380  << "Number of wanted eigenvalues " << arg1
381  << " is larger that the size of the matrix " << arg2);
382 
384  int,
385  int,
386  << "Number of wanted eigenvalues " << arg1
387  << " is larger that the size of eigenvectors " << arg2);
388 
391  int,
392  int,
393  << "To store the real and complex parts of " << arg1
394  << " eigenvectors in real-valued vectors, their size (currently set to "
395  << arg2 << ") should be greater than or equal to " << arg1 + 1);
396 
398  int,
399  int,
400  << "Number of wanted eigenvalues " << arg1
401  << " is larger that the size of eigenvalues " << arg2);
402 
404  int,
405  int,
406  << "Number of Arnoldi vectors " << arg1
407  << " is larger that the size of the matrix " << arg2);
408 
410  int,
411  int,
412  << "Number of Arnoldi vectors " << arg1
413  << " is too small to obtain " << arg2 << " eigenvalues");
414 
416  int,
417  << "This ido " << arg1
418  << " is not supported. Check documentation of ARPACK");
419 
421  int,
422  << "This mode " << arg1
423  << " is not supported. Check documentation of ARPACK");
424 
426  int,
427  << "Error with dsaupd, info " << arg1
428  << ". Check documentation of ARPACK");
429 
431  int,
432  << "Error with dnaupd, info " << arg1
433  << ". Check documentation of ARPACK");
434 
436  int,
437  << "Error with dseupd, info " << arg1
438  << ". Check documentation of ARPACK");
439 
441  int,
442  << "Error with dneupd, info " << arg1
443  << ". Check documentation of ARPACK");
444 
446  int,
447  << "Maximum number " << arg1 << " of iterations reached.");
448 
450  "No shifts could be applied during implicit"
451  " Arnoldi update, try increasing the number of"
452  " Arnoldi vectors.");
453 };
454 
455 
457  const unsigned int number_of_arnoldi_vectors,
459  const bool symmetric)
460  : number_of_arnoldi_vectors(number_of_arnoldi_vectors)
461  , eigenvalue_of_interest(eigenvalue_of_interest)
462  , symmetric(symmetric)
463 {
464  // Check for possible options for symmetric problems
465  if (symmetric)
466  {
467  Assert(
468  eigenvalue_of_interest != largest_real_part,
469  ExcMessage(
470  "'largest real part' can only be used for non-symmetric problems!"));
471  Assert(
472  eigenvalue_of_interest != smallest_real_part,
473  ExcMessage(
474  "'smallest real part' can only be used for non-symmetric problems!"));
475  Assert(
476  eigenvalue_of_interest != largest_imaginary_part,
477  ExcMessage(
478  "'largest imaginary part' can only be used for non-symmetric problems!"));
479  Assert(
480  eigenvalue_of_interest != smallest_imaginary_part,
481  ExcMessage(
482  "'smallest imaginary part' can only be used for non-symmetric problems!"));
483  }
484  // Check for possible options for asymmetric problems
485  else
486  {
487  Assert(
488  eigenvalue_of_interest != algebraically_largest,
489  ExcMessage(
490  "'largest algebraic part' can only be used for symmetric problems!"));
491  Assert(
492  eigenvalue_of_interest != algebraically_smallest,
493  ExcMessage(
494  "'smallest algebraic part' can only be used for symmetric problems!"));
495  Assert(eigenvalue_of_interest != both_ends,
496  ExcMessage(
497  "'both ends' can only be used for symmetric problems!"));
498  }
499 }
500 
501 
503  const AdditionalData &data)
504  : solver_control(control)
505  , additional_data(data)
506  , initial_vector_provided(false)
507  , sigmar(0.0)
508  , sigmai(0.0)
509 {}
510 
511 
512 
513 inline void
514 ArpackSolver::set_shift(const std::complex<double> sigma)
515 {
516  sigmar = sigma.real();
517  sigmai = sigma.imag();
518 }
519 
520 
521 
522 template <typename VectorType>
523 inline void
525 {
527  resid.resize(vec.size());
528  for (size_type i = 0; i < vec.size(); ++i)
529  resid[i] = vec[i];
530 }
531 
532 
533 template <typename VectorType,
534  typename MatrixType1,
535  typename MatrixType2,
536  typename INVERSE>
537 inline void
538 ArpackSolver::solve(const MatrixType1 & /*system_matrix*/,
539  const MatrixType2 & mass_matrix,
540  const INVERSE & inverse,
541  std::vector<std::complex<double>> &eigenvalues,
542  std::vector<VectorType> & eigenvectors,
543  const unsigned int n_eigenvalues)
544 {
545  // Problem size
546  unsigned int n = eigenvectors[0].size();
547 
548  // Number of eigenvalues
549  const unsigned int nev_const =
550  (n_eigenvalues == 0) ? eigenvalues.size() : n_eigenvalues;
551  // nev for arpack, which might change by plus one during dneupd
552  unsigned int nev = nev_const;
553 
554  // check input sizes
556  {
557  Assert(nev <= eigenvectors.size(),
558  ArpackExcInvalidEigenvectorSize(nev, eigenvectors.size()));
559  }
560  else
561  Assert(nev + 1 <= eigenvectors.size(),
563  eigenvectors.size()));
564 
565  Assert(nev <= eigenvalues.size(),
567 
568  // check large enough problem size
570 
574 
575  // check whether we have enough Arnoldi vectors
579 
580  // ARPACK mode for dsaupd/dnaupd, here only mode 3, i.e. shift-invert mode
581  int mode = 3;
582 
583  // reverse communication parameter
584  int ido = 0;
585 
586  // 'G' generalized eigenvalue problem 'I' standard eigenvalue problem
587  char bmat[2] = "G";
588 
589  // Specify the eigenvalues of interest, possible parameters "LA" algebraically
590  // largest "SA" algebraically smallest "LM" largest magnitude "SM" smallest
591  // magnitude "LR" largest real part "SR" smallest real part "LI" largest
592  // imaginary part "SI" smallest imaginary part "BE" both ends of spectrum
593  // simultaneous.
594  char which[3];
596  {
598  std::strcpy(which, "LA");
599  break;
601  std::strcpy(which, "SA");
602  break;
603  case largest_magnitude:
604  std::strcpy(which, "LM");
605  break;
606  case smallest_magnitude:
607  std::strcpy(which, "SM");
608  break;
609  case largest_real_part:
610  std::strcpy(which, "LR");
611  break;
612  case smallest_real_part:
613  std::strcpy(which, "SR");
614  break;
616  std::strcpy(which, "LI");
617  break;
619  std::strcpy(which, "SI");
620  break;
621  case both_ends:
622  std::strcpy(which, "BE");
623  break;
624  }
625 
626  // tolerance for ARPACK
627  double tol = control().tolerance();
628 
629  // if the starting vector is used it has to be in resid
630  if (!initial_vector_provided || resid.size() != n)
631  resid.resize(n, 1.);
632 
633  // number of Arnoldi basis vectors specified
634  // in additional_data
636 
637  int ldv = n;
638  std::vector<double> v(ldv * ncv, 0.0);
639 
640  // information to the routines
641  std::vector<int> iparam(11, 0);
642 
643  iparam[0] = 1; // shift strategy
644 
645  // maximum number of iterations
646  iparam[2] = control().max_steps();
647 
648  // Set the mode of dsaupd. 1 is exact shifting, 2 is user-supplied shifts,
649  // 3 is shift-invert mode, 4 is buckling mode, 5 is Cayley mode.
650 
651  iparam[6] = mode;
652  std::vector<int> ipntr(14, 0);
653 
654  // work arrays for ARPACK
655  std::vector<double> workd(3 * n, 0.);
656  int lworkl =
657  additional_data.symmetric ? ncv * ncv + 8 * ncv : 3 * ncv * ncv + 6 * ncv;
658  std::vector<double> workl(lworkl, 0.);
659 
660  // information out of the iteration
661  int info = 1;
662 
663  while (ido != 99)
664  {
665  // call of ARPACK dsaupd/dnaupd routine
667  dsaupd_(&ido,
668  bmat,
669  &n,
670  which,
671  &nev,
672  &tol,
673  resid.data(),
674  &ncv,
675  v.data(),
676  &ldv,
677  iparam.data(),
678  ipntr.data(),
679  workd.data(),
680  workl.data(),
681  &lworkl,
682  &info);
683  else
684  dnaupd_(&ido,
685  bmat,
686  &n,
687  which,
688  &nev,
689  &tol,
690  resid.data(),
691  &ncv,
692  v.data(),
693  &ldv,
694  iparam.data(),
695  ipntr.data(),
696  workd.data(),
697  workl.data(),
698  &lworkl,
699  &info);
700 
701  if (ido == 99)
702  break;
703 
704  switch (mode)
705  {
706  case 3:
707  {
708  switch (ido)
709  {
710  case -1:
711  {
712  VectorType src, dst, tmp;
713  src.reinit(eigenvectors[0]);
714  dst.reinit(src);
715  tmp.reinit(src);
716 
717 
718  for (size_type i = 0; i < src.size(); ++i)
719  src(i) = workd[ipntr[0] - 1 + i];
720 
721  // multiplication with mass matrix M
722  mass_matrix.vmult(tmp, src);
723  // solving linear system
724  inverse.vmult(dst, tmp);
725 
726  for (size_type i = 0; i < dst.size(); ++i)
727  workd[ipntr[1] - 1 + i] = dst(i);
728  }
729  break;
730 
731  case 1:
732  {
733  VectorType src, dst, tmp, tmp2;
734  src.reinit(eigenvectors[0]);
735  dst.reinit(src);
736  tmp.reinit(src);
737  tmp2.reinit(src);
738 
739  for (size_type i = 0; i < src.size(); ++i)
740  {
741  src(i) = workd[ipntr[2] - 1 + i];
742  tmp(i) = workd[ipntr[0] - 1 + i];
743  }
744  // solving linear system
745  inverse.vmult(dst, src);
746 
747  for (size_type i = 0; i < dst.size(); ++i)
748  workd[ipntr[1] - 1 + i] = dst(i);
749  }
750  break;
751 
752  case 2:
753  {
754  VectorType src, dst;
755  src.reinit(eigenvectors[0]);
756  dst.reinit(src);
757 
758  for (size_type i = 0; i < src.size(); ++i)
759  src(i) = workd[ipntr[0] - 1 + i];
760 
761  // Multiplication with mass matrix M
762  mass_matrix.vmult(dst, src);
763 
764  for (size_type i = 0; i < dst.size(); ++i)
765  workd[ipntr[1] - 1 + i] = dst(i);
766  }
767  break;
768 
769  default:
770  Assert(false, ArpackExcArpackIdo(ido));
771  break;
772  }
773  }
774  break;
775  default:
776  Assert(false, ArpackExcArpackMode(mode));
777  break;
778  }
779  }
780 
781  // Set number of used iterations in SolverControl
782  control().check(iparam[2], 0.);
783 
784  if (info < 0)
785  {
787  {
788  Assert(false, ArpackExcArpackInfodsaupd(info));
789  }
790  else
791  Assert(false, ArpackExcArpackInfodnaupd(info));
792  }
793  else
794  {
795  // 1 - compute eigenvectors, 0 - only eigenvalues
796  int rvec = 1;
797 
798  // which eigenvectors
799  char howmany = 'A';
800 
801  std::vector<int> select(ncv, 1);
802 
803  int ldz = n;
804 
805  std::vector<double> eigenvalues_real(nev + 1, 0.);
806  std::vector<double> eigenvalues_im(nev + 1, 0.);
807 
808  // call of ARPACK dseupd/dneupd routine
810  {
811  std::vector<double> z(ldz * nev, 0.);
812  dseupd_(&rvec,
813  &howmany,
814  select.data(),
815  eigenvalues_real.data(),
816  z.data(),
817  &ldz,
818  &sigmar,
819  bmat,
820  &n,
821  which,
822  &nev,
823  &tol,
824  resid.data(),
825  &ncv,
826  v.data(),
827  &ldv,
828  iparam.data(),
829  ipntr.data(),
830  workd.data(),
831  workl.data(),
832  &lworkl,
833  &info);
834  }
835  else
836  {
837  std::vector<double> workev(3 * ncv, 0.);
838  dneupd_(&rvec,
839  &howmany,
840  select.data(),
841  eigenvalues_real.data(),
842  eigenvalues_im.data(),
843  v.data(),
844  &ldz,
845  &sigmar,
846  &sigmai,
847  workev.data(),
848  bmat,
849  &n,
850  which,
851  &nev,
852  &tol,
853  resid.data(),
854  &ncv,
855  v.data(),
856  &ldv,
857  iparam.data(),
858  ipntr.data(),
859  workd.data(),
860  workl.data(),
861  &lworkl,
862  &info);
863  }
864 
865  if (info == 1)
866  {
867  Assert(false, ArpackExcArpackInfoMaxIt(control().max_steps()));
868  }
869  else if (info == 3)
870  {
872  }
873  else if (info != 0)
874  {
876  {
877  Assert(false, ArpackExcArpackInfodseupd(info));
878  }
879  else
880  Assert(false, ArpackExcArpackInfodneupd(info));
881  }
882 
883  for (unsigned int i = 0; i < nev; ++i)
884  for (unsigned int j = 0; j < n; ++j)
885  eigenvectors[i](j) = v[i * n + j];
886 
887  for (unsigned int i = 0; i < nev_const; ++i)
888  eigenvalues[i] =
889  std::complex<double>(eigenvalues_real[i], eigenvalues_im[i]);
890  }
891 }
892 
893 
894 inline SolverControl &
896 {
897  return solver_control;
898 }
899 
901 
902 
903 #endif
904 #endif
AdditionalData(const unsigned int number_of_arnoldi_vectors=15, const WhichEigenvalues eigenvalue_of_interest=largest_magnitude, const bool symmetric=false)
SolverControl & solver_control
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:538
virtual State check(const unsigned int step, const double check_value)
std::array< std::pair< Number, Tensor< 1, dim, Number > >, std::integral_constant< int, dim >::value > eigenvectors(const SymmetricTensor< 2, dim, Number > &T, const SymmetricTensorEigenvectorMethod method=SymmetricTensorEigenvectorMethod::ql_implicit_shifts)
void set_shift(const std::complex< double > sigma)
static ::ExceptionBase & ArpackExcInvalidEigenvalueSize(int arg1, int arg2)
void dsaupd_(int *ido, char *bmat, unsigned int *n, char *which, unsigned int *nev, double *tol, double *resid, int *ncv, double *v, int *ldv, int *iparam, int *ipntr, double *workd, double *workl, int *lworkl, int *info)
static ::ExceptionBase & ArpackExcInvalidEigenvectorSizeNonsymmetric(int arg1, int arg2)
static ::ExceptionBase & ArpackExcArpackMode(int arg1)
std::vector< double > resid
static ::ExceptionBase & ExcMessage(std::string arg1)
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:515
void mass_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const double factor=1.)
Definition: l2.h:58
#define Assert(cond, exc)
Definition: exceptions.h:1403
static ::ExceptionBase & ArpackExcInvalidNumberofArnoldiVectors(int arg1, int arg2)
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:493
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:363
ArpackSolver(SolverControl &control, const AdditionalData &data=AdditionalData())
void dnaupd_(int *ido, char *bmat, unsigned int *n, char *which, unsigned int *nev, const double *tol, double *resid, int *ncv, double *v, int *ldv, int *iparam, int *ipntr, double *workd, double *workl, int *lworkl, int *info)
static ::ExceptionBase & ArpackExcSmallNumberofArnoldiVectors(int arg1, int arg2)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
static ::ExceptionBase & ArpackExcInvalidEigenvectorSize(int arg1, int arg2)
double tolerance() const
void dseupd_(int *rvec, char *howmany, int *select, double *d, double *z, int *ldz, double *sigmar, char *bmat, unsigned int *n, char *which, unsigned int *nev, double *tol, double *resid, int *ncv, double *v, int *ldv, int *iparam, int *ipntr, double *workd, double *workl, int *lworkl, int *info)
const AdditionalData additional_data
static ::ExceptionBase & ArpackExcArpackInfodneupd(int arg1)
static ::ExceptionBase & ArpackExcArpackInfoMaxIt(int arg1)
static const char A
unsigned int global_dof_index
Definition: types.h:76
static ::ExceptionBase & ArpackExcInvalidNumberofEigenvalues(int arg1, int arg2)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:362
void set_initial_vector(const VectorType &vec)
void dneupd_(int *rvec, char *howmany, int *select, double *d, double *di, double *z, int *ldz, double *sigmar, double *sigmai, double *workev, char *bmat, unsigned int *n, char *which, unsigned int *nev, double *tol, double *resid, int *ncv, double *v, int *ldv, int *iparam, int *ipntr, double *workd, double *workl, int *lworkl, int *info)
unsigned int max_steps() const
static ::ExceptionBase & ArpackExcArpackInfodseupd(int arg1)
bool initial_vector_provided
static ::ExceptionBase & ArpackExcArpackInfodsaupd(int arg1)
const unsigned int number_of_arnoldi_vectors
static ::ExceptionBase & ArpackExcArpackIdo(int arg1)
Eigenvalue vector is filled.
void solve(const MatrixType1 &A, const MatrixType2 &B, const INVERSE &inverse, std::vector< std::complex< double >> &eigenvalues, std::vector< VectorType > &eigenvectors, const unsigned int n_eigenvalues=0)
SolverControl & control() const
const WhichEigenvalues eigenvalue_of_interest
static ::ExceptionBase & ArpackExcArpackInfodnaupd(int arg1)
static ::ExceptionBase & ArpackExcArpackNoShifts()