]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Final update of comments.
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Fri, 2 Oct 2009 15:21:07 +0000 (15:21 +0000)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Fri, 2 Oct 2009 15:21:07 +0000 (15:21 +0000)
git-svn-id: https://svn.dealii.org/trunk@19668 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 456a02c1a9657f1f5b51cdc1647bc2c4e3e45520..a41f5fc6d386cfb5bbded71b81561bd1f56b346f 100644 (file)
@@ -26,8 +26,8 @@ of the equation to be solved, between 40 and 95 per cent of the computational
 time is spent in performing sparse matrix-vector products.
 
 Let us briefly look at a simplified version of code for the matrix-vector
-product (the actual implementation in deal.II uses a different counter for
-the innermost loop, which avoids having one additional counter variable):
+product (the actual implementation in deal.II uses a different layout for
+the innermost loop, which avoids having the counter variable):
 @code
 // y = A * x
 std::size_t element_index = 0;
@@ -47,47 +47,47 @@ for (unsigned int row=0; row<n_rows; ++row)
 
 Looking into what actually happens in the computer when forming a
 matrix-vector product, one observes that the matrix data
-<code>matrix_values</code> is continuously traveling from main memory into
-the CPU in order to be multiplied with some vector element determined by the
-other array <code>column_indices</code> and add this product to a sum that
+<code>matrix_values</code> is continuously traveling from main memory into the
+CPU in order to be multiplied with some vector element determined by the other
+array <code>column_indices</code> and add this product to a sum that
 eventually is written into the output vector. Let us assume for simplicity
-that all the vector elements are sitting in the CPU (or some fast memory
-like caches) and that we use a compressed storage scheme for the sparse
-matrix, i.e., the format deal.II matrices (as well as PETSc and Trilinos
-matrices) use. Then each matrix element corresponds to 12 bytes of data, 8
-bytes for the respective element in <code>matrix_values</code>, and 4
-bytes for the unsigned integer position <code>column_indices</code> that
-tells which vector element we actually use for multiplication. Here we
-neglect the additional array that tells us the ranges of individual rows in
-the matrix. With that 12 bytes of data, we perform two floating point
-operations, a multiplication and an addition. If our matrix has one billion
-entries, a matrix-vector product consists of two billion floating point
-operations, 2 GFLOP. One core of a processor of 2009's standard (Intel's
-'Nehalem' processor, 3 GHz) has a peak performance of about 12 billion
-operations per second, 12 GFLOP/s. Now we might be tempted to hope for
-getting a matrix-vector product done in about one sixth of a second on such
-a machine. Remember, though, that we have to get 12 GB of data through the
-processor in order to form the matrix-vector product. Looking again at which
-hardware is available in 2009, we will hardly get more than 10 GB/s of data
-read. This means that the matrix-vector product will take more than one
-second to complete, giving a rate of 1.7 GFLOP/s at the best. This is quite
-far away from the theoretical peak performance of 12 GFLOP/s.
+that all the vector elements are sitting in the CPU (or some fast memory like
+caches) and that we use a compressed storage scheme for the sparse matrix,
+i.e., the format deal.II matrices (as well as PETSc and Trilinos matrices)
+use. Then each matrix element corresponds to 12 bytes of data, 8 bytes for the
+respective element in <code>matrix_values</code>, and 4 bytes for the unsigned
+integer position <code>column_indices</code> that tells which vector element
+we actually use for multiplication. Here we neglect the additional array that
+tells us the ranges of individual rows in the matrix. With that 12 bytes of
+data, we perform two floating point operations, a multiplication and an
+addition. If our matrix has one billion entries, a matrix-vector product
+consists of two billion floating point operations, 2 GFLOP. One core of a
+processor of 2009's standard (Intel's 'Nehalem' processor, 3 GHz) has a peak
+performance of about 12 billion operations per second, 12 GFLOP/s. Now we
+might be tempted to hope for getting a matrix-vector product done in about one
+sixth of a second on such a machine. Remember, though, that we have to get 12
+GB of data through the processor in order to form the matrix-vector
+product. Looking again at which hardware is available in 2009, we will hardly
+get more than 10 GB/s of data read. This means that the matrix-vector product
+will take more than one second to complete, giving a rate of 1.7 GFLOP/s at
+the best. This is quite far away from the theoretical peak performance of 12
+GFLOP/s. In practice, the rate is often considerably lower because, for
+example, the vectors do not fit into cache. A usual value on a 2009's machine
+is 0.5 to 1.1 GFLOP/s.
 
 What makes things worse is that today's processors have multiple cores, and
 multiple cores have to compete for memory bandwidth. Imagine we have 8 cores
