From: Denis Davydov Date: Tue, 1 Aug 2017 20:33:27 +0000 (+0200) Subject: add Utilities::LinearAlgebra::Lanczos_largest_eigenvalue() X-Git-Tag: v9.0.0-rc1~1355^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=e1e59d2fba8c07a281599c69f8869c8d305eb826;p=dealii.git add Utilities::LinearAlgebra::Lanczos_largest_eigenvalue() and Utilities::LinearAlgebra::Chebyshev_filter() --- diff --git a/doc/doxygen/images/chebyshev8.png b/doc/doxygen/images/chebyshev8.png new file mode 100644 index 0000000000..9eb23a1a0b 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 index 0000000000..8d0d164fda --- /dev/null +++ b/doc/news/changes/minor/20170801DenisDavydov-c @@ -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. +
+(Denis Davydov, 2017/08/01) diff --git a/doc/news/changes/minor/20170802DenisDavydov b/doc/news/changes/minor/20170802DenisDavydov new file mode 100644 index 0000000000..3aee291080 --- /dev/null +++ b/doc/news/changes/minor/20170802DenisDavydov @@ -0,0 +1,4 @@ +New: add Utilities::LinearAlgebra::Chebyshev_filter() to apply Chebyshev filter of a +given degree. +
+(Denis Davydov, 2017/08/02) diff --git a/include/deal.II/lac/diagonal_matrix.h b/include/deal.II/lac/diagonal_matrix.h index 621ece0709..83480f965b 100644 --- a/include/deal.II/lac/diagonal_matrix.h +++ b/include/deal.II/lac/diagonal_matrix.h @@ -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 index 0000000000..8383159838 --- /dev/null +++ b/include/deal.II/lac/utilities.h @@ -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 +#include +#include +#include +#include +#include + +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 nullptr, the eigenvalues of $T_k$ will be written there. + * + * @p vector_memory is used to allocate memory for temporary vectors. + * OperatorType has to provide vmult 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::connect_eigenvalues_slot(). + * + * @author Denis Davydov, 2017 + */ + template + double Lanczos_largest_eigenvalue(const OperatorType &H, + const VectorType &v0, + const unsigned int k, + VectorMemory &vector_memory, + std::vector *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$. + * + * + * + * + *
+ * @image html chebyshev8.png + *
+ * 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 vmult() 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 + * std::numeric_limits::infinity(), no normalization + * will be performed. + * + * @author Denis Davydov, 2017 + */ + template + void Chebyshev_filter(VectorType &x, + const OperatorType &H, + const unsigned int n, + const std::pair unwanted_spectrum, + const double tau, + VectorMemory &vector_memory); + + } + +} + + +/*------------------------- Implementation ----------------------------*/ + +#ifndef DOXYGEN + +namespace Utilities +{ + namespace LinearAlgebra + { + template + double Lanczos_largest_eigenvalue(const OperatorType &H, + const VectorType &v0_, + const unsigned int k, + VectorMemory &vector_memory, + std::vector *eigenvalues) + { + // Do k-step Lanczos: + + typename VectorMemory::Pointer v(vector_memory); + typename VectorMemory::Pointer v0(vector_memory); + typename VectorMemory::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 diagonal; + std::vector 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 Z; // unused for eigenvalues-only ("N") job + const int ldz = 1; // ^^ (>=1) + std::vector 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 + void Chebyshev_filter(VectorType &x, + const OperatorType &op, + const unsigned int degree, + const std::pair unwanted_spectrum, + const double a_L, + VectorMemory &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::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::Pointer p_y(vector_memory); + typename VectorMemory::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 index 0000000000..f747722c37 --- /dev/null +++ b/tests/lac/utilities_01.cc @@ -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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +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(Testing::rand())/static_cast(RAND_MAX); + + v0.compress(VectorOperation::insert); + GrowingVectorMemory 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::epsilon()), + 1e-10, false, false); + std::vector estimated_eigenvalues; + SolverCG solver (control); + solver.connect_eigenvalues_slot([&estimated_eigenvalues] (const std::vector &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 index 0000000000..fd5aca82d7 --- /dev/null +++ b/tests/lac/utilities_01.output @@ -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 index 0000000000..b754902120 --- /dev/null +++ b/tests/lac/utilities_02.cc @@ -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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +#include +#include +#include +#include + +// debug only: turn on calculation of eigenvalues of the operator by Arpack +// #define PARPACK + +#ifdef PARPACK +#include +#endif + +#include +#include + + +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 triangulation (mpi_communicator); + GridGenerator::hyper_cube (triangulation, -1, 1); + triangulation.refine_global (global_mesh_refinement_steps); + + + DoFHandler dof_handler(triangulation); + FE_Q 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 (), + constraints); + constraints.close (); + + std::shared_ptr > mf_data(new MatrixFree ()); + { + const QGauss<1> quad (fe_degree+1); + typename MatrixFree::AdditionalData data; + data.tasks_parallel_scheme = + MatrixFree::AdditionalData::partition_color; + data.mapping_update_flags = update_values | update_gradients | update_JxW_values; + mf_data->reinit (dof_handler, constraints, quad, data); + } + + MatrixFreeOperators::MassOperator > mass; + MatrixFreeOperators::LaplaceOperator > laplace; + mass.initialize(mf_data); + laplace.initialize(mf_data); + + // Gauss-Legendre-Lobatto mass-matrix: + DiagonalMatrix> diagonal_mass_inv; + { + LinearAlgebra::distributed::Vector inv_mass_matrix; + VectorizedArray one = make_vectorized_array (1.); + mf_data->initialize_dof_vector (inv_mass_matrix); + FEEvaluation fe_eval(*mf_data); + const unsigned int n_q_points = fe_eval.n_q_points; + for (unsigned int cell=0; celln_macro_cells(); ++cell) + { + fe_eval.reinit(cell); + for (unsigned int q=0; q1e-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>(diagonal_mass_inv); + const auto OP = invM * + linear_operator>(laplace) * + invM; + + // Do actuall work: + LinearAlgebra::distributed::Vector init_vector; + mf_data->initialize_dof_vector(init_vector); + for (auto it = init_vector.begin(); it != init_vector.end(); ++it) + *it = static_cast(Testing::rand())/static_cast(RAND_MAX); + + constraints.set_zero(init_vector); + + GrowingVectorMemory> 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 > eigenfunctions; + std::vector eigenvalues; + eigenfunctions.resize (number_of_eigenvalues); + eigenvalues.resize (number_of_eigenvalues); + for (unsigned int i=0; iinitialize_dof_vector (eigenfunctions[i]); + + std::vector > lambda(number_of_eigenvalues); + + const unsigned int num_arnoldi_vectors = 2*eigenvalues.size() + 10; + PArpackSolver >::AdditionalData + additional_data(num_arnoldi_vectors, + PArpackSolver >::largest_magnitude, + true, + 1); + + SolverControl solver_control( + dof_handler.n_dofs(), 1e-10, /*log_history*/ false, /*log_results*/ false); + + PArpackSolver > 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 init_vector; + mf_data->initialize_dof_vector(init_vector); + for (auto it = init_vector.begin(); it != init_vector.end(); ++it) + *it = static_cast(Testing::rand())/static_cast(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 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"< +#include +#include +#include + +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 ev(size), x(size), y(size), exact(size), diff(size); + GrowingVectorMemory> vector_memory; + + for (unsigned int i=0; i> 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= a && ev(i) <= b) + n_in++; + else + n_out++; + x(i) = static_cast(Testing::rand())/static_cast(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::infinity() ); + y = x; + Utilities::LinearAlgebra::Chebyshev_filter(y, mat, degree, std::make_pair(a, b), g_, vector_memory); + diff = y; + diff -=exact; + + deallog << " Filter ["< 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 index 0000000000..5c0d203ca7 --- /dev/null +++ b/tests/lac/utilities_03.output @@ -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