]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Change the parallelization scheme: Now we use WorkStream instead of parallel_for...
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Thu, 1 Oct 2009 12:00:20 +0000 (12:00 +0000)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Thu, 1 Oct 2009 12:00:20 +0000 (12:00 +0000)
git-svn-id: https://svn.dealii.org/trunk@19635 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-37/doc/intro.dox
deal.II/examples/step-37/doc/results.dox
deal.II/examples/step-37/step-37.cc

index f8e8367b4718f494962a16d5453a741c82913ef9..456a02c1a9657f1f5b51cdc1647bc2c4e3e45520 100644 (file)
@@ -431,8 +431,12 @@ matrix-vector product implementation efficient on a GPU.
 
 For our program, we choose to follow a simple strategy: We let several
 processors work together by splitting the cells into several chunks. The
-threading building blocks implemenation of a %parallel <code>for</code> loop
-implements this concept using the parallel::apply_to_subranges() function.
+threading building blocks implemenation of a %parallel pipeline implements
+this concept using the WorkStream::run() function. What the pipeline does is
+pretty much the same as a for loop. However, it can be instructed to do some
+part of the loop body by just one process at a time and in natural order. We
+need to use this for writing the local contributions into the global vector,
+in order to avoid a race condition.
 
 
 <h3>Combination with multigrid</h3>
index a444ddf67a95b10545687d8a7319ccbdbb965c87..e8bce1379073b89ce1b6b939eea66a1a8c580649 100644 (file)
@@ -9,20 +9,20 @@ output:
 @code
 Cycle 0
 Number of degrees of freedom: 337
-System matrix memory consumption: 0.0257 MBytes.
-Multigrid objects memory consumption: 0.05071 MBytes.
+System matrix memory consumption: 0.02573 MBytes.
+Multigrid objects memory consumption: 0.05083 MBytes.
 Convergence in 10 CG iterations.
 
 Cycle 1
 Number of degrees of freedom: 1313
-System matrix memory consumption: 0.0925 MBytes.
-Multigrid objects memory consumption: 0.1792 MBytes.
+System matrix memory consumption: 0.09254 MBytes.
+Multigrid objects memory consumption: 0.1794 MBytes.
 Convergence in 10 CG iterations.
 
 Cycle 2
 Number of degrees of freedom: 5185
-System matrix memory consumption: 0.3551 MBytes.
-Multigrid objects memory consumption: 0.6775 MBytes.
+System matrix memory consumption: 0.3552 MBytes.
+Multigrid objects memory consumption: 0.6777 MBytes.
 Convergence in 10 CG iterations.
 
 Cycle 3
@@ -40,7 +40,7 @@ Convergence in 10 CG iterations.
 Cycle 5
 Number of degrees of freedom: 328193
 System matrix memory consumption: 22.1 MBytes.
-Multigrid objects memory consumption: 41.63 MBytes.
+Multigrid objects memory consumption: 41.64 MBytes.
 Convergence in 10 CG iterations.
 @endcode
 
@@ -51,14 +51,14 @@ program in three dimensions:
 @code
 Cycle 0
 Number of degrees of freedom: 517
-System matrix memory consumption: 0.1 MBytes.
-Multigrid objects memory consumption: 0.1462 MBytes.
+System matrix memory consumption: 0.1001 MBytes.
+Multigrid objects memory consumption: 0.1463 MBytes.
 Convergence in 9 CG iterations.
 
 Cycle 1
 Number of degrees of freedom: 3817
-System matrix memory consumption: 0.6612 MBytes.
-Multigrid objects memory consumption: 0.8894 MBytes.
+System matrix memory consumption: 0.6613 MBytes.
+Multigrid objects memory consumption: 0.8896 MBytes.
 Convergence in 10 CG iterations.
 
 Cycle 2
@@ -75,9 +75,8 @@ Convergence in 11 CG iterations.
 
 Cycle 4
 Number of degrees of freedom: 1847617
-System matrix memory consumption: 321.9 MBytes.
+System matrix memory consumption: 322 MBytes.
 Multigrid objects memory consumption: 415.1 MBytes.
-Convergence in 11 CG iterations.
 @endcode
 
 A comparison of the memory consumption of the MatrixFree class with a 
