]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Improve comments a bit further.
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Fri, 6 Nov 2009 23:06:37 +0000 (23:06 +0000)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Fri, 6 Nov 2009 23:06:37 +0000 (23:06 +0000)
git-svn-id: https://svn.dealii.org/trunk@20057 0785d39b-7218-0410-832d-ea1e28bc413d

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

index 40dafdf9c8d8f88a3feeba3b00947f2ee3d8dcdd..09d319ddae53c9aafc5a82bc1183796b331a7bf2 100644 (file)
@@ -30,9 +30,11 @@ product when the matrix is stored in the usual sparse compressed row storage
 &mdash; in short CRS &mdash; format implemented by the SparsityPattern and
 SparseMatrix classes (the actual implementation in deal.II uses a slightly
 different structure for the innermost loop, thereby avoiding the counter
-variable):
+variable), also used by Trilinos and PETSc matrices:
 @code
 // y = A * x
+// variables: double * A_values;  unsigned int * A_column_indices;
+//            std::size_t * A_row_indices, unsigned int n_rows;
 std::size_t element_index = A_row_indices[0];
 for (unsigned int row=0; row<n_rows; ++row)
   {
@@ -48,43 +50,42 @@ for (unsigned int row=0; row<n_rows; ++row)
   }
 @endcode
 
-Looking into what actually happens in the computer when forming a
-matrix-vector product, one observes that the matrix data
-<code>A_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>A_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>A_values</code>, and 4 bytes for the unsigned
-integer position <code>A_column_indices</code> that indicates 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 those 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.
+The matrix-vector product basically goes through the @p A_values array,
+multiplies the current element <code>A_values[element_index]</code> with the
+entry in the vector that corresponds to that element, and then accumulates
+this product to the sum along the current row. Each matrix element thus
+involves two floating-point operations, a multiplication and an
+addition. Assume now we have a matrix with a billion (10<sup>9</sup>)
+entries. Then, a matrix-vector product requires 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
+floating point 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.
+
+However, that is usually not the case &mdash; because of memory. Each matrix
+element corresponds to 12 bytes of data, 8 bytes for the actual data in @p
+A_values and 4 bytes for the unsigned integer column position of that element,
+which we need to find the correct vector element. For our one-billion-sized
+matrix, these two arrays make up 12 GB of data, which we need to stream into
+the processor. 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. Here, we neglect the additional storage required by
+the additional array @p A_row_indices that tells us the range of the
+individual rows for simplicity, and assumed that the vector data is stored in
+some fast (cache) memory. 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, 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
-a nightmarish 6 per cent of the system's peak performance. And this is the
-theoretical maximum we can get!
+cores will at best 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 a nightmarish 6 per cent of the system's peak
+performance. And this is the theoretical maximum!
 
 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
@@ -110,18 +111,16 @@ can be streamed into a processor, this trade-off will be worthwhile.
 
 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:
+looking how some finite-element related matrix <i>A</i> is assembled:
 @f{eqnarray*}
 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 rectangular
+In this formula, the matrix <i>P</i><sub>cell,loc-glob</sub> is a rectangular
 matrix that defines the index mapping from local degrees of freedom in the
-current cell
-to the global degrees of freedom. The information from which this
-operator can be built is usually encoded in the
-<code>local_dof_indices</code> variable we have always used in the
-assembly of matrices.
+current cell to the global degrees of freedom. The information from which this
+operator can be built is usually encoded in the <code>local_dof_indices</code>
+variable we have always used in the assembly of matrices.
 
 If we are to perform a matrix-vector product, we can hence use that
 @f{eqnarray*}
@@ -134,8 +133,8 @@ 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}$ are the values of <i>x</i> at the degrees of freedom
-of the respective cell, and $y_\mathrm{cell}$ correspondingly for the result.
+where <i>x</i><sub>cell</sub> are the values of <i>x</i> at the degrees of freedom
+of the respective cell, and  <i>x</i><sub>cell</sub> correspondingly for the result.
 A naive attempt to implement the local action of the Laplacian would hence be
 to use the following code:
 @code
@@ -185,14 +184,16 @@ MatrixFree::vmult (Vector<double>       &dst,
 }
 @endcode
 
