From b416b3fd3e702127a4890d2a4ee0235d13b3d3d2 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Wed, 10 Aug 2005 22:02:30 +0000 Subject: [PATCH] Reimplement apply_boundary_values for PETSc matrices using clear_rows. git-svn-id: https://svn.dealii.org/trunk@11278 0785d39b-7218-0410-832d-ea1e28bc413d --- .../numerics/matrices.all_dimensions.cc | 157 ++++++------------ 1 file changed, 47 insertions(+), 110 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 71b0b21937..e2b2f853fe 100644 --- a/deal.II/deal.II/source/numerics/matrices.all_dimensions.cc +++ b/deal.II/deal.II/source/numerics/matrices.all_dimensions.cc @@ -554,26 +554,10 @@ namespace PETScWrappers const std::pair local_range = matrix.local_range(); - // make sure that there is at - // least one boundary node on - // each processor, or we get into - // trouble with synchronisation - { - bool found = false; - for (unsigned int i=local_range.first; i::const_iterator - dof = boundary_values.begin(), - endd = boundary_values.end(); - std::vector > - set_to_zero_entries; - - for (; dof != endd; ++dof) - if ((dof->first >= local_range.first) && - (dof->first < local_range.second)) - { - const unsigned int dof_number = dof->first; - - // 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.); - } - + // figure out which rows of the matrix we + // have to eliminate on this processor + std::vector constrained_rows; + for (std::map::const_iterator + dof = boundary_values.begin(); + dof != boundary_values.end(); + ++dof) + if ((dof->first >= local_range.first) && + (dof->first < local_range.second)) + constrained_rows.push_back (dof->first); + + // then eliminate these rows and set + // their diagonal entry to what we have + // determined above + matrix.clear_rows (constrained_rows, average_nonzero_diagonal_entry); + + // the next thing is to set right hand + // side to the 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. + right_hand_side.compress (); + solution.compress (); - // 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 (); - right_hand_side.compress (); - solution.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; - } - } + 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; + right_hand_side(dof_number) + = dof->second * average_nonzero_diagonal_entry; + + // preset solution vector + solution(dof_number) = dof->second; + } matrix.compress (); solution.compress (); -- 2.39.5