-available with a theoretical peak performance of 96 GFLOP/s. However, they
-can only get about 35 GB/s of data. For our matrix-vector product, we would
-get a performance of about 6 GFLOP/s, which is at nightmarish 6 per cent of
-the system's peak performance. In practice, things are likely even worse
-since, for example, a substantial amount of vector entries needs to be
-fetched from main memory, too.
-
-Things won't get better in the future, rather worse: Memory bandwidth
-will most likely continue to grow more slowly than the number of cores
-(i.e., the computing power). Consequently, we will probably not see
-that much of an improvement in the speed of our programs even though
-computers do become faster in the future if we accumulate their speed
-over all of their cores.
+available with a theoretical peak performance of 96 GFLOP/s. However, these
+cores sit on a machine with about 35 GB/s of memory bandwidth. For our
+matrix-vector product, we would get a performance of about 6 GFLOP/s, which is
+at nightmarish 6 per cent of the system's peak performance. And this is the
+theoretical maximum we can get!
+
+Things won't get better in the future, rather worse: Memory bandwidth will
+most likely continue to grow more slowly than the number of cores (i.e., the
+computing power). Consequently, we will probably not see that much of an
+improvement in the speed of our programs even though computers do become
+faster if we accumulate their speed over all of their cores.
 
 In essence, one may ask how this can be avoided. The basic fact that
 precipitates this is that matrices just consume a lot of memory &mdash;
@@ -109,7 +109,7 @@ In order to find out how we can write a code that performs a matrix-vector
 product, but does not need to store the matrix elements, let us start at
 looking how some finite-element related matrix $A$ is assembled:
 @f{eqnarray*}
-A = \sum_{\mathrm{cell}=1}^{\mathrm{n,cells}} P_\mathrm{cell,{loc-glob}}^T A_\mathrm{cell}
+A = \sum_{\mathrm{cell}=1}^{\mathrm{n\_cells}} P_\mathrm{cell,{loc-glob}}^T A_\mathrm{cell}
 P_\mathrm{cell,{loc-glob}}.
 @f}
 In this formula, the matrix $P_\mathrm{cell,{loc-glob}}$ is a permutation
@@ -121,13 +121,15 @@ assembly of matrices.
 
 If we are to perform a matrix-vector product, we can hence use that
 @f{eqnarray*}
-y &=& A\cdot x = \left(\sum_{\text{cell}=1}^{\mathrm{n,cells}} P_\mathrm{cell,{loc-glob}}^T
+y &=& A\cdot x = \left(\sum_{\text{cell}=1}^{\mathrm{n\_cells}} P_\mathrm{cell,{loc-glob}}^T
 A_\mathrm{cell} P_\mathrm{cell,{loc-glob}}\right) \cdot x\\
-&=& \sum_{\mathrm{cell}=1}^{\mathrm{n,cells}} P_\mathrm{cell,{loc-glob}}^T
+&=& \sum_{\mathrm{cell}=1}^{\mathrm{n\_cells}} P_\mathrm{cell,{loc-glob}}^T
 A_\mathrm{cell} x_\mathrm{cell},
+&=& \sum_{\mathrm{cell}=1}^{\mathrm{n\_cells}} P_\mathrm{cell,{loc-glob}}^T
+y_\mathrm{cell},
 @f}
 where $x_\mathrm{cell}$ is the values of <i>x</i> at the degrees of freedom
-of the respective cell.
+of the respective cell, and $y_\mathrm{cell}$ does so for the result.
 
 A naive attempt to implement the local action of the Laplacian would hence be
 to use the following code:
@@ -160,12 +162,10 @@ MatrixFree::vmult (Vector<double>       &dst,
 
       for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
         for (unsigned int i=0; i<dofs_per_cell; ++i)
-          {
-            for (unsigned int j=0; j<dofs_per_cell; ++j)
-              cell_matrix(i,j) += (fe_values.shape_grad(i,q_point) *
-                                   fe_values.shape_grad(j,q_point) *
-                                   fe_values.JxW(q_point));
-          }
+         for (unsigned int j=0; j<dofs_per_cell; ++j)
+            cell_matrix(i,j) += (fe_values.shape_grad(i,q_point) *
+                                 fe_values.shape_grad(j,q_point) *
+                                 fe_values.JxW(q_point));
 
       cell->get_dof_indices (local_dof_indices);
 
@@ -363,15 +363,17 @@ matrix from unit to real cell inverse_jacobian, because the transformation
 direction in deal.II is from real to unit cell).
 
 To sum up, we want to look at the additional costs introduced by this
-written-out matrix-matrix-vector product compared to a sparse matrix. To
-first approximation, we have increased the number of operations for the
-local matrix-vector product by a factor of 4 in 2D and 6 in 3D, i.e., two
+written-out matrix-matrix-vector product compared to a sparse matrix. To first
+approximation, we have increased the number of operations for the local
+matrix-vector product by a factor of 4 in 2D and 6 in 3D, i.e., two
 matrix-vector products with matrices that have <code>dim</code> as many