-Here we neglected boundary conditions as well as any hanging nodes we
-may have, though neither would be very difficult to
-include using the ConstraintMatrix class. Note how we first generate the local
-matrix in the usual way. To form the actual product as expressed in the
-above formula, we read in the values of <code>src</code> of the cell-related
-degrees of freedom, multiply by the local matrix, and finally add the result
-to the destination vector <code>dst</code>. It is not more difficult than
-that, in principle.
+Here we neglected boundary conditions as well as any hanging nodes we may
+have, though neither would be very difficult to include using the
+ConstraintMatrix class. Note how we first generate the local matrix in the
+usual way. To form the actual product as expressed in the above formula, we
+read in the values of <code>src</code> of the cell-related degrees of freedom
+(the action of <i>P</i><sub>cell,loc-glob</sub>), multiply by the local matrix
+(the action of <i>A</i><sub>cell</sub>), and finally add the result to the
+destination vector <code>dst</code> (the action of
+<i>P</i><sub>cell,loc-glob</sub><sup>T</sup>, added over all the elements). It
+is not more difficult than that, in principle.
 
 While this code is completely correct, it is very slow. For every cell, we
 generate a local matrix, which takes three nested loops with as many
@@ -205,12 +206,12 @@ matrix can be thought of as the product of three matrices,
 @f{eqnarray*}
 A_\mathrm{cell} = B_\mathrm{cell}^T D_\mathrm{cell} B_\mathrm{cell},
 @f}
-where for the example of the Laplace operator 
-the $(q*\mathrm{dim}+d,i)$-th element of <i>B</i><sub>cell</sub> is given by
+where for the example of the Laplace operator the (<i>q</i>*dim+<i>d,i</i>)-th
+element of <i>B</i><sub>cell</sub> is given by
 <code>fe_values.shape_grad(i,q)[d]</code>. The matrix consists of
-<code>dim*n_q_pointsdofs_per_cell</code> rows and <code>dofs_per_cell</code>
-columns). The matrix <i>D</i><sub>cell</sub> is diagonal and contains the values
-<code>fe_values.JxW(q)</code> (or, rather, <code>dim</code> copies of it).
+<code>dim*n_q_points</code> rows and @p dofs_per_cell columns). The matrix
+<i>D</i><sub>cell</sub> is diagonal and contains the values
+<code>fe_values.JxW(q)</code> (or, rather, @p dim copies of it).
 
 Every numerical analyst learns in one of her first classes that for
 forming a product of the form
@@ -218,9 +219,9 @@ forming a product of the form
 A_\mathrm{cell}\cdot x_\mathrm{cell} = B_\mathrm{cell} D_\mathrm{cell} 
                     B_\mathrm{cell}^T \cdot x_\mathrm{cell},
 @f}
-one should never form the matrix-matrix products, but rather multiply
-with the vector from right to left so that only matrix-vector products are
-formed. To put this into code, we can write:
+one should never form the matrix-matrix products, but rather multiply with the
+vector from right to left so that only three successive matrix-vector products
+are formed. To put this into code, we can write:
 @code
 ...
   for (; cell!=endc; ++cell)
@@ -270,24 +271,23 @@ from something like $\mathcal {O}(\mathrm{dofs\_per\_cell}^3)$ to $\mathcal
 a slightly more clever use of data in order to gain some extra speed. It does
 not change the code structure, though.
 
