From: Matthias Maier Date: Sun, 1 Feb 2015 09:40:51 +0000 (+0100) Subject: Revert parpack merge X-Git-Tag: v8.3.0-rc1~509^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=6c379894849fc783a617df783f46cbfe08090d27;p=dealii.git Revert parpack merge This commit reverts: a75f23180d6ec4b24fcf7b714d5b763dcb017efb 616a403d8f799a21b7d306ecf38f3c73d0857524 9367c1fa0371a260438dd93db74882a7e7ebe542 6ed3174c7ecb1fb61d263e7a768deb6325bca7d3 9a03bc439277d85502accef56d8f9048558ae925 --- diff --git a/cmake/configure/configure_arpack.cmake b/cmake/configure/configure_arpack.cmake index 5ac73fed24..7cffb145bd 100644 --- a/cmake/configure/configure_arpack.cmake +++ b/cmake/configure/configure_arpack.cmake @@ -20,4 +20,3 @@ SET(FEATURE_ARPACK_DEPENDS LAPACK) CONFIGURE_FEATURE(ARPACK) -SET(DEAL_II_ARPACK_WITH_PARPACK ARPACK_WITH_PARPACK) diff --git a/cmake/modules/FindARPACK.cmake b/cmake/modules/FindARPACK.cmake index b15a4c722a..0e185cbc8b 100644 --- a/cmake/modules/FindARPACK.cmake +++ b/cmake/modules/FindARPACK.cmake @@ -1,6 +1,6 @@ ## --------------------------------------------------------------------- ## -## Copyright (C) 2012 - 2015 by the deal.II authors +## Copyright (C) 2012 - 2014 by the deal.II authors ## ## This file is part of the deal.II library. ## @@ -20,7 +20,10 @@ # # ARPACK_LIBRARIES # ARPACK_LINKER_FLAGS -# ARPACK_WITH_PARPACK +# + +# +# TODO: ARPACK and mpi... # SET(ARPACK_DIR "" CACHE PATH "An optional hint to an ARPACK installation") @@ -32,37 +35,8 @@ DEAL_II_FIND_LIBRARY(ARPACK_LIBRARY PATH_SUFFIXES lib${LIB_SUFFIX} lib64 lib ) -IF(DEAL_II_WITH_MPI) - # - # Sanity check: Only search the parpack library in the same directory as - # the arpack library... - # - GET_FILENAME_COMPONENT(_path "${ARPACK_LIBRARY}" PATH) - DEAL_II_FIND_LIBRARY(PARPACK_LIBRARY - NAMES parpack - HINTS ${_path} - NO_DEFAULT_PATH - NO_CMAKE_ENVIRONMENT_PATH - NO_CMAKE_PATH - NO_SYSTEM_ENVIRONMENT_PATH - NO_CMAKE_SYSTEM_PATH - NO_CMAKE_FIND_ROOT_PATH - ) -ELSE() - SET(PARPACK_LIBRARY "PARPACK_LIBRARY-NOTFOUND") -ENDIF() - -IF(NOT PARPACK_LIBRARY MATCHES "-NOTFOUND") - SET(ARPACK_WITH_PARPACK TRUE) -ELSE() - SET(ARPACK_WITH_PARPACK FALSE) -ENDIF() - DEAL_II_PACKAGE_HANDLE(ARPACK - LIBRARIES - OPTIONAL PARPACK_LIBRARY - REQUIRED ARPACK_LIBRARY LAPACK_LIBRARIES - OPTIONAL MPI_C_LIBRARIES + LIBRARIES REQUIRED ARPACK_LIBRARY LAPACK_LIBRARIES LINKER_FLAGS OPTIONAL LAPACK_LINKER_FLAGS - CLEAR ARPACK_LIBRARY PARPACK_LIBRARY + CLEAR ARPACK_LIBRARY ) diff --git a/doc/news/changes.h b/doc/news/changes.h index e0ac7c88c2..b2eb4ad45b 100644 --- a/doc/news/changes.h +++ b/doc/news/changes.h @@ -280,12 +280,6 @@ inconvenience this causes.

