From 922d903b5051eff2cb418fad521532756f6802e8 Mon Sep 17 00:00:00 2001 From: wolf Date: Wed, 17 May 2000 12:22:49 +0000 Subject: [PATCH] When eliminating boundary values from a matrix, we only have to eliminate the respective element from the other rows for which we know that they have a nonzero element in that column. Previously, we searched all rows for this, but if the sparsity pattern is symmetric then we can do better than that by checking the nonzero elements of the present row and only visiting the respective columns. This reduces the complexity from N*log(m) to m*log(m), where N=size of the matrix and m=number of elements per row, i.e. from O(N) to O(1). In the case of non-symmetric sparsity patterns, it might not be worth to eliminate the respective row, since we have to use a solver for nonsymmetric problems anyway, and that can then eliminate the respective column iteratively. git-svn-id: https://svn.dealii.org/trunk@2874 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/dofs/dof_tools.h | 26 +++- deal.II/deal.II/include/numerics/matrices.h | 87 +++++++---- deal.II/deal.II/source/numerics/base.cc | 3 +- deal.II/deal.II/source/numerics/matrices.cc | 162 +++++++++++++------- deal.II/deal.II/source/numerics/vectors.cc | 3 +- 5 files changed, 185 insertions(+), 96 deletions(-) diff --git a/deal.II/deal.II/include/dofs/dof_tools.h b/deal.II/deal.II/include/dofs/dof_tools.h index b5165e8400..af6dc71a7a 100644 --- a/deal.II/deal.II/include/dofs/dof_tools.h +++ b/deal.II/deal.II/include/dofs/dof_tools.h @@ -83,7 +83,12 @@ class DoFTools * #make_sparsity_pattern# just * loops over all cells and * enters all couplings local to - * that cell. + * that cell. As the generation + * of the sparsity pattern is + * irrespective of the equation + * which is solved later on, the + * resulting sparsity pattern is + * symmetric. * * Since this process is purely * local, the sparsity pattern @@ -157,13 +162,18 @@ class DoFTools * than to the degrees of freedom * contained in there. * - * This function is designed to accept - * a mask, like the one shown above, - * through the #mask# parameter, which - * contains boolean values. It builds - * the matrix structure just like the - * previous function, but does not create - * elements if not specified by the mask. + * This function is designed to + * accept a mask, like the one + * shown above, through the + * #mask# parameter, which + * contains boolean values. It + * builds the matrix structure + * just like the previous + * function, but does not create + * elements if not specified by + * the mask. If the mask is + * symmetric, then so will be the + * resulting sparsity pattern. * * The actual type of the * sparsity pattern may be diff --git a/deal.II/deal.II/include/numerics/matrices.h b/deal.II/deal.II/include/numerics/matrices.h index 37c761551a..f56fae4baa 100644 --- a/deal.II/deal.II/include/numerics/matrices.h +++ b/deal.II/deal.II/include/numerics/matrices.h @@ -352,9 +352,9 @@ class MatrixCreator * * \subsection{Boundary conditions} * - * The #apply_boundar_values# function inserts boundary conditions of + * The #apply_boundary_values# function inserts boundary conditions of * into a system of equations. To actually do this you have to specify - * a list of degree of freedom indices along with the value this degree of + * a list of degree of freedom indices along with the values these degrees of * freedom shall assume. To see how to get such a list, see the discussion * of the #VectorTools::interpolate_boundary_values# function. * @@ -371,6 +371,43 @@ class MatrixCreator * eliminating all coupling between this degree of freedom and others. Now * also the column consists only of zeroes, apart from the main diagonal entry. * + * Finding which rows contain an entry in the column for which we are + * presently performing a Gauss elimination step is either difficult + * or very simple, depending on the circumstances. If the sparsity + * pattern is symmetric (whether the matrix is symmetric is irrelevant + * here), then we can infer the rows which have a nonzero entry in the + * present column by looking at which columns in the present row are + * nonempty. In this case, we only need to look into a fixed number of + * rows and need not search all rows. On the other hand, if the + * sparsity pattern is nonsymmetric, then we need to use an iterative + * solver which can handle nonsymmetric matrices in any case, so there + * may be no need to do the Gauss elimination anyway. In fact, this is + * the way the function works: it takes a parameter + * (#elininate_columns#) that specifies whether the sparsity pattern + * is symmetric; if so, then the column is eliminated and the right + * hand side is also modified accordingly. If not, then only the row + * is deleted and the column is not touched at all, and all right hand + * side values apart from the one corresponding to the present row + * remain unchanged. + * + * If the sparsity pattern for your matrix is non-symmetric, you must + * set the value of this parameter to #false# in any case, since then + * we can't eliminate the column without searching all rows, which + * would be too expensive (if #N# be the number of rows, and #m# the + * number of nonzero elements per row, then eliminating one column is + * an #O(N*log(m))# operation, since searching in each row takes + * #log(m)# operations). If your sparsity pattern is symmetric, but + * your matrix is not, then you might specify #false# as well. If your + * sparsity pattern and matrix are both symmetric, you might want to + * specify #true# (the complexity of eliminating one row is then + * #O(m*log(m))#, since we only have to search #m# rows for the + * respective element of the column). Given the fact that #m# is + * roughly constant, irrespective of the discretization, and that the + * number of boundary nodes is #sqrt(N)# in 2d, the algorithm for + * symmetric sparsity patterns is #O(sqrt(N)*m*log(m))#, while it + * would be #O(N*sqrt(N)*m*log(m))# for the general case; the latter + * is too expensive to be performed. + * * It seems as if we had to make clear not to overwrite the lines of other * boundary nodes when doing the Gauss elimination step. However, since we * reset the right hand side when passing such a node, it is not a problem @@ -381,33 +418,22 @@ class MatrixCreator * the right hand side. We need therefore not take special care of other * boundary nodes. * - * To make solving faster, we preset the solution vector with the right boundary - * values. Since boundary nodes can never be hanging nodes, and since all other - * entries of the solution vector are zero, we need not condense the solution - * vector if the condensation process is done in-place. If done by copying - * matrix and vectors to smaller ones, it would also be necessary to condense - * the solution vector to preserve the preset boundary values. - * - * It it not clear whether the deletion of coupling between the boundary degree - * of freedom and other dofs really forces the corresponding entry in the - * solution vector to have the right value when using iterative solvers, - * since their search directions may contain components in the direction - * of the boundary node. For this reason, we perform a very simple line - * balancing by not setting the main diagonal entry to unity, but rather - * to the value it had before deleting this line, or to the first nonzero - * main diagonal entry if it is zero from a previous Gauss elimination - * step. Of course we have to change - * the right hand side appropriately. This is not a very good - * strategy, but it at least should give the main diagonal entry a value - * in the right order of dimension, which makes the solving process a bit - * more stable. A refined algorithm would set the entry to the mean of the - * other diagonal entries, but this seems to be too expensive. - * - * Because of the mentioned question, whether or not a preset solution value - * which does not couple with other degrees of freedom remains its value or - * not during solving iteratively, it may or may not be necessary to set - * the correct value after solving again. This question is an open one as of - * now and may be answered by future experience. + * To make solving faster, we preset the solution vector with the + * right boundary values. It it not clear whether the deletion of + * coupling between the boundary degree of freedom and other dofs + * really forces the corresponding entry in the solution vector to + * have the right value when using iterative solvers, since their + * search directions may contain components in the direction of the + * boundary node. For this reason, we perform a very simple line + * balancing by not setting the main diagonal entry to unity, but + * rather to the value it had before deleting this line, or to the + * first nonzero main diagonal entry if it is zero for some reason. + * Of course we have to change the right hand side appropriately. This + * is not a very good strategy, but it at least should give the main + * diagonal entry a value in the right order of dimension, which makes + * the solvution process a bit more stable. A refined algorithm would + * set the entry to the mean of the other diagonal entries, but this + * seems to be too expensive. * * * @author Wolfgang Bangerth, 1998 @@ -425,7 +451,8 @@ class MatrixTools : public MatrixCreator static void apply_boundary_values (const map &boundary_values, SparseMatrix &matrix, Vector &solution, - Vector &right_hand_side); + Vector &right_hand_side, + const bool eliminate_columns = true); /** * Exception diff --git a/deal.II/deal.II/source/numerics/base.cc b/deal.II/deal.II/source/numerics/base.cc index 7a07c3711c..ddf586f3ae 100644 --- a/deal.II/deal.II/source/numerics/base.cc +++ b/deal.II/deal.II/source/numerics/base.cc @@ -140,7 +140,8 @@ void ProblemBase::assemble (const Equation &equation, boundary_value_list); MatrixTools::apply_boundary_values (boundary_value_list, system_matrix, solution, - right_hand_side); + right_hand_side, + true); }; diff --git a/deal.II/deal.II/source/numerics/matrices.cc b/deal.II/deal.II/source/numerics/matrices.cc index 8f838033f6..4d26689880 100644 --- a/deal.II/deal.II/source/numerics/matrices.cc +++ b/deal.II/deal.II/source/numerics/matrices.cc @@ -34,7 +34,8 @@ template void MatrixCreator::create_mass_matrix (const DoFHandler &dof, const Quadrature &q, SparseMatrix &matrix, - const Function * const a) { + const Function * const a) +{ Vector dummy; // no entries, should give an error if accessed UpdateFlags update_flags = UpdateFlags(update_values | update_JxW_values); if (a != 0) @@ -57,13 +58,15 @@ void MatrixCreator::create_mass_matrix (const DoFHandler &dof, }; + template void MatrixCreator::create_mass_matrix (const DoFHandler &dof, const Quadrature &q, SparseMatrix &matrix, const Function &rhs, Vector &rhs_vector, - const Function * const a) { + const Function * const a) +{ UpdateFlags update_flags = UpdateFlags(update_values | update_q_points | update_JxW_values); @@ -85,9 +88,11 @@ void MatrixCreator::create_mass_matrix (const DoFHandler &dof, }; + template void MatrixCreator::create_mass_matrix (const DoFHandler &dof, - SparseMatrix &matrix) { + SparseMatrix &matrix) +{ const FiniteElement &fe = dof.get_fe(); const unsigned int dofs_per_cell = fe.dofs_per_cell; @@ -119,7 +124,8 @@ void MatrixCreator<1>::create_boundary_mass_matrix (const DoFHandler<1> &, const FunctionMap &, Vector &, vector &, - const Function<1> *) { + const Function<1> *) +{ Assert (false, ExcNotImplemented()); }; @@ -133,7 +139,8 @@ void MatrixCreator::create_boundary_mass_matrix (const DoFHandler & const FunctionMap &rhs, Vector &rhs_vector, vector &dof_to_boundary_mapping, - const Function *a) { + const Function *a) +{ const FiniteElement &fe = dof.get_fe(); const unsigned int n_components = fe.n_components(); const bool fe_is_system = (n_components != 1); @@ -404,11 +411,13 @@ void MatrixCreator::create_boundary_mass_matrix (const DoFHandler & }; + template void MatrixCreator::create_laplace_matrix (const DoFHandler &dof, const Quadrature &q, SparseMatrix &matrix, - const Function * const a) { + const Function * const a) +{ const unsigned int n_components = dof.get_fe().n_components(); Assert ((n_components==1) || (a==0), ExcNotImplemented()); @@ -434,6 +443,8 @@ void MatrixCreator::create_laplace_matrix (const DoFHandler &dof, while ((++assembler).state() == valid); }; + + /* template @@ -468,6 +479,8 @@ void MatrixCreator::create_level_laplace_matrix (unsigned int level, */ + + template void MatrixCreator::create_laplace_matrix (const DoFHandler &dof, const Quadrature &q, @@ -500,11 +513,15 @@ void MatrixCreator::create_laplace_matrix (const DoFHandler &dof, }; + template -void MatrixTools::apply_boundary_values (const map &boundary_values, - SparseMatrix &matrix, - Vector &solution, - Vector &right_hand_side) { +void +MatrixTools::apply_boundary_values (const map &boundary_values, + SparseMatrix &matrix, + Vector &solution, + Vector &right_hand_side, + const bool preserve_symmetry) +{ Assert (matrix.n() == matrix.m(), ExcDimensionsDontMatch(matrix.n(), matrix.m())); Assert (matrix.n() == right_hand_side.size(), @@ -583,47 +600,63 @@ void MatrixTools::apply_boundary_values (const map &bo new_rhs = right_hand_side(dof_number) = dof->second * first_nonzero_diagonal_entry; }; - - // store the only nonzero entry - // of this line for the Gauss - // elimination step - const double diagonal_entry = matrix.diag_element(dof_number); - // do the Gauss step - for (unsigned int row=0; row::apply_boundary_values (const map &bo }; }; - // preset solution vector solution(dof_number) = dof->second; }; @@ -648,13 +680,16 @@ MassMatrix::MassMatrix (const Function * const rhs, const Function * const a) : Equation (1), right_hand_side (rhs), - coefficient (a) {}; + coefficient (a) +{}; + template void MassMatrix::assemble (FullMatrix &cell_matrix, const FEValues &fe_values, - const typename DoFHandler::cell_iterator &) const { + const typename DoFHandler::cell_iterator &) const +{ const unsigned int dofs_per_cell = fe_values.dofs_per_cell, n_q_points = fe_values.n_quadrature_points; const FiniteElement &fe = fe_values.get_fe(); @@ -735,11 +770,13 @@ void MassMatrix::assemble (FullMatrix &cell_matrix, }; + template void MassMatrix::assemble (FullMatrix &cell_matrix, Vector &rhs, const FEValues &fe_values, - const DoFHandler::cell_iterator &) const { + const DoFHandler::cell_iterator &) const +{ Assert (right_hand_side != 0, ExcNoRHSSelected()); const unsigned int dofs_per_cell = fe_values.dofs_per_cell, @@ -800,10 +837,12 @@ void MassMatrix::assemble (FullMatrix &cell_matrix, }; + template void MassMatrix::assemble (Vector &rhs, const FEValues &fe_values, - const DoFHandler::cell_iterator &) const { + const DoFHandler::cell_iterator &) const +{ Assert (right_hand_side != 0, ExcNoRHSSelected()); const unsigned int dofs_per_cell = fe_values.dofs_per_cell, @@ -833,6 +872,7 @@ void MassMatrix::assemble (Vector &rhs, }; + template LaplaceMatrix::LaplaceMatrix (const Function * const rhs, const Function * const a) : @@ -841,11 +881,13 @@ LaplaceMatrix::LaplaceMatrix (const Function * const rhs, coefficient (a) {}; + template void LaplaceMatrix::assemble (FullMatrix &cell_matrix, Vector &rhs, const FEValues &fe_values, - const DoFHandler::cell_iterator &) const { + const DoFHandler::cell_iterator &) const +{ Assert (right_hand_side != 0, ExcNoRHSSelected()); const unsigned int dofs_per_cell = fe_values.dofs_per_cell, @@ -909,10 +951,12 @@ void LaplaceMatrix::assemble (FullMatrix &cell_matrix, }; + template void LaplaceMatrix::assemble (FullMatrix &cell_matrix, const FEValues &fe_values, - const DoFHandler::cell_iterator &) const { + const DoFHandler::cell_iterator &) const +{ const unsigned int dofs_per_cell = fe_values.dofs_per_cell, n_q_points = fe_values.n_quadrature_points; @@ -963,10 +1007,12 @@ void LaplaceMatrix::assemble (FullMatrix &cell_matrix, }; + template void LaplaceMatrix::assemble (Vector &rhs, const FEValues &fe_values, - const DoFHandler::cell_iterator &) const { + const DoFHandler::cell_iterator &) const +{ Assert (right_hand_side != 0, ExcNoRHSSelected()); const unsigned int dofs_per_cell = fe_values.dofs_per_cell, @@ -997,6 +1043,7 @@ void LaplaceMatrix::assemble (Vector &rhs, }; + template void MatrixCreator::create_interpolation_matrix(const FiniteElement &high, @@ -1032,6 +1079,9 @@ MatrixCreator::create_interpolation_matrix(const FiniteElement &high, } + +// explicit instantiations + template class MatrixCreator; template class MatrixTools; template class MassMatrix; diff --git a/deal.II/deal.II/source/numerics/vectors.cc b/deal.II/deal.II/source/numerics/vectors.cc index c2587cb979..8c3fb2a204 100644 --- a/deal.II/deal.II/source/numerics/vectors.cc +++ b/deal.II/deal.II/source/numerics/vectors.cc @@ -437,7 +437,8 @@ void VectorTools::project (const DoFHandler &dof, constraints.condense (tmp); if (boundary_values.size() != 0) MatrixTools::apply_boundary_values (boundary_values, - mass_matrix, vec, tmp); + mass_matrix, vec, tmp, + true); SolverControl control(1000,1e-16); PrimitiveVectorMemory<> memory; -- 2.39.5