Reference documentation for deal.II version 8.1.0
arpack_solver.h
1 // ---------------------------------------------------------------------
2 // @f$Id: arpack_solver.h 31932 2013-12-08 02:15:54Z heister @f$
3 //
4 // Copyright (C) 2010 - 2013 by the deal.II authors
5 //
6 // This file is part of the deal.II library.
7 //
8 // The deal.II library is free software; you can use it, redistribute
9 // it, and/or modify it under the terms of the GNU Lesser General
10 // Public License as published by the Free Software Foundation; either
11 // version 2.1 of the License, or (at your option) any later version.
12 // The full text of the license can be found in the file LICENSE at
13 // the top level of the deal.II distribution.
14 //
15 // ---------------------------------------------------------------------
16 
17 #ifndef __deal2__arpack_solver_h
18 #define __deal2__arpack_solver_h
19 
20 #include <deal.II/base/config.h>
21 #include <deal.II/base/smartpointer.h>
22 #include <deal.II/lac/solver_control.h>
23 
24 #include <cstring>
25 
26 
27 #ifdef DEAL_II_WITH_ARPACK
28 
30 
31 
32 extern "C" void dnaupd_(int *ido, char *bmat, const unsigned int *n, char *which,
33  const unsigned int *nev, const double *tol, double *resid, int *ncv,
34  double *v, int *ldv, int *iparam, int *ipntr,
35  double *workd, double *workl, int *lworkl,
36  int *info);
37 
38 extern "C" void dneupd_(int *rvec, char *howmany, int *select, double *d,
39  double *di, double *z, int *ldz, double *sigmar,
40  double *sigmai, double *workev, char *bmat,const unsigned int *n, char *which,
41  const unsigned int *nev, const double *tol, double *resid, int *ncv,
42  double *v, int *ldv, int *iparam, int *ipntr,
43  double *workd, double *workl, int *lworkl, int *info);
44 
83 class ArpackSolver : public Subscriptor
84 {
85 public:
90 
91 
98  {
99  algebraically_largest,
100  algebraically_smallest,
101  largest_magnitude,
102  smallest_magnitude,
103  largest_real_part,
104  smallest_real_part,
105  largest_imaginary_part,
106  smallest_imaginary_part,
107  both_ends
108  };
109 
116  {
117  const unsigned int number_of_arnoldi_vectors;
118  const WhichEigenvalues eigenvalue_of_interest;
120  const unsigned int number_of_arnoldi_vectors = 15,
121  const WhichEigenvalues eigenvalue_of_interest = largest_magnitude);
122  };
123 
128  SolverControl &control () const;
129 
133  ArpackSolver(SolverControl &control,
134  const AdditionalData &data = AdditionalData());
135 
142  template <typename VECTOR, typename MATRIX1,
143  typename MATRIX2, typename INVERSE>
144  void solve(
145  const MATRIX1 &A,
146  const MATRIX2 &B,
147  const INVERSE &inverse,
148  std::vector<std::complex<double> > &eigenvalues,
149  std::vector<VECTOR> &eigenvectors,
150  const unsigned int n_eigenvalues);
151 
152 protected:
153 
160 
166 
167 private:
168 
172  DeclException2 (ExcInvalidNumberofEigenvalues, int, int,
173  << "Number of wanted eigenvalues " << arg1
174  << " is larger that the size of the matrix " << arg2);
175 
176  DeclException2 (ExcInvalidNumberofArnoldiVectors, int, int,
177  << "Number of Arnoldi vectors " << arg1
178  << " is larger that the size of the matrix " << arg2);
179 
180  DeclException2 (ExcSmallNumberofArnoldiVectors, int, int,
181  << "Number of Arnoldi vectors " << arg1
182  << " is too small to obtain " << arg2
183  << " eigenvalues");
184 
185  DeclException1 (ExcArpackIdo, int, << "This ido " << arg1
186  << " is not supported. Check documentation of ARPACK");
187 
188  DeclException1 (ExcArpackMode, int, << "This mode " << arg1
189  << " is not supported. Check documentation of ARPACK");
190 
191  DeclException1 (ExcArpackInfodsaupd, int,
192  << "Error with dsaupd, info " << arg1
193  << ". Check documentation of ARPACK");
194 
195  DeclException1 (ExcArpackInfodneupd, int,
196  << "Error with dneupd, info " << arg1
197  << ". Check documentation of ARPACK");
198 
199  DeclException1 (ExcArpackInfoMaxIt, int,
200  << "Maximum number " << arg1
201  << " of iterations reached.");
202 
203  DeclException1 (ExcArpackNoShifts, int,
204  << "No shifts could be applied during implicit"
205  << " Arnoldi update, try increasing the number of"
206  << " Arnoldi vectors.");
207 };
208 
209 
210 inline
211 ArpackSolver::AdditionalData::
212 AdditionalData (const unsigned int number_of_arnoldi_vectors,
213  const WhichEigenvalues eigenvalue_of_interest)
214  :
215  number_of_arnoldi_vectors(number_of_arnoldi_vectors),
216  eigenvalue_of_interest(eigenvalue_of_interest)
217 {}
218 
219 
220 inline
222  const AdditionalData &data)
223  :
224  solver_control (control),
225  additional_data (data)
226 
227 {}
228 
229 
230 template <typename VECTOR, typename MATRIX1,
231  typename MATRIX2, typename INVERSE>
232 inline
234  const MATRIX1 &system_matrix,
235  const MATRIX2 &mass_matrix,
236  const INVERSE &inverse,
237  std::vector<std::complex<double> > &eigenvalues,
238  std::vector<VECTOR> &eigenvectors,
239  const unsigned int n_eigenvalues)
240 {
241  //inside the routines of ARPACK the
242  //values change magically, so store
243  //them here
244 
245  const unsigned int n = system_matrix.m();
246  const unsigned int n_inside_arpack = system_matrix.m();
247 
248  /*
249  if(n < 0 || nev <0 || p < 0 || maxit < 0 )
250  std:cout << "All input parameters have to be positive.\n";
251  */
252  Assert (n_eigenvalues < n,
253  ExcInvalidNumberofEigenvalues(n_eigenvalues, n));
254 
255  Assert (additional_data.number_of_arnoldi_vectors < n,
256  ExcInvalidNumberofArnoldiVectors(
257  additional_data.number_of_arnoldi_vectors, n));
258 
259  Assert (additional_data.number_of_arnoldi_vectors > 2*n_eigenvalues+1,
260  ExcSmallNumberofArnoldiVectors(
261  additional_data.number_of_arnoldi_vectors, n_eigenvalues));
262  // ARPACK mode for dnaupd, here only mode 3
263  int mode = 3;
264 
265  // reverse communication parameter
266  int ido = 0;
267 
272  char bmat[2] = "G";
273 
286  char which[3];
287  switch (additional_data.eigenvalue_of_interest)
288  {
289  case algebraically_largest:
290  std::strcpy (which, "LA");
291  break;
292  case algebraically_smallest:
293  std::strcpy (which, "SA");
294  break;
295  case largest_magnitude:
296  std::strcpy (which, "LM");
297  break;
298  case smallest_magnitude:
299  std::strcpy (which, "SM");
300  break;
301  case largest_real_part:
302  std::strcpy (which, "LR");
303  break;
304  case smallest_real_part:
305  std::strcpy (which, "SR");
306  break;
307  case largest_imaginary_part:
308  std::strcpy (which, "LI");
309  break;
310  case smallest_imaginary_part:
311  std::strcpy (which, "SI");
312  break;
313  case both_ends:
314  std::strcpy (which, "BE");
315  break;
316  }
317 
318  // tolerance for ARPACK
319  const double tol = control().tolerance();
320 
321  // if the starting vector is used it has to be in resid
322  std::vector<double> resid(n, 1.);
323 
324  // number of Arnoldi basis vectors specified
325  // in additional_data
326  int ncv = additional_data.number_of_arnoldi_vectors;
327 
328  int ldv = n;
329  std::vector<double> v (ldv*ncv, 0.0);
330 
331  //information to the routines
332  std::vector<int> iparam (11, 0);
333 
334  iparam[0] = 1; // shift strategy
335 
336  // maximum number of iterations
337  iparam[2] = control().max_steps();
338 
347  iparam[6] = mode;
348  std::vector<int> ipntr (14, 0);
349 
350  // work arrays for ARPACK
351  double *workd;
352  workd = new double[3*n];
353 
354  for (unsigned int i=0; i<3*n; ++i)
355  workd[i] = 0.0;
356 
357  int lworkl = 3*ncv*(ncv + 6);
358  std::vector<double> workl (lworkl, 0.);
359  //information out of the iteration
360  int info = 1;
361 
362  const unsigned int nev = n_eigenvalues;
363  while (ido != 99)
364  {
365  // call of ARPACK dnaupd routine
366  dnaupd_(&ido, bmat, &n_inside_arpack, which, &nev, &tol,
367  &resid[0], &ncv, &v[0], &ldv, &iparam[0], &ipntr[0],
368  workd, &workl[0], &lworkl, &info);
369 
370  if (ido == 99)
371  break;
372 
373  switch (mode)
374  {
375  case 3:
376  {
377  switch (ido)
378  {
379  case -1:
380  {
381 
382  VECTOR src,dst,tmp;
383  src.reinit(eigenvectors[0]);
384  dst.reinit(src);
385  tmp.reinit(src);
386 
387 
388  for (size_type i=0; i<src.size(); ++i)
389  src(i) = *(workd+ipntr[0]-1+i);
390 
391  // multiplication with mass matrix M
392  mass_matrix.vmult(tmp, src);
393  // solving linear system
394  inverse.vmult(dst,tmp);
395 
396  for (size_type i=0; i<dst.size(); ++i)
397  *(workd+ipntr[1]-1+i) = dst(i);
398  }
399  break;
400 
401  case 1:
402  {
403 
404  VECTOR src,dst,tmp, tmp2;
405  src.reinit(eigenvectors[0]);
406  dst.reinit(src);
407  tmp.reinit(src);
408  tmp2.reinit(src);
409 
410  for (size_type i=0; i<src.size(); ++i)
411  {
412  src(i) = *(workd+ipntr[2]-1+i);
413  tmp(i) = *(workd+ipntr[0]-1+i);
414  }
415  // solving linear system
416  inverse.vmult(dst,src);
417 
418  for (size_type i=0; i<dst.size(); ++i)
419  *(workd+ipntr[1]-1+i) = dst(i);
420  }
421  break;
422 
423  case 2:
424  {
425 
426  VECTOR src,dst;
427  src.reinit(eigenvectors[0]);
428  dst.reinit(src);
429 
430  for (size_type i=0; i<src.size(); ++i)
431  src(i) = *(workd+ipntr[0]-1+i);
432 
433  // Multiplication with mass matrix M
434  mass_matrix.vmult(dst, src);
435 
436  for (size_type i=0; i<dst.size(); ++i)
437  *(workd+ipntr[1]-1+i) = dst(i);
438 
439  }
440  break;
441 
442  default:
443  Assert (false, ExcArpackIdo(ido));
444  break;
445  }
446  }
447  break;
448  default:
449  Assert (false, ExcArpackMode(mode));
450  break;
451  }
452  }
453 
454  if (info<0)
455  {
456  Assert (false, ExcArpackInfodsaupd(info));
457  }
458  else
459  {
463  int rvec = 1;
464 
465  // which eigenvectors
466  char howmany[4] = "All";
467 
468  std::vector<int> select (ncv, 0);
469 
470  int ldz = n;
471 
472  std::vector<double> z (ldz*ncv, 0.);
473 
474  double sigmar = 0.0; // real part of the shift
475  double sigmai = 0.0; // imaginary part of the shift
476 
477  int lworkev = 3*ncv;
478  std::vector<double> workev (lworkev, 0.);
479 
480  std::vector<double> eigenvalues_real (n_eigenvalues, 0.);
481  std::vector<double> eigenvalues_im (n_eigenvalues, 0.);
482 
483  // call of ARPACK dneupd routine
484  dneupd_(&rvec, howmany, &select[0], &eigenvalues_real[0],
485  &eigenvalues_im[0], &z[0], &ldz, &sigmar, &sigmai,
486  &workev[0], bmat, &n_inside_arpack, which, &nev, &tol,
487  &resid[0], &ncv, &v[0], &ldv,
488  &iparam[0], &ipntr[0], workd, &workl[0], &lworkl, &info);
489 
490  if (info == 1)
491  {
492  Assert (false, ExcArpackInfoMaxIt(control().max_steps()));
493  }
494  else if (info == 3)
495  {
496  Assert (false, ExcArpackNoShifts(1));
497  }
498  else if (info!=0)
499  {
500  Assert (false, ExcArpackInfodneupd(info));
501  }
502 
503  for (size_type i=0; i<eigenvectors.size(); ++i)
504  for (unsigned int j=0; j<n; ++j)
505  eigenvectors[i](j) = v[i*n+j];
506 
507  delete[] workd;
508 
509  AssertDimension (eigenvalues.size(), eigenvalues_real.size());
510  AssertDimension (eigenvalues.size(), eigenvalues_im.size());
511 
512  for (size_type i=0; i<eigenvalues.size(); ++i)
513  eigenvalues[i] = std::complex<double> (eigenvalues_real[i],
514  eigenvalues_im[i]);
515  }
516 }
517 
518 
519 inline
521 {
522  return solver_control;
523 }
524 
525 DEAL_II_NAMESPACE_CLOSE
526 
527 
528 #endif
529 #endif
SolverControl & solver_control
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:858
double tolerance() const
SolverControl & control() const
types::global_dof_index size_type
Definition: arpack_solver.h:89
void solve(const MATRIX1 &A, const MATRIX2 &B, const INVERSE &inverse, std::vector< std::complex< double > > &eigenvalues, std::vector< VECTOR > &eigenvectors, const unsigned int n_eigenvalues)
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:515
DeclException2(ExcInvalidNumberofEigenvalues, int, int,<< "Number of wanted eigenvalues "<< arg1<< " is larger that the size of the matrix "<< arg2)
unsigned int global_dof_index
Definition: types.h:100
unsigned int max_steps() const
#define Assert(cond, exc)
Definition: exceptions.h:299
ArpackSolver(SolverControl &control, const AdditionalData &data=AdditionalData())
const AdditionalData additional_data