Specific improvements

    -
  1. New: deal.II now also provides an interface to the parallel ARPACK - (PARPACK) solver with the PArpackSolver class. -
    - (Denis Davydov, 2015/01/30) -
    2. -
  2. Bugfix: deal.II now correctly links against librt in case of bundled boost being used.
    diff --git a/include/deal.II/base/config.h.in b/include/deal.II/base/config.h.in index dc5f57ccfa..fda7109f58 100644 --- a/include/deal.II/base/config.h.in +++ b/include/deal.II/base/config.h.in @@ -62,8 +62,6 @@ #cmakedefine DEAL_II_WITH_UMFPACK #cmakedefine DEAL_II_WITH_ZLIB -#cmakedefine DEAL_II_ARPACK_WITH_PARPACK - /*********************************************************************** * Compiler bugs: diff --git a/include/deal.II/lac/parpack_solver.h b/include/deal.II/lac/parpack_solver.h deleted file mode 100644 index 10861b29bd..0000000000 --- a/include/deal.II/lac/parpack_solver.h +++ /dev/null @@ -1,870 +0,0 @@ -// --------------------------------------------------------------------- -// -// 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 -#include -#include -#include -#include -#include - -#include - - -#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 pdneupd, - * pdseupd, pdnaupd, pdsaupd 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::AdditionalData - additional_data(num_arnoldi_vectors, - dealii::PArpackSolver::largest_magnitude, - true); - - PArpackSolver 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 size_of_spectrum tells PARPACK the number of - * eigenvector/eigenvalue pairs to solve for. Here, - * lambda is a vector that will contain the eigenvalues - * computed, x a vector of objects of type V - * that will contain the eigenvectors computed. OP is an - * inverse operation for the matrix A - sigma * B, where - * sigma 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 A are supported. - * - * For further information on how the PARPACK routines - * pdneupd, pdseupd, pdnaupd, - * pdsaupd work and also how to set the parameters - * appropriately please take a look into the PARPACK manual. - * - * @author Denis Davydov, 2014. - */ -template -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 A 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 A-sigma*B operator. - */ - template - class Shift : public dealii::Subscriptor - { - public: - - /** - * Constructor. - */ - Shift (const MATRIX &A, - const MATRIX &B, - const double sigma) - : - A(A), - B(B), - sigma(sigma) - {} - - /** - * Apply A-sigma * B - */ - void vmult (VECTOR &dst, const VECTOR &src) const - { - B.vmult(dst,src); - dst *= (-sigma); - A.vmult_add(dst,src); - } - - /** - * Apply A^T-sigma * B^T - */ - 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 pd(n/s)eupd and pd(n/s)aupd - * functions of PARPACK. - */ - template - void solve( - const MATRIX1 &A, - const MATRIX2 &B, - const INVERSE &inverse, - std::vector > &eigenvalues, - std::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 workl; - - /** - * Double precision work array of length 3*N - */ - std::vector 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 v; - - /** - * The initial residual vector, possibly from a previous run. On - * output, it contains the final residual vector. - */ - std::vector 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 z; - - /** - * The size of the workev array. - */ - int lworkev; - - /** - * Double precision work array of dimension 3* NCV. - */ - std::vector workev; - - /** - * A vector of dimension NCV. - */ - std::vector 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 -std::size_t -PArpackSolver::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 -PArpackSolver::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 -PArpackSolver::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 -void PArpackSolver::set_shift(const double s ) -{ - shift_value = s; -} - -template -void PArpackSolver::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 -template -void PArpackSolver::solve ( - const MATRIX1 &/*system_matrix*/, - const MATRIX2 &mass_matrix, - const INVERSE &inverse, - std::vector > &eigenvalues, - std::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 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 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 eigenvalues_real (n_eigenvalues, 0.); - std::vector 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 (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 -SolverControl &PArpackSolver::control () const -{ - return solver_control; -} - -DEAL_II_NAMESPACE_CLOSE - - -#endif -#endif diff --git a/tests/arpack/step-36_arpack.cc b/tests/arpack/step-36_ar.cc similarity index 100% rename from tests/arpack/step-36_arpack.cc rename to tests/arpack/step-36_ar.cc diff --git a/tests/arpack/step-36_arpack.output b/tests/arpack/step-36_ar.output similarity index 100% rename from tests/arpack/step-36_arpack.output rename to tests/arpack/step-36_ar.output diff --git a/tests/arpack/step-36_parpack.cc b/tests/arpack/step-36_parpack.cc deleted file mode 100644 index 7624e9cb38..0000000000 --- a/tests/arpack/step-36_parpack.cc +++ /dev/null @@ -1,342 +0,0 @@ -#include -#include -#include - -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include - -#include - - -#include -#include -#include -#include - -#include - -#include -#include - -const unsigned int dim = 2;//run in 2d to save time - -const double eps = 1e-10; - -template -std::vector -locally_owned_dofs_per_subdomain (const 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 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 triangulation; - dealii::DoFHandler dof_handler(triangulation); - dealii::FE_Q fe(1); - dealii::ConstraintMatrix constraints; - dealii::IndexSet locally_owned_dofs; - dealii::IndexSet locally_relevant_dofs; - - std::vector eigenfunctions; - std::vector 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::active_cell_iterator - cell = triangulation.begin_active(), - endc = triangulation.end(); - for (; cell!=endc; ++cell) - { - const Point ¢er = cell->center(); - const double x = center[0]; - - const unsigned int id = std::floor ( (x-x0)/dL); - cell->set_subdomain_id (id); - } - } - - dof_handler.distribute_dofs (fe); - dealii::DoFRenumbering::subdomain_wise (dof_handler); - std::vector locally_owned_dofs_per_processor - = locally_owned_dofs_per_subdomain (dof_handler); - locally_owned_dofs = locally_owned_dofs_per_processor[this_mpi_process]; - locally_relevant_dofs.clear(); - dealii::DoFTools::extract_locally_relevant_dofs (dof_handler, - locally_relevant_dofs); - - constraints.clear(); - constraints.reinit (locally_relevant_dofs); - dealii::DoFTools::make_hanging_node_constraints (dof_handler, constraints); - dealii::VectorTools::interpolate_boundary_values (dof_handler, - 0, - dealii::ZeroFunction (), - constraints); - constraints.close (); - - dealii::CompressedSimpleSparsityPattern csp (locally_relevant_dofs); - // Fill in ignoring all cells that are not locally owned - dealii::DoFTools::make_sparsity_pattern (dof_handler, csp, - constraints, - /* keep constrained dofs */ true); - std::vector n_locally_owned_dofs(n_mpi_processes); - for (unsigned int i = 0; i < n_mpi_processes; i++) - n_locally_owned_dofs[i] = locally_owned_dofs_per_processor[i].n_elements(); - - dealii::SparsityTools::distribute_sparsity_pattern - (csp, - n_locally_owned_dofs, - mpi_communicator, - locally_relevant_dofs); - - // Initialise the stiffness and mass matrices - stiffness_matrix.reinit (locally_owned_dofs, - locally_owned_dofs, - csp, - mpi_communicator); - - mass_matrix.reinit (locally_owned_dofs, - locally_owned_dofs, - csp, - mpi_communicator); - - eigenfunctions.resize (5); - for (unsigned int i=0; i quadrature_formula(2); - dealii::FEValues fe_values (fe, quadrature_formula, - dealii::update_values | - dealii::update_gradients | - dealii::update_quadrature_points | - dealii::update_JxW_values); - - const unsigned int dofs_per_cell = fe.dofs_per_cell; - const unsigned int n_q_points = quadrature_formula.size(); - - dealii::FullMatrix cell_stiffness_matrix (dofs_per_cell, dofs_per_cell); - dealii::FullMatrix cell_mass_matrix (dofs_per_cell, dofs_per_cell); - - std::vector local_dof_indices (dofs_per_cell); - - typename dealii::DoFHandler::active_cell_iterator - cell = dof_handler.begin_active (), - endc = dof_handler.end (); - for (; cell!=endc; ++cell) - { - fe_values.reinit (cell); - cell_stiffness_matrix = 0; - cell_mass_matrix = 0; - - for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); - - constraints - .distribute_local_to_global (cell_stiffness_matrix, - local_dof_indices, - stiffness_matrix); - constraints - .distribute_local_to_global (cell_mass_matrix, - local_dof_indices, - mass_matrix); - } - - stiffness_matrix.compress (dealii::VectorOperation::add); - mass_matrix.compress (dealii::VectorOperation::add); - - // test Arpack - { - std::vector > lambda(eigenfunctions.size()); - - for (unsigned int i=0; i < eigenvalues.size(); i++) - eigenfunctions[i] = PetscScalar(); - - dealii::SolverControl solver_control (dof_handler.n_dofs(), 1e-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::AdditionalData - additional_data(num_arnoldi_vectors, - dealii::PArpackSolver::largest_magnitude, - true); - - dealii::PArpackSolver 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"< -#include -#include - -#include -#include -#include -#include -#include -#include -#include -#include - -#include - -#include -#include -#include - - -#include -#include -#include -#include - -#include - -#include -#include - -const unsigned int dim = 2;//run in 2d to save time - -using namespace dealii; - -const double eps = 1e-10; - -template -std::vector -locally_owned_dofs_per_subdomain (const 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 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 triangulation; - DoFHandler dof_handler(triangulation); - FE_Q fe(1); - ConstraintMatrix constraints; - IndexSet locally_owned_dofs; - IndexSet locally_relevant_dofs; - - std::vector eigenfunctions; - std::vector 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::active_cell_iterator - cell = triangulation.begin_active(), - endc = triangulation.end(); - for (; cell!=endc; ++cell) - { - const Point ¢er = cell->center(); - const double x = center[0]; - - const unsigned int id = std::floor ( (x-x0)/dL); - cell->set_subdomain_id (id); - } - } - - - dof_handler.distribute_dofs (fe); - DoFRenumbering::subdomain_wise (dof_handler); - std::vector locally_owned_dofs_per_processor - = locally_owned_dofs_per_subdomain (dof_handler); - locally_owned_dofs = locally_owned_dofs_per_processor[this_mpi_process]; - locally_relevant_dofs.clear(); - DoFTools::extract_locally_relevant_dofs (dof_handler, - locally_relevant_dofs); - - constraints.clear(); - constraints.reinit (locally_relevant_dofs); - DoFTools::make_hanging_node_constraints (dof_handler, constraints); - VectorTools::interpolate_boundary_values (dof_handler, - 0, - ZeroFunction (), - constraints); - constraints.close (); - - CompressedSimpleSparsityPattern csp (locally_relevant_dofs); - // Fill in ignoring all cells that are not locally owned - DoFTools::make_sparsity_pattern (dof_handler, csp, - constraints, - /* keep constrained dofs */ true); - std::vector n_locally_owned_dofs(n_mpi_processes); - for (unsigned int i = 0; i < n_mpi_processes; i++) - n_locally_owned_dofs[i] = locally_owned_dofs_per_processor[i].n_elements(); - - SparsityTools::distribute_sparsity_pattern - (csp, - n_locally_owned_dofs, - mpi_communicator, - locally_relevant_dofs); - - // Initialise the stiffness and mass matrices - stiffness_matrix.reinit (locally_owned_dofs, - locally_owned_dofs, - csp, - mpi_communicator); - - mass_matrix.reinit (locally_owned_dofs, - locally_owned_dofs, - csp, - mpi_communicator); - - eigenfunctions.resize (5); - for (unsigned int i=0; i quadrature_formula(2); - FEValues fe_values (fe, quadrature_formula, - update_values | - update_gradients | - update_quadrature_points | - update_JxW_values); - - const unsigned int dofs_per_cell = fe.dofs_per_cell; - const unsigned int n_q_points = quadrature_formula.size(); - - FullMatrix cell_stiffness_matrix (dofs_per_cell, dofs_per_cell); - FullMatrix cell_mass_matrix (dofs_per_cell, dofs_per_cell); - - std::vector local_dof_indices (dofs_per_cell); - - typename DoFHandler::active_cell_iterator - cell = dof_handler.begin_active (), - endc = dof_handler.end (); - for (; cell!=endc; ++cell) - { - fe_values.reinit (cell); - cell_stiffness_matrix = 0; - cell_mass_matrix = 0; - - for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); - - constraints - .distribute_local_to_global (cell_stiffness_matrix, - local_dof_indices, - stiffness_matrix); - constraints - .distribute_local_to_global (cell_mass_matrix, - local_dof_indices, - mass_matrix); - } - - stiffness_matrix.compress (VectorOperation::add); - mass_matrix.compress (VectorOperation::add); - - // test Arpack - { - const double shift = 4.0; - std::vector > lambda(eigenfunctions.size()); - - for (unsigned int i=0; i < eigenvalues.size(); i++) - eigenfunctions[i] = 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::Shift shifted_matrix(stiffness_matrix,mass_matrix,shift); - TrilinosWrappers::PreconditionIdentity preconditioner; - IterativeInverse shift_and_invert; - shift_and_invert.initialize(shifted_matrix,preconditioner); - shift_and_invert.solver.select("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::AdditionalData - additional_data(num_arnoldi_vectors, - PArpackSolver::largest_magnitude, - true); - - PArpackSolver 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"<