This example shows how to implement a matrix-free method, that is, a
method that does not explicitly store the matrix elements, for a
second-order Poisson equation with variable coefficients on a
-hypercube. The elliptic equation will be solved with a multigrid
+hypercube. The linear system will be solved with a multigrid
method.
-The major motivation for matrix-free methods is the fact that today
-access to main memory (i.e., for objects that don't fit in the cache)
-has become the bottleneck in scientific computing: To perform a
-matrix-vector product, modern CPUs spend far more time waiting for
-data to arrive from memory than on actually doing the floating point
-multiplications and additions. Thus, if we could substitute looking up
-matrix elements in memory by re-computing them — or rather, the
-operator represented by these entries —, we may win in terms of
-overall run-time (even if this requires a significant number of
-additional floating point operations). That said, to realize this with
-a trivial implementation is not enough and one needs to really look at
-what it takes to make this happen. This tutorial program (and the
-papers referenced above) show how one can implement such a scheme and
-demonstrates the speedup that can be obtained.
+The major motivation for matrix-free methods is the fact that on today's
+processors access to main memory (i.e., for objects that do not fit in the
+caches) has become the bottleneck in scientific computing: To perform a
+matrix-vector product based on matrices, modern CPUs spend far more time
+waiting for data to arrive from memory than on actually doing the floating
+point multiplications and additions. Thus, if we could substitute looking up
+matrix elements in memory by re-computing them — or rather, the operator
+represented by these entries —, we may win in terms of overall run-time
+(even if this requires a significant number of additional floating point
+operations). That said, to realize this with a trivial implementation is not
+enough and one needs to really look at the details to gain in
+performance. This tutorial program (and the papers referenced above) show how
+one can implement such a scheme and demonstrates the speedup that can be
+obtained.
<h3>The test case</h3>
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. Moreover,
-<i>A</i><sub>cell</sub> denotes the cell-operation associated with <i>A</i>.
+variable is used in the assembly calls filling matrices in deal.II. Moreover,
+<i>A</i><sub>cell</sub> denotes the cell matrix associated with <i>A</i>.
If we are to perform a matrix-vector product, we can hence use that
@f{eqnarray*}
correspondingly for the result.
A naive attempt to implement the local action of the Laplacian would hence be
to use the following code:
-@code
-MatrixFree<dim>::vmult (Vector<double> &dst,
- const Vector<double> &src) const
+@code Matrixfree<dim>::vmult (Vector<double> &dst,
+ const Vector<double> &src) const
{
dst = 0;
QGauss<dim> quadrature_formula(fe.degree+1);
FEValues<dim> fe_values (fe, quadrature_formula,
update_gradients | update_JxW_values|
- update_quadrature_points);
+ update_quadrature_points);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
Vector<double> cell_src (dofs_per_cell),
- cell_dst (dofs_per_cell);
+ cell_dst (dofs_per_cell);
const Coefficient<dim> coefficient;
std::vector<double> coefficient_values(n_q_points);
for (unsigned int q=0; q<n_q_points; ++q)
for (unsigned int i=0; i<dofs_per_cell; ++i)
- for (unsigned int j=0; j<dofs_per_cell; ++j)
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
cell_matrix(i,j) += (fe_values.shape_grad(i,q) *
fe_values.shape_grad(j,q) *
fe_values.JxW(q)*
- coefficient_values[q]);
+ coefficient_values[q]);
cell->get_dof_indices (local_dof_indices);
@f}
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_points</code> rows and @p dofs_per_cell columns). The matrix
+<code>fe_values.shape_grad(i,q)[d]</code>. This matrix consists of
+<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) * coefficient_values[q]</code> (or, rather, @p
dim copies of each of these values). This kind of representation of
finite element matrices can often be found in the engineering literature.
-When the cell-based matrix is applied to a vector @f{eqnarray*}
+When the cell matrix is applied to a vector,
+@f{eqnarray*}
A_\mathrm{cell}\cdot u_\mathrm{cell} = B_\mathrm{cell}^T
-D_\mathrm{cell} B_\mathrm{cell} \cdot u_\mathrm{cell}, @f} one would
-then 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. This removed the three nested
-loops in the calculation of the local matrix. What happens is as
-follows: We first transform the vector of values on the local dofs to
-a vector of gradients on the quadrature points. In the second loop, we
-multiply these gradients by the integration weight. The third loop
-applies the second gradient (in transposed form), so that we get back
-to a vector of (Laplacian) values on the cell dofs. This reduces the
+D_\mathrm{cell} B_\mathrm{cell} \cdot u_\mathrm{cell},
+@f}
+one would then not form the matrix-matrix products, but rather multiply one
+matrix at a time with a vector from right to left so that only three
+successive matrix-vector products are formed. This approach removes the three
+nested loops in the calculation of the local matrix, which reduces the
complexity of the work on one cell from something like $\mathcal
{O}(\mathrm{dofs\_per\_cell}^3)$ to $\mathcal
-{O}(\mathrm{dofs\_per\_cell}^2)$.
+{O}(\mathrm{dofs\_per\_cell}^2)$. An interpretation of this algorithm is that
+we first transform the vector of values on the local DoFs to a vector of
+gradients on the quadrature points. In the second loop, we multiply these
+gradients by the integration weight. The third loop applies the second
+gradient (in transposed form), so that we get back to a vector of (Laplacian)
+values on the cell dofs.
The bottleneck in the above code is the operations done by the call to
-FEValues::reinit for every <code>cell</code>, which take about as much time as the
-other steps together (at least if the mesh is unstructured; deal.II can
+FEValues::reinit for every <code>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
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. If you only build the matrix once as we've done in all
-previous tutorial programs, there is nothing one can do about the need
-to call FEValues::reinit on every cell since this transformation has
-to be done when we want to compute the local matrix elements.
-
-However, in a matrix-free implementation, we are not interested in
-applying the matrix only once. Rather, in iterative solvers, we need
-to expect that we have to apply the matrix many times, and so we can
-think about whether we may be able to cache something between
-different applications. On the other hand, we realize that we must not
-cache too much data since otherwise we get back to the situation where
-memory access becomes the dominating factor.
-
-The trick is now to factor out the Jacobian transformation and first
-apply the gradient on the reference cell only. 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 the
-transposed Jacobian for preparing the third loop which again uses the
-gradients on the unit cell.
+general. If you only build the matrix once as we've done in all previous
+tutorial programs, there is nothing to be optimized since FEValues::reinit
+needs to be called on every cell. In this process, the transformation is
+applied while computing the local matrix elements.
+
+In a matrix-free implementation, however, we will compute those integrals much
+more often than once. For example, iterative solvers will apply the matrix
+many times, and so we can think about whether we may be able to cache
+something between different operator applications, i.e., integral
+computations. On the other hand, we realize that we must not cache too much
+data since otherwise we get back to the situation where memory access becomes
+the dominating factor.
+
+The trick is to factor out the Jacobian transformation and first apply the
+gradient on the reference cell only. 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, apply the
+weights of the quadrature, and finally apply 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
containing the values on the quadrature points. It is constructed by a
-matrix-matrix product as
-@f{eqnarray*}
-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> denotes the Jacobian transformation from unit to
-real cell (in the language of transformations, the operation represented by
-<i>J</i><sub>cell</sub> represents a covariant
-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.
+matrix-matrix product as @f{eqnarray*} B_\mathrm{cell} =
+J_\mathrm{cell}^{-\mathrm T} B_\mathrm{ref\_cell}, @f} where
+<i>B</i><sub>ref_cell</sub> denotes the gradient on the reference cell and
+<i>J</i><sup>-T</sup><sub>cell</sub> denotes the inverse transpose Jacobian of
+the transformation from unit to real cell (in the language of transformations,
+the operation represented by <i>J</i><sup>-T</sup><sub>cell</sub> represents a
+covariant transformation). <i>J</i><sup>-T</sup><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*}
A_\mathrm{cell} = B_\mathrm{cell}^T D B_\mathrm{cell}
- = B_\mathrm{ref\_cell}^T J_\mathrm{cell}^T
- D_\mathrm{cell}
- J_\mathrm{cell} B_\mathrm{ref\_cell},
+ = B_\mathrm{ref\_cell}^T J_\mathrm{cell}^{-1}
+ D_\mathrm{cell}
+ J_\mathrm{cell}^{-\mathrm T} B_\mathrm{ref\_cell},
@f}
so we calculate the product (starting the local product from the right)
@f{eqnarray*}
-v_\mathrm{cell} = B_\mathrm{ref\_cell}^T J_\mathrm{cell}^T D J_\mathrm{cell}
+v_\mathrm{cell} = B_\mathrm{ref\_cell}^T J_\mathrm{cell}^{-1} D J_\mathrm{cell}^{-\mathrm T}
B_\mathrm{ref\_cell} u_\mathrm{cell}, \quad
v = \sum_{\mathrm{cell}=1}^{\mathrm{n\_cells}} P_\mathrm{cell,{loc-glob}}^T
v_\mathrm{cell}.
@code
...
FEValues<dim> fe_values_reference (fe, quadrature_formula,
- update_gradients);
+ update_gradients);
Triangulation<dim> reference_cell;
GridGenerator::hyper_cube(reference_cell, 0., 1.);
fe_values_reference.reinit (reference_cell.begin());
FEValues<dim> fe_values (fe, quadrature_formula,
update_inverse_jacobians | update_JxW_values |
- update_quadrature_points);
+ update_quadrature_points);
for (; cell!=endc; ++cell)
{
for (unsigned int q=0; q<n_q_points; ++q)
for (unsigned int d=0; d<dim; ++d)
for (unsigned int i=0; i<dofs_per_cell; ++i)
- temp_vector(q*dim+d) +=
- fe_values_reference.shape_grad(i,q)[d] * cell_src(i);
+ temp_vector(q*dim+d) +=
+ fe_values_reference.shape_grad(i,q)[d] * cell_src(i);
for (unsigned int q=0; q<n_q_points; ++q)
{
// apply the transposed inverse Jacobian of the mapping
- Tensor<1,dim> temp;
- for (unsigned int d=0; d<dim; ++d)
- temp[d] = temp_vector(q*dim+d);
- for (unsigned int d=0; d<dim; ++d)
- {
- double sum = 0;
- for (unsigned int e=0; e<dim; ++e)
- sum += fe_values.inverse_jacobian(q)[e][d] *
- temp[e];
- temp_vector(q*dim+d) = sum;
- }
+ Tensor<1,dim> temp;
+ for (unsigned int d=0; d<dim; ++d)
+ temp[d] = temp_vector(q*dim+d);
+ for (unsigned int d=0; d<dim; ++d)
+ {
+ double sum = 0;
+ for (unsigned int e=0; e<dim; ++e)
+ sum += fe_values.inverse_jacobian(q)[e][d] *
+ temp[e];
+ temp_vector(q*dim+d) = sum;
+ }
// multiply by coefficient and integration weight
- for (unsigned int d=0; d<dim; ++d)
- temp_vector(q*dim+d) *= fe_values.JxW(q) * coefficient_values[q];
+ for (unsigned int d=0; d<dim; ++d)
+ temp_vector(q*dim+d) *= fe_values.JxW(q) * coefficient_values[q];
// apply the inverse Jacobian of the mapping
- for (unsigned int d=0; d<dim; ++d)
- temp[d] = temp_vector(q*dim+d);
- for (unsigned int d=0; d<dim; ++d)
- {
- double sum = 0;
- for (unsigned int e=0; e<dim; ++e)
- sum += fe_values.inverse_jacobian(q)[d][e] *
- temp[e];
- temp_vector(q*dim+d) = sum;
- }
+ for (unsigned int d=0; d<dim; ++d)
+ temp[d] = temp_vector(q*dim+d);
+ for (unsigned int d=0; d<dim; ++d)
+ {
+ double sum = 0;
+ for (unsigned int e=0; e<dim; ++e)
+ sum += fe_values.inverse_jacobian(q)[d][e] *
+ temp[e];
+ temp_vector(q*dim+d) = sum;
+ }
}
cell_dst = 0;
for (unsigned int i=0; i<dofs_per_cell; ++i)
for (unsigned int q=0; q<n_q_points; ++q)
for (unsigned int d=0; d<dim; ++d)
- cell_dst(i) += fe_values_reference.shape_grad(i,q)[d] *
- temp_vector(q*dim+d);
+ cell_dst(i) += fe_values_reference.shape_grad(i,q)[d] *
+ temp_vector(q*dim+d);
for (unsigned int i=0; i<dofs_per_cell; ++i)
dst(local_dof_indices(i)) += cell_dst(i);
Note how we create an additional FEValues object for the reference cell
gradients and how we initialize it to the reference cell. The actual
derivative data is then applied by the inverse, transposed Jacobians (deal.II
-calls the Jacobian matrix from unit to real cell inverse_jacobian, because the
-transformation direction in deal.II is from real to unit cell).
+calls the Jacobian matrix from real to unit cell inverse_jacobian, as the
+forward transformation is from unit to real cell). Note that the factor
+$J_\mathrm{cell}^{-1} D_\mathrm{cell} J_\mathrm{cell}^{-\mathrm T}$ is
+block-diagonal over quadrature. In this form, one realizes that variable
+coefficients (possibly expressed through a tensor) and general grid topologies
+with Jacobian transformations have a similar effect on the coefficient that
+transforms the unit-cell derivatives.
Finally, we are using tensor product basis functions and now that we have
separated out the gradient on the reference cell <i>B</i><sub>ref_cell</sub>,
\varphi_j(y)$, we get $(B_\mathrm{grad,x} \otimes
B_\mathrm{val,y})u_\mathrm{cell} = B_\mathrm{val,y} U B_\mathrm{grad,x}$. This
reduces the complexity for computing this product from $p^4$ to $2 p^3$, where
-<i>p</i>-1 is the degree of the finite element (i.e., equivalently,
-<i>p</i> is the number of shape functions in each coordinate
-direction), or $p^{2d}$ to $d p^{d+1}$ in general.
+<i>p</i>-1 is the degree of the finite element (i.e., equivalently, <i>p</i>
+is the number of shape functions in each coordinate direction), or $p^{2d}$ to
+$d p^{d+1}$ in general. Note that the concentration on the polynomial degree
+in terms of complexity is relevant as the polynomial degree is increased
+rather than the grid resolution. Good algorithms for moderate degrees like the
+one used here are linear in the polynomial degree independent on the
+dimension, as opposed to matrix-based schemes or naive evaluation through
+FEValues. These techniques have been established in the spectral element
+community since the 1980s.
Implementing a matrix-free and cell-based finite element operator requires a
somewhat different design compared to the usual matrix assembly codes shown in
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 number of
-matrix-vector products — 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 SSOR, ILU or
-algebraic multigrid (AMG) cannot 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 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. 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. 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 <a
+however, one wants more than just performing some number of matrix-vector
+products — 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. Unfortunately, most of the more
+frequently used preconditioners such as SSOR, ILU or algebraic multigrid (AMG)
+cannot be used here because their implementation requires knowledge of the
+elements of the system matrix.
+
+One solution is to use geometric multigrid methods as shown in 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 related to the a
+colletion of cells. 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. A Chebyshev iteration around the Jacobi method, eventually, is our
+choice in this tutorial program. 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 <a
href="http://www.sciencedirect.com/science/article/pii/S0021999103001943">
<i>M. Adams, M. Brezina, J. Hu, R. Tuminaro. Parallel multigrid smoothers:
polynomial versus Gauss–Seidel, J. Comput. Phys. 188:593–610,
Chebyshev smoothers that we exploit here, namely that they are easy to
parallelize, whereas SOR/Gauss–Seidel smoothing relies on substitutions,
for which a naive parallelization works on diagonal sub-blocks of the matrix,
-thereby decreases efficiency (for more detail see e.g. Y. Saad,
-Iterative Methods for Sparse Linear Systems, SIAM, 2nd edition, 2003, chapters
-11 & 12).
+thereby decreases efficiency (for more detail see e.g. Y. Saad, Iterative
+Methods for Sparse Linear Systems, SIAM, 2nd edition, 2003, chapters 11 & 12).
The implementation into the multigrid framework is then straightforward. The
-multigrid implementation in this program is based on a simplified version of
-step-16 that disregards adaptivity.
+multigrid implementation in this program is similar to step-16 and includes
+adaptivity.
<h3>Using CPU-dependent instructions (vectorization)</h3>
The computational kernels for evaluation in FEEvaluation are written in a way
-to optimally use computational resources. Indeed, they operate not on double
-data types, but something we call VectorizedArray (check e.g. the return type
-of FEEvaluationBase::get_value, which is VectorizedArray for a scalar element
-and a Tensor of VectorizedArray for a vector finite element). VectorizedArray
-is a short array of doubles or float whose length depends on the particular
-computer system in use. For example, systems based on x86-64 support the
-streaming SIMD extensions (SSE), where the processor's vector units can
-process two doubles (or four single-precision floats) by one CPU
-instruction. Newer processors with support for the so-called advanced vector
-extensions (AVX) with 256 bit operands can use four doubles and eight floats,
-respectively. Vectorization is a single-instruction/multiple-data (SIMD)
-concept, that is, one CPU instruction is used to process multiple data values
-at once. Often, finite element programs do not use vectorization explicitly as
-the benefits of this concept are only in arithmetic intensive operations. The
-bulk of typical finite element workloads are memory bandwidth limited
-(operations on sparse matrices and vectors) where the additional computational
-power is useless.
+to optimally use computational resources. To achieve this, they do not operate
+on double data types, but something we call VectorizedArray (check e.g. the
+return type of FEEvaluationBase::get_value, which is VectorizedArray for a
+scalar element and a Tensor of VectorizedArray for a vector finite
+element). VectorizedArray is a short array of doubles or float whose length
+depends on the particular computer system in use. For example, systems based
+on x86-64 support the streaming SIMD extensions (SSE), where the processor's
+vector units can process two doubles (or four single-precision floats) by one
+CPU instruction. Newer processors (from about 2012 and onwards) support the
+so-called advanced vector extensions (AVX) with 256 bit operands, which can
+use four doubles and eight floats, respectively. Vectorization is a
+single-instruction/multiple-data (SIMD) concept, that is, one CPU instruction
+is used to process multiple data values at once. Often, finite element
+programs do not use vectorization explicitly as the benefits of this concept
+are only in arithmetic intensive operations. The bulk of typical finite
+element workloads are memory bandwidth limited (operations on sparse matrices
+and vectors) where the additional computational power is useless.
Behind the scenes, optimized BLAS packages might heavily rely on
vectorization, though. Also, optimizing compilers might automatically
-transform loops involving standard code into more efficient vectorized
-form. However, the data flow must be very regular in order for compilers to
-produce efficient code. For example, already the automatic vectorization of
+transform loops involving standard code into more efficient vectorized form
+(deal.II uses OpenMP SIMD pragmas inside the regular loops of vector
+updates). However, the data flow must be very regular in order for compilers
+to produce efficient code. For example, already the automatic vectorization of
the prototype operation that benefits from vectorization, matrix-matrix
-products, fails on most compilers (as of writing this tutorial in early 2012,
-neither gcc-4.6 nor the Intel compiler v. 12 manage to produce useful
-vectorized code for the FullMatrix::mmult function, and not even on the
-simpler case where the matrix bounds are compile-time constants instead of
-run-time constants as in FullMatrix::mmult). The main reason for this is that
-the information to be processed at the innermost loop (that is where
-vectorization is applied) is not necessarily a multiple of the vector length,
-leaving parts of the resources unused. Moreover, the data that can potentially
-be processed together might not be laid out in a contiguous way in memory or
-not with the necessary alignment to address boundaries that are needed by the
-processor. Or the compiler might not be able to prove that data arrays do not
-overlap when loading several elements at once.
+products, fails on most compilers (as of writing this tutorial in early 2012
+and updating in late 2016, neither gcc nor the Intel compiler manage to
+produce useful vectorized code for the FullMatrix::mmult function, and not
+even on the simpler case where the matrix bounds are compile-time constants
+instead of run-time constants as in FullMatrix::mmult). The main reason for
+this is that the information to be processed at the innermost loop (that is
+where vectorization is applied) is not necessarily a multiple of the vector
+length, leaving parts of the resources unused. Moreover, the data that can
+potentially be processed together might not be laid out in a contiguous way in
+memory or not with the necessary alignment to address boundaries that are
+needed by the processor. Or the compiler might not be able to prove that data
+arrays do not overlap when loading several elements at once.
In the matrix-free implementation in deal.II, we have therefore chosen to
apply vectorization at the level which is most appropriate for finite element
computations: The cell-wise computations are typically exactly the same for
-all cells (except for reading from and writing to vectors), and hence SIMD can
-be used to process several cells at once. In all what follows, you can think
-of a VectorizedArray to hold data from several cells. Remember that it is not
-related to the spatial dimension and the number of elements e.g. in a Tensor
-or Point.
+all cells (except for indices in the indirect addressing used while reading
+from and writing to vectors), and hence SIMD can be used to process several
+cells at once. In all what follows, you can think of a VectorizedArray to hold
+data from several cells. Remember that it is not related to the spatial
+dimension and the number of elements e.g. in a Tensor or Point.
Note that vectorization depends on the CPU that is used for deal.II. In order
to generate the fastest kernels of FEEvaluation for your computer, you should
<tt>-DCMAKE_CXX_FLAGS="-march=native"</tt>, see the deal.II README for more
information). Similar options exist for other compilers.
+
+<h3>Running multigrid on large-scale parallel computers</h3>
+
+As mentioned above, all components in the matrix-free framework can easily be
+parallelized with MPI using domain decomposition. Thanks to the good support
+of large-scale parallel mesh framework provided by p4est (see step-40 for
+details) in deal.II, and the fact that cell-based loops with matrix-free
+evaluation <i>only</i> need a decomposition of the mesh into chunks of roughly
+equal size on each processor, there is relatively little to do for extending
+the program to MPI. While other tutorial programs using MPI have relied on
+either PETSc or Trilinos, this program uses deal.II's own parallel vector
+facilities.
+
+The deal.II parallel vector class, LinearAlgebra::distributed::Vector, holds
+the processor-local part of the solution as well as information on and data
+fields for the ghost DoFs, i.e. DoFs that are owned by a remote processor but
+needed on cells that are treated by the present processor. Moreover, it holds
+the MPI-send information for DoFs that are owned locally but needed by other
+processors. A benefit of the design of that vector class is the way ghosted
+entries are accessed. In the storage scheme of the vector, the array of data
+does extend beyond the processor-local part of the solution with further
+vector entries available for the ghost degrees of freedom by a continguous
+index range. (Note that the index range depends on the exact configuration of
+the mesh.) Moreover, it holds the MPI-send information for DoFs that are owned
+locally but needed by other processors. Since matrix-free operations can be
+thought of doing linear algebra that is performance critical, and
+performance-critical code cannot waste time on doing MPI-global to MPI-local
+index translations, the availability of such access functionality in index
+spaces local to one MPI rank is fundamental. The way things are accessed here
+is a direct array access. A user can use the query
+LinearAlgebra::distributed::Vector::local_element for that purpose, but it is
+actually rarely used because all of this happens internally in FEEvaluation.
+
+The design of LinearAlgebra::distributed::Vector is similar to the
+PETScWrappers::MPI::Vector and TrilinosWrappers::MPI::Vector data types we
+have used in step-40 and step-32 before, but since we do not need any other
+parallel functionality of these libraries, we use the
+parallel::distributed::Vector class of deal.II instead of linking in another
+large library in this tutorial program. Also note that the PETSc and Trilinos
+vectors do not provide the fine-grained control over ghost entries with direct
+array access because they abstract away the necessary implementation details.
Of more interest is to evaluate some aspects of the multigrid solver.
When we run this program in 2D for quadratic ($Q_2$) elements, we get the
-following output:
+following output (when run on one core in release mode):
@code
Cycle 0
Number of degrees of freedom: 81
-System matrix memory consumption: 0.008982 MB.
-Multigrid objects memory consumption: 0.02617 MB.
-Total setup time (wall) 0.001811s
-Time solve (5 iterations) (CPU/wall) 0s/0.0002651s
+Total setup time (wall) 0.00159788s
+Time solve (6 iterations) (CPU/wall) 0.000951s/0.000951052s
Cycle 1
Number of degrees of freedom: 289
-System matrix memory consumption: 0.01817 MB.
-Multigrid objects memory consumption: 0.05779 MB.
-Total setup time (wall) 0.001223s
-Time solve (5 iterations) (CPU/wall) 0s/0.000926s
+Total setup time (wall) 0.00114608s
+Time solve (6 iterations) (CPU/wall) 0.000935s/0.000934839s
Cycle 2
Number of degrees of freedom: 1089
-System matrix memory consumption: 0.05286 MB.
-Multigrid objects memory consumption: 0.1581 MB.
-Total setup time (wall) 0.003045s
-Time solve (6 iterations) (CPU/wall) 0.012s/0.003393s
+Total setup time (wall) 0.00244665s
+Time solve (6 iterations) (CPU/wall) 0.00207s/0.002069s
Cycle 3
Number of degrees of freedom: 4225
-System matrix memory consumption: 0.1957 MB.
-Multigrid objects memory consumption: 0.5228 MB.
-Total setup time (wall) 0.008561s
-Time solve (6 iterations) (CPU/wall) 0.02s/0.01133s
+Total setup time (wall) 0.00678205s
+Time solve (6 iterations) (CPU/wall) 0.005616s/0.00561595s
Cycle 4
Number of degrees of freedom: 16641
-System matrix memory consumption: 0.7343 MB.
-Multigrid objects memory consumption: 1.925 MB.
-Total setup time (wall) 0.02938s
-Time solve (6 iterations) (CPU/wall) 0.068s/0.03312s
+Total setup time (wall) 0.0241671s
+Time solve (6 iterations) (CPU/wall) 0.019543s/0.0195441s
Cycle 5
Number of degrees of freedom: 66049
-System matrix memory consumption: 2.856 MB.
-Multigrid objects memory consumption: 7.435 MB.
-Total setup time (wall) 0.1128s
-Time solve (6 iterations) (CPU/wall) 0.228s/0.09577s
+Total setup time (wall) 0.0967851s
+Time solve (6 iterations) (CPU/wall) 0.07457s/0.0745709s
Cycle 6
Number of degrees of freedom: 263169
-System matrix memory consumption: 11.28 MB.
-Multigrid objects memory consumption: 29.3 MB.
-Total setup time (wall) 0.4553s
-Time solve (6 iterations) (CPU/wall) 1.272s/0.3955s
+Total setup time (wall) 0.346374s
+Time solve (6 iterations) (CPU/wall) 0.260042s/0.265033s
@endcode
As in step-16, we see that the number of CG iterations remains constant with
-increasing number of degrees of freedom. We can also see that the various
-objects we have to store for the multigrid method on the individual levels of
-our mesh together make up more than twice as much as the matrix on the finest
-level. For the present example, about half the memory consumption of the
-multigrid objects are the level transfer matrices, and the other half is
-consumed by the matrix-free objects (and there, mainly the indices and the
-variable coefficient).
+increasing number of degrees of freedom. The constant number of iterations
+(together with optimal computational properties) means that the computing time
+approximately quadruples as the problem size quadruples from one cycle to the
+next. A look at the memory consumption reveals that the code is also very
+efficient in terms of storage. Around 2-4 million degrees of freedom fit into
+1 GB of memory. An interesting fact is that solving one linear system is
+cheaper than the setup, despite not building a matrix (approximately half of
+which is spent in the DoFHandler::distribute_dofs() and
+DoFHandler::distributed_mg_dofs() calls). This shows the high efficiency of
+this approach, but also the fact that most benefit will be obtained when
+solving multiple systems on the same mesh and DoFHandler structure.
Not much changes if we run the program in three spatial dimensions, with the
-exception that the multilevel objects now take up some more memory (because
-the level transfer matrices are denser) and the computing times are somewhat
-larger:
+exception that problem sizes grow by a factor eight when refining the mesh and
+the obvious increase in computing times:
@code
Cycle 0
Number of degrees of freedom: 125
-System matrix memory consumption: 0.01093 MB.
-Multigrid objects memory consumption: 0.03094 MB.
-Total setup time (wall) 0.002481s
-Time solve (5 iterations) (CPU/wall) 0s/0.000334s
+Total setup time (wall) 0.00231099s
+Time solve (6 iterations) (CPU/wall) 0.000692s/0.000922918s
Cycle 1
Number of degrees of freedom: 729
-System matrix memory consumption: 0.04105 MB.
-Multigrid objects memory consumption: 0.1274 MB.
-Total setup time (wall) 0.004471s
-Time solve (5 iterations) (CPU/wall) 0.004s/0.001979s
+Total setup time (wall) 0.00289083s
+Time solve (6 iterations) (CPU/wall) 0.001534s/0.0024128s
Cycle 2
Number of degrees of freedom: 4913
-System matrix memory consumption: 0.2821 MB.
-Multigrid objects memory consumption: 0.8048 MB.
-Total setup time (wall) 0.01651s
-Time solve (4 iterations) (CPU/wall) 0.036s/0.01295s
+Total setup time (wall) 0.0143182s
+Time solve (6 iterations) (CPU/wall) 0.010785s/0.0107841s
Cycle 3
Number of degrees of freedom: 35937
-System matrix memory consumption: 1.948 MB.
-Multigrid objects memory consumption: 5.734 MB.
-Total setup time (wall) 0.1072s
-Time solve (5 iterations) (CPU/wall) 0.16s/0.0709s
+Total setup time (wall) 0.087064s
+Time solve (6 iterations) (CPU/wall) 0.063522s/0.06545s
Cycle 4
Number of degrees of freedom: 274625
-System matrix memory consumption: 14.49 MB.
-Multigrid objects memory consumption: 44.41 MB.
-Total setup time (wall) 0.8173s
-Time solve (5 iterations) (CPU/wall) 1.52s/0.5093s
+Total setup time (wall) 0.596306s
+Time solve (6 iterations) (CPU/wall) 0.427757s/0.431765s
Cycle 5
Number of degrees of freedom: 2146689
-System matrix memory consumption: 115.9 MB.
-Multigrid objects memory consumption: 342.6 MB.
-Total setup time (wall) 6.387s
-Time solve (5 iterations) (CPU/wall) 12.45s/3.767s
+Total setup time (wall) 4.96491s
+Time solve (6 iterations) (CPU/wall) 3.53126s/3.56142s
@endcode
+Since it is so easy, we look at what happens if we increase the polynomial
+degree. When selecting the degree as four in 3D, i.e., on $\mathcal Q_4$
+elements, we get the following program output:
+
+@code
+Cycle 0
+Number of degrees of freedom: 729
+Total setup time (wall) 0.00633097s
+Time solve (6 iterations) (CPU/wall) 0.002829s/0.00379395s
+
+Cycle 1
+Number of degrees of freedom: 4913
+Total setup time (wall) 0.0174279s
+Time solve (6 iterations) (CPU/wall) 0.012255s/0.012254s
+
+Cycle 2
+Number of degrees of freedom: 35937
+Total setup time (wall) 0.082655s
+Time solve (6 iterations) (CPU/wall) 0.052362s/0.0523629s
+
+Cycle 3
+Number of degrees of freedom: 274625
+Total setup time (wall) 0.507943s
+Time solve (6 iterations) (CPU/wall) 0.341811s/0.345788s
+
+Cycle 4
+Number of degrees of freedom: 2146689
+Total setup time (wall) 3.46251s
+Time solve (7 iterations) (CPU/wall) 3.29638s/3.3265s
+
+Cycle 5
+Number of degrees of freedom: 16974593
+Total setup time (wall) 27.8989s
+Time solve (7 iterations) (CPU/wall) 26.3705s/27.1077s
+@endcode
+
+Since $\mathcal Q_4$ elements on a mesh correspond to $\mathcal Q_2$ elements
+on half the mesh size, we can start to compare the run time at cycle 4 with
+fourth degree polynomials with cycle 5 at quadratic polynomials, both at 2.1
+million degrees of freedom. The surprising effect is that the solver for
+$\mathcal Q_4$ element is actually slightly faster than for the quadratic
+case, despite using one more linear iteration. The effect that higher-degree
+polynomials are similarly fast or even faster than lower degree ones is one of
+the main strengths of matrix-free operator evaluation through sum
+factorization, see the <a
+href="http://dx.doi.org/10.1016/j.compfluid.2012.04.012">matrix-free
+paper</a>. This is fundamentally different to matrix-based methods that get
+more expensive as the degree increases and the coupling gets denser.
+
+The improved efficiency is also reflected in the initialization routines.
+
+Finally, let us look at the timings with degree 8, which corresponds to
+another round of mesh refinement in the lower order methods:
+
+@code
+Cycle 0
+Number of degrees of freedom: 4913
+Total setup time (wall) 0.0842004s
+Time solve (8 iterations) (CPU/wall) 0.019296s/0.0192959s
+
+Cycle 1
+Number of degrees of freedom: 35937
+Total setup time (wall) 0.327048s
+Time solve (8 iterations) (CPU/wall) 0.07517s/0.075999s
+
+Cycle 2
+Number of degrees of freedom: 274625
+Total setup time (wall) 2.12335s
+Time solve (8 iterations) (CPU/wall) 0.448739s/0.453698s
+
+Cycle 3
+Number of degrees of freedom: 2146689
+Total setup time (wall) 16.1743s
+Time solve (8 iterations) (CPU/wall) 3.95003s/3.97717s
+
+Cycle 4
+Number of degrees of freedom: 16974593
+Total setup time (wall) 130.8s
+Time solve (8 iterations) (CPU/wall) 31.0316s/31.767s
+@endcode
+
+Here, the initialization seems considerably slower than before, which is
+mainly due to the computation of the diagonal of the matrix, which actually
+computes a 729 x 729 matrix on each cell and throws away everything but the
+diagonal. The solver times, however, are again very close to the quartic case,
+showing that the linear increase with the polynomial degree that is
+theoretically expected is offset by better computational characteristics and
+the fact that higher order methods have a smaller share of degrees of freedom
+living on several cells that slow the method down.
<h3>Comparison with a sparse matrix</h3>
</tr>
</table>
-The table clearly shows that the matrix-free implementation is twice as fast
-for the solver, and more than six times as fast when it comes to
+The table clearly shows that the matrix-free implementation is more than twice
+as fast for the solver, and more than six times as fast when it comes to
initialization costs. As the problem size is made a factor 8 larger, we note
that the times usually go up by a factor eight, too (as the solver iterations
-are constant at 5). There are two deviations. The first is in the sparse
-matrix between 5k and 36k degrees of freedom, where the time increases by a
-factor 12. This is the threshold when the cache in the processor can no longer
-hold all data necessary for the matrix-vector products and all matrix elements
-must be fetched from main memory. The second deviation is the times for the
-matrix-free solve which increase by less than a factor 8. This is because of
-more parallelism from more cells, exploited by the (involved) dynamic task
-scheduling approach taken in the cell loop of the MatrixFree class. Note
-that about 30% of the time in the
-matrix-free solver is spent on restriction and prolongation, which use sparse
-matrices. So the speedup could be even better if all parts were done
-efficiently.
+are constant at six). The main deviation is in the sparse matrix between 5k
+and 36k degrees of freedom, where the time increases by a factor 12. This is
+the threshold when the cache in the processor can no longer hold all data
+necessary for the matrix-vector products and all matrix elements must be
+fetched from main memory.
Of course, this picture does not necessarily translate to all cases, as there
are problems where knowledge of matrix entries enables much better solvers (as
2<i>p<sup>d</sup></i> for the sparse matrix, so it will win anyway for order 4
and higher in 3d).
-<h3>Possibilities for extensions</h3>
+<h3> Results for large-scale parallel computations</h3>
-Above, we have shown figures for second-order finite elements. Our
-implementation gains more compared to sparse matrices if higher order elements
-are used. However, FE_Q elements with equidistant nodes are badly conditioned
-if the order increases. In this case, the smoother and the multigrid solver
-break down. Node clustering close to the element boundaries resolves this
-problem (and the multigrid solver converges in 5 or 6 iterations also for very
-high order). Elements with this properties are the Gauss-Lobatto FE_Q
-elements, which are presented in step-48.
+<h3> Results with adaptivity</h3>
projects / dealii.git / commitdiff