From 94b82dc7b06b975c7d82137698cb535bab22b09e Mon Sep 17 00:00:00 2001 From: Guido Kanschat Date: Fri, 10 Apr 2015 10:04:58 +0200 Subject: [PATCH] the Arpack solver in deal.II did not output the eigenvectors, but rather the Schur vectors; this pach fixes the issue --- include/deal.II/lac/arpack_solver.h | 78 +++++++++++++++++++++++------ 1 file changed, 62 insertions(+), 16 deletions(-) diff --git a/include/deal.II/lac/arpack_solver.h b/include/deal.II/lac/arpack_solver.h index 3900e1cdf9..70409a6c73 100644 --- a/include/deal.II/lac/arpack_solver.h +++ b/include/deal.II/lac/arpack_solver.h @@ -76,7 +76,16 @@ extern "C" void dneupd_(int *rvec, char *howmany, int *select, double *d, * dnaupd work and also how to set the parameters appropriately * please take a look into the ARPACK manual. * - * @author Baerbel Janssen, Agnieszka Miedlar, 2010. + * @note Whenever you eliminate degrees of freedom using + * ConstraintMatrix, you generate spurious eigenvalues and + * eigenvectors. If you make sure that the diagonals of eliminated + * matrix rows are all equal to one, you get a single additional + * eigenvalue. But beware that some functions in deal.II set these + * diagonals to rather arbitrary (from the point of view of + * eigenvalue problems) values. See also @ref step_36 "step-36" for an + * example. + * + * @author Baerbel Janssen, Agnieszka Miedlar, 2010, Guido Kanschat 2015 */ class ArpackSolver : public Subscriptor { @@ -130,7 +139,40 @@ public: /** * Solve the generalized eigensprectrum problem $A x=\lambda B x$ by calling - * the dneupd and dnaupd functions of ARPACK. + * the dneupd and dnaupd functions of + * ARPACK. + * + * The function returns a vector of eigenvalues of length n + * and a vector of eigenvectors, where the latter should be twice + * the size of the eigenvalue vector. The first n vectors in + * eigenvectors will be the real parts of the + * eigenvectors, the second n the imaginary parts. + * + * @param A The operator for which we want to compute + * eigenvalues. Actually, this parameter is entirely unused. + * + * @param B The inner product of the underlying space, typically the + * mass matrix. For constrained problems, it can be a partial mass + * matrix, like for instance the velocity mass matrix of a Stokes + * problem. Only its function vmult() is used. + * + * @param inverse This is the possibly shifted inverse that is + * actually used instead of A. Only its function + * vmult() is used. + * + * @param eigenvalues is a vector of complex numbers in which the + * eigenvalues are returned. + * + * @param eigenvectors is a real vector of eigenvectors, + * containing alternatingly the real parts and the imaginary parts of the + * eigenvectors. Therefore, its length should be twice the number of + * eigenvalues. The vectors have to be initialized to match the + * matrices. + * + * @param n_eigenvalues The purpose of this parameter is not clear, + * but it is safe to set it to the size of eigenvalues + * or greater. Leave it at its default zero, which will be reset to the size + * of eigenvalues internally. */ template @@ -140,7 +182,7 @@ public: const INVERSE &inverse, std::vector > &eigenvalues, std::vector &eigenvectors, - const unsigned int n_eigenvalues); + const unsigned int n_eigenvalues = 0); protected: @@ -233,23 +275,26 @@ void ArpackSolver::solve ( //values change magically, so store //them here - const unsigned int n = system_matrix.m(); - const unsigned int n_inside_arpack = system_matrix.m(); + const unsigned int n = eigenvectors[0].size(); + const unsigned int n_inside_arpack = eigenvectors[0].size(); + // Number of eigenvalues for arpack + const unsigned int nev = (n_eigenvalues == 0) ? eigenvalues.size() : n_eigenvalues; + AssertIndexRange(eigenvalues.size()-1, nev); /* if(n < 0 || nev <0 || p < 0 || maxit < 0 ) std:cout << "All input parameters have to be positive.\n"; */ Assert (n_eigenvalues < n, - ExcInvalidNumberofEigenvalues(n_eigenvalues, n)); + ExcInvalidNumberofEigenvalues(nev, n)); Assert (additional_data.number_of_arnoldi_vectors < n, ExcInvalidNumberofArnoldiVectors( additional_data.number_of_arnoldi_vectors, n)); - Assert (additional_data.number_of_arnoldi_vectors > 2*n_eigenvalues+1, + Assert (additional_data.number_of_arnoldi_vectors > 2*nev+1, ExcSmallNumberofArnoldiVectors( - additional_data.number_of_arnoldi_vectors, n_eigenvalues)); + additional_data.number_of_arnoldi_vectors, nev)); // ARPACK mode for dnaupd, here only mode 3 int mode = 3; @@ -341,7 +386,6 @@ void ArpackSolver::solve ( //information out of the iteration int info = 1; - const unsigned int nev = n_eigenvalues; while (ido != 99) { // call of ARPACK dnaupd routine @@ -445,9 +489,9 @@ void ArpackSolver::solve ( int rvec = 1; // which eigenvectors - char howmany[4] = "All"; + char howmany = 'A'; - std::vector select (ncv, 0); + std::vector select (ncv, 1); int ldz = n; @@ -459,11 +503,11 @@ void ArpackSolver::solve ( int lworkev = 3*ncv; std::vector workev (lworkev, 0.); - std::vector eigenvalues_real (n_eigenvalues, 0.); - std::vector eigenvalues_im (n_eigenvalues, 0.); + std::vector eigenvalues_real (nev, 0.); + std::vector eigenvalues_im (nev, 0.); // call of ARPACK dneupd routine - dneupd_(&rvec, howmany, &select[0], &eigenvalues_real[0], + dneupd_(&rvec, &howmany, &select[0], &eigenvalues_real[0], &eigenvalues_im[0], &z[0], &ldz, &sigmar, &sigmai, &workev[0], bmat, &n_inside_arpack, which, &nev, &tol, &resid[0], &ncv, &v[0], &ldv, @@ -482,9 +526,11 @@ void ArpackSolver::solve ( Assert (false, ExcArpackInfodneupd(info)); } - for (size_type i=0; i