From: Daniel Arndt Date: Tue, 6 Mar 2018 04:39:24 +0000 (+0100) Subject: Allow for setting the size of the kernel instead of a threshold X-Git-Tag: v9.0.0-rc1~345^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F6003%2Fhead;p=dealii.git Allow for setting the size of the kernel instead of a threshold --- diff --git a/doc/news/changes/minor/20180307DanielArndt b/doc/news/changes/minor/20180307DanielArndt new file mode 100644 index 0000000000..0d9af92a97 --- /dev/null +++ b/doc/news/changes/minor/20180307DanielArndt @@ -0,0 +1,4 @@ +New: LAPACKFullMatrix::compute_inverse_svd_with_kernel allows to +set an expected kernel size for an inverse singular value decomposition. +
+(Daniel Arndt, 2018/03/07) diff --git a/include/deal.II/lac/lapack_full_matrix.h b/include/deal.II/lac/lapack_full_matrix.h index 19db50e6d7..0547f29e3c 100644 --- a/include/deal.II/lac/lapack_full_matrix.h +++ b/include/deal.II/lac/lapack_full_matrix.h @@ -746,6 +746,12 @@ public: */ void compute_inverse_svd (const double threshold = 0.); + /** + * Same as above but provide the size of the kernel instead of a threshold, + * i.e. the @p kernel_size smallest eigenvalues. + */ + void compute_inverse_svd_with_kernel (const unsigned int kernel_size); + /** * Retrieve eigenvalue after compute_eigenvalues() was called. */ diff --git a/include/deal.II/lac/relaxation_block.h b/include/deal.II/lac/relaxation_block.h index 140d6a32fa..f0ae55f962 100644 --- a/include/deal.II/lac/relaxation_block.h +++ b/include/deal.II/lac/relaxation_block.h @@ -131,10 +131,23 @@ public: * If #inversion is SVD, we can compute the Penrose-Moore inverse of the * blocks. In order to do so, we can specify here the threshold below * which a singular value will be considered zero and thus not inverted. + * Setting this parameter to a value greater than zero takes precedence over + * threshold, i.e. kernel_size must be zero if you want to use threshold. * This parameter is used in the call to * LAPACKFullMatrix::compute_inverse_svd(). */ - double threshold; + double threshold = 0.; + + /** + * If #inversion is SVD, we can compute the Penrose-Moore inverse of the + * blocks. In order to do so, we can specify here the size of the kernel + * that will not be inverted but considered zero. Setting this parameter + * to a value greater than zero takes precedence over threshold, i.e. + * kernel_size must be zero if you want to use threshold. + * This parameter is used in the call to + * LAPACKFullMatrix::compute_inverse_svd(). + */ + unsigned int kernel_size = 0; /** * The order in which blocks should be traversed. This vector can initiate diff --git a/include/deal.II/lac/relaxation_block.templates.h b/include/deal.II/lac/relaxation_block.templates.h index 36b47305f8..48cdcf18f0 100644 --- a/include/deal.II/lac/relaxation_block.templates.h +++ b/include/deal.II/lac/relaxation_block.templates.h @@ -156,7 +156,10 @@ RelaxationBlock::block_kernel (const case PreconditionBlockBase::svd: this->inverse_svd(block).reinit(bs, bs); this->inverse_svd(block) = M_cell; - this->inverse_svd(block).compute_inverse_svd(this->additional_data->threshold); + if (this->additional_data->kernel_size>0) + this->inverse_svd(block).compute_inverse_svd_with_kernel(this->additional_data->kernel_size); + else + this->inverse_svd(block).compute_inverse_svd(this->additional_data->threshold); break; default: Assert(false, ExcNotImplemented()); diff --git a/source/lac/lapack_full_matrix.cc b/source/lac/lapack_full_matrix.cc index cd23a761cd..698cd94454 100644 --- a/source/lac/lapack_full_matrix.cc +++ b/source/lac/lapack_full_matrix.cc @@ -1149,6 +1149,26 @@ LAPACKFullMatrix::compute_inverse_svd(const double threshold) } + +template +void +LAPACKFullMatrix::compute_inverse_svd_with_kernel(const unsigned int kernel_size) +{ + if (state == LAPACKSupport::matrix) + compute_svd(); + + Assert (state==LAPACKSupport::svd, ExcState(state)); + + const unsigned int n_wr = wr.size(); + for (size_type i=0; i void LAPACKFullMatrix::invert() diff --git a/tests/dofs/block_relaxation_01.cc b/tests/dofs/block_relaxation_01.cc new file mode 100644 index 0000000000..8a1f27a59a --- /dev/null +++ b/tests/dofs/block_relaxation_01.cc @@ -0,0 +1,164 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2018 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 that setting the kernel size in RelaxationBlock actually works. + + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +template +void +make_stokes_matrix (const DoFHandler &dof_handler, + const Quadrature &quadrature_formula, + SparseMatrix &system_matrix) +{ + const FiniteElement &fe = dof_handler.get_fe(); + const unsigned int degree = fe.degree; + system_matrix=0; + + ConstraintMatrix constraints; + constraints.close(); + FEValues fe_values (fe, quadrature_formula, + update_values | + update_quadrature_points | + update_JxW_values | + update_gradients); + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + FullMatrix local_matrix (dofs_per_cell, dofs_per_cell); + std::vector local_dof_indices (dofs_per_cell); + const FEValuesExtractors::Vector velocities (0); + const FEValuesExtractors::Scalar pressure (dim); + std::vector > symgrad_phi_u (dofs_per_cell); + std::vector div_phi_u (dofs_per_cell); + std::vector phi_p (dofs_per_cell); + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + { + fe_values.reinit (cell); + local_matrix = 0; + for (unsigned int q=0; qget_dof_indices (local_dof_indices); + constraints.distribute_local_to_global (local_matrix, local_dof_indices, + system_matrix); + } +} + + +template +void +check () +{ + Triangulation tr(Triangulation::limit_level_difference_at_vertices); + GridGenerator::hyper_cube(tr, -1,1); + tr.refine_global (1); + + FESystem element (FESystem(FE_Q(2), dim), 1, FE_Q(1), 1); + DoFHandler dof(tr); + dof.distribute_dofs(element); + dof.distribute_mg_dofs(); + + QGauss quadrature(element.degree+1); + + SparsityPattern sparsity (dof.n_dofs(), dof.n_dofs()); + DoFTools::make_sparsity_pattern (dof, sparsity); + sparsity.compress (); + + SparseMatrix matrix; + matrix.reinit (sparsity); + + make_stokes_matrix (dof, quadrature, matrix); + + BlockMask exclude_boundary_dofs(std::vector {true, false}); + + using Smoother = RelaxationBlockJacobi>; + { + Smoother::AdditionalData smoother_data; + Smoother smoother; + + DoFTools::make_vertex_patches(smoother_data.block_list, dof, tr.n_levels()-1, exclude_boundary_dofs); + smoother_data.block_list.compress(); + smoother_data.inversion = PreconditionBlockBase::svd; + smoother_data.threshold = 1.e-8; + + smoother.initialize(matrix, smoother_data); + smoother.log_statistics(); + } + { + Smoother::AdditionalData smoother_data; + Smoother smoother; + + DoFTools::make_vertex_patches(smoother_data.block_list, dof, tr.n_levels()-1, exclude_boundary_dofs); + smoother_data.block_list.compress(); + smoother_data.inversion = PreconditionBlockBase::svd; + smoother_data.kernel_size = 1; + + smoother.initialize(matrix, smoother_data); + smoother.log_statistics(); + } +} + + + +int main () +{ + initlog(); + + deallog.push ("1d"); + check<1> (); + deallog.pop (); + deallog.push ("2d"); + check<2> (); + deallog.pop (); + deallog.push ("3d"); + check<3> (); + deallog.pop (); +} diff --git a/tests/dofs/block_relaxation_01.output b/tests/dofs/block_relaxation_01.output new file mode 100644 index 0000000000..7918499e80 --- /dev/null +++ b/tests/dofs/block_relaxation_01.output @@ -0,0 +1,7 @@ + +DEAL:1d::PreconditionBlockBase: 1 blocks; dim ker [1:1] sigma [0.0568887:8.09268] kappa [142.255:142.255] +DEAL:1d::PreconditionBlockBase: 1 blocks; dim ker [1:1] sigma [0.0568887:8.09268] kappa [142.255:142.255] +DEAL:2d::PreconditionBlockBase: 1 blocks; dim ker [1:1] sigma [0.0870473:121.120] kappa [1391.43:1391.43] +DEAL:2d::PreconditionBlockBase: 1 blocks; dim ker [1:1] sigma [0.0870473:121.120] kappa [1391.43:1391.43] +DEAL:3d::PreconditionBlockBase: 1 blocks; dim ker [1:1] sigma [0.137508:1177.71] kappa [8564.70:8564.70] +DEAL:3d::PreconditionBlockBase: 1 blocks; dim ker [1:1] sigma [0.137508:1177.71] kappa [8564.70:8564.70]