#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/block_sparsity_pattern.h>
#include <deal.II/lac/block_sparse_matrix.h>
+#include <deal.II/lac/parallel_vector.h>
+#include <deal.II/lac/parallel_block_vector.h>
#include <iomanip>
+// number of functions to select the right implementation for set_zero().
+namespace internal
+{
+ namespace ConstraintMatrix
+ {
+ namespace
+ {
+ // TODO: in general we should iterate over the constraints and not over all DoFs
+ // for performance reasons
+ template<class VEC>
+ void set_zero_parallel(const dealii::ConstraintMatrix &cm, VEC &vec, unsigned int shift = 0)
+ {
+ Assert(!vec.has_ghost_elements(), ExcInternalError());//ExcGhostsPresent());
+
+ const unsigned int
+ start = vec.local_range().first,
+ end = vec.local_range().second;
+ for (unsigned int i=start; i<end; ++i)
+ if (cm.is_constrained (shift + i))
+ vec(i) = 0;
+ }
+
+ template<class VEC>
+ void set_zero_in_parallel(const dealii::ConstraintMatrix &cm, VEC &vec, internal::bool2type<false>)
+ {
+ set_zero_parallel(cm, vec, 0);
+ }
+
+ // in parallel for BlockVectors
+ template<class VEC>
+ void set_zero_in_parallel(const dealii::ConstraintMatrix &cm, VEC &vec, internal::bool2type<true>)
+ {
+ unsigned int start_shift = 0;
+ for (unsigned int j=0; j<vec.n_blocks(); ++j)
+ {
+ set_zero_parallel(cm, vec.block(j), start_shift);
+ start_shift += vec.block(j).size();
+ }
+ }
+
+ template<class VEC>
+ void set_zero_serial(const dealii::ConstraintMatrix &cm, VEC &vec)
+ {
+ // TODO would be faster:
+ /* std::vector<dealii::ConstraintMatrix::ConstraintLine>::const_iterator constraint_line = cm.lines.begin();
+ for (; constraint_line!=cm.lines.end(); ++constraint_line)
+ vec(constraint_line->line) = 0.;*/
+ for (unsigned int i=0; i<vec.size(); ++i)
+ if (cm.is_constrained (i))
+ vec(i) = 0;
+ }
+
+ template<class VEC>
+ void set_zero_all(const dealii::ConstraintMatrix &cm, VEC &vec)
+ {
+ set_zero_in_parallel<VEC>(cm, vec, internal::bool2type<IsBlockVector<VEC>::value>());
+ }
+
+
+ template<class T>
+ void set_zero_all(const dealii::ConstraintMatrix &cm, dealii::Vector<T> &vec)
+ {
+ set_zero_serial(cm, vec);
+ }
+
+ template<class T>
+ void set_zero_all(const dealii::ConstraintMatrix &cm, dealii::BlockVector<T> &vec)
+ {
+ set_zero_serial(cm, vec);
+ }
+
+
+ template<class T>
+ void set_zero_all(const dealii::ConstraintMatrix &cm, dealii::parallel::distributed::Vector<T> &vec)
+ {
+ // should use the general template above, but requires that
+ // some member functions are added to parallel::distributed::Vector
+ Assert (false, ExcNotImplemented());
+ }
+
+ template<class T>
+ void set_zero_all(const dealii::ConstraintMatrix &cm, dealii::parallel::distributed::BlockVector<T> &vec)
+ {
+ Assert (false, ExcNotImplemented());
+ }
+ }
+ }
+}
+
+
template <class VectorType>
void
ConstraintMatrix::set_zero (VectorType &vec) const
{
- Assert (sorted == true, ExcMatrixNotClosed());
-
- std::vector<ConstraintLine>::const_iterator constraint_line = lines.begin();
- for (; constraint_line!=lines.end(); ++constraint_line)
- vec(constraint_line->line) = 0.;
+ internal::ConstraintMatrix::set_zero_all(*this, vec);
}
+
template <typename VectorType>
void
ConstraintMatrix::
# $Id$
######################################################################
#
-# Copyright (C) 2011 by the deal.II authors
+# Copyright (C) 2011, 2012 by the deal.II authors
#
# This file is subject to QPL and may not be distributed
# without copyright and license information. Please refer
template void ConstraintMatrix::condense<V1 >(const V1 &, V1&) const;
template void ConstraintMatrix::condense<V1 >(V1 &vec) const;
-template void ConstraintMatrix::set_zero<V1 >(V1&) const;
template void ConstraintMatrix::distribute_local_to_global<V1 > (
const Vector<double>&, const std::vector<unsigned int> &, V1&, const FullMatrix<double>&) const;
template void ConstraintMatrix::distribute<V1 >(const V1 &, V1&) const;
template void ConstraintMatrix::distribute<V1 >(V1 &) const;
+EOT
+ ;
+
+my $vector_functions_also_parallel = <<EOT
+
+template void ConstraintMatrix::set_zero<V1 >(V1&) const;
+
EOT
;
######################################################################
multisubst($vector_functions, ['V1'], \@sequential_vectors);
+multisubst($vector_functions_also_parallel, ['V1'], \@sequential_vectors);
+multisubst($vector_functions_also_parallel, ['V1'], \@parallel_vectors);
multisubst($scalar_functions, ['S1'], \@real_scalars);
multisubst($scalar_scalar_functions, ['S1','S2'], \@real_scalars, \@real_scalars);