From ca69963b2c70ad5557062f67d5d826f11eed7056 Mon Sep 17 00:00:00 2001 From: kronbichler Date: Fri, 6 Nov 2009 23:06:37 +0000 Subject: [PATCH] Improve comments a bit further. git-svn-id: https://svn.dealii.org/trunk@20057 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-37/doc/intro.dox | 219 ++++++------ deal.II/examples/step-37/step-37.cc | 443 ++++++++++++------------- 2 files changed, 318 insertions(+), 344 deletions(-) diff --git a/deal.II/examples/step-37/doc/intro.dox b/deal.II/examples/step-37/doc/intro.dox index 40dafdf9c8..09d319ddae 100644 --- a/deal.II/examples/step-37/doc/intro.dox +++ b/deal.II/examples/step-37/doc/intro.dox @@ -30,9 +30,11 @@ product when the matrix is stored in the usual sparse compressed row storage — in short CRS — 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; rowA_values is continuously traveling from main memory into the -CPU in order to be multiplied with some vector element determined by the other -array A_column_indices 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 A_values, and 4 bytes for the unsigned -integer position A_column_indices 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 A_values[element_index] 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 (109) +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 — 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 A 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 Pcell,loc-glob 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 -local_dof_indices 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 local_dof_indices +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 x at the degrees of freedom -of the respective cell, and $y_\mathrm{cell}$ correspondingly for the result. +where xcell are the values of x at the degrees of freedom +of the respective cell, and xcell 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 &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 src of the cell-related -degrees of freedom, multiply by the local matrix, and finally add the result -to the destination vector dst. 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 src of the cell-related degrees of freedom +(the action of Pcell,loc-glob), multiply by the local matrix +(the action of Acell), and finally add the result to the +destination vector dst (the action of +Pcell,loc-globT, 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 Bcell is given by +where for the example of the Laplace operator the (q*dim+d,i)-th +element of Bcell is given by fe_values.shape_grad(i,q)[d]. The matrix consists of -dim*n_q_pointsdofs_per_cell rows and dofs_per_cell -columns). The matrix Dcell is diagonal and contains the values -fe_values.JxW(q) (or, rather, dim copies of it). +dim*n_q_points rows and @p dofs_per_cell columns). The matrix +Dcell is diagonal and contains the values +fe_values.JxW(q) (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 fe_values.reinit(cell) 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 -temp_vector 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 +fe_values.reinit(cell), 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 temp_vector 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 Bcell 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 Bref_cell denotes the gradient on the reference cell -and JcellT denotes the -Jacobian. JcellT 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 Jcell denotes the Jacobian +transformation. Jcell 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:
  • 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 JcellT D + operations JcellT Dcell Jcell (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).
  • 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 -src and the write operation to dst 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 graphics processing unit (GP-GPU), 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.

    Combination with multigrid

    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 — 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 diff --git a/deal.II/examples/step-37/step-37.cc b/deal.II/examples/step-37/step-37.cc index ff07b2d8e3..f545b378d4 100644 --- a/deal.II/examples/step-37/step-37.cc +++ b/deal.II/examples/step-37/step-37.cc @@ -108,7 +108,7 @@ void Coefficient::value_list (const std::vector > &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 – we choose a // Table<2,Transformation> data // format) – and call a transform - // command of the Transformation + // 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 class MatrixFree : public Subscriptor @@ -253,7 +253,7 @@ class MatrixFree : public Subscriptor std::size_t memory_consumption () const; // The private member variables of the - // MatrixFree 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 - // MatrixFree class. All it does - // is to subscribe to the general deal.II - // 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. + // 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 MatrixFree::MatrixFree () : @@ -360,24 +360,20 @@ MatrixFree::get_constraints () - // The following function takes a - // vector of local dof indices on - // cell level and writes the data - // into the - // 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 - // false. + // 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 void MatrixFree:: 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 - // 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. + // 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 void MatrixFree:: 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 - // cell_range. 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 Pcell,local-global + // 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, //
      //
    1. 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. //
    2. 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 + // x1 contains the + // reference cell gradient to the local + // cell vector. //
    3. 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 - // dim. 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 - // JxW values and the + // Dcell now contains the + // @p JxW values and the // inhomogeneous coefficient. //
    4. Form $y_\mathrm{cell} = // B_\mathrm{ref\_cell}^T x_2$. This gives // the local result of the matrix-vector // product. - //
    5. Form $y += + //
    6. 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 - // copy_local_to_global + // @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 y. //
    // 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 + // Bref_cell + // xcell in the second step + // and the transpose operation in step + // 4. Note that the matrix + // JT D J is + // block-diagonal, and hence, its application + // is cheaper. Since the matrix + // Bref_cell 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 xcell. 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 - // dofs_per_cell rows and - // columns. This is an operation that can be - // much better optimized than matrix-vector - // products. The functions - // FullMatrix::mmult and - // FullMatrix::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. 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::mmult and + // @p FullMatrix::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), ©.solutions(0,0), @@ -584,12 +583,10 @@ copy_local_to_global (const WorkStreamData::CopyData ©, - // Now to the vmult function - // that is called externally: It is very - // similar to the vmult_add - // 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 template void @@ -602,10 +599,11 @@ MatrixFree::vmult (Vector &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 template void @@ -629,60 +627,47 @@ MatrixFree::Tvmult_add (Vector &dst, - // This is the vmult_add - // function that multiplies the - // matrix with vector - // src and adds the - // result to vector 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 - // matrix_sizes.chunks. What - // the rather cryptic command to - // 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 - // _1, _2, _3 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 - // local_vmult function - // are bound: to this - // and a constant reference to the - // src in the first - // case, and this and a - // reference to the output vector in - // the second. Similarly, the - // placeholder _1 in the - // copy_local_to_global - // function sets the first explicit - // argument of that function, which - // is of class - // 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. + // 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 _1, _2, _3 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 template void @@ -732,27 +717,24 @@ MatrixFree::vmult_add (Vector &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 - // B_ref_cell). The number - // of quadrature points needs to be passed - // through the last variable - // n_points_per_cell, since - // the number of columns in the small - // matrix is - // dim*n_points_per_cell for - // the Laplace problem (the Laplacian is a - // tensor and has dim - // 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::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 double, + // and with data type @p double, // about 80 per cent of the memory // consumption is due to the - // derivatives array, while in 3D + // @p derivatives array, while in 3D // this number is even 85 per cent. template std::size_t MatrixFree::memory_consumption () const @@ -973,38 +955,34 @@ std::size_t MatrixFree::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 - // transform operation needed - // by the MatrixFree - // 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 - // double %numbers. Since we - // also want to use float - // %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 - // 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. + // 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 class LaplaceOperator { @@ -1047,21 +1025,20 @@ LaplaceOperator::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 vmult + // 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 if clause at this - // position in the code at the innermost - // loop, but note once again that - // dim is known when this piece - // of code is compiled, so the compiler can - // optimize away the if statement - // (and actually even inline these few lines - // of code into the MatrixFree + // 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 void LaplaceOperator::transform (number* result) const @@ -1090,7 +1067,7 @@ void LaplaceOperator::transform (number* result) const // The final function in this group // takes the content of a rank-2 // tensor and writes it to the field - // transformation 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::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 - // 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. + // 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 void LaplaceProblem::assemble_system () { @@ -1621,7 +1596,7 @@ void LaplaceProblem::run () - // @sect3{The main function} + // @sect3{The @p main function} // This is as in all other programs: int main () -- 2.39.5