17 #ifndef __deal2__arpack_solver_h
18 #define __deal2__arpack_solver_h
20 #include <deal.II/base/config.h>
21 #include <deal.II/base/smartpointer.h>
22 #include <deal.II/lac/solver_control.h>
27 #ifdef DEAL_II_WITH_ARPACK
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,
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);
99 algebraically_largest,
100 algebraically_smallest,
105 largest_imaginary_part,
106 smallest_imaginary_part,
117 const unsigned int number_of_arnoldi_vectors;
120 const unsigned int number_of_arnoldi_vectors = 15,
142 template <
typename VECTOR,
typename MATRIX1,
143 typename MATRIX2,
typename INVERSE>
147 const INVERSE &inverse,
148 std::vector<std::complex<double> > &eigenvalues,
149 std::vector<VECTOR> &eigenvectors,
150 const unsigned int n_eigenvalues);
173 <<
"Number of wanted eigenvalues " << arg1
174 <<
" is larger that the size of the matrix " << arg2);
177 <<
"Number of Arnoldi vectors " << arg1
178 <<
" is larger that the size of the matrix " << arg2);
181 <<
"Number of Arnoldi vectors " << arg1
182 <<
" is too small to obtain " << arg2
185 DeclException1 (ExcArpackIdo,
int, <<
"This ido " << arg1
186 <<
" is not supported. Check documentation of ARPACK");
188 DeclException1 (ExcArpackMode,
int, <<
"This mode " << arg1
189 <<
" is not supported. Check documentation of ARPACK");
191 DeclException1 (ExcArpackInfodsaupd,
int,
192 <<
"Error with dsaupd, info " << arg1
193 <<
". Check documentation of ARPACK");
195 DeclException1 (ExcArpackInfodneupd,
int,
196 <<
"Error with dneupd, info " << arg1
197 <<
". Check documentation of ARPACK");
199 DeclException1 (ExcArpackInfoMaxIt,
int,
200 <<
"Maximum number " << arg1
201 <<
" of iterations reached.");
203 DeclException1 (ExcArpackNoShifts,
int,
204 <<
"No shifts could be applied during implicit"
205 <<
" Arnoldi update, try increasing the number of"
206 <<
" Arnoldi vectors.");
211 ArpackSolver::AdditionalData::
212 AdditionalData (
const unsigned int number_of_arnoldi_vectors,
213 const WhichEigenvalues eigenvalue_of_interest)
215 number_of_arnoldi_vectors(number_of_arnoldi_vectors),
216 eigenvalue_of_interest(eigenvalue_of_interest)
224 solver_control (control),
225 additional_data (data)
230 template <
typename VECTOR,
typename MATRIX1,
231 typename MATRIX2,
typename INVERSE>
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)
245 const unsigned int n = system_matrix.m();
246 const unsigned int n_inside_arpack = system_matrix.m();
252 Assert (n_eigenvalues < n,
253 ExcInvalidNumberofEigenvalues(n_eigenvalues, n));
256 ExcInvalidNumberofArnoldiVectors(
260 ExcSmallNumberofArnoldiVectors(
289 case algebraically_largest:
290 std::strcpy (which,
"LA");
292 case algebraically_smallest:
293 std::strcpy (which,
"SA");
295 case largest_magnitude:
296 std::strcpy (which,
"LM");
298 case smallest_magnitude:
299 std::strcpy (which,
"SM");
301 case largest_real_part:
302 std::strcpy (which,
"LR");
304 case smallest_real_part:
305 std::strcpy (which,
"SR");
307 case largest_imaginary_part:
308 std::strcpy (which,
"LI");
310 case smallest_imaginary_part:
311 std::strcpy (which,
"SI");
314 std::strcpy (which,
"BE");
322 std::vector<double> resid(n, 1.);
329 std::vector<double> v (ldv*ncv, 0.0);
332 std::vector<int> iparam (11, 0);
348 std::vector<int> ipntr (14, 0);
352 workd =
new double[3*n];
354 for (
unsigned int i=0; i<3*n; ++i)
357 int lworkl = 3*ncv*(ncv + 6);
358 std::vector<double> workl (lworkl, 0.);
362 const unsigned int nev = n_eigenvalues;
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);
383 src.reinit(eigenvectors[0]);
388 for (size_type i=0; i<src.size(); ++i)
389 src(i) = *(workd+ipntr[0]-1+i);
392 mass_matrix.vmult(tmp, src);
394 inverse.vmult(dst,tmp);
396 for (size_type i=0; i<dst.size(); ++i)
397 *(workd+ipntr[1]-1+i) = dst(i);
404 VECTOR src,dst,tmp, tmp2;
405 src.reinit(eigenvectors[0]);
410 for (size_type i=0; i<src.size(); ++i)
412 src(i) = *(workd+ipntr[2]-1+i);
413 tmp(i) = *(workd+ipntr[0]-1+i);
416 inverse.vmult(dst,src);
418 for (size_type i=0; i<dst.size(); ++i)
419 *(workd+ipntr[1]-1+i) = dst(i);
427 src.reinit(eigenvectors[0]);
430 for (size_type i=0; i<src.size(); ++i)
431 src(i) = *(workd+ipntr[0]-1+i);
434 mass_matrix.vmult(dst, src);
436 for (size_type i=0; i<dst.size(); ++i)
437 *(workd+ipntr[1]-1+i) = dst(i);
443 Assert (
false, ExcArpackIdo(ido));
449 Assert (
false, ExcArpackMode(mode));
456 Assert (
false, ExcArpackInfodsaupd(info));
466 char howmany[4] =
"All";
468 std::vector<int> select (ncv, 0);
472 std::vector<double> z (ldz*ncv, 0.);
478 std::vector<double> workev (lworkev, 0.);
480 std::vector<double> eigenvalues_real (n_eigenvalues, 0.);
481 std::vector<double> eigenvalues_im (n_eigenvalues, 0.);
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);
496 Assert (
false, ExcArpackNoShifts(1));
500 Assert (
false, ExcArpackInfodneupd(info));
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];
512 for (size_type i=0; i<eigenvalues.size(); ++i)
513 eigenvalues[i] = std::complex<double> (eigenvalues_real[i],
525 DEAL_II_NAMESPACE_CLOSE
SolverControl & solver_control
#define AssertDimension(dim1, dim2)
SolverControl & control() const
types::global_dof_index size_type
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)
DeclException2(ExcInvalidNumberofEigenvalues, int, int,<< "Number of wanted eigenvalues "<< arg1<< " is larger that the size of the matrix "<< arg2)
unsigned int global_dof_index
unsigned int max_steps() const
#define Assert(cond, exc)
ArpackSolver(SolverControl &control, const AdditionalData &data=AdditionalData())
const AdditionalData additional_data