]> https://gitweb.dealii.org/ - dealii.git/commitdiff
add Utilities::LinearAlgebra::Lanczos_largest_eigenvalue() 4679/head
authorDenis Davydov <davydden@gmail.com>
Tue, 1 Aug 2017 20:33:27 +0000 (22:33 +0200)
committerDenis Davydov <davydden@gmail.com>
Sat, 5 Aug 2017 17:01:38 +0000 (19:01 +0200)
and Utilities::LinearAlgebra::Chebyshev_filter()

doc/doxygen/images/chebyshev8.png [new file with mode: 0644]
doc/news/changes/minor/20170801DenisDavydov-c [new file with mode: 0644]
doc/news/changes/minor/20170802DenisDavydov [new file with mode: 0644]
include/deal.II/lac/diagonal_matrix.h
include/deal.II/lac/utilities.h [new file with mode: 0644]
tests/lac/utilities_01.cc [new file with mode: 0644]
tests/lac/utilities_01.output [new file with mode: 0644]
tests/lac/utilities_02.cc [new file with mode: 0644]
tests/lac/utilities_02.with_mpi=true.with_p4est=true.mpirun=3.output [new file with mode: 0644]
tests/lac/utilities_03.cc [new file with mode: 0644]
tests/lac/utilities_03.output [new file with mode: 0644]

diff --git a/doc/doxygen/images/chebyshev8.png b/doc/doxygen/images/chebyshev8.png
new file mode 100644 (file)
index 0000000..9eb23a1
Binary files /dev/null and b/doc/doxygen/images/chebyshev8.png differ
diff --git a/doc/news/changes/minor/20170801DenisDavydov-c b/doc/news/changes/minor/20170801DenisDavydov-c
new file mode 100644 (file)
index 0000000..8d0d164
--- /dev/null
@@ -0,0 +1,4 @@
+New: add Utilities::LinearAlgebra::Lanczos_largest_eigenvalue() to estimate the largest
+eigenvalue of a symmetric linear operator using k-steps of Lanczos algorithm.
+<br>
+(Denis Davydov, 2017/08/01)
diff --git a/doc/news/changes/minor/20170802DenisDavydov b/doc/news/changes/minor/20170802DenisDavydov
new file mode 100644 (file)
index 0000000..3aee291
--- /dev/null
@@ -0,0 +1,4 @@
+New: add Utilities::LinearAlgebra::Chebyshev_filter() to apply Chebyshev filter of a
+given degree.
+<br>
+(Denis Davydov, 2017/08/02)
index 621ece07099be4e6410e2fc6741a99187381856f..83480f965b1b9f12e654626047dd586c1f5acb86 100644 (file)
@@ -179,8 +179,11 @@ public:
                    const VectorType &src) const;
 
   /**
-   * Initialize vector @p dst. This is a part of the interface required
-   * by linear operator.
+   * Initialize vector @p dst to have the same size and partition as
+   * @p diagonal member of this class.
+   *
+   * This is a part of the interface required
+   * by linear_operator().
    */
   void initialize_dof_vector(VectorType &dst) const;
 
