From 37ac5fd25f586e7dfbd27f8b1a34f197b76ac196 Mon Sep 17 00:00:00 2001
From: Wolfgang Bangerth <bangerth@math.tamu.edu>
Date: Sun, 15 Jun 2014 15:46:13 +0000
Subject: [PATCH] Patch by Michal Wichrowski: fix the same bug in a second
 place.

git-svn-id: https://svn.dealii.org/trunk@33047 0785d39b-7218-0410-832d-ea1e28bc413d
---
 deal.II/source/numerics/matrix_tools.cc | 150 +++++++++++++-----------
 1 file changed, 79 insertions(+), 71 deletions(-)

diff --git a/deal.II/source/numerics/matrix_tools.cc b/deal.II/source/numerics/matrix_tools.cc
index dbb7d854f6..b836e73cd8 100644
--- a/deal.II/source/numerics/matrix_tools.cc
+++ b/deal.II/source/numerics/matrix_tools.cc
@@ -2553,77 +2553,85 @@ namespace MatrixTools
         Assert (matrix.n() == solution.size(),
                 ExcDimensionMismatch(matrix.m(), solution.size()));
 
-        // if no boundary values are to be applied
-        // simply return
-        if (boundary_values.size() == 0)
-          return;
-
-        const std::pair<types::global_dof_index, types::global_dof_index> local_range
-          = matrix.local_range();
-        Assert (local_range == right_hand_side.local_range(),
-                ExcInternalError());
-        Assert (local_range == solution.local_range(),
-                ExcInternalError());
-
-        // we have to read and write from this
-        // matrix (in this order). this will only
-        // work if we compress the matrix first,
-        // done here
-        matrix.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
-        TrilinosScalar average_nonzero_diagonal_entry = 1;
-        for (types::global_dof_index i=local_range.first; i<local_range.second; ++i)
-          if (matrix.diag_element(i) != 0)
-            {
-              average_nonzero_diagonal_entry = std::fabs(matrix.diag_element(i));
-              break;
-            }
-
-        // figure out which rows of the matrix we
-        // have to eliminate on this processor
-        std::vector<types::global_dof_index> constrained_rows;
-        for (std::map<types::global_dof_index,double>::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. if the value already is
-        // nonzero, it will be preserved,
-        // in accordance with the basic
-        // matrix classes in deal.II.
-        matrix.clear_rows (constrained_rows, average_nonzero_diagonal_entry);
-
-        std::vector<types::global_dof_index> indices;
-        std::vector<TrilinosScalar>  solution_values;
-        for (std::map<types::global_dof_index,double>::const_iterator
-             dof  = boundary_values.begin();
-             dof != boundary_values.end();
-             ++dof)
-          if ((dof->first >= local_range.first) &&
-              (dof->first < local_range.second))
-            {
-              indices.push_back (dof->first);
-              solution_values.push_back (dof->second);
-            }
-        solution.set (indices, solution_values);
-
-        // now also set appropriate
-        // values for the rhs
-        for (unsigned int i=0; i<solution_values.size(); ++i)
-          solution_values[i] *= matrix.diag_element(indices[i]);
-
-        right_hand_side.set (indices, solution_values);
-
+        // if no boundary values are to be applied, then
+        // jump straight to the compress() calls that we still have
+        // to perform because they are collective operations
+        if (boundary_values.size() > 0)
+	  {
+	    const std::pair<types::global_dof_index, types::global_dof_index> local_range
+	      = matrix.local_range();
+	    Assert (local_range == right_hand_side.local_range(),
+		    ExcInternalError());
+	    Assert (local_range == solution.local_range(),
+		    ExcInternalError());
+
+	    // we have to read and write from this
+	    // matrix (in this order). this will only
+	    // work if we compress the matrix first,
+	    // done here
+	    matrix.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
+	    TrilinosScalar average_nonzero_diagonal_entry = 1;
+	    for (types::global_dof_index i=local_range.first; i<local_range.second; ++i)
+	      if (matrix.diag_element(i) != 0)
+		{
+		  average_nonzero_diagonal_entry = std::fabs(matrix.diag_element(i));
+		  break;
+		}
+
+	    // figure out which rows of the matrix we
+	    // have to eliminate on this processor
+	    std::vector<types::global_dof_index> constrained_rows;
+	    for (std::map<types::global_dof_index,double>::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. if the value already is
+	    // nonzero, it will be preserved,
+	    // in accordance with the basic
+	    // matrix classes in deal.II.
+	    matrix.clear_rows (constrained_rows, average_nonzero_diagonal_entry);
+
+	    std::vector<types::global_dof_index> indices;
+	    std::vector<TrilinosScalar>  solution_values;
+	    for (std::map<types::global_dof_index,double>::const_iterator
+		   dof  = boundary_values.begin();
+		 dof != boundary_values.end();
+		 ++dof)
+	      if ((dof->first >= local_range.first) &&
+		  (dof->first < local_range.second))
+		{
+		  indices.push_back (dof->first);
+		  solution_values.push_back (dof->second);
+		}
+	    solution.set (indices, solution_values);
+
+	    // now also set appropriate
+	    // values for the rhs
+	    for (unsigned int i=0; i<solution_values.size(); ++i)
+	      solution_values[i] *= matrix.diag_element(indices[i]);
+
+	    right_hand_side.set (indices, solution_values);
+	  }
+	else
+	  {
+            // clear_rows() is a collective operation so we still have to call
+            // it:
+	    std::vector<types::global_dof_index> constrained_rows;
+	    matrix.clear_rows (constrained_rows, 1.);
+	  }
+	
         // clean up
         matrix.compress ();
         solution.compress (VectorOperation::insert);
-- 
2.39.5