-If one would implement the code above, one would soon realize that the
-operations done by the call <code>fe_values.reinit(cell)</code> take about
-as much time as the other steps together (at least if the mesh is
-unstructured; deal.II can recognize that the gradients are often unchanged
-on structured meshes). That is certainly not ideal and we would like to do
-better than this. What the reinit function does is to calculate the
-gradient in real space by transforming the gradient on the reference cell
-using the Jacobian of the transformation from real to reference cell. This is
-done for each basis function on the cell on each quadrature point. The
-Jacobian does not depend on the basis function, but it is different on
-different quadrature points in general. The trick is now to factor out the
-Jacobian transformation and first apply the operation that leads us to
-<code>temp_vector</code> only with the gradient on the reference cell. That
-transforms the vector of values on the local dofs to a vector of gradients
-on the quadrature points. There, we first apply the Jacobian that we
-factored out from the gradient, then we apply the weights of the quadrature,
-and we apply with the transposed Jacobian for preparing the third loop which
-again uses the gradients on the unit cell.
+The bottleneck in the above code is the operations done by the call
+<code>fe_values.reinit(cell)</code>, which take about as much time as the
+other steps together (at least if the mesh is unstructured; deal.II can
+recognize that the gradients are often unchanged on structured meshes). That
+is certainly not ideal and we would like to do better than this. What the
+reinit function does is to calculate the gradient in real space by
+transforming the gradient on the reference cell using the Jacobian of the
+transformation from real to reference cell. This is done for each basis
+function on the cell, for each quadrature point. The Jacobian does not depend on
+the basis function, but it is different on different quadrature points in
+general. The trick is now to factor out the Jacobian transformation and first
+apply the operation that leads us to <code>temp_vector</code> only with the
+gradient on the reference cell. That transforms the vector of values on the
+local dofs to a vector of gradients on the quadrature points. There, we first
+apply the Jacobian that we factored out from the gradient, then we apply the
+weights of the quadrature, and we apply with the transposed Jacobian for
+preparing the third loop which again uses the gradients on the unit cell.
 
 Let us again write this in terms of matrices. Let the matrix
 <i>B</i><sub>cell</sub> denote the cell-related gradient matrix, with each row
@@ -297,11 +297,10 @@ matrix-matrix product as
 B_\mathrm{cell} = J_\mathrm{cell} B_\mathrm{ref\_cell},
 @f}
 where <i>B</i><sub>ref_cell</sub> denotes the gradient on the reference cell
-and <i>J</i><sub>cell</sub><sup>T</sup> denotes the
-Jacobian. <i>J</i><sub>cell</sub><sup>T</sup> is block-diagonal, and the
-blocks size is equal to the dimension of the problem. Each diagonal block is
-the Jacobian transformation that goes from the reference cell to the real 
-cell.
+and <i>J</i><sub>cell</sub> denotes the Jacobian
+transformation. <i>J</i><sub>cell</sub> is block-diagonal, and the blocks size
+is equal to the dimension of the problem. Each diagonal block is the Jacobian
+transformation that goes from the reference cell to the real cell.
 
 Putting things together, we find that
 @f{eqnarray*}
@@ -310,10 +309,12 @@ A_\mathrm{cell} = B_\mathrm{cell}^T D B_\mathrm{cell}
                  D_\mathrm{cell} 
                  J_\mathrm{cell} B_\mathrm{ref\_cell},
 @f}
-so we calculate the product (starting from the right)
+so we calculate the product (starting the local product from the right)
 @f{eqnarray*}
-B_\mathrm{ref\_cell}^T J_\mathrm{cell}^T D J_\mathrm{cell} 
-B_\mathrm{ref\_cell}\cdot x_\mathrm{cell}.
+y_\mathrm{cell} = B_\mathrm{ref\_cell}^T J_\mathrm{cell}^T D J_\mathrm{cell} 
+B_\mathrm{ref\_cell} x_\mathrm{cell}, \quad
+y = \sum_{\mathrm{cell}=1}^{\mathrm{n\_cells}} P_\mathrm{cell,{loc-glob}}^T
+y_\mathrm{cell}.
 @f}
 @code
 ...
@@ -414,11 +415,11 @@ there are a few more points done to be even more efficient, namely:
 <ul>
   <li>We pre-compute the inverse of the Jacobian of the transformation and
       store it in an extra array. This allows us to fuse the three
-      operations <i>J</i><sub>cell</sub><sup>T</sup><i> D </i>
+      operations <i>J</i><sub>cell</sub><sup>T</sup><i> D</i><sub>cell</sub>
       <i>J</i><sub>cell</sub> (apply Jacobian, multiply by 
       weights, apply transposed
       Jacobian) into one second-rank tensor that is also symmetric (so we
-      can save only half the tensor).
+      only need to store half the tensor).
   <li>We work on several cells at once when we apply the gradients of the
       unit cell (it is always the same matrix with the reference cell
       data). This allows us to replace the matrix-vector product by a
@@ -446,20 +447,19 @@ 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 source 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
-require random access to a vector. This kind of rather uniform data access
-should make it not too difficult to implement a matrix-free
+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 @p src and the write operation to @p dst in
+the end perform more random access to a vector. This kind of rather uniform
+data access should make it not too difficult to implement a matrix-free
 matrix-vector product on a <a
 href="http://en.wikipedia.org/wiki/GPGPU">graphics processing unit
 (GP-GPU)</a>, for example. On the contrary, it would be quite complex to make