diff --git a/include/deal.II/lac/utilities.h b/include/deal.II/lac/utilities.h
new file mode 100644 (file)
index 0000000..8383159
--- /dev/null
@@ -0,0 +1,318 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 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 dealii__lac_utilities_h
+#define dealii__lac_utilities_h
+
+#include <deal.II/base/config.h>
+#include <deal.II/base/exceptions.h>
+#include <deal.II/base/signaling_nan.h>
+#include <deal.II/lac/lapack_templates.h>
+#include <deal.II/lac/lapack_support.h>
+#include <deal.II/lac/vector_memory.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace Utilities
+{
+  /**
+   * A collection of linear-algebra utilities.
+   */
+  namespace LinearAlgebra
+  {
+
+    /**
+     * Estimate an upper bound for the largest eigenvalue of @p H by a @p k -step
+     * Lanczos process starting from the initial vector @p v0. Typical
+     * values of @p k are below 10. This estimator computes a k-step Lanczos
+     * decomposition $H V_k=V_k T_k+f_k e_k^T$ where $V_k$ contains k Lanczos basis,
+     * $V_k^TV_k=I_k$, $T_k$ is the tridiagonal Lanczos matrix, $f_k$ is a residual
+     * vector $f_k^TV_k=0$, and $e_k$ is the k-th canonical basis of $R^k$.
+     * The returned value is $ ||T_k||_2 + ||f_k||_2$.
+     * If @p eigenvalues is not <code>nullptr</code>, the eigenvalues of $T_k$ will be written there.
+     *
+     * @p vector_memory is used to allocate memory for temporary vectors.
+     * OperatorType has to provide <code>vmult</code> operation with
+     * VectorType.
+     *
+     * This function implements the algorithm from
+     * @code{.bib}
+     * @Article{Zhou2006,
+     *  Title     = {Self-consistent-field Calculations Using Chebyshev-filtered Subspace Iteration},
+     *  Author    = {Zhou, Yunkai and Saad, Yousef and Tiago, Murilo L. and Chelikowsky, James R.},
+     *  Journal   = {Journal of Computational Physics},
+     *  Year      = {2006},
+     *  Volume    = {219},
+     *  Pages     = {172--184},
+     * }
+     * @endcode
+     *
+     * @note This function uses Lapack routines to compute the largest
+     * eigenvalue of $T_k$.
+     *
+     * @note This function provides an alternate estimate to that obtained from
+     * several steps of SolverCG with SolverCG<VectorType>::connect_eigenvalues_slot().
+     *
+     * @author Denis Davydov, 2017
+     */
+    template <typename OperatorType, typename VectorType>
+    double Lanczos_largest_eigenvalue(const OperatorType &H,
+                                      const VectorType &v0,
+                                      const unsigned int k,
+                                      VectorMemory<VectorType> &vector_memory,
+                                      std::vector<double> *eigenvalues = nullptr);
+
+    /**
+     * Apply Chebyshev polynomial of the operator @p H to @p x. For a
+     * non-defective operator $H$ with a complete set of eigenpairs
+     * $H \psi_i = \lambda_i \psi_i$, the action of a polynomial filter $p$ is given by
+     * $p(H)x =\sum_i a_i p(\lambda_i) \psi_i$, where $x=: \sum_i a_i \psi_i$. Thus
+     * by appropriately choosing the polynomial filter, one can alter
+     * the eigenmodes contained in $x$.
+     *
+     * This function uses Chebyshev polynomials of first kind. Below is an
+     * example of polynomial $T_n(x)$ of degree $n=8$ normalized to unity at $-1.2$.
+     * <table>
+     *   <tr>
+     *     <td align="center">
+     *       @image html chebyshev8.png
+     *     </td>
+     *   </tr>
+     * </table>
+     * By introducing a linear mapping $L$ from @p unwanted_spectrum to
+     * $[-1,1]$, we can dump the corresponding modes in @p x. The higher
+     * the polynomial degree $n$, the more rapid it grows outside of the
+     * $[-1,1]$. In order to avoid numerical overflow, we normalize
+     * polynomial filter to unity at @p tau. Thus, the filtered operator
+     * is $p(H) = T_n(L(H))/T_n(L(\tau))$.
+     *
+     * The action of the Chebyshev filter only requires
+     * evaluation of <code>vmult()</code> of @p H and is based on the
+     * recursion equation for Chebyshev polynomial of degree $n$:
+     * $T_{n}(x) = 2x T_{n-1}(x) - T_{n-2}(x)$ with $T_0(x)=1$ and $T_1(x)=x$.
+     *
+     * @p vector_memory is used to allocate memory for temporary objects.
+     *
+     * This function implements the algorithm (with a minor fix of sign of $\sigma_1$) from
+     * @code{.bib}
+     * @Article{Zhou2014,
+     *  Title     = {Chebyshev-filtered subspace iteration method free of sparse diagonalization for solving the Kohn--Sham equation},
+     *  Author    = {Zhou, Yunkai and Chelikowsky, James R and Saad, Yousef},
+     *  Journal   = {Journal of Computational Physics},
+     *  Year      = {2014},
+     *  Volume    = {274},
+     *  Pages     = {770--782},
+     * }
+     * @endcode
+     *
+     * @note If @p tau is equal to
+     * <code>std::numeric_limits<double>::infinity()</code>, no normalization
+     * will be performed.
+     *
+     * @author Denis Davydov, 2017
+     */
+    template <typename OperatorType, typename VectorType>
+    void Chebyshev_filter(VectorType &x,
+                          const OperatorType &H,
+                          const unsigned int n,
+                          const std::pair<double,double> unwanted_spectrum,
+                          const double tau,
+                          VectorMemory<VectorType> &vector_memory);
+
+  }
+
+}
+
+
+/*------------------------- Implementation ----------------------------*/
+
+#ifndef DOXYGEN
+
+namespace Utilities
+{
+  namespace LinearAlgebra
+  {
+    template <typename OperatorType, typename VectorType>
+    double Lanczos_largest_eigenvalue(const OperatorType &H,
+                                      const VectorType &v0_,
+                                      const unsigned int k,
+                                      VectorMemory<VectorType> &vector_memory,
+                                      std::vector<double> *eigenvalues)
+    {
+      // Do k-step Lanczos:
+
+      typename VectorMemory<VectorType>::Pointer v(vector_memory);
+      typename VectorMemory<VectorType>::Pointer v0(vector_memory);
+      typename VectorMemory<VectorType>::Pointer f(vector_memory);
+
+      v->reinit(v0_);
+      v0->reinit(v0_);
+      f->reinit(v0_);
+
+      // two vectors to store diagonal and subdiagonal of the Lanczos
+      // matrix
+      std::vector<double> diagonal;
+      std::vector<double> subdiagonal;
+
+      // scalars to store norms and inner products
+      double a = 0, b = 0;
+
+      // 1. Normalize input vector
+      (*v) = v0_;
+      a = v->l2_norm();
+      Assert (a!=0, ExcDivideByZero());
+      (*v) *= 1./a;
+
+      // 2. Compute f = Hv; a = f*v; f <- f - av; T(0,0)=a;
+      H.vmult(*f,*v);
+      a = (*f)*(*v);
+      f->add(-a,*v);
+      diagonal.push_back(a);
+
+      // 3. Loop over steps
+      for (unsigned int i = 1; i < k; ++i)
+        {
+          // 4. L2 norm of f
+          b = f->l2_norm();
+          Assert (b!=0, ExcDivideByZero());
+          // 5. v0 <- v; v <- f/b
+          *v0 = *v;
+          *v = *f;
+          (*v) *= 1./b;
+          // 6. f = Hv; f <- f - b v0;
+          H.vmult(*f,*v);
+          f->add(-b, *v0);
+          // 7. a = f*v; f <- f - a v;
+          a = (*f) * (*v);
+          f->add(-a, *v);
+          // 8. T(i,i-1) = T(i-1,i) = b;  T(i,i) = a;
+          diagonal.push_back(a);
+          subdiagonal.push_back(b);
+        }
+
+      Assert (diagonal.size() == k,
+              ExcInternalError());
+      Assert (subdiagonal.size() == k-1,
+              ExcInternalError());
+
+      // Use Lapack dstev to get ||T||_2 norm, i.e. the largest eigenvalue
+      // of T
+      const int n = k;
+      std::vector<double> Z;    // unused for eigenvalues-only ("N") job
+      const int ldz = 1;        // ^^   (>=1)
+      std::vector<double> work; // ^^
+      int info;
+      // call lapack_templates.h wrapper:
+      stev ("N", &n,
+            &diagonal[0], &subdiagonal[0],
+            &Z[0], &ldz, &work[0],
+            &info);
+
+      Assert (info == 0,
+              LAPACKSupport::ExcErrorCode("dstev", info));
+
+      if (eigenvalues != nullptr)
+        {
+          eigenvalues->resize(diagonal.size());
+          std::copy(diagonal.begin(), diagonal.end(),
+                    eigenvalues->begin());
+        }
+
+      // note that the largest eigenvalue of T is below the largest
+      // eigenvalue of the operator.
+      // return ||T||_2 + ||f||_2, although it is not guaranteed to be an upper bound.
+      return diagonal[k-1] + f->l2_norm();
+    }
+
+
+    template <typename OperatorType, typename VectorType>
+    void Chebyshev_filter(VectorType &x,
+                          const OperatorType &op,
+                          const unsigned int degree,
+                          const std::pair<double,double> unwanted_spectrum,
+                          const double a_L,
+                          VectorMemory<VectorType> &vector_memory)
+    {
+      const double a = unwanted_spectrum.first;
+      const double b = unwanted_spectrum.second;
+      Assert (degree > 0,
+              ExcMessage ("Only positive degrees make sense."));
+
+      const bool scale = (a_L < std::numeric_limits<double>::infinity());
+      Assert (a < b,
+              ExcMessage("Lower bound of the unwanted spectrum should be smaller than the upper bound."));
+
+      Assert (a_L <= a || a_L >= b || !scale,
+              ExcMessage("Scaling point should be outside of the unwanted spectrum."));
+
+      // Setup auxiliary vectors:
+      typename VectorMemory<VectorType>::Pointer p_y(vector_memory);
+      typename VectorMemory<VectorType>::Pointer p_yn(vector_memory);
+
+      p_y->reinit(x);
+      p_yn->reinit(x);
+
+      // convenience to avoid pointers
+      VectorType &y   = *p_y;
+      VectorType &yn  = *p_yn;
+
+      // Below is an implementation of
+      // Algorithm 3.2 in Zhou et al, Journal of Computational Physics 274 (2014) 770-782
+      // with **a bugfix for sigma1**. Here is the original algorithm verbatim:
+      //
+      // [Y]=Chebyshev_filter_scaled(X, m, a, b, aL).
+      // e=(b−a)/2; c=(a+b)/2; σ=e/(c−aL); τ=2/σ;
+      // Y=(H∗X−c∗X)∗(σ/e);
+      // for i=2 to m do
+      //   σnew =1/(τ −σ);
+      //   Yt =(H∗Y−c∗Y)∗(2∗σnew/e)−(σ∗σnew)∗X;
+      //   X =Y; Y =Yt; σ =σnew;
+
+      const double e = (b-a)/2.;
+      const double c = (a+b)/2.;
+      const double alpha = 1./e;
+      const double beta = - c/e;
+
+      const double sigma1 = e/(a_L - c); // BUGFIX which is relevant for odd degrees
+      double sigma = scale ? sigma1 : 1.;
+      const double tau = 2./sigma;
+      op.vmult(y,x);
+      y.sadd(alpha*sigma, beta*sigma, x);
+
+      for (unsigned int i = 2; i <= degree; ++i)
+        {
+          const double sigma_new = scale ? 1./(tau-sigma) : 1.;
+          op.vmult(yn,y);
+          yn.sadd(2.*alpha*sigma_new, 2.*beta*sigma_new, y);
+          yn.add(-sigma*sigma_new, x);
+          x.swap(y);
+          y.swap(yn);
+          sigma = sigma_new;
+        }
+
+      x.swap(y);
+    }
+
+  }
+}
+
+#endif
+
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+
+#endif
diff --git a/tests/lac/utilities_01.cc b/tests/lac/utilities_01.cc
new file mode 100644 (file)
index 0000000..f747722
--- /dev/null
@@ -0,0 +1,97 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 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.
+//
+// ---------------------------------------------------------------------
+
+
+// test estimate of largest eigenvalue of a matrix.
+// The matrix is the same as in slepc/solve_04 test which has eingenvalues:
+// 3.98974 > 3.95906 > 3.90828 > 3.83792
+
+#include "../tests.h"
+#include "../testmatrix.h"
+#include "../slepc/testmatrix.h"
+#include <cmath>
+#include <fstream>
+#include <iostream>
+#include <iomanip>
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/petsc_compatibility.h>
+#include <deal.II/lac/petsc_sparse_matrix.h>
+#include <deal.II/lac/petsc_parallel_vector.h>
+#include <deal.II/lac/vector_memory.h>
+#include <deal.II/lac/utilities.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_cg.h>
+#include <typeinfo>
+
+int main(int argc, char **argv)
+{
+  initlog();
+  deallog << std::setprecision(6);
+  deallog.threshold_double(1.e-10);
+
+  Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1);
+  {
+
+    const unsigned int size = 31;
+    unsigned int dim = (size-1);
+
+    deallog << "Size " << size << " Unknowns " << dim << std::endl << std::endl;
+
+    // Make matrix
+    FD1DLaplaceMatrix testproblem(size);
+    PETScWrappers::SparseMatrix  A(dim, dim, 3);
+    testproblem.three_point(A);
+    A.compress (VectorOperation::insert);
+
+    PETScWrappers::MPI::Vector v0(MPI_COMM_WORLD, dim, dim);
+    PETScWrappers::MPI::Vector y(MPI_COMM_WORLD, dim, dim);
+    PETScWrappers::MPI::Vector x(MPI_COMM_WORLD, dim, dim);
+    for (unsigned int j=0; j<v0.size(); ++j)
+      v0[j] = static_cast<double>(Testing::rand())/static_cast<double>(RAND_MAX);
+
+    v0.compress(VectorOperation::insert);
+    GrowingVectorMemory<PETScWrappers::MPI::Vector> vector_memory;
+
+    for (unsigned int k = 4; k < 10; ++k)
+      {
+        const double est = Utilities::LinearAlgebra::Lanczos_largest_eigenvalue(A,v0,k,vector_memory);
+        Assert (est > 3.98974, ExcInternalError());
+        deallog << k << std::endl
+                << "Lanczos " << est << std::endl;
+
+        // estimate from CG
+        {
+          ReductionControl control (k,
+                                    std::sqrt(std::numeric_limits<double>::epsilon()),
+                                    1e-10, false, false);
+          std::vector<double> estimated_eigenvalues;
+          SolverCG<PETScWrappers::MPI::Vector> solver (control);
+          solver.connect_eigenvalues_slot([&estimated_eigenvalues] (const std::vector<double> &ev) -> void {estimated_eigenvalues = ev;});
+          y = v0;
+          PreconditionIdentity preconditioner;
+          try
+            {
+              solver.solve(A, x, y, preconditioner);
+            }
+          catch (SolverControl::NoConvergence &)
+            {
+            }
+
+          deallog << "CG " <<  estimated_eigenvalues.back() << std::endl;
+        }
+      }
+  }
+
+}
diff --git a/tests/lac/utilities_01.output b/tests/lac/utilities_01.output
new file mode 100644 (file)
index 0000000..fd5aca8
--- /dev/null
@@ -0,0 +1,21 @@
+
+DEAL::Size 31 Unknowns 30
+DEAL::
+DEAL::4
+DEAL::Lanczos 4.83338
+DEAL::CG 3.72123
+DEAL::5
+DEAL::Lanczos 4.95881
+DEAL::CG 3.90336
+DEAL::6
+DEAL::Lanczos 4.81862
+DEAL::CG 3.82808
+DEAL::7
+DEAL::Lanczos 4.98166
+DEAL::CG 3.94033
+DEAL::8
+DEAL::Lanczos 4.70266
+DEAL::CG 3.89841
+DEAL::9
+DEAL::Lanczos 4.93451
+DEAL::CG 3.97020
diff --git a/tests/lac/utilities_02.cc b/tests/lac/utilities_02.cc
new file mode 100644 (file)
index 0000000..b754902
--- /dev/null
@@ -0,0 +1,299 @@
+/* ---------------------------------------------------------------------
+ *
+ * Copyright (C) 2017 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.
+ *
+ * ---------------------------------------------------------------------
+
+ *
+ * Test estimated largest eigenvalue of M^{-1/2} L M^{-1/2} using
+ * k-steps of Lanczos algorithm. Here M is diagoal mass matrix obtained
+ * from Gauss-Legendre-Lobatto quadrature and L is Laplace operator.
+ *
+ * Largest eigenvalues from pArpack are:
+ *
+ * DEAL::1014.26
+ * DEAL::1018.29
+ * DEAL::1018.29
+ * DEAL::1020.72
+ * DEAL::1020.72
+ *
+ */
+
+#include "../tests.h"
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/index_set.h>
+#include <deal.II/distributed/tria.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/matrix_free/operators.h>
+#include <deal.II/lac/linear_operator.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/lac/diagonal_matrix.h>
+#include <deal.II/lac/linear_operator.h>
+#include <deal.II/lac/utilities.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>
+
+// debug only: turn on calculation of eigenvalues of the operator by Arpack
+// #define PARPACK
+
+#ifdef PARPACK
+#include <deal.II/lac/parpack_solver.h>
+#endif
+
+#include <fstream>
+#include <iostream>
+
+
+const unsigned int dim = 2;
+
+using namespace dealii;
+
+const double eps = 1e-10;
+
+const unsigned int fe_degree = 1;
+
+void test ()
+{
+  const unsigned int global_mesh_refinement_steps = 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);
+
+  parallel::distributed::Triangulation<dim> triangulation (mpi_communicator);
+  GridGenerator::hyper_cube (triangulation, -1, 1);
+  triangulation.refine_global (global_mesh_refinement_steps);
+
+
+  DoFHandler<dim> dof_handler(triangulation);
+  FE_Q<dim> fe(fe_degree);
+  dof_handler.distribute_dofs (fe);
+
+
+  IndexSet locally_relevant_dofs;
+  DoFTools::extract_locally_relevant_dofs (dof_handler,
+                                           locally_relevant_dofs);
+  ConstraintMatrix constraints;
+  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 ();
+
+  std::shared_ptr<MatrixFree<dim,double> > mf_data(new MatrixFree<dim,double> ());
+  {
+    const QGauss<1> quad (fe_degree+1);
+    typename MatrixFree<dim,double>::AdditionalData data;
+    data.tasks_parallel_scheme =
+      MatrixFree<dim,double>::AdditionalData::partition_color;
+    data.mapping_update_flags = update_values | update_gradients | update_JxW_values;
+    mf_data->reinit (dof_handler, constraints, quad, data);
+  }
+
+  MatrixFreeOperators::MassOperator<dim,fe_degree, fe_degree+1, 1, LinearAlgebra::distributed::Vector<double> > mass;
+  MatrixFreeOperators::LaplaceOperator<dim,fe_degree, fe_degree+1, 1, LinearAlgebra::distributed::Vector<double> > laplace;
+  mass.initialize(mf_data);
+  laplace.initialize(mf_data);
+
+  // Gauss-Legendre-Lobatto mass-matrix:
+  DiagonalMatrix<LinearAlgebra::distributed::Vector<double>> diagonal_mass_inv;
+  {
+    LinearAlgebra::distributed::Vector<double> inv_mass_matrix;
+    VectorizedArray<double> one = make_vectorized_array (1.);
+    mf_data->initialize_dof_vector (inv_mass_matrix);
+    FEEvaluation<dim,fe_degree>   fe_eval(*mf_data);
+    const unsigned int            n_q_points = fe_eval.n_q_points;
+    for (unsigned int cell=0; cell<mf_data->n_macro_cells(); ++cell)
+      {
+        fe_eval.reinit(cell);
+        for (unsigned int q=0; q<n_q_points; ++q)
+          fe_eval.submit_value(one,q);
+        fe_eval.integrate (true,false);
+        fe_eval.distribute_local_to_global (inv_mass_matrix);
+      }
+    inv_mass_matrix.compress(VectorOperation::add);
+    for (unsigned int k=0; k<inv_mass_matrix.local_size(); ++k)
+      if (inv_mass_matrix.local_element(k)>1e-15)
+        {
+          inv_mass_matrix.local_element(k) = std::sqrt(1./inv_mass_matrix.local_element(k));
+        }
+      else
+        inv_mass_matrix.local_element(k) = 0;
+
+    diagonal_mass_inv.reinit(inv_mass_matrix);
+  }
+
+  const auto invM = linear_operator<LinearAlgebra::distributed::Vector<double>>(diagonal_mass_inv);
+  const auto OP =  invM *
+                   linear_operator<LinearAlgebra::distributed::Vector<double>>(laplace) *
+                   invM;
+
+  // Do actuall work:
+  LinearAlgebra::distributed::Vector<double> init_vector;
+  mf_data->initialize_dof_vector(init_vector);
+  for (auto it = init_vector.begin(); it != init_vector.end(); ++it)
+    *it = static_cast<double>(Testing::rand())/static_cast<double>(RAND_MAX);
+
+  constraints.set_zero(init_vector);
+
+  GrowingVectorMemory<LinearAlgebra::distributed::Vector<double>> vector_memory;
+  for (unsigned int k = 4; k < 10; ++k)
+    {
+      const double est = Utilities::LinearAlgebra::Lanczos_largest_eigenvalue(OP,init_vector,k,vector_memory);
+      deallog << k << " " << est << std::endl;
+    }
+
+  // exact eigenvectors via PArpack
+#ifdef PARPACK
+  {
+    const unsigned int number_of_eigenvalues        = 5;
+
+    std::vector<LinearAlgebra::distributed::Vector<double> > eigenfunctions;
+    std::vector<double>                                      eigenvalues;
+    eigenfunctions.resize (number_of_eigenvalues);
+    eigenvalues.resize (number_of_eigenvalues);
+    for (unsigned int i=0; i<eigenfunctions.size (); ++i)
+      mf_data->initialize_dof_vector (eigenfunctions[i]);
+
+    std::vector<std::complex<double> > lambda(number_of_eigenvalues);
+
+    const unsigned int num_arnoldi_vectors = 2*eigenvalues.size() + 10;
+    PArpackSolver<LinearAlgebra::distributed::Vector<double> >::AdditionalData
+    additional_data(num_arnoldi_vectors,
+                    PArpackSolver<LinearAlgebra::distributed::Vector<double> >::largest_magnitude,
+                    true,
+                    1);
+
+    SolverControl solver_control(
+      dof_handler.n_dofs(), 1e-10, /*log_history*/ false, /*log_results*/ false);
+
+    PArpackSolver<LinearAlgebra::distributed::Vector<double> > eigensolver(
+      solver_control, mpi_communicator, additional_data);
+
+    eigensolver.reinit(eigenfunctions[0]);
+    // make sure initial vector is orthogonal to the space due to constraints
+    {
+      LinearAlgebra::distributed::Vector<double> init_vector;
+      mf_data->initialize_dof_vector(init_vector);
+      for (auto it = init_vector.begin(); it != init_vector.end(); ++it)
+        *it = static_cast<double>(Testing::rand())/static_cast<double>(RAND_MAX);
+
+      constraints.set_zero(init_vector);
+      eigensolver.set_initial_vector(init_vector);
+    }
+    // avoid output of iterative solver:
+    const unsigned int previous_depth = deallog.depth_file(0);
+    eigensolver.solve (OP,
+                       mass,
+                       OP,
+                       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;
+
+    // make sure that we have eigenvectors and they are mass-orthonormal:
+    // a) (A*x_i-\lambda*x_i).L2() == 0
+    // b) x_j*x_i=\delta_{ij}
+    {
+      const double precision = 1e-7;
+      LinearAlgebra::distributed::Vector<double> Ax(eigenfunctions[0]);
+      for (unsigned int i=0; i < eigenfunctions.size(); ++i)
+        {
+          for (unsigned int j=0; j < eigenfunctions.size(); j++)
+            {
+              const double err = std::abs( eigenfunctions[j] * eigenfunctions[i] - (i==j));
+              Assert( err< precision,
+                      ExcMessage("Eigenvectors " +
+                                 Utilities::int_to_string(i) +
+                                 " and " +
+                                 Utilities::int_to_string(j) +
+                                 " are not orthonormal: " +
+                                 std::to_string(err)));
+            }
+
+          OP.vmult(Ax,eigenfunctions[i]);
+          Ax.add(-1.0*eigenvalues[i],eigenfunctions[i]);
+          const double err = Ax.l2_norm();
+          Assert (err < precision,
+                  ExcMessage("Returned vector " +
+                             Utilities::int_to_string(i) +
+                             " is not an eigenvector: " +
+                             std::to_string(err)));
+        }
+    }
+  }
+#endif
+
+  dof_handler.clear ();
+  deallog << "Ok"<<std::endl;
+}
+
+
+int main (int argc,char **argv)
+{
+  std::ofstream logfile("output");
+  deallog.attach(logfile,/*do not print job id*/false);
+  deallog.threshold_double(eps);
+
+  try
+    {
+      Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+      {
+        test ();
+      }
+
+    }
+  catch (std::exception &exc)
+    {
+      std::cerr << std::endl << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+                << exc.what() << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+
+      return 1;
+    }
+  catch (...)
+    {
+      std::cerr << std::endl << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      return 1;
+    };
+}
diff --git a/tests/lac/utilities_02.with_mpi=true.with_p4est=true.mpirun=3.output b/tests/lac/utilities_02.with_mpi=true.with_p4est=true.mpirun=3.output
new file mode 100644 (file)
index 0000000..47e2344
--- /dev/null
@@ -0,0 +1,7 @@
+DEAL::4 1188.92
+DEAL::5 1221.68
+DEAL::6 1243.58
+DEAL::7 1234.77
+DEAL::8 1258.59
+DEAL::9 1256.07
+DEAL::Ok
diff --git a/tests/lac/utilities_03.cc b/tests/lac/utilities_03.cc
new file mode 100644 (file)
index 0000000..ab2c834
--- /dev/null
@@ -0,0 +1,194 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 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.
+//
+// ---------------------------------------------------------------------
+
+
+// Test Chebyshev filter on Diagonal matrix with equidistance eigenvalues in (-1,1).
+// Example is taken from Seciton 2.2. in
+// Pieper et al, Journal of Computational Physics 325 (2016), 226-243
+//
+// The test exploits the fact that:
+// x_0 =: a_i \psi_i
+// H \psi_i = \lambda_i \psi_i
+// p(H)x_0 = a_i p(\lambda_i) \psi_i
+
+// Note that linfty_norm() norm will show non-zero difference for non-scaled
+// case of high degree. This is ok, as the difference is small compared
+// to the absolute value on the order 10^12
+
+/* MWE in Maxima for degree 3 with scaling. The calculations are
+ * done for the mode which has the maximum difference in the resulting value
+
+Cheb2(n, x) :=
+block ( [],
+    if n = 0
+       then 1
+       else
+        if n = 1
+           then x
+           else expand(2*x*Cheb2 (n - 1, x)
+                        - Cheb2 (n - 2, x))
+);
+min:-0.01;
+max:0.01;
+c:(max+min)/2;
+h:(max-min)/2;
+L(x):= (x-c)/h;
+ev: -0.0989011;
+x:  0.904187;
+aL: -0.1;
+Cheb2(3,L(aL));
+Cheb2(3,L(ev));
+x * Cheb2(3,L(ev)) / Cheb2(3,L(aL));
+
+// the last three produce:
+
+-3970.0
+-3839.905459408033
+0.8745573293767684
+
+ */
+
+//#define EXTRA_OUTPUT
+
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/diagonal_matrix.h>
+#include <deal.II/lac/la_parallel_vector.h>
+#include <deal.II/lac/utilities.h>
+
+double cheb2(const unsigned int d, const double x)
+{
+  if (d == 0)
+    {
+      return 1.;
+    }
+  else if (d == 1)
+    {
+      return x;
+    }
+
+  return 2.*x*cheb2(d-1,x) - cheb2(d-2,x);
+}
+
+
+
+void
+check(const int degree, const bool scale = false, const double a_L = -0.1, const double a = -0.01, const double b = 0.01, const unsigned int size = 1000)
+{
+  deallog << "Degree " << degree << std::endl;
+  LinearAlgebra::distributed::Vector<double> ev(size), x(size), y(size), exact(size), diff(size);
+  GrowingVectorMemory<LinearAlgebra::distributed::Vector<double>> vector_memory;
+
+  for (unsigned int i=0; i<size; ++i)
+    ev(i) = -1.+2.*(1.+i)/(1.+size);
+  DiagonalMatrix<LinearAlgebra::distributed::Vector<double>> mat;
+  mat.reinit(ev);
+
+  x = 0.;
+  // prevent overflow by not perturbing modes far away from the region
+  // to be filtered
+  unsigned int n_in = 0;
+  unsigned int n_out = 0;
+  for (unsigned int i=0; i<size; ++i)
+    if (std::abs(ev(i)) <= std::abs(a_L))
+      {
+        if ( ev(i) >= a && ev(i) <= b)
+          n_in++;
+        else
+          n_out++;
+        x(i) = static_cast<double>(Testing::rand())/static_cast<double>(RAND_MAX);
+      }
+
+  deallog << " Modes inside/outside: " << n_in << " " << n_out << std::endl;
+
+  // for x = x_i v_i , where v_i are eigenvectors
+  // p[H]x = \sum_i x_i p(\lambda_i) v_i
+  const double c = (a+b)/2.;
+  const double e = (b-a)/2.;
+  auto L = [&](const double &x)
+  {
+    return (x-c)/e;
+  };
+
+  const double scaling = scale ? cheb2(degree,L(a_L)) : 1.; // p(L(a_L))
+  deallog << " Scaling: " << scaling << " @ " << a_L << std::endl;
+  exact = 0.;
+  for (unsigned int i=0; i<size; ++i)
+    exact(i) = x(i) * cheb2(degree,L(ev(i))) / scaling;
+  deallog << " Input norm: " << x.l2_norm() << std::endl;
+  deallog << " Exact norm: " << exact.l2_norm() << std::endl;
+
+  const double g_ = (scale ? a_L :  std::numeric_limits<double>::infinity() );
+  y = x;
+  Utilities::LinearAlgebra::Chebyshev_filter(y, mat, degree, std::make_pair(a, b), g_, vector_memory);
+  diff = y;
+  diff -=exact;
+
+  deallog << " Filter ["<<a<<","<<b<<"]" << std::endl;
+  deallog << " Error: " << diff.linfty_norm() << std::endl;
+
+#ifdef EXTRA_OUTPUT
+  // extra output for debugging:
+  unsigned int max_i = 0;
+  for (unsigned int i = 1; i < size; ++i)
+    if ( std::abs(diff(i)) > std::abs(diff(max_i)))
+      max_i = i;
+
+  deallog << " i =" << max_i << std::endl
+          << " d =" << diff(max_i) << std::endl
+          << " ev=" << ev(max_i) << std::endl
+          << " x =" << x(max_i) << std::endl
+          << " y =" << y(max_i) << std::endl
+          << " ex=" << exact(max_i) << std::endl;
+#endif
+}
+
+
+int main()
+{
+  std::ofstream logfile("output");
+  deallog << std::setprecision(6);
+  deallog.attach(logfile);
+  deallog.threshold_double(1.e-10);
+
+  deallog << "No scaling:" << std::endl;
+  // no scaling:
+  check(1);
+  check(2);
+  check(3);
+  check(4);
+  check(10);
+
+  deallog << "Lower scaling:" << std::endl;
+  // scaling at the lower end
+  check(1,true);
+  check(2,true);
+  check(3,true);
+  check(4,true);
+  check(10,true);
+  check(30,true);
+
+  deallog << "Upper scaling:" << std::endl;
+  // scaling at the upper end
+  check(1,true,  0.1);
+  check(2,true,  0.1);
+  check(3,true,  0.1);
+  check(4,true,  0.1);
+  check(10,true, 0.1);
+  check(30,true, 0.1);
+
+  return 0;
+}
diff --git a/tests/lac/utilities_03.output b/tests/lac/utilities_03.output
new file mode 100644 (file)
index 0000000..5c0d203
--- /dev/null
@@ -0,0 +1,123 @@
+
+DEAL::No scaling:
+DEAL::Degree 1
+DEAL:: Modes inside/outside: 10 90
+DEAL:: Scaling: 1.00000 @ -0.100000
+DEAL:: Input norm: 6.16127
+DEAL:: Exact norm: 37.4711
+DEAL:: Filter [-0.0100000,0.0100000]
+DEAL:: Error: 0
+DEAL::Degree 2
+DEAL:: Modes inside/outside: 10 90
+DEAL:: Scaling: 1.00000 @ -0.100000
+DEAL:: Input norm: 5.92779
+DEAL:: Exact norm: 586.339
+DEAL:: Filter [-0.0100000,0.0100000]
+DEAL:: Error: 0
+DEAL::Degree 3
+DEAL:: Modes inside/outside: 10 90
+DEAL:: Scaling: 1.00000 @ -0.100000
+DEAL:: Input norm: 5.89974
+DEAL:: Exact norm: 10657.5
+DEAL:: Filter [-0.0100000,0.0100000]
+DEAL:: Error: 0
+DEAL::Degree 4
+DEAL:: Modes inside/outside: 10 90
+DEAL:: Scaling: 1.00000 @ -0.100000
+DEAL:: Input norm: 5.60617
+DEAL:: Exact norm: 157052.
+DEAL:: Filter [-0.0100000,0.0100000]
+DEAL:: Error: 0
+DEAL::Degree 10
+DEAL:: Modes inside/outside: 10 90
+DEAL:: Scaling: 1.00000 @ -0.100000
+DEAL:: Input norm: 5.36039
+DEAL:: Exact norm: 5.05476e+12
+DEAL:: Filter [-0.0100000,0.0100000]
+DEAL:: Error: 0.00146484
+DEAL::Lower scaling:
+DEAL::Degree 1
+DEAL:: Modes inside/outside: 10 90
+DEAL:: Scaling: -10.0000 @ -0.100000
+DEAL:: Input norm: 6.43748
+DEAL:: Exact norm: 3.99850
+DEAL:: Filter [-0.0100000,0.0100000]
+DEAL:: Error: 0
+DEAL::Degree 2
+DEAL:: Modes inside/outside: 10 90
+DEAL:: Scaling: 199.000 @ -0.100000
+DEAL:: Input norm: 5.74697
+DEAL:: Exact norm: 2.70598
+DEAL:: Filter [-0.0100000,0.0100000]
+DEAL:: Error: 0
+DEAL::Degree 3
+DEAL:: Modes inside/outside: 10 90
+DEAL:: Scaling: -3970.00 @ -0.100000
+DEAL:: Input norm: 5.64201
+DEAL:: Exact norm: 2.34171
+DEAL:: Filter [-0.0100000,0.0100000]
+DEAL:: Error: 0
+DEAL::Degree 4
+DEAL:: Modes inside/outside: 10 90
+DEAL:: Scaling: 79201.0 @ -0.100000
+DEAL:: Input norm: 5.70416
+DEAL:: Exact norm: 1.92446
+DEAL:: Filter [-0.0100000,0.0100000]
+DEAL:: Error: 0
+DEAL::Degree 10
+DEAL:: Modes inside/outside: 10 90
+DEAL:: Scaling: 4.99312e+12 @ -0.100000
+DEAL:: Input norm: 5.64307
+DEAL:: Exact norm: 0.989370
+DEAL:: Filter [-0.0100000,0.0100000]
+DEAL:: Error: 0
+DEAL::Degree 30
+DEAL:: Modes inside/outside: 10 90
+DEAL:: Scaling: 4.97938e+38 @ -0.100000
+DEAL:: Input norm: 5.82031
+DEAL:: Exact norm: 0.745548
+DEAL:: Filter [-0.0100000,0.0100000]
+DEAL:: Error: 0
+DEAL::Upper scaling:
+DEAL::Degree 1
+DEAL:: Modes inside/outside: 10 90
+DEAL:: Scaling: 10.0000 @ 0.100000
+DEAL:: Input norm: 5.51882
+DEAL:: Exact norm: 3.18201
+DEAL:: Filter [-0.0100000,0.0100000]
+DEAL:: Error: 0
+DEAL::Degree 2
+DEAL:: Modes inside/outside: 10 90
+DEAL:: Scaling: 199.000 @ 0.100000
+DEAL:: Input norm: 5.83880
+DEAL:: Exact norm: 2.78386
+DEAL:: Filter [-0.0100000,0.0100000]
+DEAL:: Error: 0
+DEAL::Degree 3
+DEAL:: Modes inside/outside: 10 90
+DEAL:: Scaling: 3970.00 @ 0.100000
+DEAL:: Input norm: 5.49421
+DEAL:: Exact norm: 2.11010
+DEAL:: Filter [-0.0100000,0.0100000]
+DEAL:: Error: 0
+DEAL::Degree 4
+DEAL:: Modes inside/outside: 10 90
+DEAL:: Scaling: 79201.0 @ 0.100000
+DEAL:: Input norm: 5.55782
+DEAL:: Exact norm: 1.56672
+DEAL:: Filter [-0.0100000,0.0100000]
+DEAL:: Error: 0
+DEAL::Degree 10
+DEAL:: Modes inside/outside: 10 90
+DEAL:: Scaling: 4.99312e+12 @ 0.100000
+DEAL:: Input norm: 5.85002
+DEAL:: Exact norm: 1.24742
+DEAL:: Filter [-0.0100000,0.0100000]
+DEAL:: Error: 0
+DEAL::Degree 30
+DEAL:: Modes inside/outside: 10 90
+DEAL:: Scaling: 4.97938e+38 @ 0.100000
+DEAL:: Input norm: 6.25978
+DEAL:: Exact norm: 0.896899
+DEAL:: Filter [-0.0100000,0.0100000]
+DEAL:: Error: 0

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.