]> https://gitweb.dealii.org/ - dealii.git/commitdiff
use DiagonalMatrix in MatrixFreeOperators::Base 3340/head
authorDenis Davydov <davydden@gmail.com>
Wed, 2 Nov 2016 07:34:37 +0000 (08:34 +0100)
committerDenis Davydov <davydden@gmail.com>
Wed, 2 Nov 2016 16:01:59 +0000 (17:01 +0100)
To that end add a few missing function to DiagonalMatrix and also
switch to Jacobi preconditioner in VectorTools::project() from
quadrature points.

include/deal.II/lac/diagonal_matrix.h
include/deal.II/matrix_free/operators.h
include/deal.II/numerics/vector_tools.templates.h
tests/matrix_free/mass_operator_02.cc

index 3238b340e4622977e779d03f326a8fa443dfe887..ff9e27b8fdad8b04ab928897a6036b7a3af81314 100644 (file)
@@ -53,6 +53,11 @@ public:
    */
   VectorType &get_vector();
 
+  /**
+   * Clear content of this object and reset to the state of default constructor.
+   */
+  void clear();
+
   /**
    * Returns a read-only reference to the underlying vector.
    */
@@ -143,6 +148,11 @@ public:
   void Tvmult_add (VectorType       &dst,
                    const VectorType &src) const;
 
+  /**
+   * Return the memory consumption of this object.
+   */
+  std::size_t memory_consumption () const;
+
 private:
   /**
    * The stored vector.
@@ -162,6 +172,24 @@ DiagonalMatrix<VectorType>::DiagonalMatrix()
 
 
 
+template <typename VectorType>
+void
+DiagonalMatrix<VectorType>::clear()
+{
+  diagonal.reinit(0);
+}
+
+
+
+template <typename VectorType>
+std::size_t
+DiagonalMatrix<VectorType>::memory_consumption () const
+{
+  return diagonal.memory_consumption();
+}
+
+
+
 template <typename VectorType>
 void
 DiagonalMatrix<VectorType>::reinit(const VectorType &vec)
index 4e5a78c13252a13eb3784f59b898e294e32c3de2..3b178053382e6d285aa9bc8c919df516aac39d41 100644 (file)
@@ -22,6 +22,7 @@
 #include <deal.II/base/vectorization.h>
 #include <deal.II/base/subscriptor.h>
 
+#include <deal.II/lac/diagonal_matrix.h>
 #include <deal.II/lac/la_parallel_vector.h>
 #include <deal.II/multigrid/mg_constrained_dofs.h>
 #include <deal.II/matrix_free/matrix_free.h>
@@ -38,7 +39,7 @@ namespace MatrixFreeOperators
    * the finest mesh or at a certain level in geometric multigrid.
    *
    * A derived class has to implement apply_add() method as well as
-   * compute_diagonal() to fill the protected member inverse_diagonal_entries.
+   * compute_diagonal() to initialize the protected member inverse_diagonal_entries.
    * In case of a non-symmetric operator, Tapply_add() should be additionally
    * implemented.
    *
@@ -160,7 +161,7 @@ namespace MatrixFreeOperators
     /**
      * Get read access to the inverse diagonal of this operator.
      */
-    const LinearAlgebra::distributed::Vector<Number> &get_matrix_diagonal_inverse() const;
+    const DiagonalMatrix<LinearAlgebra::distributed::Vector<Number> > &get_matrix_diagonal_inverse() const;
 
     /**
      * Apply the Jacobi preconditioner, which multiplies every element of the
@@ -201,7 +202,7 @@ namespace MatrixFreeOperators
     /**
      * A vector to store inverse of diagonal elements.
      */
-    LinearAlgebra::distributed::Vector<Number> inverse_diagonal_entries;
+    DiagonalMatrix<LinearAlgebra::distributed::Vector<Number> > inverse_diagonal_entries;
 
   private:
 
@@ -611,7 +612,7 @@ namespace MatrixFreeOperators
   Base<dim,Number>::clear ()
   {
     data = NULL;
-    inverse_diagonal_entries.reinit(0);
+    inverse_diagonal_entries.clear();
   }
 
 
@@ -623,8 +624,8 @@ namespace MatrixFreeOperators
   {
     (void) col;
     Assert (row == col, ExcNotImplemented());
-    Assert (inverse_diagonal_entries.size() > 0, ExcNotInitialized());
-    return 1.0/inverse_diagonal_entries(row);
+    Assert (inverse_diagonal_entries.m() > 0, ExcNotInitialized());
+    return 1.0/inverse_diagonal_entries(row,row);
   }
 
 
@@ -866,10 +867,10 @@ namespace MatrixFreeOperators
 
 
   template <int dim, typename Number>
-  const LinearAlgebra::distributed::Vector<Number> &
+  const DiagonalMatrix<LinearAlgebra::distributed::Vector<Number> > &
   Base<dim,Number>::get_matrix_diagonal_inverse() const
   {
-    Assert(inverse_diagonal_entries.size() > 0, ExcNotInitialized());
+    Assert(inverse_diagonal_entries.m() > 0, ExcNotInitialized());
     return inverse_diagonal_entries;
   }
 
@@ -891,10 +892,8 @@ namespace MatrixFreeOperators
                                         const LinearAlgebra::distributed::Vector<Number> &src,
                                         const Number omega) const
   {
-    Assert(inverse_diagonal_entries.size() > 0, ExcNotInitialized());
-
-    dst = src;
-    dst.scale(inverse_diagonal_entries);
+    Assert(inverse_diagonal_entries.m() > 0, ExcNotInitialized());
+    inverse_diagonal_entries.vmult(dst,src);
     dst*= omega;
   }
 
