From 105d6c7bddc1435d80a19bd249572cf947b626a6 Mon Sep 17 00:00:00 2001 From: wolf Date: Tue, 3 May 2005 03:44:46 +0000 Subject: [PATCH] Unify some code, and make it more efficient. git-svn-id: https://svn.dealii.org/trunk@10608 0785d39b-7218-0410-832d-ea1e28bc413d --- .../numerics/matrices.all_dimensions.cc | 375 +++++++----------- 1 file changed, 152 insertions(+), 223 deletions(-) diff --git a/deal.II/deal.II/source/numerics/matrices.all_dimensions.cc b/deal.II/deal.II/source/numerics/matrices.all_dimensions.cc index 5891d805a0..66cd52eca1 100644 --- a/deal.II/deal.II/source/numerics/matrices.all_dimensions.cc +++ b/deal.II/deal.II/source/numerics/matrices.all_dimensions.cc @@ -519,125 +519,165 @@ MatrixTools::apply_boundary_values (const std::map &boundar #ifdef DEAL_II_USE_PETSC -void -MatrixTools:: -apply_boundary_values (const std::map &boundary_values, - PETScWrappers::SparseMatrix &matrix, - PETScWrappers::Vector &solution, - PETScWrappers::Vector &right_hand_side, - const bool preserve_symmetry) +namespace PETScWrappers { - Assert (matrix.n() == right_hand_side.size(), - ExcDimensionMismatch(matrix.n(), right_hand_side.size())); - Assert (matrix.n() == solution.size(), - ExcDimensionMismatch(matrix.n(), solution.size())); - - // if no boundary values are to be applied - // simply return - if (boundary_values.size() == 0) - return; - - - // we have to read and write from this - // matrix (in this order). this will only - // work if we compress the matrix first, - // done here. do the same with the other - // objects just to be on the safe side: - matrix.compress (); - solution.compress (); - right_hand_side.compress (); - - const unsigned int n_dofs = matrix.m(); - - // if a diagonal entry is zero - // later, then we use another - // number instead. take it to be - // the first nonzero diagonal - // element of the matrix, or 1 if - // there is no such thing - PetscScalar first_nonzero_diagonal_entry = 1; - for (unsigned int i=0; i + void + apply_boundary_values (const std::map &boundary_values, + PETScMatrix &matrix, + PETScVector &solution, + PETScVector &right_hand_side, + const bool preserve_symmetry) + { + Assert (preserve_symmetry == false, ExcNotImplemented()); + + Assert (matrix.n() == right_hand_side.size(), + ExcDimensionMismatch(matrix.n(), right_hand_side.size())); + Assert (matrix.n() == solution.size(), + ExcDimensionMismatch(matrix.n(), solution.size())); + + // if no boundary values are to be applied + // simply return + if (boundary_values.size() == 0) + return; + + + // we have to read and write from this + // matrix (in this order). this will only + // work if we compress the matrix first, + // done here. do the same with the other + // objects just to be on the safe side: + matrix.compress (); + solution.compress (); + right_hand_side.compress (); + + const std::pair local_range + = matrix.local_range(); + + // determine the first nonzero diagonal + // entry from within the part of the matrix + // that we can see. if we can't find such + // an entry, take one + PetscScalar average_nonzero_diagonal_entry = 1; + for (unsigned int i=local_range.first; + i::const_iterator - dof = boundary_values.begin(), - endd = boundary_values.end(); - for (; dof != endd; ++dof) + // iterate over all fixed degrees of + // freedom that are within the local + // range of this matrix. the function is + // pretty awkward to implement, since we + // can't freely mix reading and writing + // from the matrix without global + // synchronisation, so we first determine + // which are the entries we have to work + // on in a first step, and then come back + // and do that all at once { - Assert (dof->first < n_dofs, ExcInternalError()); + matrix.compress (); - const unsigned int dof_number = dof->first; + std::map::const_iterator + dof = boundary_values.begin(), + endd = boundary_values.end(); + std::vector > + set_to_zero_entries; - // for each constrained dof: + for (; dof != endd; ++dof) + if ((dof->first >= local_range.first) && + (dof->first < local_range.second)) + { + const unsigned int dof_number = dof->first; - // set entries of this line - // to zero except for the diagonal - // entry. - { - PETScWrappers::MatrixBase::const_iterator - p = matrix.begin(dof_number), - e = matrix.end(dof_number); - - // iterate over all elements of this - // row and set elements to zero - // except for the diagonal - // element. note that this is not - // exactly clean programming, since - // we change the matrix underneath, - // while we still keep working with - // the iterators into it - for (; p!=e; ++p) - if (p->column() != dof_number) - matrix.set (dof_number, p->column(), 0.); - } + // for each constrained dof: + + // store which entries of this line + // to set to zero except for the + // diagonal entry. + PETScWrappers::MatrixBase::const_iterator + p = matrix.begin(dof_number), + e = matrix.end(dof_number); + + for (; p!=e; ++p) + { + Assert (p->row() == dof_number, ExcInternalError()); + + if (p->column() != dof_number) + set_to_zero_entries.push_back (std::make_pair (dof_number, + p->column())); + } + } + + // now set all these entries to zero in + // one bulk operation that requires + // only a single synchronisation: + matrix.compress (); + for (std::vector >::const_iterator + i = set_to_zero_entries.begin(); + i != set_to_zero_entries.end(); ++i) + matrix.set (i->first, i->second, 0.); + } - // set right hand side to wanted value; - // also store the new rhs entry to make - // the gauss step (when preserving the - // symmetry of the matrix) more - // efficient - // - // note that for petsc matrices - // interleaving read with write - // operations is very expensive. thus, - // we here always replace the diagonal - // element, rather than first checking - // whether it is nonzero and in that - // case preserving it. this is - // different from the case of deal.II - // sparse matrices treated in the other - // functions. - PetscScalar new_rhs; - matrix.set (dof_number, dof_number, - first_nonzero_diagonal_entry); - new_rhs = dof->second * first_nonzero_diagonal_entry; - right_hand_side(dof_number) = new_rhs; + // the next thing is to set right + // hand side to wanted value. note + // that for petsc matrices + // interleaving read with write + // operations is very + // expensive. thus, we here always + // replace the diagonal element, + // rather than first checking + // whether it is nonzero and in + // that case preserving it. this is + // different from the case of + // deal.II sparse matrices treated + // in the other functions. + { + matrix.compress (); + + std::map::const_iterator + dof = boundary_values.begin(), + endd = boundary_values.end(); + + for (; dof != endd; ++dof) + if ((dof->first >= local_range.first) && + (dof->first < local_range.second)) + { + const unsigned int dof_number = dof->first; + + matrix.set (dof_number, dof_number, + average_nonzero_diagonal_entry); + right_hand_side(dof_number) + = dof->second * average_nonzero_diagonal_entry; + + // preset solution vector + solution(dof_number) = dof->second; + } + } + + matrix.compress (); + solution.compress (); + right_hand_side.compress (); + } +} - // if the user wants to have - // the symmetry of the matrix - // preserved, and if the - // sparsity pattern is - // symmetric, then do a Gauss - // elimination step with the - // present row - if (preserve_symmetry) - { - Assert (false, ExcNotImplemented()); - } - // preset solution vector - solution(dof_number) = dof->second; - } - matrix.compress (); - solution.compress (); - right_hand_side.compress (); +void +MatrixTools:: +apply_boundary_values (const std::map &boundary_values, + PETScWrappers::SparseMatrix &matrix, + PETScWrappers::Vector &solution, + PETScWrappers::Vector &right_hand_side, + const bool preserve_symmetry) +{ + // simply redirect to the generic function + // used for both petsc matrix types + PETScWrappers::apply_boundary_values (boundary_values, matrix, solution, + right_hand_side, preserve_symmetry); } @@ -650,121 +690,10 @@ apply_boundary_values (const std::map &boundary_values, PETScWrappers::MPI::Vector &right_hand_side, const bool preserve_symmetry) { - // this function works almost exactly as - // the one above, with the only exception - // that we have to make sure that we only - // access elements that belong to the slice - // of the matrix that is available locally - Assert (matrix.n() == right_hand_side.size(), - ExcDimensionMismatch(matrix.n(), right_hand_side.size())); - Assert (matrix.n() == solution.size(), - ExcDimensionMismatch(matrix.n(), solution.size())); - - // if no boundary values are to be applied - // simply return - if (boundary_values.size() == 0) - return; - - - // we have to read and write from this - // matrix (in this order). this will only - // work if we compress the matrix first, - // done here. do the same with the other - // objects just to be on the safe side: - matrix.compress (); - solution.compress (); - right_hand_side.compress (); - - // determine the first nonzero diagonal - // entry from within the part of the matrix - // that we can see. if we can't find such - // an entry, take one - PetscScalar first_nonzero_diagonal_entry = 1; - for (unsigned int i=matrix.local_range().first; - i::const_iterator - dof = boundary_values.begin(), - endd = boundary_values.end(); - for (; dof != endd; ++dof) - if ((dof->first >= matrix.local_range().first) && - (dof->first < matrix.local_range().second)) - { - const unsigned int dof_number = dof->first; - - // for each constrained dof: - - // set entries of this line - // to zero except for the diagonal - // entry. - { - PETScWrappers::MatrixBase::const_iterator - p = matrix.begin(dof_number), - e = matrix.end(dof_number); - - // iterate over all elements of this - // row and set elements to zero - // except for the diagonal - // element. note that this is not - // exactly clean programming, since - // we change the matrix underneath, - // while we still keep working with - // the iterators into it - for (; p!=e; ++p) - if (p->column() != dof_number) - matrix.set (dof_number, p->column(), 0.); - } - - - // set right hand side to wanted value; - // also store the new rhs entry to make - // the gauss step (when preserving the - // symmetry of the matrix) more - // efficient - // - // note that for petsc matrices - // interleaving read with write - // operations is very expensive. thus, - // we here always replace the diagonal - // element, rather than first checking - // whether it is nonzero and in that - // case preserving it. this is - // different from the case of deal.II - // sparse matrices treated in the other - // functions. - PetscScalar new_rhs; - matrix.set (dof_number, dof_number, - first_nonzero_diagonal_entry); - new_rhs = dof->second * first_nonzero_diagonal_entry; - right_hand_side(dof_number) = new_rhs; - - // if the user wants to have - // the symmetry of the matrix - // preserved, and if the - // sparsity pattern is - // symmetric, then do a Gauss - // elimination step with the - // present row - if (preserve_symmetry) - { - Assert (false, ExcNotImplemented()); - } - - // preset solution vector - solution(dof_number) = dof->second; - } - - matrix.compress (); - solution.compress (); - right_hand_side.compress (); + // simply redirect to the generic function + // used for both petsc matrix types + PETScWrappers::apply_boundary_values (boundary_values, matrix, solution, + right_hand_side, preserve_symmetry); } -- 2.39.5