-columns than before. Then, we also need to think that we touch some degrees
-of freedom several times because they belong to several cells. This also
-increases computational costs. A realistic value compared to a sparse matrix
-is that we now have to perform about 10 times as many operations (a bit less
-in 2D, a bit more in 3D).
+columns than before (here we assume that the number of quadrature points is
+the same as the number of degrees of freedom per cell, which is usual for
+scalar problems). Then, we also need to keep in mind that we touch some
+degrees of freedom several times because they belong to several cells. This
+also increases computational costs. A realistic value compared to a sparse
+matrix is that we now have to perform about 10 times as many operations (a bit
+less in 2D, a bit more in 3D).
 
 The above is, in essence, what happens in the code below and if you have
 difficulties in understanding the implementation, you should try to first
@@ -392,36 +394,36 @@ there are a few more points done to be even more efficient, namely:
       is nice is that matrix-matrix products are close to the processor's
       peak performance if the matrices are neither too small nor too large
       &mdash; and these operations are the most expensive part in the
-      implementation shown here. This implies that the matrix-free
-      matrix-vector product is slower than a matrix-vector product using a
-      sparse matrix for linear and quadratic elements, but on par with third
-      order elements and faster for even higher order elements.
+      implementation shown here. 
 </ul>
 
-An additional gain with this implementation is that we do not have to build
-the sparse matrix itself, which can be quite expensive depending on the
-underlying differential equation.
+The implementation of the matrix-free matrix-vector product shown in this
+tutorial is slower than a matrix-vector product using a sparse matrix for
+linear and quadratic elements, but on par with third order elements and faster
+for even higher order elements. An additional gain with this implementation is
+that we do not have to build the sparse matrix itself, which can also be quite
+expensive depending on the underlying differential equation.
 
 
 <h4>Parallelization issues</h4>
 
 We mentioned in the philosophical section above that parallelization with
 multiple cores in a shared memory machine is an issue where sparse
-matrix-vector products are not particularly well suited because processors
-are memory bandwidth limited. There is a lot of data traffic involved, and
-the access patterns in the vector are not very regular. Also, different rows
-might have different %numbers of nonzero elements. The matrix-free
-implementation, however, is more favorable in this respect. It does not need
-to save all the elements (only the product of transposed Jacobian, weights,
-and Jacobian, for all quadrature points on all cells, which is about 4 times
-the size of the solution vector in 2D and 9 times the size of the solution
-vector in 3D), whereas the number of nonzeros grows with the element
-order. Moreover, most of the work is done on a very regular pattern with
-stride-one access to data: Perform matrix-vector products with the same
-matrix, perform (equally many) transformations on the vector related
-quadrature points, and do one more matrix-vector product. Only the read
-operation from the global vector <code>src</code> and the write operation to
-<code>dst</code> in the end need to have random access to a vector.
+matrix-vector products are not particularly well suited because processors are
+memory bandwidth limited. There is a lot of data traffic involved, and the
+access patterns in the vector are not very regular. Also, different rows might
+have different %numbers of nonzero elements. The matrix-free implementation,
+however, is more favorable in this respect. It does not need to save all the
+elements (only the product of transposed Jacobian, weights, and Jacobian, for
+all quadrature points on all cells, which is about 4 times the size of the
+solution vector in 2D and 9 times the size of the solution vector in 3D),
+whereas the number of nonzeros grows with the element order. Moreover, most of
+the work is done on a very regular pattern with stride-one access to data:
+Performing matrix-vector products with the same matrix, performing (equally
+many) transformations on the vector related quadrature points, and doing one
+more matrix-vector product. Only the read operation from the global vector
+<code>src</code> and the write operation to <code>dst</code> in the end
+request for random access to a vector.
 
 For example, it should not be too difficult to implement a matrix-free
 matrix-vector product on a <a
@@ -429,14 +431,14 @@ href="http://en.wikipedia.org/wiki/GPGPU">graphics processing unit
 (GP-GPU)</a>. However, it would be quite complex to make a sparse
 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 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.
+For our program, we choose to follow a simple strategy to make the code
+%parallel: We let several processors work together by splitting the cells into
+several chunks. The threading building blocks implemenation of a %parallel
+pipeline implements this concept using the WorkStream::run() function. What
+the pipeline does closely resembles the work done by 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>
@@ -455,27 +457,27 @@ system matrix.
 
 One solution is to use multigrid methods as shown in @ref step_16
 "step-16". They are known to be very fast, and they are suitable for our
