--- /dev/null
+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)
--- /dev/null
+New: add Utilities::LinearAlgebra::Chebyshev_filter() to apply Chebyshev filter of a
+given degree.
+<br>
+(Denis Davydov, 2017/08/02)
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;
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+ }
+ }
+ }
+
+}
--- /dev/null
+
+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
--- /dev/null
+/* ---------------------------------------------------------------------
+ *
+ * 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;
+ };
+}
--- /dev/null
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+}
--- /dev/null
+
+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
projects / dealii.git / commitdiff
