Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
arpack_solver.h
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 
21 #include <deal.II/base/smartpointer.h>
22 
23 #include <deal.II/lac/solver_control.h>
24 
25 #include <cstring>
26 
27 
28 #ifdef DEAL_II_WITH_ARPACK
29 
30 DEAL_II_NAMESPACE_OPEN
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,
460  const WhichEigenvalues eigenvalue_of_interest,
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(
471  ExcMessage(
472  "'largest real part' can only be used for non-symmetric problems!"));
473  Assert(
475  ExcMessage(
476  "'smallest real part' can only be used for non-symmetric problems!"));
477  Assert(
479  ExcMessage(
480  "'largest imaginary part' can only be used for non-symmetric problems!"));
481  Assert(
483  ExcMessage(
484  "'smallest imaginary part' can only be used for non-symmetric problems!"));
485  }
486 }
487 
488 
490  const AdditionalData &data)
492  , additional_data(data)
493  , initial_vector_provided(false)
494  , sigmar(0.0)
495  , sigmai(0.0)
496 {}
497 
498 
499 
500 inline void
501 ArpackSolver::set_shift(const std::complex<double> sigma)
502 {
503  sigmar = sigma.real();
504  sigmai = sigma.imag();
505 }
506 
507 
508 
509 template <typename VectorType>
510 inline void
511 ArpackSolver::set_initial_vector(const VectorType &vec)
512 {
514  resid.resize(vec.size());
515  for (size_type i = 0; i < vec.size(); ++i)
516  resid[i] = vec[i];
517 }
518 
519 
520 template <typename VectorType,
521  typename MatrixType1,
522  typename MatrixType2,
523  typename INVERSE>
524 inline void
525 ArpackSolver::solve(const MatrixType1 & /*system_matrix*/,
526  const MatrixType2 & mass_matrix,
527  const INVERSE & inverse,
528  std::vector<std::complex<double>> &eigenvalues,
529  std::vector<VectorType> & eigenvectors,
530  const unsigned int n_eigenvalues)
531 {
532  // Problem size
533  unsigned int n = eigenvectors[0].size();
534 
535  // Number of eigenvalues
536  const unsigned int nev_const =
537  (n_eigenvalues == 0) ? eigenvalues.size() : n_eigenvalues;
538  // nev for arpack, which might change by plus one during dneupd
539  unsigned int nev = nev_const;
540 
541  // check input sizes
543  {
544  Assert(nev <= eigenvectors.size(),
546  }
547  else
548  Assert(nev + 1 <= eigenvectors.size(),
550  eigenvectors.size()));
551 
552  Assert(nev <= eigenvalues.size(),
554 
555  // check large enough problem size
557 
561 
562  // check whether we have enough Arnoldi vectors
566 
567  // ARPACK mode for dsaupd/dnaupd, here only mode 3, i.e. shift-invert mode
568  int mode = 3;
569 
570  // reverse communication parameter
571  int ido = 0;
572 
573  // 'G' generalized eigenvalue problem 'I' standard eigenvalue problem
574  char bmat[2] = "G";
575 
576  // Specify the eigenvalues of interest, possible parameters "LA" algebraically
577  // largest "SA" algebraically smallest "LM" largest magnitude "SM" smallest
578  // magnitude "LR" largest real part "SR" smallest real part "LI" largest
579  // imaginary part "SI" smallest imaginary part "BE" both ends of spectrum
580  // simultaneous.
581  char which[3];
583  {
585  std::strcpy(which, "LA");
586  break;
588  std::strcpy(which, "SA");
589  break;
590  case largest_magnitude:
591  std::strcpy(which, "LM");
592  break;
593  case smallest_magnitude:
594  std::strcpy(which, "SM");
595  break;
596  case largest_real_part:
597  std::strcpy(which, "LR");
598  break;
599  case smallest_real_part:
600  std::strcpy(which, "SR");
601  break;
603  std::strcpy(which, "LI");
604  break;
606  std::strcpy(which, "SI");
607  break;
608  case both_ends:
609  std::strcpy(which, "BE");
610  break;
611  }
612 
613  // tolerance for ARPACK
614  double tol = control().tolerance();
615 
616  // if the starting vector is used it has to be in resid
617  if (!initial_vector_provided || resid.size() != n)
618  resid.resize(n, 1.);
619 
620  // number of Arnoldi basis vectors specified
621  // in additional_data
623 
624  int ldv = n;
625  std::vector<double> v(ldv * ncv, 0.0);
626 
627  // information to the routines
628  std::vector<int> iparam(11, 0);
629 
630  iparam[0] = 1; // shift strategy
631 
632  // maximum number of iterations
633  iparam[2] = control().max_steps();
634 
635  // Set the mode of dsaupd. 1 is exact shifting, 2 is user-supplied shifts,
636  // 3 is shift-invert mode, 4 is buckling mode, 5 is Cayley mode.
637 
638  iparam[6] = mode;
639  std::vector<int> ipntr(14, 0);
640 
641  // work arrays for ARPACK
642  std::vector<double> workd(3 * n, 0.);
643  int lworkl =
644  additional_data.symmetric ? ncv * ncv + 8 * ncv : 3 * ncv * ncv + 6 * ncv;
645  std::vector<double> workl(lworkl, 0.);
646 
647  // information out of the iteration
648  int info = 1;
649 
650  while (ido != 99)
651  {
652  // call of ARPACK dsaupd/dnaupd routine
654  dsaupd_(&ido,
655  bmat,
656  &n,
657  which,
658  &nev,
659  &tol,
660  resid.data(),
661  &ncv,
662  v.data(),
663  &ldv,
664  iparam.data(),
665  ipntr.data(),
666  workd.data(),
667  workl.data(),
668  &lworkl,
669  &info);
670  else
671  dnaupd_(&ido,
672  bmat,
673  &n,
674  which,
675  &nev,
676  &tol,
677  resid.data(),
678  &ncv,
679  v.data(),
680  &ldv,
681  iparam.data(),
682  ipntr.data(),
683  workd.data(),
684  workl.data(),
685  &lworkl,
686  &info);
687 
688  if (ido == 99)
689  break;
690 
691  switch (mode)
692  {
693  case 3:
694  {
695  switch (ido)
696  {
697  case -1:
698  {
699  VectorType src, dst, tmp;
700  src.reinit(eigenvectors[0]);
701  dst.reinit(src);
702  tmp.reinit(src);
703 
704 
705  for (size_type i = 0; i < src.size(); ++i)
706  src(i) = workd[ipntr[0] - 1 + i];
707 
708  // multiplication with mass matrix M
709  mass_matrix.vmult(tmp, src);
710  // solving linear system
711  inverse.vmult(dst, tmp);
712 
713  for (size_type i = 0; i < dst.size(); ++i)
714  workd[ipntr[1] - 1 + i] = dst(i);
715  }
716  break;
717 
718  case 1:
719  {
720  VectorType src, dst, tmp, tmp2;
721  src.reinit(eigenvectors[0]);
722  dst.reinit(src);
723  tmp.reinit(src);
724  tmp2.reinit(src);
725 
726  for (size_type i = 0; i < src.size(); ++i)
727  {
728  src(i) = workd[ipntr[2] - 1 + i];
729  tmp(i) = workd[ipntr[0] - 1 + i];
730  }
731  // solving linear system
732  inverse.vmult(dst, src);
733 
734  for (size_type i = 0; i < dst.size(); ++i)
735  workd[ipntr[1] - 1 + i] = dst(i);
736  }
737  break;
738 
739  case 2:
740  {
741  VectorType src, dst;
742  src.reinit(eigenvectors[0]);
743  dst.reinit(src);
744 
745  for (size_type i = 0; i < src.size(); ++i)
746  src(i) = workd[ipntr[0] - 1 + i];
747 
748  // Multiplication with mass matrix M
749  mass_matrix.vmult(dst, src);
750 
751  for (size_type i = 0; i < dst.size(); ++i)
752  workd[ipntr[1] - 1 + i] = dst(i);
753  }
754  break;
755 
756  default:
757  Assert(false, ArpackExcArpackIdo(ido));
758  break;
759  }
760  }
761  break;
762  default:
763  Assert(false, ArpackExcArpackMode(mode));
764  break;
765  }
766  }
767 
768  if (info < 0)
769  {
771  {
772  Assert(false, ArpackExcArpackInfodsaupd(info));
773  }
774  else
775  Assert(false, ArpackExcArpackInfodnaupd(info));
776  }
777  else
778  {
779  // 1 - compute eigenvectors, 0 - only eigenvalues
780  int rvec = 1;
781 
782  // which eigenvectors
783  char howmany = 'A';
784 
785  std::vector<int> select(ncv, 1);
786 
787  int ldz = n;
788 
789  std::vector<double> eigenvalues_real(nev + 1, 0.);
790  std::vector<double> eigenvalues_im(nev + 1, 0.);
791 
792  // call of ARPACK dseupd/dneupd routine
794  {
795  std::vector<double> z(ldz * nev, 0.);
796  dseupd_(&rvec,
797  &howmany,
798  select.data(),
799  eigenvalues_real.data(),
800  z.data(),
801  &ldz,
802  &sigmar,
803  bmat,
804  &n,
805  which,
806  &nev,
807  &tol,
808  resid.data(),
809  &ncv,
810  v.data(),
811  &ldv,
812  iparam.data(),
813  ipntr.data(),
814  workd.data(),
815  workl.data(),
816  &lworkl,
817  &info);
818  }
819  else
820  {
821  std::vector<double> workev(3 * ncv, 0.);
822  dneupd_(&rvec,
823  &howmany,
824  select.data(),
825  eigenvalues_real.data(),
826  eigenvalues_im.data(),
827  v.data(),
828  &ldz,
829  &sigmar,
830  &sigmai,
831  workev.data(),
832  bmat,
833  &n,
834  which,
835  &nev,
836  &tol,
837  resid.data(),
838  &ncv,
839  v.data(),
840  &ldv,
841  iparam.data(),
842  ipntr.data(),
843  workd.data(),
844  workl.data(),
845  &lworkl,
846  &info);
847  }
848 
849  if (info == 1)
850  {
851  Assert(false, ArpackExcArpackInfoMaxIt(control().max_steps()));
852  }
853  else if (info == 3)
854  {
856  }
857  else if (info != 0)
858  {
860  {
861  Assert(false, ArpackExcArpackInfodseupd(info));
862  }
863  else
864  Assert(false, ArpackExcArpackInfodneupd(info));
865  }
866 
867  for (unsigned int i = 0; i < nev; ++i)
868  for (unsigned int j = 0; j < n; ++j)
869  eigenvectors[i](j) = v[i * n + j];
870 
871  for (unsigned int i = 0; i < nev_const; ++i)
872  eigenvalues[i] =
873  std::complex<double>(eigenvalues_real[i], eigenvalues_im[i]);
874  }
875 }
876 
877 
878 inline SolverControl &
880 {
881  return solver_control;
882 }
883 
884 DEAL_II_NAMESPACE_CLOSE
885 
886 
887 #endif
888 #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
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)
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
void set_shift(const std::complex< double > sigma)
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)
static ::ExceptionBase & ArpackExcInvalidEigenvalueSize(int arg1, int arg2)
static ::ExceptionBase & ArpackExcInvalidEigenvectorSizeNonsymmetric(int arg1, int arg2)
static ::ExceptionBase & ArpackExcArpackMode(int arg1)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:518
#define Assert(cond, exc)
Definition: exceptions.h:1407
static ::ExceptionBase & ArpackExcInvalidNumberofArnoldiVectors(int arg1, int arg2)
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:496
types::global_dof_index size_type
ArpackSolver(SolverControl &control, const AdditionalData &data=AdditionalData())
static ::ExceptionBase & ArpackExcSmallNumberofArnoldiVectors(int arg1, int arg2)
static ::ExceptionBase & ArpackExcInvalidEigenvectorSize(int arg1, int arg2)
double tolerance() const
const AdditionalData additional_data
static ::ExceptionBase & ArpackExcArpackInfodneupd(int arg1)
static ::ExceptionBase & ArpackExcArpackInfoMaxIt(int arg1)
unsigned int global_dof_index
Definition: types.h:89
static ::ExceptionBase & ArpackExcInvalidNumberofEigenvalues(int arg1, int arg2)
void set_initial_vector(const VectorType &vec)
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
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
static ::ExceptionBase & ArpackExcArpackIdo(int arg1)
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()