-purpose since they can be designed based purely on matrix-vector products
-(and knowledge of the diagonal entries of the matrix, though). All one needs
-to do is to find a smoother that works with matrix-vector products only. One
-such candidate would be a damped Jacobi iteration, but that one is often not
-sufficiently good in damping high-frequency errors than one would like. A
-Chebyshev preconditioner, eventually, is what we use here. It can be seen as
-an extension of the Jacobi method by using Chebyshev polynomials. With
-degree zero, the Jacobi method is retrieved, whereas higher order
-corrections improve the smoothing properties if some parameters are suitably
-chosen. The effectiveness of Chebyshev smoothing in multigrid has been
-demonstrated, e.g., in the article <i>M. Adams, M. Brezina, J. Hu,
-R. Tuminaro. Parallel multigrid smoothers: polynomial versus
+purpose since they can be designed based purely on matrix-vector products. All
+one needs to do is to find a smoother that works with matrix-vector products
+only (our choice requires knowledge of the diagonal entries of the matrix,
+though). One such candidate would be a damped Jacobi iteration, but that is
+often not sufficiently good in damping high-frequency errors than one would
+like. A Chebyshev preconditioner, eventually, is what we use here. It can be
+seen as an extension of the Jacobi method by using Chebyshev polynomials. With
+degree zero, the Jacobi method with optimal damping parameter is retrieved,
+whereas higher order corrections improve the smoothing properties if some
+parameters are suitably chosen. The effectiveness of Chebyshev smoothing in
+multigrid has been demonstrated, e.g., in the article <i>M. Adams, M. Brezina,
+J. Hu, R. Tuminaro. Parallel multigrid smoothers: polynomial versus
 Gauss&ndash;Seidel, J. Comput. Phys. 188:593&ndash;610, 2003</i>. This
-publication also identifies one more advantage of Chebyshev smoothers that
-we exploit here, namely that they are easy to parallelize, whereas
-SOR/Gauss&ndash;Seidel smoothing relies on substitutions, which can often
-only be parallelized by working on diagonal subblocks of the matrix, which
+publication also identifies one more advantage of Chebyshev smoothers that we
+exploit here, namely that they are easy to parallelize, whereas
+SOR/Gauss&ndash;Seidel smoothing relies on substitutions, which can often only
+be parallelized by working on diagonal subblocks of the matrix, which
 decreases efficiency.
 
 The implementation into the multigrid framework is then very
-straight-forward. We will only need some few modifications compared to @ref
+straight-forward. We will only need some minor modifications compared to @ref
 step_16 "step-16".
 
 <h3>The test case</h3>
index e8bce1379073b89ce1b6b939eea66a1a8c580649..e67a97ec4810e3d8b0ea1fe3cdb90f232d4b1bcf 100644 (file)
@@ -44,9 +44,9 @@ Multigrid objects memory consumption: 41.64 MBytes.
 Convergence in 10 CG iterations.
 @endcode
 
-As in step-16, we see that the number of CG iterations stays constant with 
-increasing number of degrees of freedom. Not much changes if we run the 
-program in three dimensions:
+As in step-16, we see that the number of CG iterations remains constant with
+increasing number of degrees of freedom. Not much changes if we run the
+program in three spatial dimensions:
 
 @code
 Cycle 0
@@ -84,4 +84,4 @@ standard SparseMatrix object shows that, for the largest problem considered
 in 3D, the memory consumption for SparseMatrix and its sparsity pattern 
 would be about 1.3 GBytes, which is a factor of 4 more memory. If combined
 with cubic elements in 3D, the MatrixFree class would decrease memory 
-consumption by a factor of 12.
+consumption even by a factor of 12.
index 4843001592df22941a20e55a96f9797e527247d4..92da263056feefd93059f252b78964e0ad7fd1fb 100644 (file)
@@ -23,9 +23,6 @@
 #include <lac/solver_cg.h>
 #include <lac/precondition.h>
 
-#include <dofs/dof_handler.h>
-#include <dofs/dof_accessor.h>
-
 #include <fe/fe_q.h>
 #include <fe/fe_values.h>
 
@@ -56,15 +53,14 @@ using namespace dealii;
 
                                 // @sect3{Equation data.}
 