@@ -479,16 +479,15 @@ contributions into the global vector, in order to avoid a race condition.
 <h3>Combination with multigrid</h3>
 
 Above, we have gone to significant lengths to implement a matrix-vector
-product that does not actually store the matrix elements. In many user
-codes, however, one wants more than just performing some uncertain number of
+product that does not actually store the matrix elements. In many user codes,
+however, one wants more than just performing some uncertain number of
 matrix-vector products &mdash; one wants to do as little of these operations
-as possible when solving linear equation systems. In theory, we could use
-the CG method without preconditioning; however, that would not be very
-efficient. Rather, one uses preconditioners for improving speed. On the
-other hand, most of the more frequently used preconditioners such as Jacobi,
-SSOR, ILU or algebraic multigrid (AMG) can now no longer be used here
-because their implementation requires knowledge of the elements of the
-system matrix.
+as possible when solving linear equation systems. In theory, we could use the
+CG method without preconditioning; however, that would not be very
+efficient. Rather, one uses preconditioners for improving speed. On the other
+hand, most of the more frequently used preconditioners such as SSOR, ILU or
+algebraic multigrid (AMG) can now no longer be used here because their
+implementation requires knowledge of the elements of the 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
index ff07b2d8e3590271ef2caf1cce766778fcbcad09..f545b378d40b8a84d0084f7eecf901b4331f1bcd 100644 (file)
@@ -108,7 +108,7 @@ void Coefficient<dim>::value_list (const std::vector<Point<dim> > &points,
 
                                 // In this program, we want to make
                                 // use of the ability of deal.II to
-                                // runs things in parallel if compute
+                                // runs things in %parallel if compute
                                 // resources are available. We will
                                 // follow the general framework laid
                                 // out in the @ref threads module and
@@ -205,11 +205,11 @@ namespace WorkStreamData
                                 // 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>
+                                // command of the @p Transformation
                                 // 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
+                                // operation without the need for substantial
                                 // changes.
 template <typename number, class Transformation>
 class MatrixFree : public Subscriptor
@@ -253,7 +253,7 @@ class MatrixFree : public Subscriptor
     std::size_t memory_consumption () const;
 
                                     // The private member variables of the
-                                    // <code>MatrixFree</code> class are a
+                                    // @p MatrixFree class are a
                                     // small matrix that does the
                                     // transformation from solution values to
                                     // quadrature points, a list with the
@@ -302,13 +302,13 @@ class MatrixFree : public Subscriptor
 
 
 
-                                // This is the constructor of the
-                                // <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.
+                                // This is the constructor of the @p
+                                // MatrixFree class. All it does is to
+                                // subscribe to the general deal.II @p
+                                // Subscriptor 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 ()
                :
@@ -360,24 +360,20 @@ MatrixFree<number,Transformation>::get_constraints ()
 
 
 
-                                // The following 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 have
-                                // computed to have fast access to
-                                // those elements, is invalidated. We
-                                // set the respective flag to
-                                // <code>false</code>.
+                                // The following function takes a vector of
+                                // local dof indices on cell level and writes
+                                // the data into the
+                                // @p indices_local_to_global 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 have computed to have fast access to
+                                // those elements, is invalidated. We set the
+                                // respective flag to @p false.
 template <typename number, class Transformation>
 void MatrixFree<number,Transformation>::
 set_local_dof_indices (const unsigned int               cell_no,
@@ -396,18 +392,16 @@ set_local_dof_indices (const unsigned int               cell_no,
 
 
 
-                                // Next a function that writes the
-                                // derivative data on a certain cell
-                                // and a certain quadrature point to
-                                // the array that keeps the data
-                                // around. Even though the array
-                                // <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.
+                                // Next a function that writes the derivative
+                                // data on a certain cell and a certain
+                                // quadrature point to the array that keeps
+                                // the data around. Even though the array @p
+                                // derivatives 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,
@@ -425,9 +419,9 @@ set_derivative_data (const unsigned int cell_no,
                                 // matrix-free class, implementing the
                                 // multiplication of the matrix with a
                                 // vector. This function does not actually
-                                // work on all cells of a mesh, but only the subset
-                                // of cells specified by the first argument
-                                // <code>cell_range</code>. Since this
+                                // work on all cells of a mesh, but only the
+                                // subset of cells specified by the first
+                                // argument @p cell_range. Since this
                                 // function operates similarly irrespective
                                 // on which cell chunk we are sitting, we can
                                 // call it simultaneously on many processors,
@@ -450,8 +444,8 @@ set_derivative_data (const unsigned int cell_no,
                                 // D_\mathrm{cell} 
                                 // J_\mathrm{cell} B_\mathrm{ref\_cell}
                                 // @f}
-                                // and $P_\mathrm{cell,local-global}$ is the
-                                // transformation from local to global
+                                // and <i>P</i><sub>cell,local-global</sub>
+                                // is the transformation from local to global
                                 // indices.
                                 //
                                 // To do this, we would have to do the
@@ -459,31 +453,32 @@ set_derivative_data (const unsigned int cell_no,
                                 // <ol>
                                 //   <li> Form $x_\mathrm{cell} =
                                 //   P_\mathrm{cell,local-global} x$. This is
-                                //   done using the command
+                                //   done by using the command
                                 //   ConstraintMatrix::get_dof_values.
                                 //   <li> Form $x_1 = B_\mathrm{ref\_cell}
-                                //   x_\mathrm{cell}$. The vector $x_1$
-                                //   contains the reference cell gradient to
-                                //   the local cell vector.
+                                //   x_\mathrm{cell}$. The vector
+                                //   <i>x</i><sub>1</sub> contains the
+                                //   reference cell gradient to the local
+                                //   cell vector.
                                 //   <li> Form $x_2 = J_\mathrm{cell}^T
                                 //   D_\mathrm{cell} J_\mathrm{cell}
                                 //   x_1$. This is a block-diagonal
                                 //   operation, with the block size equal to
-                                //   <code>dim</code>. The blocks just
+                                //   @p dim. The blocks just
                                 //   correspond to the individual quadrature
                                 //   points. The operation on each quadrature
                                 //   point is implemented by the
                                 //   Transformation class object that this
                                 //   class is equipped with. Compared to the
                                 //   introduction, the matrix
-                                //   $D_\mathrm{cell}$ now contains the
-                                //   <code>JxW</code> values and the
+                                //   <i>D</i><sub>cell</sub> now contains the
+                                //   @p JxW values and the
                                 //   inhomogeneous coefficient.
                                 //   <li> Form $y_\mathrm{cell} =
                                 //   B_\mathrm{ref\_cell}^T x_2$. This gives
                                 //   the local result of the matrix-vector
                                 //   product.
-                                //   <li> Form $y +=
+                                //   <li> Form $y \leftarrow y +
                                 //   P_\mathrm{cell,local-global}^T
                                 //   y_\mathrm{cell}$. This adds the local
                                 //   result to the global vector, which is
@@ -491,49 +486,53 @@ set_derivative_data (const unsigned int cell_no,
                                 //   ConstraintMatrix::distribute_local_to_global.
                                 //   Note that we do this in an extra
                                 //   function called
-                                //   <code>copy_local_to_global</code>
+                                //   @p copy_local_to_global
                                 //   because that operation must not be done
                                 //   in %parallel, in order to avoid two or
                                 //   more processes trying to add to the same
-                                //   positions in the result vector $y$.
+                                //   positions in the result vector <i>y</i>.
                                 //   </ol>
                                 // The steps 1 to 4 can be done in %parallel
                                 // by multiple processes.
 
                                 // Now, it turns out that the most expensive
                                 // part of the above is the multiplication
-                                // $B_\mathrm{ref\_cell} x_\mathrm{cell}$ in
-                                // the second step and the transpose
-                                // operation in step 4. Note that the matrix
-                                // $J^T D J$ is block-diagonal, and hence,
-                                // its application is cheaper. Since the
-                                // matrix $B_\mathrm{ref\_cell}$ is the same
+                                // <i>B</i><sub>ref_cell</sub>
+                                // <i>x</i><sub>cell</sub> in the second step
+                                // and the transpose operation in step
+                                // 4. Note that the matrix
+                                // <i>J</i><sup>T</sup><i> D J</i> is
+                                // block-diagonal, and hence, its application
+                                // is cheaper. Since the matrix
+                                // <i>B</i><sub>ref_cell</sub> is the same
                                 // for all cells, all that changes is the
-                                // vector $x_\mathrm{cell}$. Hence, nothing
-                                // prevents us from collecting several cell
-                                // vectors to a (rectangular) matrix, and
-                                // then perform a matrix-matrix
+                                // vector <i>x</i><sub>cell</sub>. Hence,
+                                // nothing prevents us from collecting
+                                // several cell vectors to a (rectangular)
+                                // matrix, and then perform a matrix-matrix
                                 // product. These matrices are both full, but
-                                // not very large, having of the order
-                                // <code>dofs_per_cell</code> rows and
-                                // columns. This is an operation that can be
-                                // much better optimized than matrix-vector
-                                // products. The functions
-                                // <code>FullMatrix<number>::mmult</code> and
-                                // <code>FullMatrix<number>::mTmult</code>
-                                // use the BLAS dgemm function (as long as
-                                // 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
-                                // three and five times faster than doing the
-                                // matrix-vector product on one cell after
-                                // the other. The variables that hold the
-                                // solution on the respective cell's support
-                                // points and the quadrature points are thus
-                                // full matrices. The number of rows is given
-                                // by the number of cells they work on, and
-                                // the number of columns is the number of
+                                // not very large, having of the order @p
+                                // dofs_per_cell rows and columns. This is an
+                                // operation that can be much better
+                                // optimized than matrix-vector products. The
+                                // functions @p FullMatrix<number>::mmult and
+                                // @p FullMatrix<number>::mTmult use the BLAS
+                                // dgemm function (as long as 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 three and five times
+                                // faster than doing the matrix-vector
+                                // product on one cell after the other. The
+                                // variables that hold the solution on the
+                                // respective cell's support points and the
+                                // quadrature points are thus full matrices,
+                                // which we set to the correct size as a
+                                // first action in this function. The number
+                                // of rows in the two matrices @p
+                                // scratch.solutions and @p copy.solutions is
+                                // given by the number of cells they work on,
+                                // and the number of columns is the number of
                                 // degrees of freedom per cell for the first
                                 // and the number of quadrature points times
                                 // the number of components per point for the
@@ -549,10 +548,10 @@ local_vmult (CellChunkIterator                    cell_range,
 {
   const unsigned int chunk_size = cell_range->second - cell_range->first;
 
-  copy.solutions.reinit    (chunk_size,matrix_sizes.m, true);
+  scratch.solutions.reinit (chunk_size, matrix_sizes.n, true);
+  copy.solutions.reinit    (chunk_size, matrix_sizes.m, true);
   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_cell,0),
                             &copy.solutions(0,0),
@@ -584,12 +583,10 @@ copy_local_to_global (const WorkStreamData::CopyData<number> &copy,
 
 
 
-                                // Now to the <code>vmult</code> function
-                                // that is called externally: It is very
-                                // similar to the <code>vmult_add</code>
-                                // function, so just set the destination to
-                                // zero first, and then go to the other
-                                // function.
+                                // Now to the @p vmult function that is
+                                // called externally: In addition to what we
+                                // do in a @p vmult_add function, we set the
+                                // destination to zero first.
 template <typename number, class Transformation>
 template <typename number2>
 void
@@ -602,10 +599,11 @@ MatrixFree<number,Transformation>::vmult (Vector<number2>       &dst,
 
 
 
-                                // Transposed matrix-vector products: do
-                                // the same. Since we implement a symmetric
-                                // operation, we can refer to the vmult_add
-                                // operation.
+                                // Transposed matrix-vector products (needed
+                                // for the multigrid operations to be
+                                // well-defined): do the same. Since we
+                                // implement a symmetric operation, we can
+                                // refer to the @ vmult_add operation.
 template <typename number, class Transformation>
 template <typename number2>
 void
@@ -629,60 +627,47 @@ 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 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>. 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 three arguments
-                                // (in the first case) or one
-                                // argument (in the second), 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 that the WorkStream::run
-                                // function will provide, whereas the
-                                // other arguments to the
-                                // <code>local_vmult</code> function
-                                // are bound: to <code>this</code>
-                                // and a constant reference to the
-                                // <code>src</code> in the first
-                                // case, and <code>this</code> and a
-                                // reference to the output vector in
-                                // the second. Similarly, the
-                                // placeholder <code>_1</code> in the
-                                // <code>copy_local_to_global</code>
-                                // function sets the first explicit
-                                // 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.
+                                // This is the @p vmult_add function that
+                                // multiplies the matrix with vector @p src
+                                // and adds the result to vector @p dst.  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 @p
+                                // matrix_sizes.chunks. What the rather
+                                // cryptic command to @p std_cxx1x::bind does
+                                // is to transform a function that has
+                                // several arguments (source vector, chunk
+                                // information) into a function which has
+                                // three arguments (in the first case) or one
+                                // argument (in the second), 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 that the
+                                // WorkStream::run function will provide,
+                                // whereas the other arguments to the @p
+                                // local_vmult function are bound: to @p this
+                                // and a constant reference to the @p src in
+                                // the first case, and @p this and a
+                                // reference to the output vector in the
+                                // second. Similarly, the placeholder
+                                // @p _1 argument in the
+                                // @p copy_local_to_global function
+                                // sets the first explicit argument of that
+                                // function, which is of class
+                                // @p CopyData. 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
@@ -732,27 +717,24 @@ MatrixFree<number,Transformation>::vmult_add (Vector<number2>       &dst,
 
 
 
-                                // The next 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>B_ref_cell</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
+                                // The next 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 @p
+                                // B_ref_cell). The number of quadrature
+                                // points needs to be passed through the last
+                                // variable @p n_points_per_cell, since the
+                                // number of columns in the small matrix is
+                                // @p dim*n_points_per_cell for the Laplace
+                                // problem (the Laplacian is a tensor and has
+                                // @p dim 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
@@ -951,10 +933,10 @@ MatrixFree<number,Transformation>::calculate_diagonal() const
                                 // consumption of the arrays, the
                                 // constraints, the small matrix and of the
                                 // local variables. Just as a remark: In 2D
-                                // and with data type <code>double</code>,
+                                // and with data type @p double,
                                 // about 80 per cent of the memory
                                 // consumption is due to the
-                                // <code>derivatives</code> array, while in 3D
+                                // @p derivatives array, while in 3D
                                 // this number is even 85 per cent.
 template <typename number, class Transformation>
 std::size_t MatrixFree<number,Transformation>::memory_consumption () const
@@ -973,38 +955,34 @@ std::size_t MatrixFree<number,Transformation>::memory_consumption () const
 
                                 // @sect3{Laplace operator implementation}
 
-                                // This class implements the local action
-                                // of a Laplace operator on a quadrature
+                                // This class implements the local action of
+                                // a Laplace operator on a quadrature
                                 // point. This is a very basic class
                                 // implementation, providing functions for
-                                // initialization with a Tensor of rank 2
-                                // and implementing the
-                                // <code>transform</code> operation needed
-                                // by the <code>MatrixFree</code>
-                                // class. There is one point worth noting:
-                                // The quadrature-point related action of
-                                // the Laplace operator is a tensor of rank
+                                // initialization with a Tensor of rank 2 and
+                                // implementing the @p transform operation
+                                // needed by the @p MatrixFree class. There
+                                // is one point worth noting: The
+                                // quadrature-point related action of the
+                                // Laplace operator is a tensor of rank
                                 // two. It is symmetric since it is the
                                 // product of the inverse Jacobian
-                                // transformation between unit and real
-                                // cell with its transpose (times
-                                // quadrature weights and a coefficient,
-                                // which are scalar), so we can just save
-                                // the diagonal and upper diagonal part. We
-                                // could use the SymmetricTensor<2,dim>
-                                // class for doing this, however, that
-                                // class is only based on
-                                // <code>double</code> %numbers. Since we
-                                // also want to use <code>float</code>
-                                // %numbers for the multigrid
+                                // transformation between unit and real cell
+                                // with its transpose (times quadrature
+                                // weights and a coefficient, which are
+                                // scalar), so we can just save the diagonal
+                                // and upper diagonal part. We could use the
+                                // SymmetricTensor<2,dim> class for doing
+                                // this, however, that class is only based on
+                                // @p double %numbers. Since we also want to
+                                // use @p float %numbers for the multigrid
                                 // preconditioner (in order to save memory
-                                // and computing time), we manually
-                                // implement this operator. Note that
-                                // <code>dim</code> is a template argument
-                                // and hence known at compile-time, so the
-                                // compiler knows that this symmetric
-                                // rank-2 tensor has 3 entries if used in
-                                // 2D and 6 entries if used in 3D.
+                                // and computing time), we manually implement
+                                // this operator. Note that @p dim is a
+                                // template argument and hence known at
+                                // compile-time, so the compiler knows that
+                                // this symmetric rank-2 tensor has 3 entries
+                                // if used in 2D and 6 entries if used in 3D.
 template <int dim,typename number>
 class LaplaceOperator
 {
@@ -1047,21 +1025,20 @@ LaplaceOperator<dim,number>::LaplaceOperator(const Tensor<2,dim> &tensor)
                                 // 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>
+                                // the loop in the @p vmult
                                 // operation of the MatrixFree class. We need
                                 // to pay attention to the fact that we only
                                 // saved half of the (symmetric) rank-two
                                 // tensor.
                                 //
                                 // 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
-                                // optimize away the <code>if</code> statement
-                                // (and actually even inline these few lines
-                                // of code into the <code>MatrixFree</code>
+                                // we have an @p if clause at this position
+                                // in the code at the innermost loop, but
+                                // note once again that @p dim is known when
+                                // this piece of code is compiled, so the
+                                // compiler can optimize away the @p if
+                                // statement (and actually even inline these
+                                // few lines of code into the @p MatrixFree
                                 // class).
 template <int dim, typename number>
 void LaplaceOperator<dim,number>::transform (number* result) const
@@ -1090,7 +1067,7 @@ void LaplaceOperator<dim,number>::transform (number* result) const
                                 // The final function in this group
                                 // takes the content of a rank-2
                                 // tensor and writes it to the field
-                                // <code>transformation</code> of
+                                // @p transformation of
                                 // this class. We save the upper part
                                 // of the symmetric tensor row-wise:
                                 // we first take the (0,0)-entry,
@@ -1296,27 +1273,25 @@ void LaplaceProblem<dim>::setup_system ()
                                 // reduced compared to step-16. All we need
                                 // to do is to assemble the right hand side
                                 // and to calculate the cell-dependent part
-                                // of the Laplace operator. The first task
-                                // is standard. The second is also not too
-                                // hard given the discussion in the
-                                // introduction: We need to take the
-                                // inverse of the Jacobian of the
-                                // transformation from unit to real cell,
-                                // multiply it with its transpose and
-                                // multiply the resulting rank-2 tensor
-                                // with the quadrature weights and the
-                                // coefficient values at the quadrature
-                                // points. To make this work, we add the
-                                // update flag
-                                // <code>update_inverse_jacobians</code> to
-                                // the FEValues constructor, and query the
-                                // inverse of the Jacobian in a loop over
-                                // the quadrature points (note that the
-                                // Jacobian is not related to any kind of
-                                // degrees of freedom directly). In the
-                                // end, we condense the constraints from
-                                // Dirichlet boundary conditions away from
-                                // the right hand side.
+                                // of the Laplace operator. The first task is
+                                // standard. The second is also not too hard
+                                // given the discussion in the introduction:
+                                // We need to take the inverse of the
+                                // Jacobian of the transformation from unit
+                                // to real cell, multiply it with its
+                                // transpose and multiply the resulting
+                                // rank-2 tensor with the quadrature weights
+                                // and the coefficient values at the
+                                // quadrature points. To make this work, we
+                                // add the update flag @p
+                                // update_inverse_jacobians to the FEValues
+                                // constructor, and query the inverse of the
+                                // Jacobian in a loop over the quadrature
+                                // points (note that the Jacobian is not
+                                // related to any kind of degrees of freedom
+                                // directly). In the end, we condense the
+                                // constraints from Dirichlet boundary
+                                // conditions away from the right hand side.
 template <int dim>
 void LaplaceProblem<dim>::assemble_system ()
 {
@@ -1621,7 +1596,7 @@ void LaplaceProblem<dim>::run ()
 
 
 
-                                // @sect3{The <code>main</code> function}
+                                // @sect3{The @p main function}
 
                                 // This is as in all other programs:
 int main ()

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.