]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Revert parpack merge 483/head
authorMatthias Maier <tamiko@43-1.org>
Sun, 1 Feb 2015 09:40:51 +0000 (10:40 +0100)
committerMatthias Maier <tamiko@43-1.org>
Sun, 1 Feb 2015 09:43:01 +0000 (10:43 +0100)
This commit reverts:
  a75f23180d6ec4b24fcf7b714d5b763dcb017efb
  616a403d8f799a21b7d306ecf38f3c73d0857524
  9367c1fa0371a260438dd93db74882a7e7ebe542
  6ed3174c7ecb1fb61d263e7a768deb6325bca7d3
  9a03bc439277d85502accef56d8f9048558ae925

cmake/configure/configure_arpack.cmake
cmake/modules/FindARPACK.cmake
doc/news/changes.h
include/deal.II/base/config.h.in
include/deal.II/lac/parpack_solver.h [deleted file]
tests/arpack/step-36_ar.cc [moved from tests/arpack/step-36_arpack.cc with 100% similarity]
tests/arpack/step-36_ar.output [moved from tests/arpack/step-36_arpack.output with 100% similarity]
tests/arpack/step-36_parpack.cc [deleted file]
tests/arpack/step-36_parpack.with_mpi=true.mpirun=3.output [deleted file]
tests/arpack/step-36_parpack_trilinos.cc [deleted file]
tests/arpack/step-36_parpack_trilinos.with_mpi=true.with_trilinos=true.mpirun=3.output [deleted file]

index 5ac73fed24ed042c9d2639160957b33a66f61e5c..7cffb145bd0d0064532c59cf2312338b95ddd62b 100644 (file)
@@ -20,4 +20,3 @@
 SET(FEATURE_ARPACK_DEPENDS LAPACK)
 
 CONFIGURE_FEATURE(ARPACK)
-SET(DEAL_II_ARPACK_WITH_PARPACK ARPACK_WITH_PARPACK)
index b15a4c722afbf6d3884f0ed4b8a96ceff84a5650..0e185cbc8bcd82ca934c8d6fd26fe6d6414cb98e 100644 (file)
@@ -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.
 ##
 #
 #   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
   )
index e0ac7c88c2d61d341bd2cc24a70c8a6498a35cad..b2eb4ad45b7c9397f7f9dd681e85660b8abf77a8 100644 (file)
@@ -280,12 +280,6 @@ inconvenience this causes.
 <h3>Specific improvements</h3>
 
 <ol>
-  <li> New: deal.II now also provides an interface to the parallel ARPACK
-  (PARPACK) solver with the <code>PArpackSolver</code> class.
-  <br>
-  (Denis Davydov, 2015/01/30)
-  </li>
-
   <li> Bugfix: deal.II now correctly links against librt in case of bundled
   boost being used.
   <br>
index dc5f57ccfa1febc5a9f44d6779f86bd2ac92874b..fda7109f589ca120403d1a6b6689579ec8486946 100644 (file)
@@ -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 (file)
index 10861b2..0000000
+++ /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 <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
diff --git a/tests/arpack/step-36_parpack.cc b/tests/arpack/step-36_parpack.cc
deleted file mode 100644 (file)
index 7624e9c..0000000
+++ /dev/null
@@ -1,342 +0,0 @@
-#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> &center = 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;
-  };
-}
diff --git a/tests/arpack/step-36_parpack.with_mpi=true.mpirun=3.output b/tests/arpack/step-36_parpack.with_mpi=true.mpirun=3.output
deleted file mode 100644 (file)
index f0a5f87..0000000
+++ /dev/null
@@ -1,6 +0,0 @@
-DEAL::4.93877
-DEAL::12.3707
-DEAL::12.3707
-DEAL::19.8027
-DEAL::24.8370
-DEAL::Ok
diff --git a/tests/arpack/step-36_parpack_trilinos.cc b/tests/arpack/step-36_parpack_trilinos.cc
deleted file mode 100644 (file)
index 03a17a8..0000000
+++ /dev/null
@@ -1,339 +0,0 @@
-#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> &center = 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;
-  };
-}
diff --git a/tests/arpack/step-36_parpack_trilinos.with_mpi=true.with_trilinos=true.mpirun=3.output b/tests/arpack/step-36_parpack_trilinos.with_mpi=true.with_trilinos=true.mpirun=3.output
deleted file mode 100644 (file)
index f0a5f87..0000000
+++ /dev/null
@@ -1,6 +0,0 @@
-DEAL::4.93877
-DEAL::12.3707
-DEAL::12.3707
-DEAL::19.8027
-DEAL::24.8370
-DEAL::Ok

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.