-                                // We define a variable coefficient
-                                // function for the Poisson problem. It is
-                                // similar to the function in step-5. As a
-                                // difference, we use the formulation
-                                // $\frac{1}{0.1 + \|\bf x\|^2}$ instead of
-                                // a discontinuous one. It is merely to
-                                // demonstrate the possibilities of this
-                                // implemenation, rather than being
-                                // physically reasonable.
+                                // We define a variable coefficient function
+                                // for the Poisson problem. It is similar to
+                                // the function in step-5. As a difference,
+                                // we use the formulation $\frac{1}{0.1 +
+                                // \|\bf x\|^2}$ instead of a discontinuous
+                                // one. It is merely to demonstrate the
+                                // possibilities of this implemenation,
+                                // rather than making much sense physically.
 template <int dim>
 class Coefficient : public Function<dim>
 {
@@ -110,10 +106,21 @@ 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.
+                               // First com a few declarations that we use
+                               // for defining the %parallel layout of the
+                               // vector multiplication function with the
+                               // WorkStream concept in the Matrix-free
+                               // class. These comprise so-called scratch
+                               // data that we use for calculating
+                               // cell-related information, and copy data
+                               // that is eventually used in a separate
+                               // function for writing local data into the
+                               // global vector. The reason for this split-up
+                               // definition is that many threads at a time
+                               // can execute the local multiplications (and
+                               // filling up the copy data), but than that
+                               // copy data needs to be worked on by one
+                               // process at a time.
 namespace WorkStreamData
 {
   template <typename number>
@@ -141,7 +148,7 @@ namespace WorkStreamData
   {
     CopyData ();
     CopyData (const CopyData &scratch);
-    unsigned int first_dof;
+    unsigned int first_cell;
     unsigned int n_dofs;
   };
 
@@ -173,24 +180,23 @@ namespace WorkStreamData
                                 // We choose to make this class generic,
                                 // i.e., we do not implement the actual
                                 // differential operator (here: Laplace
-                                // operator) directly in this class. What
-                                // we do is to let the actual
-                                // transformation (which happens on the
-                                // level of quadrature points, see the
-                                // discussion in the introduction) be a
-                                // template parameter that is implemented
-                                // by another class. We then only have to
-                                // store a list of these transformation for
-                                // each quadrature point on each cell in a
-                                // big list &ndash; we choose a
-                                // <code>Table<2,Transformation></code>
-                                // data format) &ndash; and call a
-                                // transform command of the
-                                // <code>Transformation</code> class. This
-                                // template magic makes it easy to reuse
-                                // this MatrixFree class for other problems
-                                // that are based on a symmetric operation
-                                // without the need for further changes.
+                                // operator) directly in this class.  We
+                                // instead let the actual transformation
+                                // (which happens on the level of quadrature
+                                // points, see the discussion in the
+                                // introduction) be a template parameter that
+                                // is implemented by another class. We then
+                                // only have to store a list of these objects
+                                // for each quadrature point on each cell in
+                                // a big list &ndash; we choose a
+                                // <code>Table<2,Transformation></code> data
+                                // format) &ndash; and call a transform
+                                // command of the <code>Transformation</code>
+                                // class. This template magic makes it easy
+                                // to reuse this MatrixFree class for other
+                                // problems that are based on a symmetric
+                                // operation without the need for further
+                                // changes.
 template <typename number, class Transformation>
 class MatrixFree : public Subscriptor
 {
@@ -282,8 +288,12 @@ class MatrixFree : public Subscriptor
 
 
                                 // This is the constructor of the
-                                // <code>MatrixFree</code> class. It does
-                                // nothing.
+                                // <code>MatrixFree</code> class. All it does
+                                // is to subscribe to the general deal.II
+                                // <code>Subscriptor</code> scheme that makes
+                                // sure that we do not delete an object of
+                                // this class as long as it used somewhere
+                                // else, e.g. in a preconditioner.
 template <typename number, class Transformation>
 MatrixFree<number,Transformation>::MatrixFree ()
     :
@@ -329,21 +339,20 @@ MatrixFree<number,Transformation>::get_constraints ()
 
 
 
-                                // This function takes a vector of local
-                                // dof indices on cell level and writes the
-                                // data into the
-                                // <code>indices_local_to_global</code>
-                                // field in order to have fast access to
-                                // it. It performs a few sanity checks like
-                                // whether the sizes in the matrix are set
+                                // This function takes a vector of local dof
+                                // indices on cell level and writes the data
+                                // into the
+                                // <code>indices_local_to_global</code> field
+                                // in order to have fast access to it. It
+                                // performs a few sanity checks like whether
+                                // the sizes in the matrix are set
                                 // correctly. One tiny thing: Whenever we
-                                // enter this function, we probably make
-                                // some modification to the matrix. This
-                                // means that the diagonal of the matrix,
-                                // which we might compute to have access to
-                                // the matrix diagonal, is invalidated. We
-                                // set the respective flag to
-                                // <code>false</code>.
+                                // enter this function, we probably make some
+                                // modification to the matrix. This means
+                                // that the diagonal of the matrix, which we
+                                // might have computed to have access to the
+                                // matrix diagonal, is invalidated. We set
+                                // the respective flag to <code>false</code>.
 template <typename number, class Transformation>
 void MatrixFree<number,Transformation>::
 set_local_dof_indices (const unsigned int               cell_no,
@@ -366,12 +375,12 @@ set_local_dof_indices (const unsigned int               cell_no,
                                 // certain cell and a certain quadrature
                                 // point to the array that keeps the data
                                 // around. Even though the array
-                                // <code>derivatives</code> takes the
-                                // majority of the matrix memory
-                                // consumptions, it still pays off to have
-                                // that data around since it would be quite
-                                // expensive to manually compute it every
-                                // time we make a matrix-vector product.
+                                // <code>derivatives</code> stands for the
+                                // majority of the matrix memory consumption,
+                                // it still pays off to have that data around
+                                // since it would be quite expensive to
+                                // manually compute it every time we make a
+                                // matrix-vector product.
 template <typename number, class Transformation>
 void MatrixFree<number,Transformation>::
 set_derivative_data (const unsigned int cell_no,
@@ -394,20 +403,20 @@ set_derivative_data (const unsigned int cell_no,
                                 // <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.
+                                // call it simultaneously on many processors,
+                                // but with different cell range data.
                                 //
                                 // Following the discussion in the
                                 // introduction, we try to work on multiple
                                 // cells at a time. This is possible
-                                // because the small matrix stays the same
+                                // because the small matrix remains the same
                                 // on all the cells, and only the
                                 // derivative information from the Jacobian
                                 // is different. That way, the operation
                                 // that is actually the multiplication of
                                 // the small matrix with a vector (on the
                                 // local dofs) becomes a multiplication of
-                                // two full (small) matrices with each
+                                // two full but small matrices with each
                                 // other. This is an operation that can be
                                 // much better optimized than matrix-vector
                                 // products. The functions
@@ -415,7 +424,7 @@ set_derivative_data (const unsigned int cell_no,
                                 // and
                                 // <code>FullMatrix<number>::mTmult</code>
                                 // use the BLAS dgemm function (as long as
-                                // it is detected in deal.II
+                                // BLAS has been detected in deal.II
                                 // configuration), which provides optimized
                                 // kernels for doing this product. In our
                                 // case, a matrix-matrix product is between
@@ -440,8 +449,7 @@ MatrixFree<number,Transformation>::
               WorkStreamData::CopyData<number>    &copy,
               const Vector<number2>               &src) const
 {
-  const unsigned int first_cell = cell_range->first,
-    chunk_size = cell_range->second - cell_range->first;
+  const unsigned int 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
@@ -497,30 +505,30 @@ MatrixFree<number,Transformation>::
                                 // are related to each other. Since we
                                 // simultaneously apply the constraints, we
                                 // hand this task off to the ConstraintMatrix
-                                // object. We do this in an extra function
+                                // object. 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 feed the
+                                // function with data from all the cells at
+                                // once. 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.
+                                // function).
   copy.solutions.reinit    (chunk_size,matrix_sizes.m, true);
-  copy.first_dof          = first_cell;
+  copy.first_cell         = cell_range->first;
   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),
+  constraints.get_dof_values(src, &indices_local_to_global(copy.first_cell,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 i=0, k = copy.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));
 
@@ -538,7 +546,7 @@ MatrixFree<number,Transformation>::
 {
   constraints.distribute_local_to_global (&copy.solutions(0,0),
                                          &copy.solutions(0,0)+copy.n_dofs,
-                                         &indices_local_to_global(copy.first_dof,0),
+                                         &indices_local_to_global(copy.first_cell,0),
                                          dst);
 }
 
@@ -592,20 +600,46 @@ MatrixFree<number,Transformation>::Tvmult_add (Vector<number2>       &dst,
                                 // This is the <code>vmult_add</code>
                                 // function that multiplies the matrix with
                                 // vector <code>src</code> and adds the
-                                // result to vector <code>dst</code>. We call
-                                // a %parallel function that applies the
+                                // result to vector <code>dst</code>.  We
+                                // include a few sanity checks to make sure
+                                // that the size of the vectors is the same
+                                // as the dimension of the matrix. We call a
+                                // %parallel function that applies the
                                 // 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>.
+                                // <code>matrix_sizes.chunks</code>. What the
+                                // rather cryptic command to
+                                // <code>std_cxx1x::bind</code> does is to
+                                // transform a function that has several
+                                // arguments (source vector, chunk
+                                // information) into a function which has no
+                                // arguments, which is what the
+                                // WorkStream::run function expects. The
+                                // placeholders <code>_1, _2, _3</code> in
+                                // the local vmult specify variable input
+                                // values, given by the chunk information,
+                                // scratch data and copy data. Similarly, the
+                                // placeholder <code>_1</code> in the
+                                // <code>copy_local_to_global</code> function
+                                // sets the first argument of that function,
+                                // which is of class
+                                // <code>CopyData</code>. We need to
+                                // abstractly specify these arguments because
+                                // the tasks defined by different cell chunks
+                                // will be scheduled by the WorkStream class,
+                                // and we will reuse available scratch and
+                                // copy data.
 template <typename number, class Transformation>
 template <typename number2>
 void
 MatrixFree<number,Transformation>::vmult_add (Vector<number2>       &dst,
                                              const Vector<number2> &src) const
 {
+  Assert (src.size() == n(), ExcDimensionMismatch(src.size(), n()));
+  Assert (dst.size() == m(), ExcDimensionMismatch(dst.size(), m()));
 
   WorkStream::run (matrix_sizes.chunks.begin(), matrix_sizes.chunks.end(),
                   std_cxx1x::bind(&MatrixFree<number,Transformation>::
@@ -626,12 +660,13 @@ MatrixFree<number,Transformation>::vmult_add (Vector<number2>       &dst,
                                 // want). Since the
                                 // <code>distribute_local_to_global</code>
                                 // command of the constraint matrix which we
-                                // used for add the local elements into the
-                                // global matrix does not do write anything
+                                // used for adding the local elements into
+                                // the global vector does not do anything
                                 // with constrained elements, we have to
-                                // circumvent that problem by setting the
-                                // diagonal to some non-zero value. We simply
-                                // set it to one.
+                                // circumvent that problem by artificially
+                                // setting the diagonal to some non-zero
+                                // value and adding the source values. We
+                                // simply set it to one.
   for (unsigned int i=0; i<matrix_sizes.n_dofs; ++i)
     if (constraints.is_constrained(i) == true)
       dst(i) += 1.0 * src(i);
@@ -663,7 +698,7 @@ MatrixFree<number,Transformation>::vmult_add (Vector<number2>       &dst,
                                 // derivative information and the local dof
                                 // indices the correct sizes. They will be
                                 // filled by calling the respective set
-                                // function.
+                                // function defined above.
 template <typename number, class Transformation>
 void MatrixFree<number,Transformation>::
 reinit (const unsigned int        n_dofs_in,
@@ -689,7 +724,7 @@ reinit (const unsigned int        n_dofs_in,
 
                                 // One thing to make the matrix-vector
                                 // product with this class efficient is to
-                                // decide how many cells should be summarized
+                                // decide how many cells should be combined
                                 // to one chunk, which will determine the
                                 // size of the full matrix that we work
                                 // on. If we choose too few cells, then the
@@ -698,26 +733,27 @@ reinit (const unsigned int        n_dofs_in,
                                 // 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
+                                // parallelization (we need to have
+                                // sufficiently independent tasks), 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
+                                // still be present in processor cache and
+                                // not needed to be fetched from main
+                                // memory. The total memory usage of the data
+                                // on quadrature points should not be more
+                                // than about a third of 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 150 kB as a size
+                                // to leave some room for other things to be
+                                // stored in the CPU. Clearly, this is an
                                 // architecture-dependent value and the
                                 // interested user can squeeze out some extra
                                 // performance by hand-tuning this
@@ -727,7 +763,7 @@ reinit (const unsigned int        n_dofs_in,
                                 // 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));
+  const unsigned int divisor = 150000/(matrix_sizes.n*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;
@@ -777,18 +813,16 @@ MatrixFree<number,Transformation>::clear ()
 
                                 // This function returns the entries of the
                                 // matrix. Since this class is intended not
-                                // to store the matrix entries, it would
-                                // not make sense to provide all those
+                                // to store the matrix entries, it would make
+                                // no sense to provide all those
                                 // elements. However, diagonal entries are
-                                // explicitly needed in some places, like
-                                // handling the matrix-vector product on
-                                // constrained degrees of freedom or for
-                                // the implementation of the Chebyshev
-                                // smoother that we intend to use in the
-                                // multigrid preconditioner. This matrix is
-                                // equipped with a vector that stores the
-                                // diagonal, and we compute it when this
-                                // function is called for the first time.
+                                // explicitly needed for the implementation
+                                // of the Chebyshev smoother that we intend
+                                // to use in the multigrid
+                                // preconditioner. This matrix is equipped
+                                // with a vector that stores the diagonal,
+                                // and we compute it when this function is
+                                // called for the first time.
 template <typename number, class Transformation>
 number
 MatrixFree<number,Transformation>::el (const unsigned int row,
@@ -803,27 +837,26 @@ MatrixFree<number,Transformation>::el (const unsigned int row,
 
 
 
-                                // Regarding the calculation of the
-                                // diagonal, remember that this is as
-                                // simple (or complicated) as assembling a
-                                // right hand side in deal.II. Well, it is
-                                // a bit easier to do this within this
-                                // class since have all the derivative
-                                // information available. What we do is to
-                                // go through all the cells (now in serial,
-                                // since this function should not be called
-                                // very often anyway), then all the degrees
-                                // of freedom. At this place, we first copy
-                                // the first basis functions in all the
+                                // Regarding the calculation of the diagonal,
+                                // remember that this is as simple (or
+                                // complicated) as assembling a right hand
+                                // side in deal.II. Well, it is a bit easier
+                                // to do this within this class since have
+                                // all the derivative information
+                                // available. What we do is to go through all
+                                // the cells (now in serial, since this
+                                // function should not be called very often
+                                // anyway), then all the degrees of
+                                // freedom. At this place, we first copy the
+                                // first basis functions in all the
                                 // quadrature points to a temporary array,
                                 // apply the derivatives from the Jacobian
                                 // matrix, and finally multiply with the
-                                // second basis function. This is exactly
-                                // the value that would be written into the
-                                // diagonal of a sparse matrix. Note that
-                                // we need to condense hanging node
-                                // constraints and set the diagonals to
-                                // one.
+                                // second basis function. This is exactly the
+                                // value that would be written into the
+                                // diagonal of a sparse matrix. Note that we
+                                // need to condense hanging node constraints
+                                // and set the constrained diagonals to one.
 template <typename number, class Transformation>
 void
 MatrixFree<number,Transformation>::calculate_diagonal() const
@@ -940,35 +973,35 @@ LaplaceOperator<dim,number>::LaplaceOperator(const Tensor<2,dim> &tensor)
   *this = tensor;
 }
 
-                                // Now implement the transformation, which
-                                // is nothing else than a so-called
-                                // contract operation of a tensor of second
-                                // rank on a tensor of first
-                                // rank. Unfortunately, we need to
-                                // implement this by hand, since we don't
-                                // have tensors (note that the result
-                                // values are entries in a full matrix that
-                                // consists of doubles or floats). It might
-                                // feel a bit unsafe to operate on a
-                                // pointer to the data, but that is the
-                                // only possibility if we do not want to
-                                // copy data back and forth, which is
-                                // expensive since this is the innermost
-                                // position of the loop in the
-                                // <code>vmult</code> operation of the
-                                // MatrixFree class. We need to remember
-                                // that we only saved half the (symmetric)
-                                // rank-two tensor.
+                                // Now implement the transformation, which is
+                                // nothing else than a so-called contract
+                                // operation of a tensor of second rank on a
+                                // tensor of first rank. Unfortunately, we
+                                // need to implement this by hand, since we
+                                // chose not to use the
+                                // SymmetricTensor<2,dim> class (note that
+                                // the resulting values are entries in a full
+                                // matrix that consists of doubles or
+                                // floats). It feels a bit unsafe to operate
+                                // on a pointer to the data, but that is the
+                                // only possibility if we do not want to copy
+                                // data back and forth, which is expensive
+                                // since this is the innermost position of
+                                // the loop in the <code>vmult</code>
+                                // operation of the MatrixFree class. We need
+                                // to pay attention to the fact that we only
+                                // saved half the (symmetric) rank-two
+                                // tensor.
                                 //
-                                // It might seem inefficient that we have
-                                // an <code>if</code> clause at this place
-                                // (which is the innermost loop), but note
-                                // once again that <code>dim</code> is
-                                // known when this piece of code is
-                                // compiled, so the compiler can optize
-                                // away the <code>if</code> statement (and
-                                // actually even inline these few lines of
-                                // code into the <code>MatrixFree</code>
+                                // At first sight, it seems inefficient that
+                                // we have an <code>if</code> clause at this
+                                // position in the code at the innermost
+                                // loop, but note once again that
+                                // <code>dim</code> is known when this piece
+                                // of code is compiled, so the compiler can
+                                // optize away the <code>if</code> statement
+                                // (and actually even inline these few lines
+                                // of code into the <code>MatrixFree</code>
                                 // class).
 template <int dim, typename number>
 void LaplaceOperator<dim,number>::transform (number* result) const
@@ -1365,14 +1398,15 @@ void LaplaceProblem<dim>::assemble_multigrid ()
 
                                 // The solution process again looks like
                                 // step-16. We now use a Chebyshev smoother
-                                // instead of SOR (SOR would very
+                                // instead of SOR (SOR would be very
                                 // difficult to implement because we do not
                                 // have the matrix elements explicitly
-                                // available, and it is difficult to make
-                                // it work efficiently in %parallel). The
+                                // available, and it is difficult to make it
+                                // work efficiently in %parallel). The
                                 // multigrid classes provide a simple
-                                // interface for using the Chebyshev
-                                // smoother: MGSmootherPrecondition.
+                                // interface for using the Chebyshev smoother
+                                // which is defined in a preconditioner
+                                // class: MGSmootherPrecondition.
 template <int dim>
 void LaplaceProblem<dim>::solve ()
 {
@@ -1393,9 +1427,15 @@ void LaplaceProblem<dim>::solve ()
                                   // data for the Chebyshev smoother. Use a
                                   // higher polynomial degree for higher
                                   // order elements, since smoothing gets
-                                  // more difficult then. Smooth out a
-                                  // range of
-                                  // $[\lambda_{\max}/10,\lambda_{\max}]$.
+                                  // more difficult then. Smooth out a range
+                                  // of
+                                  // $[\lambda_{\max}/10,\lambda_{\max}]$. In
+                                  // order to compute the maximum eigenvalue
+                                  // of the corresponding matrix, the
+                                  // Chebyshev initializations performs a few
+                                  // steps of a CG algorithm. Since all we
+                                  // need is a rough estimate, we choose some
+                                  // eight iterations.
   typename SMOOTHER::AdditionalData smoother_data;
   smoother_data.smoothing_range = 10.;
   smoother_data.degree = fe.degree;

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.