+++ /dev/null
-// ---------------------------------------------------------------------
-//
-// Copyright (C) 2010 - 2013 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.
-//
-// ---------------------------------------------------------------------
-
-#ifndef __deal2__parpack_solver_h
-#define __deal2__parpack_solver_h
-
-#include <deal.II/base/config.h>
-#include <deal.II/base/smartpointer.h>
-#include <deal.II/base/memory_consumption.h>
-#include <deal.II/lac/solver_control.h>
-#include <deal.II/lac/vector.h>
-#include <deal.II/base/index_set.h>
-
-#include <cstring>
-
-
-#ifdef DEAL_II_ARPACK_WITH_PARPACK
-
-DEAL_II_NAMESPACE_OPEN
-
-extern "C" {
-
- // http://www.mathkeisan.com/usersguide/man/pdnaupd.html
- void pdnaupd_(MPI_Fint *comm, int *ido, char *bmat, int *n, char *which,
- int *nev, double *tol, double *resid, int *ncv,
- double *v, int *nloc, int *iparam, int *ipntr,
- double *workd, double *workl, int *lworkl,
- int *info);
-
- // http://www.mathkeisan.com/usersguide/man/pdsaupd.html
- void pdsaupd_(MPI_Fint *comm, int *ido, char *bmat, int *n, char *which,
- int *nev, double *tol, double *resid, int *ncv,
- double *v, int *nloc, int *iparam, int *ipntr,
- double *workd, double *workl, int *lworkl,
- int *info);
-
- // http://www.mathkeisan.com/usersguide/man/pdneupd.html
- void pdneupd_(MPI_Fint *comm, int *rvec, char *howmany, int *select, double *d,
- double *di, double *z, int *ldz, double *sigmar,
- double *sigmai, double *workev, char *bmat, int *n, char *which,
- int *nev, double *tol, double *resid, int *ncv,
- double *v, int *nloc, int *iparam, int *ipntr,
- double *workd, double *workl, int *lworkl, int *info);
-
- // http://www.mathkeisan.com/usersguide/man/pdseupd.html
- void pdseupd_(MPI_Fint *comm, int *rvec, char *howmany, int *select, double *d,
- double *z, int *ldz, double *sigmar,
- char *bmat, int *n, char *which,
- int *nev, double *tol, double *resid, int *ncv,
- double *v, int *nloc, int *iparam, int *ipntr,
- double *workd, double *workl, int *lworkl, int *info);
-
- // other resources:
- // http://acts.nersc.gov/superlu/example5/pnslac.c.html
- // https://github.com/phpisciuneri/tijo/blob/master/dvr_parpack.cpp
-
-}
-
-/**
- * Interface for using PARPACK. PARPACK is a collection of Fortran77
- * subroutines designed to solve large scale eigenvalue problems.
- * Here we interface to the routines <code>pdneupd</code>,
- * <code>pdseupd</code>, <code>pdnaupd</code>, <code>pdsaupd</code> of
- * PARPACK. 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.
- *
- * In this class we make use of the method applied to the generalized
- * eigenspectrum problem $(A-\lambda B)x=0$, for $x\neq0$; where $A$
- * is a system matrix, $B$ is a mass matrix, and $\lambda, x$ are a
- * set of eigenvalues and eigenvectors respectively.
- *
- * The ArpackSolver can be used in application codes in the
- * following way:
- @code
- SolverControl solver_control (1000, 1e-9);
- const unsigned int num_arnoldi_vectors = 2*size_of_spectrum + 2;
- PArpackSolver<V>::AdditionalData
- additional_data(num_arnoldi_vectors,
- dealii::PArpackSolver<V>::largest_magnitude,
- true);
-
- PArpackSolver<V> eigensolver (solver_control,
- mpi_communicator,
- additional_data);
- eigensolver.set_shift(sigma);
- eigensolver.reinit(locally_owned_dofs);
- eigensolver.solve (A,
- B,
- OP,
- lambda,
- x,
- size_of_spectrum);
- @endcode
- * for the generalized eigenvalue problem $Ax=B\lambda x$, where the
- * variable <code>size_of_spectrum</code> tells PARPACK the number of
- * eigenvector/eigenvalue pairs to solve for. Here,
- * <code>lambda</code> is a vector that will contain the eigenvalues
- * computed, <code>x</code> a vector of objects of type <code>V</code>
- * that will contain the eigenvectors computed. <code>OP</code> is an
- * inverse operation for the matrix <code>A - sigma * B</code>, where
- * <code> sigma </code> is a shift value, set to zero by default.
- *
- * Through the AdditionalData the user can specify some of the
- * parameters to be set.
- *
- * The class is intended to be used with MPI and can work on arbitrary
- * vector and matrix distributed classes. Both symmetric and
- * non-symmetric <code>A</code> are supported.
- *
- * For further information on how the PARPACK routines
- * <code>pdneupd</code>, <code>pdseupd</code>, <code>pdnaupd</code>,
- * <code>pdsaupd</code> work and also how to set the parameters
- * appropriately please take a look into the PARPACK manual.
- *
- * @author Denis Davydov, 2014.
- */
-template <typename VECTOR>
-class PArpackSolver : public Subscriptor
-{
-public:
- /**
- * Declare the type for container size.
- */
- typedef types::global_dof_index size_type;
-
- /**
- * An enum that lists the possible choices for which eigenvalues to
- * compute in the solve() function.
- *
- * A particular choice is limited based on symmetric or
- * non-symmetric matrix <code>A</code> considered.
- */
- enum WhichEigenvalues
- {
- algebraically_largest,
- algebraically_smallest,
- largest_magnitude,
- smallest_magnitude,
- largest_real_part,
- smallest_real_part,
- largest_imaginary_part,
- smallest_imaginary_part,
- both_ends
- };
-
- /**
- * Auxiliary class to represent <code>A-sigma*B</code> operator.
- */
- template <typename MATRIX>
- class Shift : public dealii::Subscriptor
- {
- public:
-
- /**
- * Constructor.
- */
- Shift (const MATRIX &A,
- const MATRIX &B,
- const double sigma)
- :
- A(A),
- B(B),
- sigma(sigma)
- {}
-
- /**
- * Apply <code>A-sigma * B</code>
- */
- void vmult (VECTOR &dst, const VECTOR &src) const
- {
- B.vmult(dst,src);
- dst *= (-sigma);
- A.vmult_add(dst,src);
- }
-
- /**
- * Apply <code>A^T-sigma * B^T</code>
- */
- void Tvmult (VECTOR &dst, const VECTOR &src) const
- {
- B.Tvmult(dst,src);
- dst *= (-sigma);
- A.Tvmult_add(dst,src);
- }
-
- private:
- const MATRIX &A;
- const MATRIX &B;
- const double sigma;
- };
-
- /**
- * Standardized data struct to pipe additional data to the solver,
- * should it be needed.
- */
- struct AdditionalData
- {
- 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 bool symmetric = false);
- };
-
- /**
- * Access to the object that controls convergence.
- */
- SolverControl &control () const;
-
- /**
- * Constructor.
- */
- PArpackSolver(SolverControl &control,
- const MPI_Comm &mpi_communicator,
- const AdditionalData &data = AdditionalData());
-
- /**
- * Initialise internal variables.
- */
- void reinit(const dealii::IndexSet &locally_owned_dofs );
-
- /**
- * Set desired shift value.
- */
- void set_shift(const double s );
-
- /**
- * Solve the generalized eigensprectrum problem $A x=\lambda B x$ by
- * calling the <code>pd(n/s)eupd</code> and <code>pd(n/s)aupd</code>
- * functions of PARPACK.
- */
- template <typename MATRIX1,
- typename MATRIX2, typename INVERSE>
- void solve(
- const MATRIX1 &A,
- const MATRIX2 &B,
- const INVERSE &inverse,
- std::vector<std::complex<double> > &eigenvalues,
- std::vector<VECTOR> &eigenvectors,
- const unsigned int n_eigenvalues);
-
- std::size_t memory_consumption() const;
-
-protected:
-
- /**
- * Reference to the object that controls convergence of the
- * iterative solver.
- */
- SolverControl &solver_control;
-
- /**
- * Store a copy of the flags for this particular solver.
- */
- const AdditionalData additional_data;
-
- // keep MPI communicator non-const as Arpack functions are not const either:
-
- /**
- * C++ MPI communicator.
- */
- MPI_Comm mpi_communicator;
-
- /**
- * Fortran MPI communicator.
- */
- MPI_Fint mpi_communicator_fortran;
-
- // C++98 guarantees that the elements of a vector are stored contiguously
-
- /**
- * Length of the work array workl.
- */
- int lworkl;
-
- /**
- * Double precision work array of length lworkl
- */
- std::vector<double> workl;
-
- /**
- * Double precision work array of length 3*N
- */
- std::vector<double> workd;
-
- /**
- * Number of local degrees of freedom.
- */
- int nloc;
-
- /**
- * Number of Arnoldi basis vectors specified in additional_data
- */
- int ncv;
-
-
- /**
- * The leading dimension of the array v
- */
- int ldv;
-
- /**
- * Double precision vector of size ldv by NCV. Will contains the
- * final set of Arnoldi basis vectors.
- */
- std::vector<double> v;
-
- /**
- * The initial residual vector, possibly from a previous run. On
- * output, it contains the final residual vector.
- */
- std::vector<double> resid;
-
- /**
- * The leading dimension of the array Z equal to nloc.
- */
- int ldz;
-
- /**
- * A vector of minimum size of nloc by NEV+1. Z contains the
- * B-orthonormal Ritz vectors of the eigensystem A*z = lambda*B*z
- * corresponding to the Ritz value approximations.
- */
- std::vector<double> z;
-
- /**
- * The size of the workev array.
- */
- int lworkev;
-
- /**
- * Double precision work array of dimension 3* NCV.
- */
- std::vector<double> workev;
-
- /**
- * A vector of dimension NCV.
- */
- std::vector<int> select;
-
- /**
- * Temporary vectors used between Arpack and deal.II
- */
- VECTOR src,dst,tmp;
-
- /**
- * Indices of local degrees of freedom.
- */
- std::vector< types::global_dof_index > local_indices;
-
- /**
- * The shift value to be applied during solution
- */
- double shift_value;
-
-private:
-
- /**
- * PArpackExcInfoPdnaupds.
- */
- DeclException2 (PArpackExcConvergedEigenvectors, int, int,
- << arg1 << "eigenpairs were requested, but only"
- << arg2 << " converged");
-
- DeclException2 (PArpackExcInvalidNumberofEigenvalues, int, int,
- << "Number of wanted eigenvalues " << arg1
- << " is larger that the size of the matrix " << arg2);
-
- DeclException2 (PArpackExcInvalidEigenvectorSize, int, int,
- << "Number of wanted eigenvalues " << arg1
- << " is larger that the size of eigenvectors " << arg2);
-
- DeclException2 (PArpackExcInvalidEigenvalueSize, int, int,
- << "Number of wanted eigenvalues " << arg1
- << " is larger that the size of eigenvalues " << arg2);
-
- DeclException2 (PArpackExcInvalidNumberofArnoldiVectors, int, int,
- << "Number of Arnoldi vectors " << arg1
- << " is larger that the size of the matrix " << arg2);
-
- DeclException2 (PArpackExcSmallNumberofArnoldiVectors, int, int,
- << "Number of Arnoldi vectors " << arg1
- << " is too small to obtain " << arg2
- << " eigenvalues");
-
- DeclException1 (PArpackExcIdo, int, << "This ido " << arg1
- << " is not supported. Check documentation of ARPACK");
-
- DeclException1 (PArpackExcMode, int, << "This mode " << arg1
- << " is not supported. Check documentation of ARPACK");
-
- DeclException1 (PArpackExcInfoPdnaupd, int,
- << "Error with Pdnaupd, info " << arg1
- << ". Check documentation of ARPACK");
-
- DeclException1 (PArpackExcInfoPdneupd, int,
- << "Error with Pdneupd, info " << arg1
- << ". Check documentation of ARPACK");
-
- DeclException1 (PArpackExcInfoMaxIt, int,
- << "Maximum number " << arg1
- << " of iterations reached.");
-
- DeclException1 (PArpackExcNoShifts, int,
- << "No shifts could be applied during implicit"
- << " Arnoldi update, try increasing the number of"
- << " Arnoldi vectors.");
-};
-
-template <typename VECTOR>
-std::size_t
-PArpackSolver<VECTOR>::memory_consumption() const
-{
- return MemoryConsumption::memory_consumption (double()) *
- (workl.size() +
- workd.size() +
- v.size() +
- resid.size() +
- z.size() +
- workev.size() ) +
- src.memory_consumption() +
- dst.memory_consumption() +
- tmp.memory_consumption() +
- MemoryConsumption::memory_consumption (types::global_dof_index()) * local_indices.size();
-}
-
-template <typename VECTOR>
-PArpackSolver<VECTOR>::AdditionalData::
-AdditionalData (const unsigned int number_of_arnoldi_vectors,
- const WhichEigenvalues eigenvalue_of_interest,
- const bool symmetric)
- :
- number_of_arnoldi_vectors(number_of_arnoldi_vectors),
- eigenvalue_of_interest(eigenvalue_of_interest),
- symmetric(symmetric)
-{}
-
-template <typename VECTOR>
-PArpackSolver<VECTOR>::PArpackSolver (SolverControl &control,
- const MPI_Comm &mpi_communicator,
- const AdditionalData &data)
- :
- solver_control (control),
- additional_data (data),
- mpi_communicator( mpi_communicator ),
- mpi_communicator_fortran ( MPI_Comm_c2f( mpi_communicator ) ),
- shift_value(0.0)
-
-{}
-
-template <typename VECTOR>
-void PArpackSolver<VECTOR>::set_shift(const double s )
-{
- shift_value = s;
-}
-
-template <typename VECTOR>
-void PArpackSolver<VECTOR>::reinit(const dealii::IndexSet &locally_owned_dofs)
-{
- // store local indices to write to vectors
- locally_owned_dofs.fill_index_vector(local_indices);
-
- // scalars
- nloc = locally_owned_dofs.n_elements ();
- ncv = additional_data.number_of_arnoldi_vectors;
-
- Assert (local_indices.size() == nloc, ExcInternalError() );
-
- // vectors
- ldv = nloc;
- v.resize (ldv*ncv, 0.0);
-
- // TODO: add optional input for resid
- resid.resize(nloc, 1.0);
-
- // work arrays for ARPACK
- workd.resize(3*nloc,0.0);
-
- lworkl = additional_data.symmetric ?
- ncv*ncv + 8*ncv
- :
- 3*ncv*ncv+6*ncv;
- workl.resize (lworkl, 0.);
-
- ldz = nloc;
- z.resize (ldz*ncv, 0.); // TODO we actually need only ldz*nev
-
- // WORKEV Double precision work array of dimension 3*NCV.
- lworkev = additional_data.symmetric ?
- 0 /*not used in symmetric case*/
- :
- 3*ncv;
- workev.resize (lworkev, 0.);
-
- select.resize (ncv, 0);
-
- // deal.II vectors:
- src.reinit (locally_owned_dofs,mpi_communicator);
- dst.reinit (locally_owned_dofs,mpi_communicator);
- tmp.reinit (locally_owned_dofs,mpi_communicator);
-
-}
-
-template <typename VECTOR>
-template <typename MATRIX1,typename MATRIX2, typename INVERSE>
-void PArpackSolver<VECTOR>::solve (
- const MATRIX1 &/*system_matrix*/,
- const MATRIX2 &mass_matrix,
- const INVERSE &inverse,
- std::vector<std::complex<double> > &eigenvalues,
- std::vector<VECTOR> &eigenvectors,
- const unsigned int n_eigenvalues)
-{
-
- Assert (n_eigenvalues <= eigenvectors.size(),
- PArpackExcInvalidEigenvectorSize(n_eigenvalues, eigenvectors.size()));
-
- Assert (n_eigenvalues <= eigenvalues.size(),
- PArpackExcInvalidEigenvalueSize(n_eigenvalues, eigenvalues.size()));
-
-
- Assert (n_eigenvalues < mass_matrix.m(),
- PArpackExcInvalidNumberofEigenvalues(n_eigenvalues, mass_matrix.m()));
-
- Assert (additional_data.number_of_arnoldi_vectors < mass_matrix.m(),
- PArpackExcInvalidNumberofArnoldiVectors(
- additional_data.number_of_arnoldi_vectors, mass_matrix.m()));
-
- Assert (additional_data.number_of_arnoldi_vectors > 2*n_eigenvalues+1,
- PArpackExcSmallNumberofArnoldiVectors(
- additional_data.number_of_arnoldi_vectors, n_eigenvalues));
- // ARPACK mode for dnaupd, here only mode 3
- int mode = 3;
-
- // reverse communication parameter
- // must be zero on the first call to pdnaupd
- int ido = 0;
-
- /**
- * 'G' generalized eigenvalue problem
- * 'I' standard eigenvalue problem
- */
- char bmat[2] = "G";
-
- /** Specify the eigenvalues of interest,
- * possible parameters
- * "LA" algebraically largest
- * "SA" algebraically smallest
- * "LM" largest magnitude
- * "SM" smallest magnitude
- * "LR" largest real part
- * "SR" smallest real part
- * "LI" largest imaginary part
- * "SI" smallest imaginary part
- * "BE" both ends of spectrum simultaneous
- */
- char which[3];
- switch (additional_data.eigenvalue_of_interest)
- {
- case algebraically_largest:
- std::strcpy (which, "LA");
- break;
- case algebraically_smallest:
- std::strcpy (which, "SA");
- break;
- case largest_magnitude:
- std::strcpy (which, "LM");
- break;
- case smallest_magnitude:
- std::strcpy (which, "SM");
- break;
- case largest_real_part:
- std::strcpy (which, "LR");
- break;
- case smallest_real_part:
- std::strcpy (which, "SR");
- break;
- case largest_imaginary_part:
- std::strcpy (which, "LI");
- break;
- case smallest_imaginary_part:
- std::strcpy (which, "SI");
- break;
- case both_ends:
- std::strcpy (which, "BE");
- break;
- }
-
- // tolerance for ARPACK
- double tol = control().tolerance();
-
- //information to the routines
- std::vector<int> iparam (11, 0);
-
- iparam[0] = 1;
- // shift strategy: exact shifts with respect to the current Hessenberg matrix H.
-
- // maximum number of iterations
- iparam[2] = control().max_steps();
-
- /** Sets the mode of dsaupd.
- * 1 is exact shifting,
- * 2 is user-supplied shifts,
- * 3 is shift-invert mode,
- * 4 is buckling mode,
- * 5 is Cayley mode.
- */
-
- iparam[6] = mode;
- std::vector<int> ipntr (14, 0);
-
- //information out of the iteration
- // If INFO .EQ. 0, a random initial residual vector is used.
- // If INFO .NE. 0, RESID contains the initial residual vector,
- // possibly from a previous run.
- // Typical choices in this situation might be to use the final value
- // of the starting vector from the previous eigenvalue calculation
- int info = 1;
-
- // Number of eigenvalues of OP to be computed. 0 < NEV < N.
- int nev = n_eigenvalues;
- int n_inside_arpack = nloc;
-
- while (ido != 99)
- {
- // call of ARPACK pdnaupd routine
- if (additional_data.symmetric)
- pdsaupd_(&mpi_communicator_fortran,&ido, bmat, &n_inside_arpack, which, &nev, &tol,
- &resid[0], &ncv, &v[0], &ldv, &iparam[0], &ipntr[0],
- &workd[0], &workl[0], &lworkl, &info);
- else
- pdnaupd_(&mpi_communicator_fortran,&ido, bmat, &n_inside_arpack, which, &nev, &tol,
- &resid[0], &ncv, &v[0], &ldv, &iparam[0], &ipntr[0],
- &workd[0], &workl[0], &lworkl, &info);
-
- if (ido == 99)
- break;
-
- switch (mode)
- {
-// OP = (inv[K - sigma*M])*M
- case 3:
- {
- switch (ido)
- {
-// compute Y = OP * X where
-// IPNTR(1) is the pointer into WORKD for X,
-// IPNTR(2) is the pointer into WORKD for Y.
- case -1:
- {
- const int shift_x = ipntr[0]-1;
- const int shift_y = ipntr[1]-1;
- Assert (shift_x>=0, dealii::ExcInternalError() );
- Assert (shift_x+nloc <= workd.size(), dealii::ExcInternalError() );
- Assert (shift_y>=0, dealii::ExcInternalError() );
- Assert (shift_y+nloc <= workd.size(), dealii::ExcInternalError() );
-
- src = 0.0;
- src.add (nloc,
- &local_indices[0],
- &workd[0]+shift_x );
- src.compress (VectorOperation::add);
-
- // multiplication with mass matrix M
- mass_matrix.vmult(tmp, src);
- // solving linear system
- inverse.vmult(dst,tmp);
-
- // store the result
- dst.extract_subvector_to (local_indices.begin(),
- local_indices.end(),
- &workd[0]+shift_y );
- }
- break;
-
-// compute Y = OP * X where
-// IPNTR(1) is the pointer into WORKD for X,
-// IPNTR(2) is the pointer into WORKD for Y.
-// In mode 3,4 and 5, the vector B * X is already
-// available in WORKD(ipntr(3)). It does not
-// need to be recomputed in forming OP * X.
- case 1:
- {
- const int shift_x = ipntr[0]-1;
- const int shift_y = ipntr[1]-1;
- const int shift_b_x = ipntr[2]-1;
-
- Assert (shift_x>=0, dealii::ExcInternalError() );
- Assert (shift_x+nloc <= workd.size(), dealii::ExcInternalError() );
- Assert (shift_y>=0, dealii::ExcInternalError() );
- Assert (shift_y+nloc <= workd.size(), dealii::ExcInternalError() );
- Assert (shift_b_x>=0, dealii::ExcInternalError() );
- Assert (shift_b_x+nloc <= workd.size(), dealii::ExcInternalError() );
- Assert (shift_y>=0, dealii::ExcInternalError() );
- Assert (shift_y+nloc <= workd.size(), dealii::ExcInternalError() );
-
- src = 0.0; // B*X
- src.add (nloc,
- &local_indices[0],
- &workd[0]+shift_b_x );
-
- tmp = 0.0; // X
- tmp.add (nloc,
- &local_indices[0],
- &workd[0]+shift_x);
-
- src.compress (VectorOperation::add);
- tmp.compress (VectorOperation::add);
-
- // solving linear system
- inverse.vmult(dst,src);
-
- // store the result
- dst.extract_subvector_to (local_indices.begin(),
- local_indices.end(),
- &workd[0]+shift_y );
-
- }
- break;
-
-// compute Y = B * X where
-// IPNTR(1) is the pointer into WORKD for X,
-// IPNTR(2) is the pointer into WORKD for Y.
- case 2:
- {
-
- const int shift_x = ipntr[0]-1;
- const int shift_y = ipntr[1]-1;
- Assert (shift_x>=0, dealii::ExcInternalError() );
- Assert (shift_x+nloc <= workd.size(), dealii::ExcInternalError() );
- Assert (shift_y>=0, dealii::ExcInternalError() );
- Assert (shift_y+nloc <= workd.size(), dealii::ExcInternalError() );
-
- src = 0.0;
- src.add (nloc,
- &local_indices[0],
- &workd[0]+shift_x );
- src.compress (VectorOperation::add);
-
- // Multiplication with mass matrix M
- mass_matrix.vmult(dst, src);
-
- // store the result
- dst.extract_subvector_to (local_indices.begin(),
- local_indices.end(),
- &workd[0]+shift_y);
-
- }
- break;
-
- default:
- Assert (false, PArpackExcIdo(ido));
- break;
- }
- }
- break;
- default:
- Assert (false, PArpackExcMode(mode));
- break;
- }
- }
-
- if (info<0)
- {
- Assert (false, PArpackExcInfoPdnaupd(info));
- }
- else
- {
- /** 1 - compute eigenvectors,
- * 0 - only eigenvalues
- */
- int rvec = 1;
-
- // which eigenvectors
- char howmany[4] = "All";
-
- 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.);
-
- // call of ARPACK pdneupd routine
- if (additional_data.symmetric)
- pdseupd_(&mpi_communicator_fortran, &rvec, howmany, &select[0], &eigenvalues_real[0],
- &z[0], &ldz, &sigmar,
- bmat, &n_inside_arpack, which, &nev, &tol,
- &resid[0], &ncv, &v[0], &ldv,
- &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,
- &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);
-
- if (info == 1)
- {
- Assert (false, PArpackExcInfoMaxIt(control().max_steps()));
- }
- else if (info == 3)
- {
- Assert (false, PArpackExcNoShifts(1));
- }
- else if (info!=0)
- {
- Assert (false, PArpackExcInfoPdneupd(info));
- }
-
- for (size_type i=0; i<n_eigenvalues; ++i)
- {
- eigenvectors[i] = 0.0;
- Assert (i*nloc + nloc <= v.size(), dealii::ExcInternalError() );
-
- eigenvectors[i].add (nloc,
- &local_indices[0],
- &v[i*nloc] );
- eigenvectors[i].compress (VectorOperation::add);
- }
-
- for (size_type i=0; i<n_eigenvalues; ++i)
- eigenvalues[i] = std::complex<double> (eigenvalues_real[i],
- eigenvalues_im[i]);
- }
-
- Assert (iparam[4] == n_eigenvalues,
- PArpackExcConvergedEigenvectors(iparam[4], n_eigenvalues));
-
- // both PDNAUPD and PDSAUPD compute eigenpairs of inv[A - sigma*M]*M
- // with respect to a semi-inner product defined by M.
-
- // resid likely contains residual with respect to M-norm.
- {
-
- tmp = 0.0;
- tmp.add (nloc,
- &local_indices[0],
- &resid[0]);
- solver_control.check ( iparam[2], tmp.l2_norm() );
- }
-
-
-}
-
-template <typename VECTOR>
-SolverControl &PArpackSolver<VECTOR>::control () const
-{
- return solver_control;
-}
-
-DEAL_II_NAMESPACE_CLOSE
-
-
-#endif
-#endif
+++ /dev/null
-#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 <class DH>
-std::vector<dealii::IndexSet>
-locally_owned_dofs_per_subdomain (const DH &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)
- {
- solver.set_symmetric_mode(true);
- }
-
- void vmult ( dealii::PETScWrappers::MPI::Vector & dst,
- const dealii::PETScWrappers::MPI::Vector & src) const
- {
- solver.solve(matrix, dst, src);
- }
-
-
-private:
- mutable dealii::PETScWrappers::SparseDirectMUMPS solver;
- const dealii::PETScWrappers::MatrixBase & matrix;
-
-};
-
-void test ()
-{
- const unsigned int global_mesh_refinement_steps = 5;
- const unsigned int number_of_eigenvalues = 5;
-
- 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<PetscScalar> eigenvalues;
- dealii::PETScWrappers::MPI::SparseMatrix stiffness_matrix, mass_matrix;
-
- dealii::GridGenerator::hyper_cube (triangulation, -1, 1);
- triangulation.refine_global (global_mesh_refinement_steps);
-
- // METIS:
- //dealii::GridTools::partition_triangulation (n_mpi_processes, triangulation);
-
- // partition by hand
- {
- 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);
- 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);
-
- eigenfunctions.resize (5);
- for (unsigned int i=0; i<eigenfunctions.size (); ++i)
- eigenfunctions[i].reinit (locally_owned_dofs, mpi_communicator);//without ghost dofs
-
- eigenvalues.resize (eigenfunctions.size ());
-
-
- // 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)
- {
- 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)
- 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)
- ) * 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-9,/*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_magnitude,
- true);
-
- dealii::PArpackSolver<dealii::PETScWrappers::MPI::Vector> eigensolver (solver_control,
- mpi_communicator,
- additional_data);
- eigensolver.reinit(locally_owned_dofs);
- eigensolver.solve (stiffness_matrix,
- mass_matrix,
- inverse,
- lambda,
- eigenfunctions,
- eigenvalues.size());
-
- for (unsigned int i = 0; i < lambda.size(); i++)
- eigenvalues[i] = lambda[i].real();
-
- for (unsigned int i=0; i < eigenvalues.size(); i++)
- dealii::deallog << eigenvalues[i] << 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.depth_console(0);
- 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
-#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 <class DH>
-std::vector<IndexSet>
-locally_owned_dofs_per_subdomain (const DH &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 = 5;
-
- 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<double> eigenvalues;
- TrilinosWrappers::SparseMatrix stiffness_matrix, mass_matrix;
-
- GridGenerator::hyper_cube (triangulation, -1, 1);
- triangulation.refine_global (global_mesh_refinement_steps);
-
- //METIS:
- //GridTools::partition_triangulation (n_mpi_processes, triangulation);
-
- // partition by hand
- {
- 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);
-
- eigenfunctions.resize (5);
- for (unsigned int i=0; i<eigenfunctions.size (); ++i)
- eigenfunctions[i].reinit (locally_owned_dofs, mpi_communicator);//without ghost dofs
-
- eigenvalues.resize (eigenfunctions.size ());
-
-
- // 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)
- {
- 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)
- 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)
- ) * 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] = PetscScalar();
-
- 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("cg");
- 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_magnitude,
- true);
-
- PArpackSolver<TrilinosWrappers::MPI::Vector> eigensolver (solver_control,
- mpi_communicator,
- additional_data);
- eigensolver.reinit(locally_owned_dofs);
- eigensolver.set_shift(shift);
-
- // avoid output of iterative solver:
- const unsigned int previous_depth = deallog.depth_file(0);
- eigensolver.solve (stiffness_matrix,
- mass_matrix,
- shift_and_invert,
- lambda,
- eigenfunctions,
- eigenvalues.size());
- deallog.depth_file(previous_depth);
-
- for (unsigned int i = 0; i < lambda.size(); i++)
- eigenvalues[i] = lambda[i].real();
-
- for (unsigned int i=0; i < eigenvalues.size(); i++)
- deallog << eigenvalues[i] << 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.depth_console(0);
- 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;
- };
-}