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);