Reference documentation for deal.II version 9.2.0
\(\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 
169 class ArpackSolver : public Subscriptor
170 {
171 public:
176 
177 
183  {
223  };
224 
229  {
235  explicit AdditionalData(
236  const unsigned int number_of_arnoldi_vectors = 15,
238  const bool symmetric = false);
239 
245  const unsigned int number_of_arnoldi_vectors;
246 
251 
255  const bool symmetric;
256  };
257 
261  SolverControl &
262  control() const;
263 
268  const AdditionalData &data = AdditionalData());
269 
273  template <typename VectorType>
274  void
275  set_initial_vector(const VectorType &vec);
276 
282  void
283  set_shift(const std::complex<double> sigma);
284 
334  template <typename VectorType,
335  typename MatrixType1,
336  typename MatrixType2,
337  typename INVERSE>
338  void
339  solve(const MatrixType1 & A,
340  const MatrixType2 & B,
341  const INVERSE & inverse,
342  std::vector<std::complex<double>> &eigenvalues,
343  std::vector<VectorType> & eigenvectors,
344  const unsigned int n_eigenvalues = 0);
345 
346 protected:
352 
357 
362  std::vector<double> resid;
363 
367  double sigmar;
368 
372  double sigmai;
373 
374 
375 private:
380  int,
381  int,
382  << "Number of wanted eigenvalues " << arg1
383  << " is larger that the size of the matrix " << arg2);
384 
386  int,
387  int,
388  << "Number of wanted eigenvalues " << arg1
389  << " is larger that the size of eigenvectors " << arg2);
390 
393  int,
394  int,
395  << "To store the real and complex parts of " << arg1
396  << " eigenvectors in real-valued vectors, their size (currently set to "
397  << arg2 << ") should be greater than or equal to " << arg1 + 1);
398 
400  int,
401  int,
402  << "Number of wanted eigenvalues " << arg1
403  << " is larger that the size of eigenvalues " << arg2);
404 
406  int,
407  int,
408  << "Number of Arnoldi vectors " << arg1
409  << " is larger that the size of the matrix " << arg2);
410 
412  int,
413  int,
414  << "Number of Arnoldi vectors " << arg1
415  << " is too small to obtain " << arg2 << " eigenvalues");
416 
418  int,
419  << "This ido " << arg1
420  << " is not supported. Check documentation of ARPACK");
421 
423  int,
424  << "This mode " << arg1
425  << " is not supported. Check documentation of ARPACK");
426 
428  int,
429  << "Error with dsaupd, info " << arg1
430  << ". Check documentation of ARPACK");
431 
433  int,
434  << "Error with dnaupd, info " << arg1
435  << ". Check documentation of ARPACK");
436 
438  int,
439  << "Error with dseupd, info " << arg1
440  << ". Check documentation of ARPACK");
441 
443  int,
444  << "Error with dneupd, info " << arg1
445  << ". Check documentation of ARPACK");
446 
448  int,
449  << "Maximum number " << arg1 << " of iterations reached.");
450 
452  "No shifts could be applied during implicit"
453  " Arnoldi update, try increasing the number of"
454  " Arnoldi vectors.");
455 };
456 
457 
459  const unsigned int number_of_arnoldi_vectors,
461  const bool symmetric)
462  : number_of_arnoldi_vectors(number_of_arnoldi_vectors)
463  , eigenvalue_of_interest(eigenvalue_of_interest)
464  , symmetric(symmetric)
465 {
466  // Check for possible options for symmetric problems
467  if (symmetric)
468  {
469  Assert(
470  eigenvalue_of_interest != largest_real_part,
471  ExcMessage(
472  "'largest real part' can only be used for non-symmetric problems!"));
473  Assert(
474  eigenvalue_of_interest != smallest_real_part,
475  ExcMessage(
476  "'smallest real part' can only be used for non-symmetric problems!"));
477  Assert(
478  eigenvalue_of_interest != largest_imaginary_part,
479  ExcMessage(
480  "'largest imaginary part' can only be used for non-symmetric problems!"));
481  Assert(
482  eigenvalue_of_interest != smallest_imaginary_part,
483  ExcMessage(
484  "'smallest imaginary part' can only be used for non-symmetric problems!"));
485  }
486  // Check for possible options for asymmetric problems
487  else
488  {
489  Assert(
490  eigenvalue_of_interest != algebraically_largest,
491  ExcMessage(
492  "'largest algebraic part' can only be used for symmetric problems!"));
493  Assert(
494  eigenvalue_of_interest != algebraically_smallest,
495  ExcMessage(
496  "'smallest algebraic part' can only be used for symmetric problems!"));
497  Assert(eigenvalue_of_interest != both_ends,
498  ExcMessage(
499  "'both ends' can only be used for symmetric problems!"));
500  }
501 }
502 
503 
505  const AdditionalData &data)
506  : solver_control(control)
507  , additional_data(data)
508  , initial_vector_provided(false)
509  , sigmar(0.0)
510  , sigmai(0.0)
511 {}
512 
513 
514 
515 inline void
516 ArpackSolver::set_shift(const std::complex<double> sigma)
517 {
518  sigmar = sigma.real();
519  sigmai = sigma.imag();
520 }
521 
522 
523 
524 template <typename VectorType>
525 inline void
527 {
529  resid.resize(vec.size());
530  for (size_type i = 0; i < vec.size(); ++i)
531  resid[i] = vec[i];
532 }
533 
534 
535 template <typename VectorType,
536  typename MatrixType1,
537  typename MatrixType2,
538  typename INVERSE>
539 inline void
540 ArpackSolver::solve(const MatrixType1 & /*system_matrix*/,
541  const MatrixType2 & mass_matrix,
542  const INVERSE & inverse,
543  std::vector<std::complex<double>> &eigenvalues,
544  std::vector<VectorType> & eigenvectors,
545  const unsigned int n_eigenvalues)
546 {
547  // Problem size
548  unsigned int n = eigenvectors[0].size();
549 
550  // Number of eigenvalues
551  const unsigned int nev_const =
552  (n_eigenvalues == 0) ? eigenvalues.size() : n_eigenvalues;
553  // nev for arpack, which might change by plus one during dneupd
554  unsigned int nev = nev_const;
555 
556  // check input sizes
558  {
559  Assert(nev <= eigenvectors.size(),
560  ArpackExcInvalidEigenvectorSize(nev, eigenvectors.size()));
561  }
562  else
563  Assert(nev + 1 <= eigenvectors.size(),
565  eigenvectors.size()));
566 
567  Assert(nev <= eigenvalues.size(),
569 
570  // check large enough problem size
572 
576 
577  // check whether we have enough Arnoldi vectors
581 
582  // ARPACK mode for dsaupd/dnaupd, here only mode 3, i.e. shift-invert mode
583  int mode = 3;
584 
585  // reverse communication parameter
586  int ido = 0;
587 
588  // 'G' generalized eigenvalue problem 'I' standard eigenvalue problem
589  char bmat[2] = "G";
590 
591  // Specify the eigenvalues of interest, possible parameters "LA" algebraically
592  // largest "SA" algebraically smallest "LM" largest magnitude "SM" smallest
593  // magnitude "LR" largest real part "SR" smallest real part "LI" largest
594  // imaginary part "SI" smallest imaginary part "BE" both ends of spectrum
595  // simultaneous.
596  char which[3];
598  {
600  std::strcpy(which, "LA");
601  break;
603  std::strcpy(which, "SA");
604  break;
605  case largest_magnitude:
606  std::strcpy(which, "LM");
607  break;
608  case smallest_magnitude:
609  std::strcpy(which, "SM");
610  break;
611  case largest_real_part:
612  std::strcpy(which, "LR");
613  break;
614  case smallest_real_part:
615  std::strcpy(which, "SR");
616  break;
618  std::strcpy(which, "LI");
619  break;
621  std::strcpy(which, "SI");
622  break;
623  case both_ends:
624  std::strcpy(which, "BE");
625  break;
626  }
627 
628  // tolerance for ARPACK
629  double tol = control().tolerance();
630 
631  // if the starting vector is used it has to be in resid
632  if (!initial_vector_provided || resid.size() != n)
633  resid.resize(n, 1.);
634 
635  // number of Arnoldi basis vectors specified
636  // in additional_data
638 
639  int ldv = n;
640  std::vector<double> v(ldv * ncv, 0.0);
641 
642  // information to the routines
643  std::vector<int> iparam(11, 0);
644 
645  iparam[0] = 1; // shift strategy
646 
647  // maximum number of iterations
648  iparam[2] = control().max_steps();
649 
650  // Set the mode of dsaupd. 1 is exact shifting, 2 is user-supplied shifts,
651  // 3 is shift-invert mode, 4 is buckling mode, 5 is Cayley mode.
652 
653  iparam[6] = mode;
654  std::vector<int> ipntr(14, 0);
655 
656  // work arrays for ARPACK
657  std::vector<double> workd(3 * n, 0.);
658  int lworkl =
659  additional_data.symmetric ? ncv * ncv + 8 * ncv : 3 * ncv * ncv + 6 * ncv;
660  std::vector<double> workl(lworkl, 0.);
661 
662  // information out of the iteration
663  int info = 1;
664 
665  while (ido != 99)
666  {
667  // call of ARPACK dsaupd/dnaupd routine
669  dsaupd_(&ido,
670  bmat,
671  &n,
672  which,
673  &nev,
674  &tol,
675  resid.data(),
676  &ncv,
677  v.data(),
678  &ldv,
679  iparam.data(),
680  ipntr.data(),
681  workd.data(),
682  workl.data(),
683  &lworkl,
684  &info);
685  else
686  dnaupd_(&ido,
687  bmat,
688  &n,
689  which,
690  &nev,
691  &tol,
692  resid.data(),
693  &ncv,
694  v.data(),
695  &ldv,
696  iparam.data(),
697  ipntr.data(),
698  workd.data(),
699  workl.data(),
700  &lworkl,
701  &info);
702 
703  if (ido == 99)
704  break;
705 
706  switch (mode)
707  {
708  case 3:
709  {
710  switch (ido)
711  {
712  case -1:
713  {
714  VectorType src, dst, tmp;
715  src.reinit(eigenvectors[0]);
716  dst.reinit(src);
717  tmp.reinit(src);
718 
719 
720  for (size_type i = 0; i < src.size(); ++i)
721  src(i) = workd[ipntr[0] - 1 + i];
722 
723  // multiplication with mass matrix M
724  mass_matrix.vmult(tmp, src);
725  // solving linear system
726  inverse.vmult(dst, tmp);
727 
728  for (size_type i = 0; i < dst.size(); ++i)
729  workd[ipntr[1] - 1 + i] = dst(i);
730  }
731  break;
732 
733  case 1:
734  {
735  VectorType src, dst, tmp, tmp2;
736  src.reinit(eigenvectors[0]);
737  dst.reinit(src);
738  tmp.reinit(src);
739  tmp2.reinit(src);
740 
741  for (size_type i = 0; i < src.size(); ++i)
742  {
743  src(i) = workd[ipntr[2] - 1 + i];
744  tmp(i) = workd[ipntr[0] - 1 + i];
745  }
746  // solving linear system
747  inverse.vmult(dst, src);
748 
749  for (size_type i = 0; i < dst.size(); ++i)
750  workd[ipntr[1] - 1 + i] = dst(i);
751  }
752  break;
753 
754  case 2:
755  {
756  VectorType src, dst;
757  src.reinit(eigenvectors[0]);
758  dst.reinit(src);
759 
760  for (size_type i = 0; i < src.size(); ++i)
761  src(i) = workd[ipntr[0] - 1 + i];
762 
763  // Multiplication with mass matrix M
764  mass_matrix.vmult(dst, src);
765 
766  for (size_type i = 0; i < dst.size(); ++i)
767  workd[ipntr[1] - 1 + i] = dst(i);
768  }
769  break;
770 
771  default:
772  Assert(false, ArpackExcArpackIdo(ido));
773  break;
774  }
775  }
776  break;
777  default:
778  Assert(false, ArpackExcArpackMode(mode));
779  break;
780  }
781  }
782 
783  // Set number of used iterations in SolverControl
784  control().check(iparam[2], 0.);
785 
786  if (info < 0)
787  {
789  {
790  Assert(false, ArpackExcArpackInfodsaupd(info));
791  }
792  else
793  Assert(false, ArpackExcArpackInfodnaupd(info));
794  }
795  else
796  {
797  // 1 - compute eigenvectors, 0 - only eigenvalues
798  int rvec = 1;
799 
800  // which eigenvectors
801  char howmany = 'A';
802 
803  std::vector<int> select(ncv, 1);
804 
805  int ldz = n;
806 
807  std::vector<double> eigenvalues_real(nev + 1, 0.);
808  std::vector<double> eigenvalues_im(nev + 1, 0.);
809 
810  // call of ARPACK dseupd/dneupd routine
812  {
813  std::vector<double> z(ldz * nev, 0.);
814  dseupd_(&rvec,
815  &howmany,
816  select.data(),
817  eigenvalues_real.data(),
818  z.data(),
819  &ldz,
820  &sigmar,
821  bmat,
822  &n,
823  which,
824  &nev,
825  &tol,
826  resid.data(),
827  &ncv,
828  v.data(),
829  &ldv,
830  iparam.data(),
831  ipntr.data(),
832  workd.data(),
833  workl.data(),
834  &lworkl,
835  &info);
836  }
837  else
838  {
839  std::vector<double> workev(3 * ncv, 0.);
840  dneupd_(&rvec,
841  &howmany,
842  select.data(),
843  eigenvalues_real.data(),
844  eigenvalues_im.data(),
845  v.data(),
846  &ldz,
847  &sigmar,
848  &sigmai,
849  workev.data(),
850  bmat,
851  &n,
852  which,
853  &nev,
854  &tol,
855  resid.data(),
856  &ncv,
857  v.data(),
858  &ldv,
859  iparam.data(),
860  ipntr.data(),
861  workd.data(),
862  workl.data(),
863  &lworkl,
864  &info);
865  }
866 
867  if (info == 1)
868  {
869  Assert(false, ArpackExcArpackInfoMaxIt(control().max_steps()));
870  }
871  else if (info == 3)
872  {
874  }
875  else if (info != 0)
876  {
878  {
879  Assert(false, ArpackExcArpackInfodseupd(info));
880  }
881  else
882  Assert(false, ArpackExcArpackInfodneupd(info));
883  }
884 
885  for (unsigned int i = 0; i < nev; ++i)
886  for (unsigned int j = 0; j < n; ++j)
887  eigenvectors[i](j) = v[i * n + j];
888 
889  for (unsigned int i = 0; i < nev_const; ++i)
890  eigenvalues[i] =
891  std::complex<double>(eigenvalues_real[i], eigenvalues_im[i]);
892  }
893 }
894 
895 
896 inline SolverControl &
898 {
899  return solver_control;
900 }
901 
903 
904 
905 #endif
906 #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:541
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:518
void mass_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const double factor=1.)
Definition: l2.h:63
#define Assert(cond, exc)
Definition: exceptions.h:1419
static ::ExceptionBase & ArpackExcInvalidNumberofArnoldiVectors(int arg1, int arg2)
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:496
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
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:358
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()