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 (file)
index 0000000..0d9af92
--- /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.
+<br>
+(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 19db50e6d786c9c611422489cd39d39c600c22ae..0547f29e3c3e8c9ee289928a864159b656146ef8 100644 (file)
--- 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 140d6a32fa3e4ec6a8ba20e922b8dc99e950c312..f0ae55f962d4fefe47a118ec7bb389a4c41b634d 100644 (file)
--- 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 36b47305f83bef4212f428c08dd7da4cc7d6d28a..48cdcf18f004ef7426dfc1e5e7160696b68eb332 100644 (file)
--- a/include/deal.II/lac/relaxation_block.templates.h
+++ b/include/deal.II/lac/relaxation_block.templates.h
@@ -156,7 +156,10 @@ RelaxationBlock<MatrixType, InverseNumberType, VectorType>::block_kernel (const
         case PreconditionBlockBase<InverseNumberType>::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 cd23a761cd0dd59ca68e9e2511765c0a2f20154c..698cd9445452179496342616de73e829930ffe7a 100644 (file)
--- a/source/lac/lapack_full_matrix.cc
+++ b/source/lac/lapack_full_matrix.cc
@@ -1149,6 +1149,26 @@ LAPACKFullMatrix<number>::compute_inverse_svd(const double threshold)
 }
 
 
+
+template <typename number>
+void
+LAPACKFullMatrix<number>::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<n_wr-kernel_size; ++i)
+    wr[i] = number(1.)/wr[i];
+  for (size_type i=n_wr-kernel_size; i<n_wr; ++i)
+    wr[i] = 0.;
+  state = LAPACKSupport::inverse_svd;
+}
+
+
+
 template <typename number>
 void
 LAPACKFullMatrix<number>::invert()

diff --git a/tests/dofs/block_relaxation_01.cc b/tests/dofs/block_relaxation_01.cc
new file mode 100644 (file)
index 0000000..8a1f27a
--- /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 <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/function_lib.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/lac/relaxation_block.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/numerics/matrix_tools.h>
+
+
+template <int dim>
+void
+make_stokes_matrix (const DoFHandler<dim> &dof_handler,
+                    const Quadrature<dim> &quadrature_formula,
+                    SparseMatrix<double> &system_matrix)
+{
+  const FiniteElement<dim> &fe = dof_handler.get_fe();
+  const unsigned int degree = fe.degree;
+  system_matrix=0;
+
+  ConstraintMatrix constraints;
+  constraints.close();
+  FEValues<dim> 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<double>   local_matrix (dofs_per_cell, dofs_per_cell);
+  std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
+  const FEValuesExtractors::Vector velocities (0);
+  const FEValuesExtractors::Scalar pressure (dim);
+  std::vector<SymmetricTensor<2,dim> > symgrad_phi_u (dofs_per_cell);
+  std::vector<double>                  div_phi_u   (dofs_per_cell);
+  std::vector<double>                  phi_p       (dofs_per_cell);
+
+  typename DoFHandler<dim>::active_cell_iterator
+  cell = dof_handler.begin_active(),
+  endc = dof_handler.end();
+  for (; cell!=endc; ++cell)
+    {
+      fe_values.reinit (cell);
+      local_matrix = 0;
+      for (unsigned int q=0; q<n_q_points; ++q)
+        {
+          for (unsigned int k=0; k<dofs_per_cell; ++k)
+            {
+              symgrad_phi_u[k] = fe_values[velocities].symmetric_gradient (k, q);
+              div_phi_u[k]     = fe_values[velocities].divergence (k, q);
+              phi_p[k]         = fe_values[pressure].value (k, q);
+            }
+          for (unsigned int i=0; i<dofs_per_cell; ++i)
+            for (unsigned int j=0; j<dofs_per_cell; ++j)
+              local_matrix(i,j) += (2 * (symgrad_phi_u[i] * symgrad_phi_u[j])
+                                    - div_phi_u[i] * phi_p[j]
+                                    + phi_p[i] * div_phi_u[j])
+                                   * fe_values.JxW(q);
+        }
+      cell->get_dof_indices (local_dof_indices);
+      constraints.distribute_local_to_global (local_matrix, local_dof_indices,
+                                              system_matrix);
+    }
+}
+
+
+template <int dim>
+void
+check ()
+{
+  Triangulation<dim> tr(Triangulation<dim>::limit_level_difference_at_vertices);
+  GridGenerator::hyper_cube(tr, -1,1);
+  tr.refine_global (1);
+
+  FESystem<dim> element (FESystem<dim>(FE_Q<dim>(2), dim), 1, FE_Q<dim>(1), 1);
+  DoFHandler<dim> dof(tr);
+  dof.distribute_dofs(element);
+  dof.distribute_mg_dofs();
+
+  QGauss<dim> quadrature(element.degree+1);
+
+  SparsityPattern sparsity (dof.n_dofs(), dof.n_dofs());
+  DoFTools::make_sparsity_pattern (dof, sparsity);
+  sparsity.compress ();
+
+  SparseMatrix<double> matrix;
+  matrix.reinit (sparsity);
+
+  make_stokes_matrix (dof, quadrature, matrix);
+
+  BlockMask exclude_boundary_dofs(std::vector<bool> {true, false});
+
+  using Smoother = RelaxationBlockJacobi<SparseMatrix<double>>;
+  {
+    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<double>::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<double>::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 (file)
index 0000000..7918499
--- /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]