-> 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
<h3>Specific improvements</h3>
<ol>
+ <li> Fixed: (P)ARPACK interface for non-symmetric matrices.
+ <br>
+ (Joscha Gedicke, 2016/08/01)
+
<li> Fixed: The TrilinosWrappers::SparsityPattern::print() and
TrilinosWrappers::SparsityPattern::print_gnuplot() methods did not produce
correct output on distributed computations. This is now fixed.
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 <code>dneupd</code> and <code>dnaupd</code> of ARPACK. The
+ * the routines <code>dnaupd</code> and <code>dneupd</code> of ARPACK.
+ * If the operator is specified to be symmetric we use the symmetric interface
+ * <code>dsaupd</code> and <code>dseupd</code> 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.
* Through the AdditionalData the user can specify some of the parameters to
* be set.
*
- * For further information on how the ARPACK routines <code>dneupd</code> and
- * <code>dnaupd</code> work and also how to set the parameters appropriately
+ * For further information on how the ARPACK routines <code>dsaupd</code>,
+ * <code>dseupd</code>, <code>dnaupd</code> and <code>dneupd</code> 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
* @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
{
{
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);
};
/**
ArpackSolver(SolverControl &control,
const AdditionalData &data = AdditionalData());
+ /**
+ * Set initial vector for building Krylov space.
+ */
+ template <typename VectorType>
+ void set_initial_vector(const VectorType &vec);
+
/**
* Solve the generalized eigensprectrum problem $A x=\lambda B x$ by calling
- * the <code>dneupd</code> and <code>dnaupd</code> functions of ARPACK.
+ * the <code>dsaupd</code> and <code>dseupd</code> or
+ * <code>dnaupd</code> and <code>dneupd</code> functions of ARPACK.
*
* The function returns a vector of eigenvalues of length <i>n</i> and a
- * vector of eigenvectors, where the latter should be twice the size of the
- * eigenvalue vector. The first <i>n</i> vectors in
- * <code>eigenvectors</code> will be the real parts of the eigenvectors, the
- * second <i>n</i> the imaginary parts.
+ * vector of eigenvectors of length <i>n</i> in the symmetric case
+ * and of length <i>n+1</i> 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.
* eigenvalues are returned.
*
* @param eigenvectors is a <b>real</b> 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 <i>n</i> in the symmetric case and <i>n+1</i>
+ * 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 <i>[real(v1),real(v2),
+ * real(v3),imag(v3)]</i>. 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 <code>eigenvalues</code> or greater.
*/
const AdditionalData additional_data;
+ /**
+ * Store an initial vector
+ */
+ bool initial_vector_provided;
+ std::vector<double> 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.");
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)
{}
const AdditionalData &data)
:
solver_control (control),
- additional_data (data)
-
+ additional_data (data),
+ initial_vector_provided(false)
{}
+template <typename VectorType>
+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 <typename VectorType, typename MatrixType1,
typename MatrixType2, typename INVERSE>
std::vector<VectorType> &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
}
// 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<double> resid(n, 1.);
+ if (!initial_vector_provided || resid.size() != n)
+ resid.resize(n, 1.);
// number of Arnoldi basis vectors specified
// in additional_data
std::vector<int> 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<double> workd (3*n, 0.);
+ int lworkl = additional_data.symmetric ? ncv*ncv + 8*ncv : 3*ncv*ncv+6*ncv;
std::vector<double> 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;
for (size_type i=0; i<src.size(); ++i)
- src(i) = *(workd+ipntr[0]-1+i);
+ src(i) = workd[ipntr[0]-1+i];
// multiplication with mass matrix M
mass_matrix.vmult(tmp, src);
inverse.vmult(dst,tmp);
for (size_type i=0; i<dst.size(); ++i)
- *(workd+ipntr[1]-1+i) = dst(i);
+ workd[ipntr[1]-1+i] = dst(i);
}
break;
for (size_type i=0; i<src.size(); ++i)
{
- src(i) = *(workd+ipntr[2]-1+i);
- tmp(i) = *(workd+ipntr[0]-1+i);
+ src(i) = workd[ipntr[2]-1+i];
+ tmp(i) = workd[ipntr[0]-1+i];
}
// solving linear system
inverse.vmult(dst,src);
for (size_type i=0; i<dst.size(); ++i)
- *(workd+ipntr[1]-1+i) = dst(i);
+ workd[ipntr[1]-1+i] = dst(i);
}
break;
dst.reinit(src);
for (size_type i=0; i<src.size(); ++i)
- src(i) = *(workd+ipntr[0]-1+i);
+ src(i) = workd[ipntr[0]-1+i];
// Multiplication with mass matrix M
mass_matrix.vmult(dst, src);
for (size_type i=0; i<dst.size(); ++i)
- *(workd+ipntr[1]-1+i) = dst(i);
+ workd[ipntr[1]-1+i] = dst(i);
}
break;
default:
- Assert (false, ExcArpackIdo(ido));
+ Assert (false, ArpackExcArpackIdo(ido));
break;
}
}
break;
default:
- Assert (false, ExcArpackMode(mode));
+ Assert (false, ArpackExcArpackMode(mode));
break;
}
}
if (info<0)
{
- Assert (false, ExcArpackInfodsaupd(info));
+ if (additional_data.symmetric)
+ {
+ Assert (false, ArpackExcArpackInfodsaupd(info));
+ }
+ else
+ Assert (false, ArpackExcArpackInfodnaupd(info));
}
else
{
int ldz = n;
- std::vector<double> 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<double> workev (lworkev, 0.);
-
- std::vector<double> eigenvalues_real (nev, 0.);
- std::vector<double> eigenvalues_im (nev, 0.);
+ std::vector<double> eigenvalues_real (nev+1, 0.);
+ std::vector<double> 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<double> 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<double> 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<n_eigenvecs; ++i)
+ for (unsigned int i=0; i<nev; ++i)
for (unsigned int j=0; j<n; ++j)
eigenvectors[i](j) = v[i*n+j];
- delete[] workd;
-
- AssertDimension (eigenvalues.size(), eigenvalues_real.size());
- AssertDimension (eigenvalues.size(), eigenvalues_im.size());
-
- for (size_type i=0; i<eigenvalues.size(); ++i)
+ for (unsigned int i=0; i<nev_const; ++i)
eigenvalues[i] = std::complex<double> (eigenvalues_real[i],
eigenvalues_im[i]);
}
<< "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);
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()));
double sigmar = shift_value; // real part of the shift
double sigmai = 0.0; // imaginary part of the shift
- std::vector<double> eigenvalues_real (n_eigenvalues, 0.);
- std::vector<double> eigenvalues_im (n_eigenvalues, 0.);
+ std::vector<double> eigenvalues_real (n_eigenvalues+1, 0.);
+ std::vector<double> eigenvalues_im (n_eigenvalues+1, 0.);
// call of ARPACK pdneupd routine
if (additional_data.symmetric)
&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);
AssertThrow (false, PArpackExcInfoPdneupd(info));
}
- for (size_type i=0; i<n_eigenvalues; ++i)
+ for (size_type i=0; i<nev; ++i)
{
eigenvectors[i] = 0.0;
Assert (i*nloc + nloc <= v.size(), dealii::ExcInternalError() );
}
// 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
--- /dev/null
+/* ---------------------------------------------------------------------
+ *
+ * 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 <deal.II/base/logstream.h>
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/function_parser.h>
+#include <deal.II/base/parameter_handler.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/lac/full_matrix.h>
+
+#include <deal.II/lac/arpack_solver.h>
+#include <deal.II/lac/sparse_direct.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+
+#include <complex>
+
+#include <fstream>
+#include <iostream>
+#include <algorithm>
+
+using namespace dealii;
+
+
+template <int dim>
+class EigenvalueProblem
+{
+public:
+ EigenvalueProblem (unsigned int n_eigenvalues);
+ void run ();
+
+private:
+ void make_grid_and_dofs ();
+ void assemble_system ();
+ void solve ();
+
+ Triangulation<dim> triangulation;
+ FE_Q<dim> fe;
+ DoFHandler<dim> dof_handler;
+
+ SparsityPattern sparsity_pattern;
+ SparseMatrix<double> stiffness_matrix, mass_matrix;
+ std::vector<Vector<double> > arpack_vectors;
+ std::vector<Vector<double> > eigenvectors;
+ std::vector<std::complex<double> > eigenvalues;
+
+ ConstraintMatrix constraints;
+
+ unsigned int n_eigenvalues;
+};
+
+
+
+template <int dim>
+EigenvalueProblem<dim>::EigenvalueProblem (unsigned int n_eigenvalues)
+ :
+ fe (1),
+ dof_handler (triangulation),
+ n_eigenvalues(n_eigenvalues)
+{
+}
+
+
+
+template <int dim>
+void EigenvalueProblem<dim>::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<arpack_vectors.size (); ++i)
+ arpack_vectors[i].reinit (dof_handler.n_dofs ());
+
+ eigenvectors.resize(2*n_eigenvalues);
+ for (unsigned int i = 0; i < eigenvectors.size(); ++i)
+ eigenvectors[i].reinit (dof_handler.n_dofs ());
+
+}
+
+
+
+template <int dim>
+void EigenvalueProblem<dim>::assemble_system ()
+{
+ QGauss<dim> quadrature_formula(2);
+
+ FEValues<dim> 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<double> cell_stiffness_matrix (dofs_per_cell, dofs_per_cell);
+ FullMatrix<double> cell_mass_matrix (dofs_per_cell, dofs_per_cell);
+
+ std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
+
+ typename DoFHandler<dim>::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<n_q_points; ++q_point)
+ {
+ const Point<dim> 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; i<dofs_per_cell; ++i)
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ {
+ cell_stiffness_matrix (i, j)
+ += (fe_values.shape_grad (i, q_point) *
+ fe_values.shape_grad (j, q_point)
+ +
+ (advection * fe_values.shape_grad (i, q_point)) *
+ fe_values.shape_value (j, q_point)
+ ) * fe_values.JxW (q_point);
+
+ cell_mass_matrix (i, j)
+ += (fe_values.shape_value (i, q_point) *
+ fe_values.shape_value (j, q_point)
+ ) * fe_values.JxW (q_point);
+ }
+ }
+
+ cell->get_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 <int dim>
+void EigenvalueProblem<dim>::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<double> Ax(eigenvectors[0]), Bx(eigenvectors[0]);
+ Vector<double> 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<double> 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<double> a, std::complex<double> b)
+{
+ if (a.imag()==0.)
+ return a.imag() < a.imag();
+ else
+ return a.real() < b.real();
+}
+
+template <int dim>
+void EigenvalueProblem<dim>::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;
+}
--- /dev/null
+
+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)
--- /dev/null
+/* ---------------------------------------------------------------------
+ *
+ * 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 <deal.II/base/logstream.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/index_set.h>
+
+#include <deal.II/dofs/dof_renumbering.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/fe_tools.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/lac/sparsity_tools.h>
+#include <deal.II/lac/compressed_simple_sparsity_pattern.h>
+
+#include <deal.II/lac/petsc_parallel_vector.h>
+#include <deal.II/lac/petsc_parallel_sparse_matrix.h>
+#include <deal.II/lac/petsc_solver.h>
+#include <deal.II/lac/petsc_precondition.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_in.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/lac/parpack_solver.h>
+
+#include <fstream>
+#include <iostream>
+
+const unsigned int dim = 2;//run in 2d to save time
+
+const double eps = 1e-10;
+
+template <typename DoFHandlerType>
+std::vector<dealii::IndexSet>
+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<dealii::IndexSet> 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<dim> triangulation;
+ dealii::DoFHandler<dim> dof_handler(triangulation);
+ dealii::FE_Q<dim> fe(1);
+ dealii::ConstraintMatrix constraints;
+ dealii::IndexSet locally_owned_dofs;
+ dealii::IndexSet locally_relevant_dofs;
+
+ std::vector<dealii::PETScWrappers::MPI::Vector> eigenfunctions;
+ std::vector<dealii::PETScWrappers::MPI::Vector> arpack_vectors;
+ std::vector<std::complex<PetscScalar>> 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<dim>::active_cell_iterator
+ cell = triangulation.begin_active(),
+ endc = triangulation.end();
+ for (; cell!=endc; ++cell)
+ {
+ const dealii::Point<dim> ¢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<dealii::IndexSet> 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<dim> (),
+ 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<dealii::types::global_dof_index> 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<arpack_vectors.size (); ++i)
+ arpack_vectors[i].reinit (locally_owned_dofs, mpi_communicator);//without ghost dofs
+
+ eigenfunctions.resize (2*number_of_eigenvalues);
+ for (unsigned int i=0; i<eigenfunctions.size (); ++i)
+ eigenfunctions[i].reinit (locally_owned_dofs, mpi_communicator);//without ghost dofs
+
+
+ // ready for assembly
+ stiffness_matrix = 0;
+ mass_matrix = 0;
+
+ dealii::QGauss<dim> quadrature_formula(2);
+ dealii::FEValues<dim> 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<double> cell_stiffness_matrix (dofs_per_cell, dofs_per_cell);
+ dealii::FullMatrix<double> cell_mass_matrix (dofs_per_cell, dofs_per_cell);
+
+ std::vector<dealii::types::global_dof_index> local_dof_indices (dofs_per_cell);
+
+ typename dealii::DoFHandler<dim>::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<n_q_points; ++q_point)
+ {
+ const Point<dim> 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; i<dofs_per_cell; ++i)
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ {
+ cell_stiffness_matrix (i, j)
+ += (fe_values.shape_grad (i, q_point) *
+ fe_values.shape_grad (j, q_point)
+ +
+ (advection * fe_values.shape_grad (i, q_point)) *
+ fe_values.shape_value (j, q_point)
+ ) * fe_values.JxW (q_point);
+
+ cell_mass_matrix (i, j)
+ += (fe_values.shape_value (i, q_point) *
+ fe_values.shape_value (j, q_point)
+ ) * fe_values.JxW (q_point);
+ }
+
+ cell->get_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<std::complex<double> > 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<dealii::PETScWrappers::MPI::Vector>::AdditionalData
+ additional_data(num_arnoldi_vectors,
+ dealii::PArpackSolver<dealii::PETScWrappers::MPI::Vector>::largest_real_part,
+ false);
+
+ dealii::PArpackSolver<dealii::PETScWrappers::MPI::Vector> 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"<<std::endl;
+}
+
+
+int main (int argc,char **argv)
+{
+ std::ofstream logfile("output");
+ dealii::deallog.attach(logfile,/*do not print job id*/false);
+ dealii::deallog.threshold_double(eps);
+
+ try
+ {
+ dealii::Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+ {
+ test ();
+ }
+
+ }
+ 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;
+ };
+}
--- /dev/null
+DEAL::(32.2359,0.00000)
+DEAL::(40.3938,-4.48216)
+DEAL::(40.3938,4.48216)
+DEAL::(50.5706,-9.07560)
+DEAL::Ok
--- /dev/null
+/* ---------------------------------------------------------------------
+ *
+ * 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 Trilinos 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 <deal.II/base/logstream.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/index_set.h>
+
+#include <deal.II/dofs/dof_renumbering.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/fe_tools.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/lac/sparsity_tools.h>
+#include <deal.II/lac/compressed_simple_sparsity_pattern.h>
+
+#include <deal.II/lac/iterative_inverse.h>
+
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/lac/trilinos_precondition.h>
+#include <deal.II/numerics/vector_tools.h>
+
+
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_in.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/lac/parpack_solver.h>
+
+#include <fstream>
+#include <iostream>
+
+const unsigned int dim = 2;//run in 2d to save time
+
+using namespace dealii;
+
+const double eps = 1e-10;
+
+template <typename DoFHandlerType>
+std::vector<IndexSet>
+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<IndexSet> 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<dim> triangulation;
+ DoFHandler<dim> dof_handler(triangulation);
+ FE_Q<dim> fe(1);
+ ConstraintMatrix constraints;
+ IndexSet locally_owned_dofs;
+ IndexSet locally_relevant_dofs;
+
+ std::vector<TrilinosWrappers::MPI::Vector> eigenfunctions;
+ std::vector<TrilinosWrappers::MPI::Vector> arpack_vectors;
+ std::vector<std::complex<double>> 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<dim>::active_cell_iterator
+ cell = triangulation.begin_active(),
+ endc = triangulation.end();
+ for (; cell!=endc; ++cell)
+ {
+ const Point<dim> ¢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<IndexSet> 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<dim> (),
+ 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<types::global_dof_index> 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<arpack_vectors.size (); ++i)
+ arpack_vectors[i].reinit (locally_owned_dofs, mpi_communicator);//without ghost dofs
+
+ eigenfunctions.resize (2*number_of_eigenvalues);
+ for (unsigned int i=0; i<eigenfunctions.size (); ++i)
+ eigenfunctions[i].reinit (locally_owned_dofs, mpi_communicator);//without ghost dofs
+
+ // ready for assembly
+ stiffness_matrix = 0;
+ mass_matrix = 0;
+
+ QGauss<dim> quadrature_formula(2);
+ FEValues<dim> 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<double> cell_stiffness_matrix (dofs_per_cell, dofs_per_cell);
+ FullMatrix<double> cell_mass_matrix (dofs_per_cell, dofs_per_cell);
+
+ std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
+
+ typename DoFHandler<dim>::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<n_q_points; ++q_point)
+ {
+ const Point<dim> 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; i<dofs_per_cell; ++i)
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ {
+ cell_stiffness_matrix (i, j)
+ += (fe_values.shape_grad (i, q_point) *
+ fe_values.shape_grad (j, q_point)
+ +
+ (advection * fe_values.shape_grad (i, q_point)) *
+ fe_values.shape_value (j, q_point)
+ ) * fe_values.JxW (q_point);
+
+ cell_mass_matrix (i, j)
+ += (fe_values.shape_value (i, q_point) *
+ fe_values.shape_value (j, q_point)
+ ) * fe_values.JxW (q_point);
+ }
+
+ cell->get_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<std::complex<double> > 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<TrilinosWrappers::MPI::Vector>::Shift<TrilinosWrappers::SparseMatrix> shifted_matrix(stiffness_matrix,mass_matrix,shift);
+ TrilinosWrappers::PreconditionIdentity preconditioner;
+ IterativeInverse<TrilinosWrappers::MPI::Vector > 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<TrilinosWrappers::MPI::Vector>::AdditionalData
+ additional_data(num_arnoldi_vectors,
+ PArpackSolver<TrilinosWrappers::MPI::Vector>::largest_real_part,
+ false);
+
+ PArpackSolver<TrilinosWrappers::MPI::Vector> 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"<<std::endl;
+}
+
+
+int main (int argc,char **argv)
+{
+ std::ofstream logfile("output");
+ deallog.attach(logfile,/*do not print job id*/false);
+ deallog.threshold_double(eps);
+
+ try
+ {
+ Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+ {
+ test ();
+ }
+
+ }
+ 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;
+ };
+}
--- /dev/null
+DEAL::(32.2359,0.00000)
+DEAL::(40.3938,-4.48216)
+DEAL::(40.3938,4.48216)
+DEAL::(50.5706,-9.07560)
+DEAL::Ok
/* ---------------------------------------------------------------------
*
- * Copyright (C) 2009 - 2015 by the deal.II authors
+ * Copyright (C) 2009 - 2016 by the deal.II authors
*
* This file is part of the deal.II library.
*
*
* Authors: Toby D. Young, Polish Academy of Sciences,
* Wolfgang Bangerth, Texas A&M University
+ *
+ * This file tests the ARPACK interface for a symmetric operator taken from step-36.
+ *
+ * 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"
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::largest_magnitude,
+ true);
ArpackSolver eigensolver (solver_control, additional_data);
eigensolver.solve (stiffness_matrix,
mass_matrix,
mass_matrix.vmult(Bx,eigenfunctions[i]);
for (unsigned int j=0; j < eigenfunctions.size(); j++)
- Assert( std::abs( eigenfunctions[j] * Bx - (i==j))< 1e-8,
- ExcMessage("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) )
- )
- ));
+ if ( std::abs( eigenfunctions[j] * Bx - (i==j))> 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<eigenfunctions.size(); ++i)
+/* ---------------------------------------------------------------------
+ *
+ * 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 PETSc 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 <deal.II/base/logstream.h>
+/* ---------------------------------------------------------------------
+ *
+ * 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 <deal.II/base/logstream.h>