From: Joscha Gedicke Date: Thu, 21 Jul 2016 13:19:58 +0000 (+0200) Subject: fix (P)ARPACK interface for non-symmetric matrices X-Git-Tag: v8.5.0-rc1~821^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=9e487bc1f593e230438d99b66586cc5ac6495100;p=dealii.git fix (P)ARPACK interface for non-symmetric matrices -> test for length n+1 of the eigenvectors in the non-symmetric case, which is related to the storage scheme of complex components -> allow ARPACK to store n+1 eigenvalues in the work arrays in case the eigenvalues n and n+1 form a complex conjugate pair -> copy nev eigenvectors from the work arrays as nev is increased by one by dneupd in some cases -> add set_initial_vector function -> add tests for advection-diffusion FE matrices with complex eigenvalues -> improve description of the tests --- diff --git a/doc/news/changes.h b/doc/news/changes.h index d917083a7d..761e56ba88 100644 --- a/doc/news/changes.h +++ b/doc/news/changes.h @@ -395,6 +395,10 @@ inconvenience this causes.

Specific improvements

    +
  1. Fixed: (P)ARPACK interface for non-symmetric matrices. +
    + (Joscha Gedicke, 2016/08/01) +
  2. Fixed: The TrilinosWrappers::SparsityPattern::print() and TrilinosWrappers::SparsityPattern::print_gnuplot() methods did not produce correct output on distributed computations. This is now fixed. diff --git a/include/deal.II/lac/arpack_solver.h b/include/deal.II/lac/arpack_solver.h index 63751a8250..dee70d45c7 100644 --- a/include/deal.II/lac/arpack_solver.h +++ b/include/deal.II/lac/arpack_solver.h @@ -28,23 +28,38 @@ DEAL_II_NAMESPACE_OPEN -extern "C" void dnaupd_(int *ido, char *bmat, const unsigned int *n, char *which, - const unsigned int *nev, const double *tol, double *resid, int *ncv, +extern "C" 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); + +extern "C" 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); extern "C" void dneupd_(int *rvec, char *howmany, int *select, double *d, double *di, double *z, int *ldz, double *sigmar, - double *sigmai, double *workev, char *bmat,const unsigned int *n, char *which, - const unsigned int *nev, const double *tol, double *resid, int *ncv, + 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); + +extern "C" 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); /** * Interface for using ARPACK. ARPACK is a collection of Fortran77 subroutines * designed to solve large scale eigenvalue problems. Here we interface to - * the routines dneupd and dnaupd of ARPACK. The + * the routines dnaupd and dneupd of ARPACK. + * If the operator is specified to be symmetric we use the symmetric interface + * dsaupd and dseupd of ARPACK instead. The * package is designed to compute a few eigenvalues and corresponding * eigenvectors of a general n by n matrix A. It is most appropriate for large * sparse matrices A. @@ -72,8 +87,9 @@ extern "C" void dneupd_(int *rvec, char *howmany, int *select, double *d, * Through the AdditionalData the user can specify some of the parameters to * be set. * - * For further information on how the ARPACK routines dneupd and - * dnaupd work and also how to set the parameters appropriately + * For further information on how the ARPACK routines dsaupd, + * dseupd, dnaupd and dneupd work + * and also how to set the parameters appropriately * please take a look into the ARPACK manual. * * @note Whenever you eliminate degrees of freedom using ConstraintMatrix, you @@ -85,7 +101,7 @@ extern "C" void dneupd_(int *rvec, char *howmany, int *select, double *d, * @ref step_36 "step-36" * for an example. * - * @author Baerbel Janssen, Agnieszka Miedlar, 2010, Guido Kanschat 2015 + * @author Baerbel Janssen, Agnieszka Miedlar, 2010, Guido Kanschat 2015, Joscha Gedicke 2016 */ class ArpackSolver : public Subscriptor { @@ -121,9 +137,11 @@ public: { const unsigned int number_of_arnoldi_vectors; const WhichEigenvalues eigenvalue_of_interest; + const bool symmetric; AdditionalData( const unsigned int number_of_arnoldi_vectors = 15, - const WhichEigenvalues eigenvalue_of_interest = largest_magnitude); + const WhichEigenvalues eigenvalue_of_interest = largest_magnitude, + const bool symmetric = false); }; /** @@ -137,15 +155,28 @@ public: ArpackSolver(SolverControl &control, const AdditionalData &data = AdditionalData()); + /** + * Set initial vector for building Krylov space. + */ + template + void set_initial_vector(const VectorType &vec); + /** * Solve the generalized eigensprectrum problem $A x=\lambda B x$ by calling - * the dneupd and dnaupd functions of ARPACK. + * the dsaupd and dseupd or + * dnaupd and dneupd 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. + * vector of eigenvectors of length n in the symmetric case + * and of length n+1 in the non-symmetric case. In the symmetric case + * all eigenvectors are real. In the non-symmetric case complex eigenvalues + * always occur as complex conjugate pairs. Therefore the eigenvector for an + * eigenvalue with nonzero complex part is stored by putting the real and + * the imaginary parts in consecutive real-valued vectors. The eigenvector + * of the complex conjugate eigenvalue does not need to be stored, since it + * is just the complex conjugate of the stored eigenvector. Thus, if the last + * n-th eigenvalue has a nonzero imaginary part, Arpack needs in total n+1 + * real-valued vectors to store real and imaginary parts of the eigenvectors. * * @param A The operator for which we want to compute eigenvalues. Actually, * this parameter is entirely unused. @@ -163,9 +194,16 @@ public: * 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. + * the real parts of all eigenvectors and the imaginary parts of the eigenvectors + * corresponding to complex conjugate eigenvalue pairs. + * Therefore, its length should be n in the symmetric case and n+1 + * in the non-symmetric case. In the non-symmetric case the storage scheme + * leads for example to the following pattern. Suppose that the first two + * eigenvalues are real and the third and fourth are a complex conjugate + * pair. Asking for three eigenpairs results in [real(v1),real(v2), + * real(v3),imag(v3)]. Note that we get the same pattern if we ask for four + * eigenpairs in this example, since the fourth eigenvector is simply the + * complex conjugate of the third one. * * @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. @@ -194,43 +232,70 @@ protected: */ const AdditionalData additional_data; + /** + * Store an initial vector + */ + bool initial_vector_provided; + std::vector resid; + private: /** * Exceptions. */ - DeclException2 (ExcInvalidNumberofEigenvalues, int, int, + DeclException2 (ArpackExcInvalidNumberofEigenvalues, int, int, << "Number of wanted eigenvalues " << arg1 << " is larger that the size of the matrix " << arg2); - DeclException2 (ExcInvalidNumberofArnoldiVectors, int, int, + DeclException2 (ArpackExcInvalidEigenvectorSize, int, int, + << "Number of wanted eigenvalues " << arg1 + << " is larger that the size of eigenvectors " << arg2); + + DeclException2 (ArpackExcInvalidEigenvectorSizeNonsymmetric, int, int, + << "To store the real and complex parts of " << arg1 + << " eigenvectors in real-valued vectors, their size (currently set to " << arg2 + << ") should be greater than or equal to " << arg1+1); + + DeclException2 (ArpackExcInvalidEigenvalueSize, int, int, + << "Number of wanted eigenvalues " << arg1 + << " is larger that the size of eigenvalues " << arg2); + + DeclException2 (ArpackExcInvalidNumberofArnoldiVectors, int, int, << "Number of Arnoldi vectors " << arg1 << " is larger that the size of the matrix " << arg2); - DeclException2 (ExcSmallNumberofArnoldiVectors, int, int, + DeclException2 (ArpackExcSmallNumberofArnoldiVectors, int, int, << "Number of Arnoldi vectors " << arg1 << " is too small to obtain " << arg2 << " eigenvalues"); - DeclException1 (ExcArpackIdo, int, << "This ido " << arg1 + DeclException1 (ArpackExcArpackIdo, int, << "This ido " << arg1 << " is not supported. Check documentation of ARPACK"); - DeclException1 (ExcArpackMode, int, << "This mode " << arg1 + DeclException1 (ArpackExcArpackMode, int, << "This mode " << arg1 << " is not supported. Check documentation of ARPACK"); - DeclException1 (ExcArpackInfodsaupd, int, + DeclException1 (ArpackExcArpackInfodsaupd, int, << "Error with dsaupd, info " << arg1 << ". Check documentation of ARPACK"); - DeclException1 (ExcArpackInfodneupd, int, + DeclException1 (ArpackExcArpackInfodnaupd, int, + << "Error with dnaupd, info " << arg1 + << ". Check documentation of ARPACK"); + + DeclException1 (ArpackExcArpackInfodseupd, int, + << "Error with dseupd, info " << arg1 + << ". Check documentation of ARPACK"); + + DeclException1 (ArpackExcArpackInfodneupd, int, << "Error with dneupd, info " << arg1 << ". Check documentation of ARPACK"); - DeclException1 (ExcArpackInfoMaxIt, int, + DeclException1 (ArpackExcArpackInfoMaxIt, int, << "Maximum number " << arg1 << " of iterations reached."); - DeclExceptionMsg (ExcArpackNoShifts, + DeclExceptionMsg (ArpackExcArpackNoShifts, "No shifts could be applied during implicit" " Arnoldi update, try increasing the number of" " Arnoldi vectors."); @@ -240,10 +305,12 @@ private: inline ArpackSolver::AdditionalData:: AdditionalData (const unsigned int number_of_arnoldi_vectors, - const WhichEigenvalues eigenvalue_of_interest) + const WhichEigenvalues eigenvalue_of_interest, + const bool symmetric) : number_of_arnoldi_vectors(number_of_arnoldi_vectors), - eigenvalue_of_interest(eigenvalue_of_interest) + eigenvalue_of_interest(eigenvalue_of_interest), + symmetric(symmetric) {} @@ -252,10 +319,21 @@ ArpackSolver::ArpackSolver (SolverControl &control, const AdditionalData &data) : solver_control (control), - additional_data (data) - + additional_data (data), + initial_vector_provided(false) {} +template +inline +void ArpackSolver:: +set_initial_vector(const VectorType &vec) +{ + initial_vector_provided = true; + resid.resize(vec.size()); + for (size_type i = 0; i < vec.size(); ++i) + resid[i] = vec[i]; +} + template @@ -267,31 +345,43 @@ void ArpackSolver::solve (const MatrixType1 &/*system_matrix*/, std::vector &eigenvectors, const unsigned int n_eigenvalues) { - //inside the routines of ARPACK the - //values change magically, so store - //them here - - 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(nev, n)); + // Problem size + unsigned int n = eigenvectors[0].size(); + + // Number of eigenvalues + const unsigned int nev_const = (n_eigenvalues == 0) ? eigenvalues.size() : n_eigenvalues; + // nev for arpack, which might change by plus one during dneupd + unsigned int nev = nev_const; + + // check input sizes + if (additional_data.symmetric) + { + Assert (nev <= eigenvectors.size(), + ArpackExcInvalidEigenvectorSize(nev, eigenvectors.size())); + } + else + Assert (nev+1 <= eigenvectors.size(), + ArpackExcInvalidEigenvectorSizeNonsymmetric(nev, eigenvectors.size())); + + Assert (nev <= eigenvalues.size(), + ArpackExcInvalidEigenvalueSize(nev, eigenvalues.size())); + + // check large enough problem size + Assert (nev < n, + ArpackExcInvalidNumberofEigenvalues(nev, n)); Assert (additional_data.number_of_arnoldi_vectors < n, - ExcInvalidNumberofArnoldiVectors( + ArpackExcInvalidNumberofArnoldiVectors( additional_data.number_of_arnoldi_vectors, n)); + // check whether we have enough Arnoldi vectors Assert (additional_data.number_of_arnoldi_vectors > 2*nev+1, - ExcSmallNumberofArnoldiVectors( + ArpackExcSmallNumberofArnoldiVectors( additional_data.number_of_arnoldi_vectors, nev)); - // ARPACK mode for dnaupd, here only mode 3 + + /* ARPACK mode for dsaupd/dnaupd, here only mode 3, + * i.e. shift-invert mode + */ int mode = 3; // reverse communication parameter @@ -342,10 +432,11 @@ void ArpackSolver::solve (const MatrixType1 &/*system_matrix*/, } // tolerance for ARPACK - const double tol = control().tolerance(); + double tol = control().tolerance(); // if the starting vector is used it has to be in resid - std::vector resid(n, 1.); + if (!initial_vector_provided || resid.size() != n) + resid.resize(n, 1.); // number of Arnoldi basis vectors specified // in additional_data @@ -371,23 +462,24 @@ void ArpackSolver::solve (const MatrixType1 &/*system_matrix*/, std::vector ipntr (14, 0); // work arrays for ARPACK - double *workd; - workd = new double[3*n]; - - for (unsigned int i=0; i<3*n; ++i) - workd[i] = 0.0; - - int lworkl = 3*ncv*(ncv + 6); + std::vector workd (3*n, 0.); + int lworkl = additional_data.symmetric ? ncv*ncv + 8*ncv : 3*ncv*ncv+6*ncv; std::vector workl (lworkl, 0.); + //information out of the iteration int info = 1; while (ido != 99) { - // call of ARPACK dnaupd routine - dnaupd_(&ido, bmat, &n_inside_arpack, which, &nev, &tol, - &resid[0], &ncv, &v[0], &ldv, &iparam[0], &ipntr[0], - workd, &workl[0], &lworkl, &info); + // call of ARPACK dsaupd/dnaupd routine + if (additional_data.symmetric) + dsaupd_(&ido, bmat, &n, which, &nev, &tol, + &resid[0], &ncv, &v[0], &ldv, &iparam[0], &ipntr[0], + &workd[0], &workl[0], &lworkl, &info); + else + dnaupd_(&ido, bmat, &n, which, &nev, &tol, + &resid[0], &ncv, &v[0], &ldv, &iparam[0], &ipntr[0], + &workd[0], &workl[0], &lworkl, &info); if (ido == 99) break; @@ -408,7 +500,7 @@ void ArpackSolver::solve (const MatrixType1 &/*system_matrix*/, for (size_type i=0; i z (ldz*ncv, 0.); - double sigmar = 0.0; // real part of the shift double sigmai = 0.0; // imaginary part of the shift - int lworkev = 3*ncv; - std::vector workev (lworkev, 0.); - - std::vector eigenvalues_real (nev, 0.); - std::vector eigenvalues_im (nev, 0.); + std::vector eigenvalues_real (nev+1, 0.); + std::vector eigenvalues_im (nev+1, 0.); - // call of ARPACK dneupd routine - 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, - &iparam[0], &ipntr[0], workd, &workl[0], &lworkl, &info); + // call of ARPACK dseupd/dneupd routine + if (additional_data.symmetric) + { + std::vector z (ldz*nev, 0.); + dseupd_(&rvec, &howmany, &select[0], &eigenvalues_real[0], + &z[0], &ldz, &sigmar, bmat, &n, which, &nev, &tol, + &resid[0], &ncv, &v[0], &ldv, + &iparam[0], &ipntr[0], &workd[0], &workl[0], &lworkl, &info); + } + else + { + std::vector workev (3*ncv, 0.); + dneupd_(&rvec, &howmany, &select[0], &eigenvalues_real[0], + &eigenvalues_im[0], &v[0], &ldz, &sigmar, &sigmai, + &workev[0], bmat, &n, which, &nev, &tol, + &resid[0], &ncv, &v[0], &ldv, + &iparam[0], &ipntr[0], &workd[0], &workl[0], &lworkl, &info); + } if (info == 1) { - Assert (false, ExcArpackInfoMaxIt(control().max_steps())); + Assert (false, ArpackExcArpackInfoMaxIt(control().max_steps())); } else if (info == 3) { - Assert (false, ExcArpackNoShifts()); + Assert (false, ArpackExcArpackNoShifts()); } else if (info!=0) { - Assert (false, ExcArpackInfodneupd(info)); + if (additional_data.symmetric) + { + Assert (false, ArpackExcArpackInfodseupd(info)); + } + else + Assert (false, ArpackExcArpackInfodneupd(info)); } - - const unsigned int n_eigenvecs = eigenvectors.size(); - for (size_type i=0; i (eigenvalues_real[i], eigenvalues_im[i]); } diff --git a/include/deal.II/lac/parpack_solver.h b/include/deal.II/lac/parpack_solver.h index de67f19d17..ba1a80db23 100644 --- a/include/deal.II/lac/parpack_solver.h +++ b/include/deal.II/lac/parpack_solver.h @@ -428,6 +428,11 @@ private: << "Number of wanted eigenvalues " << arg1 << " is larger that the size of eigenvectors " << arg2); + DeclException2 (PArpackExcInvalidEigenvectorSizeNonsymmetric, int, int, + << "To store the real and complex parts of " << arg1 + << " eigenvectors in real-valued vectors, their size (currently set to " << arg2 + << ") should be greater than or equal to " << arg1+1); + DeclException2 (PArpackExcInvalidEigenvalueSize, int, int, << "Number of wanted eigenvalues " << arg1 << " is larger that the size of eigenvalues " << arg2); @@ -603,8 +608,14 @@ void PArpackSolver::solve const unsigned int n_eigenvalues) { - Assert (n_eigenvalues <= eigenvectors.size(), - PArpackExcInvalidEigenvectorSize(n_eigenvalues, eigenvectors.size())); + if (additional_data.symmetric) + { + Assert (n_eigenvalues <= eigenvectors.size(), + PArpackExcInvalidEigenvectorSize(n_eigenvalues, eigenvectors.size())); + } + else + Assert (n_eigenvalues+1 <= eigenvectors.size(), + PArpackExcInvalidEigenvectorSizeNonsymmetric(n_eigenvalues, eigenvectors.size())); Assert (n_eigenvalues <= eigenvalues.size(), PArpackExcInvalidEigenvalueSize(n_eigenvalues, eigenvalues.size())); @@ -872,8 +883,8 @@ void PArpackSolver::solve double sigmar = shift_value; // real part of the shift double sigmai = 0.0; // imaginary part of the shift - std::vector eigenvalues_real (n_eigenvalues, 0.); - std::vector eigenvalues_im (n_eigenvalues, 0.); + std::vector eigenvalues_real (n_eigenvalues+1, 0.); + std::vector eigenvalues_im (n_eigenvalues+1, 0.); // call of ARPACK pdneupd routine if (additional_data.symmetric) @@ -884,7 +895,7 @@ void PArpackSolver::solve &iparam[0], &ipntr[0], &workd[0], &workl[0], &lworkl, &info); else pdneupd_(&mpi_communicator_fortran, &rvec, howmany, &select[0], &eigenvalues_real[0], - &eigenvalues_im[0], &z[0], &ldz, &sigmar, &sigmai, + &eigenvalues_im[0], &v[0], &ldz, &sigmar, &sigmai, &workev[0], bmat, &n_inside_arpack, which, &nev, &tol, &resid[0], &ncv, &v[0], &ldv, &iparam[0], &ipntr[0], &workd[0], &workl[0], &lworkl, &info); @@ -902,7 +913,7 @@ void PArpackSolver::solve AssertThrow (false, PArpackExcInfoPdneupd(info)); } - for (size_type i=0; i::solve } // Throw an error if the solver did not converge. - AssertThrow (iparam[4] == (int)n_eigenvalues, + AssertThrow (iparam[4] >= (int)n_eigenvalues, PArpackExcConvergedEigenvectors(n_eigenvalues,iparam[4])); // both PDNAUPD and PDSAUPD compute eigenpairs of inv[A - sigma*M]*M diff --git a/tests/arpack/arpack_advection_diffusion.cc b/tests/arpack/arpack_advection_diffusion.cc new file mode 100644 index 0000000000..c08966e7cb --- /dev/null +++ b/tests/arpack/arpack_advection_diffusion.cc @@ -0,0 +1,356 @@ +/* --------------------------------------------------------------------- + * + * Copyright (C) 2009 - 2016 by the deal.II authors + * + * This file is part of the deal.II library. + * + * The deal.II library is free software; you can use it, redistribute + * it, and/or modify it under the terms of the GNU Lesser General + * Public License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * The full text of the license can be found in the file LICENSE at + * the top level of the deal.II distribution. + * + * --------------------------------------------------------------------- + + * + * Authors: Toby D. Young, Polish Academy of Sciences, + * Wolfgang Bangerth, Texas A&M University + * Joscha Gedicke, U Heidelberg + * + * This file tests the non-symmetric interface to ARPACK for an advection-diffussion + * operator. The advection is chosen in such a way that we compute complex eigenvalues. + * The most critical case when we cut a complex conjugate pair is tested, i.e. the + * last eigenvalue is complex but the conjugate pair is not included in the range + * of eigenvalues we asked for. This is a particular case in the nonsymmetric arpack + * interface that needs to be taken care of. + * + * We test that the computed vectors are eigenvectors and mass-normal, i.e. + * a) (A*x_i-\lambda*B*x_i).L2() == 0 + * b) x_i*B*x_i = 1 + * + */ + +#include "../tests.h" + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include + +#include +#include +#include + +using namespace dealii; + + +template +class EigenvalueProblem +{ +public: + EigenvalueProblem (unsigned int n_eigenvalues); + void run (); + +private: + void make_grid_and_dofs (); + void assemble_system (); + void solve (); + + Triangulation triangulation; + FE_Q fe; + DoFHandler dof_handler; + + SparsityPattern sparsity_pattern; + SparseMatrix stiffness_matrix, mass_matrix; + std::vector > arpack_vectors; + std::vector > eigenvectors; + std::vector > eigenvalues; + + ConstraintMatrix constraints; + + unsigned int n_eigenvalues; +}; + + + +template +EigenvalueProblem::EigenvalueProblem (unsigned int n_eigenvalues) + : + fe (1), + dof_handler (triangulation), + n_eigenvalues(n_eigenvalues) +{ +} + + + +template +void EigenvalueProblem::make_grid_and_dofs () +{ + GridGenerator::hyper_cube (triangulation, -1, 1); + triangulation.refine_global (5); + dof_handler.distribute_dofs (fe); + + DoFTools::make_zero_boundary_constraints (dof_handler, constraints); + constraints.close (); + + sparsity_pattern.reinit (dof_handler.n_dofs(), + dof_handler.n_dofs(), + dof_handler.max_couplings_between_dofs()); + DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern); + constraints.condense (sparsity_pattern); + sparsity_pattern.compress(); + stiffness_matrix.reinit (sparsity_pattern); + mass_matrix.reinit (sparsity_pattern); + + eigenvalues.resize (n_eigenvalues); + + arpack_vectors.resize (n_eigenvalues + 1); + for (unsigned int i=0; i +void EigenvalueProblem::assemble_system () +{ + QGauss quadrature_formula(2); + + FEValues fe_values (fe, quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + FullMatrix cell_stiffness_matrix (dofs_per_cell, dofs_per_cell); + FullMatrix cell_mass_matrix (dofs_per_cell, dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active (), + endc = dof_handler.end (); + for (; cell!=endc; ++cell) + { + fe_values.reinit (cell); + cell_stiffness_matrix = 0; + cell_mass_matrix = 0; + + for (unsigned int q_point=0; q_point cur_point = fe_values.quadrature_point(q_point); + Tensor<1,dim> advection; + advection[0] = 10.; + advection[1] = 10.*cur_point[0]; + for (unsigned int i=0; iget_dof_indices (local_dof_indices); + + constraints + .distribute_local_to_global (cell_stiffness_matrix, + local_dof_indices, + stiffness_matrix); + constraints + .distribute_local_to_global (cell_mass_matrix, + local_dof_indices, + mass_matrix); + } + + stiffness_matrix.compress (VectorOperation::add); + mass_matrix.compress (VectorOperation::add); + +} + + + +template +void EigenvalueProblem::solve () +{ + SolverControl solver_control (dof_handler.n_dofs(), 1e-10); + SparseDirectUMFPACK inverse; + inverse.initialize (stiffness_matrix); + const unsigned int num_arnoldi_vectors = 2*eigenvalues.size() + 2; + ArpackSolver::AdditionalData additional_data(num_arnoldi_vectors, + ArpackSolver::largest_real_part); + ArpackSolver eigensolver (solver_control, additional_data); + arpack_vectors[0] = 1.; + eigensolver.set_initial_vector(arpack_vectors[0]); + eigensolver.solve (stiffness_matrix, + mass_matrix, + inverse, + eigenvalues, + arpack_vectors, + eigenvalues.size()); + + // extract real and complex components of eigenvectors + for (unsigned int i = 0; i < n_eigenvalues; ++i) + { + eigenvectors[i] = arpack_vectors[i]; + if (eigenvalues[i].imag() != 0.) + { + eigenvectors[i + eigenvalues.size()] = arpack_vectors[i+1]; + if ( i+1 < eigenvalues.size()) + { + eigenvectors[i+1] = arpack_vectors[i]; + eigenvectors[i+1 + eigenvalues.size()] = arpack_vectors[i+1]; + eigenvectors[i+1 + eigenvalues.size()] *= -1; + ++i; + } + } + } + + // make sure that we have eigenvectors and they are mass-normal: + // a) (A*x_i-\lambda*B*x_i).L2() == 0 + // b) x_i*B*x_i = 1 + { + Vector Ax(eigenvectors[0]), Bx(eigenvectors[0]); + Vector Ay(eigenvectors[0]), By(eigenvectors[0]); + for (unsigned int i=0; i < n_eigenvalues; ++i) + { + stiffness_matrix.vmult(Ax,eigenvectors[i]); + stiffness_matrix.vmult(Ay,eigenvectors[i + n_eigenvalues]); + mass_matrix.vmult(Bx,eigenvectors[i]); + mass_matrix.vmult(By,eigenvectors[i + n_eigenvalues]); + + Ax.add(-1.0*std::real(eigenvalues[i]),Bx); + Ax.add(std::imag(eigenvalues[i]),By); + Ay.add(-1.0*std::real(eigenvalues[i]),By); + Ay.add(-1.0*std::imag(eigenvalues[i]),Bx); + Vector tmpx(Ax), tmpy(Ay); + tmpx.scale(Ax); + tmpy.scale(Ay); + tmpx+=tmpy; + if (std::sqrt(tmpx.l1_norm()) > 1e-8) + deallog << "Returned vector " << i << " is not an eigenvector!" + << " L2 norm of the residual is " << std::sqrt(tmpx.l1_norm()) + << std::endl; + + const double tmp = + std::abs( eigenvectors[i] * Bx + eigenvectors[i+n_eigenvalues] * By - 1.) + + std::abs( eigenvectors[i+n_eigenvalues] * Bx - eigenvectors[i] * By); + if ( tmp > 1e-8) + deallog << "Eigenvector " << i << " is not normal! failing norm is" + << tmp << std::endl; + + + } + } +} + +bool my_compare(std::complex a, std::complex b) +{ + if (a.imag()==0.) + return a.imag() < a.imag(); + else + return a.real() < b.real(); +} + +template +void EigenvalueProblem::run () +{ + make_grid_and_dofs (); + + assemble_system (); + + solve (); + + std::sort(eigenvalues.begin(), eigenvalues.end(), my_compare); + + for (unsigned int i = 0; i < n_eigenvalues; ++i) + deallog << " Eigenvalue " << i + << " : " << eigenvalues[i] + << std::endl; +} + + +int main (int argc, char **argv) +{ + + try + { + using namespace dealii; + + std::ofstream logfile("output"); + deallog.attach(logfile); + deallog.threshold_double(1.e-10); + + + EigenvalueProblem<2> problem(4); + problem.run (); + } + + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + + return 0; +} diff --git a/tests/arpack/arpack_advection_diffusion.output b/tests/arpack/arpack_advection_diffusion.output new file mode 100644 index 0000000000..492d63b38c --- /dev/null +++ b/tests/arpack/arpack_advection_diffusion.output @@ -0,0 +1,5 @@ + +DEAL:: Eigenvalue 0 : (31.3555,0.00000) +DEAL:: Eigenvalue 1 : (39.6418,-4.53224) +DEAL:: Eigenvalue 2 : (39.6418,4.53224) +DEAL:: Eigenvalue 3 : (49.9700,-9.19446) diff --git a/tests/arpack/parpack_advection_diffusion_petsc.cc b/tests/arpack/parpack_advection_diffusion_petsc.cc new file mode 100644 index 0000000000..1a9a79b260 --- /dev/null +++ b/tests/arpack/parpack_advection_diffusion_petsc.cc @@ -0,0 +1,434 @@ +/* --------------------------------------------------------------------- + * + * Copyright (C) 2009 - 2016 by the deal.II authors + * + * This file is part of the deal.II library. + * + * The deal.II library is free software; you can use it, redistribute + * it, and/or modify it under the terms of the GNU Lesser General + * Public License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * The full text of the license can be found in the file LICENSE at + * the top level of the deal.II distribution. + * + * --------------------------------------------------------------------- + + * This file tests the non-symmetric interface to PARPACK for an advection-diffussion + * operator with PETSc mpi vectors. + * + * We test that the computed vectors are eigenvectors and mass-normal, i.e. + * a) (A*x_i-\lambda*B*x_i).L2() == 0 + * b) x_i*B*x_i = 1 + * + */ + +#include "../tests.h" + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include + + +#include +#include +#include +#include + +#include + +#include +#include + +const unsigned int dim = 2;//run in 2d to save time + +const double eps = 1e-10; + +template +std::vector +locally_owned_dofs_per_subdomain (const DoFHandlerType &dof_handler) +{ + std::vector< dealii::types::subdomain_id > subdomain_association (dof_handler.n_dofs ()); + dealii::DoFTools::get_subdomain_association (dof_handler, subdomain_association); + + const unsigned int n_subdomains = 1 + (*max_element (subdomain_association.begin (), + subdomain_association.end () )); + + std::vector index_sets (n_subdomains,dealii::IndexSet(dof_handler.n_dofs())); + + // loop over subdomain_association and populate IndexSet when a + // change in subdomain ID is found + dealii::types::global_dof_index i_min = 0; + dealii::types::global_dof_index this_subdomain = subdomain_association[0]; + + for (dealii::types::global_dof_index index = 1; + index < subdomain_association.size (); ++index) + { + //found index different from the current one + if (subdomain_association[index] != this_subdomain) + { + index_sets[this_subdomain].add_range (i_min, index); + i_min = index; + this_subdomain = subdomain_association[index]; + } + } + + // the very last element is of different index + if (i_min == subdomain_association.size () - 1) + { + index_sets[this_subdomain].add_index (i_min); + } + + // otherwise there are at least two different indices + else + { + index_sets[this_subdomain].add_range ( + i_min, subdomain_association.size ()); + } + + for (unsigned int i = 0; i < n_subdomains; i++) + index_sets[i].compress (); + + return index_sets; +} //locally_owned_dofs_per_subdomain + +class PETScInverse +{ +public: + PETScInverse(const dealii::PETScWrappers::MatrixBase &A, dealii::SolverControl &cn,const MPI_Comm &mpi_communicator = PETSC_COMM_SELF): + solver(cn,mpi_communicator), + matrix(A), + preconditioner(matrix) + { + + } + + void vmult ( dealii::PETScWrappers::MPI::Vector &dst, + const dealii::PETScWrappers::MPI::Vector &src) const + { + ; + solver.solve(matrix, dst, src,preconditioner); + } + + +private: + mutable dealii::PETScWrappers::SolverGMRES solver; + const dealii::PETScWrappers::MatrixBase &matrix; + PETScWrappers::PreconditionBlockJacobi preconditioner; + +}; + +void test () +{ + const unsigned int global_mesh_refinement_steps = 5; + const unsigned int number_of_eigenvalues = 4; + + MPI_Comm mpi_communicator = MPI_COMM_WORLD; + const unsigned int n_mpi_processes = dealii::Utilities::MPI::n_mpi_processes(mpi_communicator); + const unsigned int this_mpi_process = dealii::Utilities::MPI::this_mpi_process(mpi_communicator); + + + dealii::Triangulation triangulation; + dealii::DoFHandler dof_handler(triangulation); + dealii::FE_Q fe(1); + dealii::ConstraintMatrix constraints; + dealii::IndexSet locally_owned_dofs; + dealii::IndexSet locally_relevant_dofs; + + std::vector eigenfunctions; + std::vector arpack_vectors; + std::vector> eigenvalues; + dealii::PETScWrappers::MPI::SparseMatrix stiffness_matrix, mass_matrix; + + dealii::GridGenerator::hyper_cube (triangulation, -1, 1); + triangulation.refine_global (global_mesh_refinement_steps); + + // we do not use metis but rather partition by hand below. + //dealii::GridTools::partition_triangulation (n_mpi_processes, triangulation); + { + const double x0 = -1.0; + const double x1 = 1.0; + const double dL = (x1-x0) / n_mpi_processes; + + dealii::Triangulation::active_cell_iterator + cell = triangulation.begin_active(), + endc = triangulation.end(); + for (; cell!=endc; ++cell) + { + const dealii::Point ¢er = cell->center(); + const double x = center[0]; + + const unsigned int id = std::floor ( (x-x0)/dL); + cell->set_subdomain_id (id); + } + } + + dof_handler.distribute_dofs (fe); + dealii::DoFRenumbering::subdomain_wise (dof_handler); + std::vector locally_owned_dofs_per_processor + = locally_owned_dofs_per_subdomain (dof_handler); + locally_owned_dofs = locally_owned_dofs_per_processor[this_mpi_process]; + locally_relevant_dofs.clear(); + dealii::DoFTools::extract_locally_relevant_dofs (dof_handler, + locally_relevant_dofs); + + constraints.clear(); + constraints.reinit (locally_relevant_dofs); + dealii::DoFTools::make_hanging_node_constraints (dof_handler, constraints); + dealii::VectorTools::interpolate_boundary_values (dof_handler, + 0, + dealii::ZeroFunction (), + constraints); + constraints.close (); + + dealii::CompressedSimpleSparsityPattern csp (locally_relevant_dofs); + // Fill in ignoring all cells that are not locally owned + dealii::DoFTools::make_sparsity_pattern (dof_handler, csp, + constraints, + /* keep constrained dofs */ true); + std::vector n_locally_owned_dofs(n_mpi_processes); + for (unsigned int i = 0; i < n_mpi_processes; i++) + n_locally_owned_dofs[i] = locally_owned_dofs_per_processor[i].n_elements(); + + dealii::SparsityTools::distribute_sparsity_pattern + (csp, + n_locally_owned_dofs, + mpi_communicator, + locally_relevant_dofs); + + // Initialise the stiffness and mass matrices + stiffness_matrix.reinit (locally_owned_dofs, + locally_owned_dofs, + csp, + mpi_communicator); + + mass_matrix.reinit (locally_owned_dofs, + locally_owned_dofs, + csp, + mpi_communicator); + + eigenvalues.resize (number_of_eigenvalues); + + arpack_vectors.resize (number_of_eigenvalues+1); + for (unsigned int i=0; i quadrature_formula(2); + dealii::FEValues fe_values (fe, quadrature_formula, + dealii::update_values | + dealii::update_gradients | + dealii::update_quadrature_points | + dealii::update_JxW_values); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + dealii::FullMatrix cell_stiffness_matrix (dofs_per_cell, dofs_per_cell); + dealii::FullMatrix cell_mass_matrix (dofs_per_cell, dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + typename dealii::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active (), + endc = dof_handler.end (); + for (; cell!=endc; ++cell) + if (cell->subdomain_id() == this_mpi_process) + { + fe_values.reinit (cell); + cell_stiffness_matrix = 0; + cell_mass_matrix = 0; + + for (unsigned int q_point=0; q_point cur_point = fe_values.quadrature_point(q_point); + Tensor<1,dim> advection; + advection[0] = 10.; + advection[1] = 10.*cur_point[0]; + for (unsigned int i=0; iget_dof_indices (local_dof_indices); + + constraints + .distribute_local_to_global (cell_stiffness_matrix, + local_dof_indices, + stiffness_matrix); + constraints + .distribute_local_to_global (cell_mass_matrix, + local_dof_indices, + mass_matrix); + } + } + + stiffness_matrix.compress (dealii::VectorOperation::add); + mass_matrix.compress (dealii::VectorOperation::add); + + // test Arpack + { + std::vector > lambda(eigenfunctions.size()); + + for (unsigned int i=0; i < eigenvalues.size(); i++) + eigenfunctions[i] = PetscScalar(); + + dealii::SolverControl solver_control (dof_handler.n_dofs(), 1e-10,/*log_history*/false,/*log_results*/false); + PETScInverse inverse(stiffness_matrix,solver_control,mpi_communicator); + const unsigned int num_arnoldi_vectors = 2*eigenvalues.size() + 2; + + dealii::PArpackSolver::AdditionalData + additional_data(num_arnoldi_vectors, + dealii::PArpackSolver::largest_real_part, + false); + + dealii::PArpackSolver eigensolver (solver_control, + mpi_communicator, + additional_data); + eigensolver.reinit(locally_owned_dofs); + arpack_vectors[0] = 1.; + eigensolver.set_initial_vector(arpack_vectors[0]); + eigensolver.solve (stiffness_matrix, + mass_matrix, + inverse, + eigenvalues, + arpack_vectors, + eigenvalues.size()); + + // extract real and complex components of eigenvectors + for (unsigned int i = 0; i < eigenvalues.size(); ++i) + { + eigenfunctions[i] = arpack_vectors[i]; + if (eigenvalues[i].imag() != 0.) + { + eigenfunctions[i + eigenvalues.size()] = arpack_vectors[i+1]; + if ( i+1 < eigenvalues.size()) + { + eigenfunctions[i+1] = arpack_vectors[i]; + eigenfunctions[i+1 + eigenvalues.size()] = arpack_vectors[i+1]; + eigenfunctions[i+1 + eigenvalues.size()] *= -1; + ++i; + } + } + } + + for (unsigned int i=0; i < eigenvalues.size(); i++) + dealii::deallog << eigenvalues[i] << std::endl; + + // make sure that we have eigenvectors and they are mass-normal: + // a) (A*x_i-\lambda*B*x_i).L2() == 0 + // b) x_i*B*x_i=1 + { + const double precision = 1e-7; + PETScWrappers::MPI::Vector Ax(eigenfunctions[0]), Bx(eigenfunctions[0]); + PETScWrappers::MPI::Vector Ay(eigenfunctions[0]), By(eigenfunctions[0]); + for (unsigned int i=0; i < eigenvalues.size(); ++i) + { + stiffness_matrix.vmult(Ax,eigenfunctions[i]); + stiffness_matrix.vmult(Ay,eigenfunctions[i + eigenvalues.size()]); + mass_matrix.vmult(Bx,eigenfunctions[i]); + mass_matrix.vmult(By,eigenfunctions[i + eigenvalues.size()]); + + Ax.add(-1.0*std::real(eigenvalues[i]),Bx); + Ax.add(std::imag(eigenvalues[i]),By); + Ay.add(-1.0*std::real(eigenvalues[i]),By); + Ay.add(-1.0*std::imag(eigenvalues[i]),Bx); + PETScWrappers::MPI::Vector tmpx(Ax), tmpy(Ay); + tmpx.scale(Ax); + tmpy.scale(Ay); + tmpx+=tmpy; + if (std::sqrt(tmpx.l1_norm()) > precision) + deallog << "Returned vector " << i << " is not an eigenvector!" + << " L2 norm of the residual is " << std::sqrt(tmpx.l1_norm()) + << std::endl; + + const double tmp = + std::abs( eigenfunctions[i] * Bx + eigenfunctions[i+eigenvalues.size()] * By - 1.) + + std::abs( eigenfunctions[i+eigenvalues.size()] * Bx - eigenfunctions[i] * By); + if ( tmp > precision) + deallog << "Eigenvector " << i << " is not normal! failing norm is " + << tmp << std::endl; + } + } + } + + + dof_handler.clear (); + dealii::deallog << "Ok"< +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include + + +#include +#include +#include +#include + +#include + +#include +#include + +const unsigned int dim = 2;//run in 2d to save time + +using namespace dealii; + +const double eps = 1e-10; + +template +std::vector +locally_owned_dofs_per_subdomain (const DoFHandlerType &dof_handler) +{ + std::vector< types::subdomain_id > subdomain_association (dof_handler.n_dofs ()); + DoFTools::get_subdomain_association (dof_handler, subdomain_association); + + const unsigned int n_subdomains = 1 + (*max_element (subdomain_association.begin (), + subdomain_association.end () )); + + std::vector index_sets (n_subdomains,IndexSet(dof_handler.n_dofs())); + + // loop over subdomain_association and populate IndexSet when a + // change in subdomain ID is found + types::global_dof_index i_min = 0; + types::global_dof_index this_subdomain = subdomain_association[0]; + + for (types::global_dof_index index = 1; + index < subdomain_association.size (); ++index) + { + //found index different from the current one + if (subdomain_association[index] != this_subdomain) + { + index_sets[this_subdomain].add_range (i_min, index); + i_min = index; + this_subdomain = subdomain_association[index]; + } + } + + // the very last element is of different index + if (i_min == subdomain_association.size () - 1) + { + index_sets[this_subdomain].add_index (i_min); + } + + // otherwise there are at least two different indices + else + { + index_sets[this_subdomain].add_range ( + i_min, subdomain_association.size ()); + } + + for (unsigned int i = 0; i < n_subdomains; i++) + index_sets[i].compress (); + + return index_sets; +} //locally_owned_dofs_per_subdomain + + +void test () +{ + const unsigned int global_mesh_refinement_steps = 5; + const unsigned int number_of_eigenvalues = 4; + + MPI_Comm mpi_communicator = MPI_COMM_WORLD; + const unsigned int n_mpi_processes = Utilities::MPI::n_mpi_processes(mpi_communicator); + const unsigned int this_mpi_process = Utilities::MPI::this_mpi_process(mpi_communicator); + + + Triangulation triangulation; + DoFHandler dof_handler(triangulation); + FE_Q fe(1); + ConstraintMatrix constraints; + IndexSet locally_owned_dofs; + IndexSet locally_relevant_dofs; + + std::vector eigenfunctions; + std::vector arpack_vectors; + std::vector> eigenvalues; + TrilinosWrappers::SparseMatrix stiffness_matrix, mass_matrix; + + GridGenerator::hyper_cube (triangulation, -1, 1); + triangulation.refine_global (global_mesh_refinement_steps); + + // we do not use metis but rather partition by hand below. + //dealii::GridTools::partition_triangulation (n_mpi_processes, triangulation); + { + const double x0 = -1.0; + const double x1 = 1.0; + const double dL = (x1-x0) / n_mpi_processes; + + Triangulation::active_cell_iterator + cell = triangulation.begin_active(), + endc = triangulation.end(); + for (; cell!=endc; ++cell) + { + const Point ¢er = cell->center(); + const double x = center[0]; + + const unsigned int id = std::floor ( (x-x0)/dL); + cell->set_subdomain_id (id); + } + } + + + dof_handler.distribute_dofs (fe); + DoFRenumbering::subdomain_wise (dof_handler); + std::vector locally_owned_dofs_per_processor + = locally_owned_dofs_per_subdomain (dof_handler); + locally_owned_dofs = locally_owned_dofs_per_processor[this_mpi_process]; + locally_relevant_dofs.clear(); + DoFTools::extract_locally_relevant_dofs (dof_handler, + locally_relevant_dofs); + + constraints.clear(); + constraints.reinit (locally_relevant_dofs); + DoFTools::make_hanging_node_constraints (dof_handler, constraints); + VectorTools::interpolate_boundary_values (dof_handler, + 0, + ZeroFunction (), + constraints); + constraints.close (); + + CompressedSimpleSparsityPattern csp (locally_relevant_dofs); + // Fill in ignoring all cells that are not locally owned + DoFTools::make_sparsity_pattern (dof_handler, csp, + constraints, + /* keep constrained dofs */ true); + std::vector n_locally_owned_dofs(n_mpi_processes); + for (unsigned int i = 0; i < n_mpi_processes; i++) + n_locally_owned_dofs[i] = locally_owned_dofs_per_processor[i].n_elements(); + + SparsityTools::distribute_sparsity_pattern + (csp, + n_locally_owned_dofs, + mpi_communicator, + locally_relevant_dofs); + + // Initialise the stiffness and mass matrices + stiffness_matrix.reinit (locally_owned_dofs, + locally_owned_dofs, + csp, + mpi_communicator); + + mass_matrix.reinit (locally_owned_dofs, + locally_owned_dofs, + csp, + mpi_communicator); + + eigenvalues.resize (number_of_eigenvalues); + + arpack_vectors.resize (number_of_eigenvalues+1); + for (unsigned int i=0; i quadrature_formula(2); + FEValues fe_values (fe, quadrature_formula, + update_values | + update_gradients | + update_quadrature_points | + update_JxW_values); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + FullMatrix cell_stiffness_matrix (dofs_per_cell, dofs_per_cell); + FullMatrix cell_mass_matrix (dofs_per_cell, dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active (), + endc = dof_handler.end (); + for (; cell!=endc; ++cell) + if (cell->subdomain_id() == this_mpi_process) + { + fe_values.reinit (cell); + cell_stiffness_matrix = 0; + cell_mass_matrix = 0; + + for (unsigned int q_point=0; q_point cur_point = fe_values.quadrature_point(q_point); + Tensor<1,dim> advection; + advection[0] = 10.; + advection[1] = 10.*cur_point[0]; + for (unsigned int i=0; iget_dof_indices (local_dof_indices); + + constraints + .distribute_local_to_global (cell_stiffness_matrix, + local_dof_indices, + stiffness_matrix); + constraints + .distribute_local_to_global (cell_mass_matrix, + local_dof_indices, + mass_matrix); + } + } + + stiffness_matrix.compress (VectorOperation::add); + mass_matrix.compress (VectorOperation::add); + + // test Arpack + { + const double shift = 4.0; + std::vector > lambda(eigenfunctions.size()); + + for (unsigned int i=0; i < eigenvalues.size(); i++) + eigenfunctions[i] = 0.; + + SolverControl solver_control (dof_handler.n_dofs(), 1e-9,/*log_history*/false,/*log_results*/false); + SolverControl solver_control_lin (dof_handler.n_dofs(), 1e-10,/*log_history*/false,/*log_results*/false); + + PArpackSolver::Shift shifted_matrix(stiffness_matrix,mass_matrix,shift); + TrilinosWrappers::PreconditionIdentity preconditioner; + IterativeInverse shift_and_invert; + shift_and_invert.initialize(shifted_matrix,preconditioner); + shift_and_invert.solver.select("gmres"); + static ReductionControl inner_control_c(/*maxiter*/stiffness_matrix.m(), + /*tolerance (global)*/ 0.0, + /*reduce (w.r.t. initial)*/ 1.e-13); + shift_and_invert.solver.set_control(inner_control_c); + + const unsigned int num_arnoldi_vectors = 2*eigenvalues.size() + 2; + + PArpackSolver::AdditionalData + additional_data(num_arnoldi_vectors, + PArpackSolver::largest_real_part, + false); + + PArpackSolver eigensolver (solver_control, + mpi_communicator, + additional_data); + eigensolver.reinit(locally_owned_dofs); + eigensolver.set_shift(shift); + arpack_vectors[0] = 1.; + eigensolver.set_initial_vector(arpack_vectors[0]); + // avoid output of iterative solver: + const unsigned int previous_depth = deallog.depth_file(0); + eigensolver.solve (stiffness_matrix, + mass_matrix, + shift_and_invert, + eigenvalues, + arpack_vectors, + eigenvalues.size()); + deallog.depth_file(previous_depth); + + // extract real and complex components of eigenvectors + for (unsigned int i = 0; i < eigenvalues.size(); ++i) + { + eigenfunctions[i] = arpack_vectors[i]; + if (eigenvalues[i].imag() != 0.) + { + eigenfunctions[i + eigenvalues.size()] = arpack_vectors[i+1]; + if ( i+1 < eigenvalues.size()) + { + eigenfunctions[i+1] = arpack_vectors[i]; + eigenfunctions[i+1 + eigenvalues.size()] = arpack_vectors[i+1]; + eigenfunctions[i+1 + eigenvalues.size()] *= -1; + ++i; + } + } + } + + for (unsigned int i=0; i < eigenvalues.size(); i++) + deallog << eigenvalues[i] << std::endl; + + // make sure that we have eigenvectors and they are mass-normal: + // a) (A*x_i-\lambda*B*x_i).L2() == 0 + // b) x_i*B*x_i=1 + { + const double precision = 1e-7; + TrilinosWrappers::MPI::Vector Ax(eigenfunctions[0]), Bx(eigenfunctions[0]); + TrilinosWrappers::MPI::Vector Ay(eigenfunctions[0]), By(eigenfunctions[0]); + for (unsigned int i=0; i < eigenvalues.size(); ++i) + { + stiffness_matrix.vmult(Ax,eigenfunctions[i]); + stiffness_matrix.vmult(Ay,eigenfunctions[i + eigenvalues.size()]); + mass_matrix.vmult(Bx,eigenfunctions[i]); + mass_matrix.vmult(By,eigenfunctions[i + eigenvalues.size()]); + + Ax.add(-1.0*std::real(eigenvalues[i]),Bx); + Ax.add(std::imag(eigenvalues[i]),By); + Ay.add(-1.0*std::real(eigenvalues[i]),By); + Ay.add(-1.0*std::imag(eigenvalues[i]),Bx); + TrilinosWrappers::MPI::Vector tmpx(Ax), tmpy(Ay); + tmpx.scale(Ax); + tmpy.scale(Ay); + tmpx+=tmpy; + if (std::sqrt(tmpx.l1_norm()) > precision) + deallog << "Returned vector " << i << " is not an eigenvector!" + << " L2 norm of the residual is " << std::sqrt(tmpx.l1_norm()) + << std::endl; + + const double tmp = + std::abs( eigenfunctions[i] * Bx + eigenfunctions[i+eigenvalues.size()] * By - 1.) + + std::abs( eigenfunctions[i+eigenvalues.size()] * Bx - eigenfunctions[i] * By); + if ( tmp > precision) + deallog << "Eigenvector " << i << " is not normal! failing norm is " + << tmp << std::endl; + } + + } + } + + + dof_handler.clear (); + deallog << "Ok"< 1e-8) + deallog << "Eigenvectors " + + Utilities::int_to_string(i) + + " and " + + Utilities::int_to_string(j) + + " are not orthonormal!" + " failing norm is " + + Utilities::to_string(std::abs( eigenfunctions[j] * Bx - (i==j))) + << std::endl; stiffness_matrix.vmult(Ax,eigenfunctions[i]); Ax.add(-1.0*std::real(eigenvalues[i]),Bx); - Assert (Ax.l2_norm() < 1e-8, - ExcMessage("Returned vector " + - Utilities::int_to_string(i) + - " is not an eigenvector!" - " L2 norm of the residual is " + - Utilities::to_string(Ax.l2_norm()) - )); + if (Ax.l2_norm() > 1e-8) + deallog << "Returned vector " + + Utilities::int_to_string(i) + + " is not an eigenvector!" + " L2 norm of the residual is " + + Utilities::to_string(Ax.l2_norm()) + << std::endl; } } for (unsigned int i=0; i diff --git a/tests/arpack/step-36_parpack_trilinos.cc b/tests/arpack/step-36_parpack_trilinos.cc index b063152c53..c02481b66a 100644 --- a/tests/arpack/step-36_parpack_trilinos.cc +++ b/tests/arpack/step-36_parpack_trilinos.cc @@ -1,3 +1,28 @@ +/* --------------------------------------------------------------------- + * + * Copyright (C) 2009 - 2016 by the deal.II authors + * + * This file is part of the deal.II library. + * + * The deal.II library is free software; you can use it, redistribute + * it, and/or modify it under the terms of the GNU Lesser General + * Public License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * The full text of the license can be found in the file LICENSE at + * the top level of the deal.II distribution. + * + * --------------------------------------------------------------------- + + * + * This file tests the PARPACK interface for a symmetric operator taken from step-36 + * using Trilinos mpi vectors. + * + * We test that the computed vectors are eigenvectors and mass-orthonormal, i.e. + * a) (A*x_i-\lambda*B*x_i).L2() == 0 + * b) x_j*B*x_i = \delta_{i,j} + * + */ + #include "../tests.h" #include