index 1d3a3edd8455a9ef8712ed7dbb17e4c5aea89271..c6e7d51f82b1e51d00ce1cc887aae95355f22838 100644 (file)
@@ -16,6 +16,7 @@
 #include <base/quadrature_lib.h>
 #include <base/function.h>
 #include <base/logstream.h>
+#include <base/work_stream.h>
 
 #include <lac/vector.h>
 #include <lac/full_matrix.h>
@@ -109,6 +110,57 @@ void Coefficient<dim>::value_list (const std::vector<Point<dim> > &points,
 
                                 // @sect3{Matrix-free implementation.}
 
+                               // First com a few variables that we use for
+                               // defining the %parallel layout of the vector
+                               // multiplication function with the WorkStream
+                               // concept in the Matrix-free class.
+namespace WorkStreamData
+{
+  template <typename number>
+  struct ScratchData
+  {
+    ScratchData ();
+    ScratchData (const ScratchData &scratch);
+    FullMatrix<number> solutions;
+  };
+
+  template<typename number>
+  ScratchData<number>::ScratchData ()
+    :
+    solutions ()
+  {}
+
+  template<typename number>
+  ScratchData<number>::ScratchData (const ScratchData &scratch)
+    :
+    solutions ()
+  {}
+
+  template <typename number>
+  struct CopyData : public ScratchData<number>
+  {
+    CopyData ();
+    CopyData (const CopyData &scratch);
+    unsigned int first_dof;
+    unsigned int n_dofs;
+  };
+
+  template <typename number>
+  CopyData<number>::CopyData ()
+    :
+    ScratchData<number> ()
+  {}
+
+  template <typename number>
+  CopyData<number>::CopyData (const CopyData &scratch)
+    :
+    ScratchData<number> ()
+  {}
+
+}
+
+
+
                                 // Next comes the implemenation of the
                                 // matrix-free class. It provides some
                                 // standard information we expect for
@@ -196,11 +248,18 @@ class MatrixFree : public Subscriptor
                                 // as well as a few other variables that
                                 // store matrix properties.
   private:
+    typedef std::vector<std::pair<unsigned int,unsigned int> >::const_iterator
+      CellChunkIterator;
     template <typename number2>
-    void vmult_on_subrange (const unsigned int first_cell,
-                           const unsigned int last_cell,
-                           Vector<number2> &dst,
-                           const Vector<number2> &src) const;
+    void local_vmult (CellChunkIterator                    cell_range,
+                     WorkStreamData::ScratchData<number> &scratch,
+                     WorkStreamData::CopyData<number>    &copy,
+                     const Vector<number2>               &src) const;
+
+    template <typename number2>
+    void
+    copy_local_to_global (const WorkStreamData::CopyData<number> &copy,
+                         Vector<number2>                        &dst) const;
 
     FullMatrix<number>      small_matrix;
     Table<2,unsigned int>   indices_local_to_global;
@@ -216,6 +275,7 @@ class MatrixFree : public Subscriptor
       unsigned int n_dofs, n_cells;
       unsigned int m, n;
       unsigned int n_points, n_comp;
+      std::vector<std::pair<unsigned int,unsigned int> > chunks;
     }  matrix_sizes;
 };
 
@@ -232,81 +292,6 @@ MatrixFree<number,Transformation>::MatrixFree ()
 
 
 