@@ -919,24 +918,22 @@ namespace MatrixFreeOperators
     Assert((Base<dim, Number>::data != NULL), ExcNotInitialized());
 
     LinearAlgebra::distributed::Vector<Number> ones;
-    Base<dim, Number>::initialize_dof_vector(Base<dim, Number>::inverse_diagonal_entries);
-    Base<dim, Number>::initialize_dof_vector(ones);
+    LinearAlgebra::distributed::Vector<Number> &inverse_diagonal_entries = Base<dim, Number>::inverse_diagonal_entries.get_vector();
+    this->initialize_dof_vector(ones);
+    this->initialize_dof_vector(inverse_diagonal_entries);
     ones = 1.;
     ones.update_ghost_values();
-    apply_add(Base<dim, Number>::inverse_diagonal_entries, ones);
+    apply_add(inverse_diagonal_entries, ones);
 
-    const std::vector<unsigned int> &constrained_dofs
-      = Base<dim, Number>::data->get_constrained_dofs();
-    for (unsigned int i=0; i< constrained_dofs.size(); ++i)
-      Base<dim, Number>::inverse_diagonal_entries.local_element(constrained_dofs[i]) = 1.;
+    this->set_constrained_entries_to_one(inverse_diagonal_entries);
 
-    const unsigned int local_size = Base<dim, Number>::inverse_diagonal_entries.local_size();
+    const unsigned int local_size = inverse_diagonal_entries.local_size();
     for (unsigned int i=0; i<local_size; ++i)
-      Base<dim, Number>::inverse_diagonal_entries.local_element(i)
-        =1./Base<dim, Number>::inverse_diagonal_entries.local_element(i);
+      inverse_diagonal_entries.local_element(i)
+        =1./inverse_diagonal_entries.local_element(i);
 
-    Base<dim, Number>::inverse_diagonal_entries.compress(VectorOperation::insert);
-    Base<dim, Number>::inverse_diagonal_entries.update_ghost_values();
+    inverse_diagonal_entries.compress(VectorOperation::insert);
+    inverse_diagonal_entries.update_ghost_values();
   }
 
 
@@ -1028,7 +1025,7 @@ namespace MatrixFreeOperators
     Assert((Base<dim, Number>::data != NULL), ExcNotInitialized());
 
     unsigned int dummy = 0;
-    LinearAlgebra::distributed::Vector<Number> &inverse_diagonal_entries = Base<dim,Number>::inverse_diagonal_entries;
+    LinearAlgebra::distributed::Vector<Number> &inverse_diagonal_entries = Base<dim,Number>::inverse_diagonal_entries.get_vector();
     this->initialize_dof_vector(inverse_diagonal_entries);
     Base<dim,Number>::
     data->cell_loop (&LaplaceOperator::local_diagonal_cell,
@@ -1042,8 +1039,8 @@ namespace MatrixFreeOperators
       else
         inverse_diagonal_entries.local_element(i) = 1.;
 
-    Base<dim, Number>::inverse_diagonal_entries.compress(VectorOperation::insert);
-    Base<dim, Number>::inverse_diagonal_entries.update_ghost_values();
+    inverse_diagonal_entries.compress(VectorOperation::insert);
+    inverse_diagonal_entries.update_ghost_values();
   }
 
 
index 453e4aa351d7353a9b85951f2a56358910f532bf..1decd4bbfd413168972b69ce42352c2d20050ba3 100644 (file)
@@ -970,9 +970,8 @@ namespace VectorTools
       //now invert the matrix
       ReductionControl     control(rhs.size(), 0., 1e-12, false, false);
       SolverCG<LinearAlgebra::distributed::Vector<Number> > cg(control);
-      typename PreconditionChebyshev<MatrixType, LocalVectorType>::AdditionalData data;
-      data.matrix_diagonal_inverse = mass_matrix.get_matrix_diagonal_inverse();
-      PreconditionChebyshev<MatrixType, LocalVectorType> preconditioner;
+      typename PreconditionJacobi<MatrixType>::AdditionalData data(0.8);
+      PreconditionJacobi<MatrixType> preconditioner;
       preconditioner.initialize(mass_matrix, data);
       cg.solve (mass_matrix, vec, rhs, preconditioner);
       vec+=inhomogeneities;
@@ -1041,9 +1040,8 @@ namespace VectorTools
       //now invert the matrix
       ReductionControl     control(rhs.size(), 0., 1e-12, false, false);
       SolverCG<LinearAlgebra::distributed::Vector<Number> > cg(control);
-      typename PreconditionChebyshev<MatrixType, LocalVectorType>::AdditionalData data;
-      data.matrix_diagonal_inverse = mass_matrix.get_matrix_diagonal_inverse();
-      PreconditionChebyshev<MatrixType, LocalVectorType> preconditioner;
+      typename PreconditionJacobi<MatrixType>::AdditionalData data(0.8);
+      PreconditionJacobi<MatrixType> preconditioner;
       preconditioner.initialize(mass_matrix, data);
       cg.solve (mass_matrix, vec, rhs, preconditioner);
       vec+=inhomogeneities;
index 1b83bb0c1b3a18fd215443197bab5bac427d576a..4bcc628f4e67c6d0660de92671fd35f5d2d5931a 100644 (file)
@@ -83,8 +83,8 @@ void test ()
   MatrixFreeOperators::MassOperator<dim,fe_degree, 1, number> mf;
   mf.initialize(mf_data);
   mf.compute_diagonal();
-  const LinearAlgebra::distributed::Vector<double> diagonal
-    = mf.get_matrix_diagonal_inverse();
+  const LinearAlgebra::distributed::Vector<double> &diagonal
+    = mf.get_matrix_diagonal_inverse().get_vector();
 
   LinearAlgebra::distributed::Vector<number> in, out, ref;
   mf_data.initialize_dof_vector (in);

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.