-                                // This function initializes the structures
-                                // of the matrix. It writes the number of
-                                // total degrees of freedom in the problem
-                                // as well as the number of cells to the
-                                // MatrixSizes struct and copies the small
-                                // matrix that transforms the solution from
-                                // support points to quadrature points. It
-                                // uses the small matrix for determining
-                                // the number of degrees of freedom per
-                                // cell (number of rows in
-                                // <code>small_matrix</code>). The number
-                                // of quadrature points needs to be passed
-                                // through the last variable
-                                // <code>n_points_per_cell</code>, since
-                                // the number of columns in the small
-                                // matrix is
-                                // <code>dim*n_points_per_cell</code> for
-                                // the Laplace problem (the Laplacian is a
-                                // tensor and has <code>dim</code>
-                                // components). In this function, we also
-                                // give the fields containing the
-                                // derivative information and the local dof
-                                // indices the correct sizes. They will be
-                                // filled by calling the respective set
-                                // function.
-template <typename number, class Transformation>
-void MatrixFree<number,Transformation>::
-reinit (const unsigned int        n_dofs_in,
-       const unsigned int        n_cells_in,
-       const FullMatrix<double> &small_matrix_in,
-       const unsigned int        n_points_per_cell)
-{
-  small_matrix = small_matrix_in;
-
-  derivatives.reinit (n_cells_in, n_points_per_cell);
-  indices_local_to_global.reinit (n_cells_in, small_matrix.m());
-
-  diagonal_is_calculated = false;
-
-  matrix_sizes.n_dofs = n_dofs_in;
-  matrix_sizes.n_cells = n_cells_in;
-  matrix_sizes.m = small_matrix.m();
-  matrix_sizes.n = small_matrix.n();
-  matrix_sizes.n_points = n_points_per_cell;
-  matrix_sizes.n_comp   = small_matrix.n()/matrix_sizes.n_points;
-  Assert(matrix_sizes.n_comp * n_points_per_cell == small_matrix.n(),
-        ExcInternalError());
-}
-
-
-
-                                // This function we need if we want to
-                                // delete the content of the matrix,
-                                // e.g. when we are finished with one grid
-                                // level and continue to the next one. Just
-                                // put all the field sizes to 0.
-template <typename number, class Transformation>
-void
-MatrixFree<number,Transformation>::clear ()
-{
-  small_matrix.reinit(0,0);
-  derivatives.reinit (0,0);
-  indices_local_to_global.reinit(0,0);
-
-  constraints.clear();
-
-  diagonal_values.reinit (0);
-  diagonal_is_calculated = false;
-
-  matrix_sizes.n_dofs = 0;
-  matrix_sizes.n_cells = 0;
-}
-
-
-
                                 // This function returns the number of rows
                                 // of the global matrix, and the next one
                                 // the number of columns (which is the
@@ -403,13 +388,14 @@ set_derivative_data (const unsigned int cell_no,
                                 // This is the central function of the
                                 // matrix-free class, implementing the
                                 // multiplication of the matrix with a
-                                // vector. This function actually not work
-                                // on all the cells, but only a subset of
-                                // cells. Since this function operates
-                                // similarly irrespective on which cell
-                                // chunk we are sitting, we can parallelize
-                                // it and get very regular operation
-                                // patterns.
+                                // vector. This function does actually not
+                                // work on all the cells, but only a subset
+                                // of cells, specified by the first argument
+                                // <code>cell_range</code>. Since this
+                                // function operates similarly irrespective
+                                // on which cell chunk we are sitting, we can
+                                // parallelize it and get very regular
+                                // operation patterns.
                                 //
                                 // Following the discussion in the
                                 // introduction, we try to work on multiple
@@ -445,65 +431,17 @@ set_derivative_data (const unsigned int cell_no,
                                 // cell for the first and the number of
                                 // quadrature points times the number of
                                 // components per point for the latter.
-                                //
-                                // One more thing to make this work
-                                // efficiently is to decide how many cells
-                                // should be included in the matrix that
-                                // contains the solution values at local
-                                // dofs for several cells. If we choose too
-                                // few cells, then the gains from using the
-                                // matrix-matrix product will not be fully
-                                // utilized (dgemm tends to provide more
-                                // efficiency the larger the matrix
-                                // dimensions get). If we choose too many,
-                                // we will firstly degrade parallelization
-                                // (which is based on some these chunks),
-                                // and secondly introduce an inefficiency
-                                // that comes from the computer
-                                // architecture: Right after the first
-                                // matrix-matrix multiplication, we
-                                // transform the solution on quadrature
-                                // points by using derivatives. Obviously,
-                                // we want to have fast access to that
-                                // data, so it should still be present in
-                                // L2 cache and not to be fetched from main
-                                // memory. The total memory usage of the
-                                // data on quadrature points should be not
-                                // more than about half the cache size of
-                                // the processor in order to be on the safe
-                                // side. Since most today's processors
-                                // provide 512 kBytes or more cache memory
-                                // per core, we choose about 250 kB as a
-                                // size. Clearly, this is an
-                                // architecture-dependent value and the
-                                // interested user can squeeze out some
-                                // extra performance by hand-tuning this
-                                // parameter. Once we have chosen the
-                                // number of cells we collect in one chunk,
-                                // we determine how many chunks we have on
-                                // the given cell range and recalculate the
-                                // actual chunk size in order to evenly
-                                // distribute the chunks.
 template <typename number, class Transformation>
 template <typename number2>
 void
 MatrixFree<number,Transformation>::
-vmult_on_subrange (const unsigned int    first_cell,
-                  const unsigned int    last_cell,
-                  Vector<number2>       &dst,
-                  const Vector<number2> &src) const
+  local_vmult (CellChunkIterator                    cell_range,
+              WorkStreamData::ScratchData<number> &scratch,
+              WorkStreamData::CopyData<number>    &copy,
+              const Vector<number2>               &src) const
 {
-  FullMatrix<number> solution_cells, solution_points;
-
-  const unsigned int divisor = 250000/(matrix_sizes.n*sizeof(number));
-  const unsigned int n_chunks = (last_cell-first_cell)/divisor + 1;
-  const unsigned int chunk_size =
-    (last_cell-first_cell)/n_chunks + ((last_cell-first_cell)%n_chunks>0);
-
-  for (unsigned int k=first_cell; k<last_cell; k+=chunk_size)
-    {
-      const unsigned int current_chunk_size =
-       k+chunk_size>last_cell ? last_cell-k : chunk_size;
+  const unsigned int first_cell = cell_range->first,
+    chunk_size = cell_range->second - cell_range->first;
 
                                 // OK, now we are sitting in the loop that
                                 // goes over our chunks of cells. What we
@@ -559,33 +497,49 @@ vmult_on_subrange (const unsigned int    first_cell,
                                 // are related to each other. Since we
                                 // simultaneously apply the constraints, we
                                 // hand this task off to the ConstraintMatrix
-                                // object. Most often, itis used to work on
-                                // one cell at a time, but since we work on a
-                                // whole chunk of dofs, we can do that just
-                                // as easily for all the cells at once.
-      solution_cells.reinit (current_chunk_size,matrix_sizes.m, true);
-      solution_points.reinit (current_chunk_size,matrix_sizes.n, true);
-
-      const unsigned int n_cell_entries = current_chunk_size*matrix_sizes.m;
-      constraints.get_dof_values(src, &indices_local_to_global(k,0),
-                                &solution_cells(0,0),
-                                &solution_cells(0,0)+n_cell_entries);
-
-      solution_cells.mmult (solution_points, small_matrix);
-
-      for (unsigned int i=0; i<current_chunk_size; ++i)
-       for (unsigned int j=0; j<matrix_sizes.n_points; ++j)
-         derivatives(i+k,j).transform(&solution_points(i, j*matrix_sizes.n_comp));
-
-      solution_points.mTmult (solution_cells, small_matrix);
-
-      static Threads::Mutex mutex;
-      Threads::Mutex::ScopedLock lock (mutex);
-      constraints.distribute_local_to_global (&solution_cells(0,0),
-                                             &solution_cells(0,0)+n_cell_entries,
-                                             &indices_local_to_global(k,0),
-                                             dst);
-    }
+                                // object. We do this in an extra function
+                                // since we split between %parallel code that
+                                // can be run independently (this function)
+                                // and code that needs to be synchronized
+                                // between threads
+                                // (<code>copy_local_to_global</code>
+                                // function). Most often, the
+                                // ConstraintMatrix function is used to be
+                                // applied to data from one cell at a time,
+                                // but since we work on a whole chunk of
+                                // dofs, we can do that just as easily for
+                                // all the cells at once.
+  copy.solutions.reinit    (chunk_size,matrix_sizes.m, true);
+  copy.first_dof          = first_cell;
+  copy.n_dofs             = chunk_size*matrix_sizes.m;
+  scratch.solutions.reinit (chunk_size,matrix_sizes.n, true);
+
+  constraints.get_dof_values(src, &indices_local_to_global(copy.first_dof,0),
+                            &copy.solutions(0,0),
+                            &copy.solutions(0,0)+copy.n_dofs);
+
+  copy.solutions.mmult (scratch.solutions, small_matrix);
+
+  for (unsigned int i=0, k = first_cell; i<chunk_size; ++i, ++k)
+    for (unsigned int j=0; j<matrix_sizes.n_points; ++j)
+      derivatives(k,j).transform(&scratch.solutions(i, j*matrix_sizes.n_comp));
+
+  scratch.solutions.mTmult (copy.solutions, small_matrix);
+}
+
+
+
+template <typename number, class Transformation>
+template <typename number2>
+void
+MatrixFree<number,Transformation>::
+  copy_local_to_global (const WorkStreamData::CopyData<number> &copy,
+                       Vector<number2>                        &dst) const
+{
+  constraints.distribute_local_to_global (&copy.solutions(0,0),
+                                         &copy.solutions(0,0)+copy.n_dofs,
+                                         &indices_local_to_global(copy.first_dof,0),
+                                         dst);
 }
 
 
@@ -640,27 +594,29 @@ MatrixFree<number,Transformation>::Tvmult_add (Vector<number2>       &dst,
                                 // vector <code>src</code> and adds the
                                 // result to vector <code>dst</code>. We call
                                 // a %parallel function that applies the
-                                // multiplication on a subrange of cells
-                                // (cf. the @ref threads module).
-                                //
-                                // TODO: Use WorkStream for parallelization
-                                // instead of apply_to_subranges, once we
-                                // have realized the best way for doing
-                                // that.
+                                // multiplication on a chunk of cells at once
+                                // using the WorkStream module (cf. also the
+                                // @ref threads module). The subdivision into
+                                // chunks will be performed in the reinit
+                                // function and is stored in the field
+                                // <code>matrix_sizes.chunks</code>.
 template <typename number, class Transformation>
 template <typename number2>
 void
 MatrixFree<number,Transformation>::vmult_add (Vector<number2>       &dst,
                                              const Vector<number2> &src) const
 {
-  parallel::apply_to_subranges (0, matrix_sizes.n_cells,
-                               std_cxx1x::bind(&MatrixFree<number,Transformation>::
-                                               template vmult_on_subrange<number2>,
-                                               this,
-                                               _1,_2,
-                                               boost::ref(dst),
-                                               boost::cref(src)),
-                               200);
+
+  WorkStream::run (matrix_sizes.chunks.begin(), matrix_sizes.chunks.end(),
+                  std_cxx1x::bind(&MatrixFree<number,Transformation>::
+                                  template local_vmult<number2>,
+                                  this, _1, _2, _3, boost::cref(src)),
+                  std_cxx1x::bind(&MatrixFree<number,Transformation>::
+                                  template copy_local_to_global<number2>,
+                                  this, _1, boost::ref(dst)),
+                  WorkStreamData::ScratchData<number>(),
+                  WorkStreamData::CopyData<number>(),
+                  4,1);
 
                                 // One thing to be cautious about: The
                                 // deal.II classes expect that the matrix
@@ -683,6 +639,142 @@ MatrixFree<number,Transformation>::vmult_add (Vector<number2>       &dst,
 
 
 
+                                // This function initializes the structures
+                                // of the matrix. It writes the number of
+                                // total degrees of freedom in the problem
+                                // as well as the number of cells to the
+                                // MatrixSizes struct and copies the small
+                                // matrix that transforms the solution from
+                                // support points to quadrature points. It
+                                // uses the small matrix for determining
+                                // the number of degrees of freedom per
+                                // cell (number of rows in
+                                // <code>small_matrix</code>). The number
+                                // of quadrature points needs to be passed
+                                // through the last variable
+                                // <code>n_points_per_cell</code>, since
+                                // the number of columns in the small
+                                // matrix is
+                                // <code>dim*n_points_per_cell</code> for
+                                // the Laplace problem (the Laplacian is a
+                                // tensor and has <code>dim</code>
+                                // components). In this function, we also
+                                // give the fields containing the
+                                // derivative information and the local dof
+                                // indices the correct sizes. They will be
+                                // filled by calling the respective set
+                                // function.
+template <typename number, class Transformation>
+void MatrixFree<number,Transformation>::
+reinit (const unsigned int        n_dofs_in,
+       const unsigned int        n_cells_in,
+       const FullMatrix<double> &small_matrix_in,
+       const unsigned int        n_points_per_cell)
+{
+  small_matrix = small_matrix_in;
+
+  derivatives.reinit (n_cells_in, n_points_per_cell);
+  indices_local_to_global.reinit (n_cells_in, small_matrix.m());
+
+  diagonal_is_calculated = false;
+
+  matrix_sizes.n_dofs = n_dofs_in;
+  matrix_sizes.n_cells = n_cells_in;
+  matrix_sizes.m = small_matrix.m();
+  matrix_sizes.n = small_matrix.n();
+  matrix_sizes.n_points = n_points_per_cell;
+  matrix_sizes.n_comp   = small_matrix.n()/matrix_sizes.n_points;
+  Assert(matrix_sizes.n_comp * n_points_per_cell == small_matrix.n(),
+        ExcInternalError());
+
+                                // One thing to make the matrix-vector
+                                // product with this class efficient is to
+                                // decide how many cells should be summarized
+                                // to one chunk, which will determine the
+                                // size of the full matrix that we work
+                                // on. If we choose too few cells, then the
+                                // gains from using the matrix-matrix product
+                                // will not be fully utilized (dgemm tends to
+                                // provide more efficiency the larger the
+                                // matrix dimensions get). If we choose too
+                                // many, we will firstly degrade
+                                // parallelization (which is based on some
+                                // these chunks), and secondly introduce an
+                                // inefficiency that comes from the computer
+                                // architecture: In the actual working
+                                // function above, right after the first
+                                // matrix-matrix multiplication, we transform
+                                // the solution on quadrature points by using
+                                // derivatives. Obviously, we want to have
+                                // fast access to that data, so it should
+                                // still be present in L2 cache and not
+                                // needed to be fetched from main memory. The
+                                // total memory usage of the data on
+                                // quadrature points should be not more than
+                                // about half the cache size of the processor
+                                // in order to be on the safe side. Since
+                                // most today's processors provide 512 kBytes
+                                // or more cache memory per core, we choose
+                                // about 50 kB as a size to be on the safe
+                                // side (other things need to be stored in
+                                // the CPU as well). Clearly, this is an
+                                // architecture-dependent value and the
+                                // interested user can squeeze out some extra
+                                // performance by hand-tuning this
+                                // parameter. Once we have chosen the number
+                                // of cells we collect in one chunk, we
+                                // determine how many chunks we have on the
+                                // given cell range and recalculate the
+                                // actual chunk size in order to evenly
+                                // distribute the chunks.
+  const unsigned int divisor = 50000/(matrix_sizes.m*sizeof(double));
+  unsigned int n_chunks = matrix_sizes.n_cells/divisor + 1;
+  if (n_chunks<2*multithread_info.n_default_threads)
+    n_chunks = 2*multithread_info.n_default_threads;
+
+  const unsigned int chunk_size = (matrix_sizes.n_cells/n_chunks +
+                                  (matrix_sizes.n_cells%n_chunks>0));
+  matrix_sizes.chunks.resize (n_chunks);
+  for (unsigned int i=0; i<n_chunks; ++i)
+    {
+      matrix_sizes.chunks[i].first = i*chunk_size;
+      if ((i+1)*chunk_size > matrix_sizes.n_cells)
+       {
+         matrix_sizes.chunks[i].second = matrix_sizes.n_cells;
+         break;
+       }
+      else
+       matrix_sizes.chunks[i].second = (i+1)*chunk_size;
+    }
+}
+
+
+
+                                // This function we need if we want to
+                                // delete the content of the matrix,
+                                // e.g. when we are finished with one grid
+                                // level and continue to the next one. Just
+                                // put all the field sizes to 0.
+template <typename number, class Transformation>
+void
+MatrixFree<number,Transformation>::clear ()
+{
+  small_matrix.reinit(0,0);
+  derivatives.reinit (0,0);
+  indices_local_to_global.reinit(0,0);
+
+  constraints.clear();
+
+  diagonal_values.reinit (0);
+  diagonal_is_calculated = false;
+
+  matrix_sizes.n_dofs = 0;
+  matrix_sizes.n_cells = 0;
+  matrix_sizes.chunks.clear();
+}
+
+
+
                                 // This function returns the entries of the
                                 // matrix. Since this class is intended not
                                 // to store the matrix entries, it would
@@ -778,7 +870,9 @@ std::size_t MatrixFree<number,Transformation>::memory_consumption () const
     indices_local_to_global.memory_consumption() +
     constraints.memory_consumption() +
     small_matrix.memory_consumption() +
-    diagonal_values.memory_consumption() + sizeof(*this);
+    diagonal_values.memory_consumption() +
+    matrix_sizes.chunks.size()*2*sizeof(unsigned int) +
+    sizeof(*this);
   return glob_size;
 }
 

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.