Vector<Number, MemorySpace>::local_element(
const size_type local_index) const
{
- Assert(false, ExcNotImplemented());
+ Assert((std::is_same<MemorySpace, ::dealii::MemorySpace::Host>::value),
+ ExcMessage(
+ "This function is only implemented for the Host memory space"));
+ AssertIndexRange(local_index,
+ partitioner->local_size() +
+ partitioner->n_ghost_indices());
+ // do not allow reading a vector which is not in ghost mode
+ Assert(local_index < local_size() || vector_is_ghosted == true,
+ ExcMessage("You tried to read a ghost element of this vector, "
+ "but it has not imported its ghost values."));
return data.values[local_index];
}
inline Number &
Vector<Number, MemorySpace>::local_element(const size_type local_index)
{
- Assert(false, ExcNotImplemented());
+ Assert((std::is_same<MemorySpace, ::dealii::MemorySpace::Host>::value),
+ ExcMessage(
+ "This function is only implemented for the Host memory space"));
+
+ AssertIndexRange(local_index,
+ partitioner->local_size() +
+ partitioner->n_ghost_indices());
return data.values[local_index];
}
const Vector<Number2, MemorySpaceType> &v,
const bool omit_zeroing_entries)
{
- Assert(false, ExcNotImplemented());
- (void)v;
- (void)omit_zeroing_entries;
+ clear_mpi_requests();
+ Assert(v.partitioner.get() != nullptr, ExcNotInitialized());
+
+ // check whether the partitioners are
+ // different (check only if the are allocated
+ // differently, not if the actual data is
+ // different)
+ if (partitioner.get() != v.partitioner.get())
+ {
+ partitioner = v.partitioner;
+ partitioner_sm = v.partitioner_sm;
+ const size_type new_allocated_size =
+ partitioner->local_size() + partitioner->n_ghost_indices();
+ resize_val(new_allocated_size,
+ partitioner_sm->get_sm_mpi_communicator());
+ }
+
+ if (omit_zeroing_entries == false)
+ this->operator=(Number());
+ else
+ zero_out_ghosts();
+
+ // do not reallocate import_data directly, but only upon request. It
+ // is only used as temporary storage for compress() and
+ // update_ghost_values, and we might have vectors where we never
+ // call these methods and hence do not need to have the storage.
+ import_data.values.reset();
+ import_data.values_dev.reset();
}
, allocated_size(0)
, vector_is_ghosted(false)
{
- Assert(false, ExcNotImplemented());
- (void)v;
+ reinit(v, true);
+
+ const size_type this_size = local_size();
+ if (this_size > 0)
+ std::memcpy(this->begin(), v.begin(), this_size * sizeof(Number));
}
Vector<Number, MemorySpaceType>::
operator-=(const VectorSpaceVector<Number> &vv)
{
- Assert(false, ExcNotImplemented());
- (void)vv;
+ // Downcast. Throws an exception if invalid.
+ using VectorType = Vector<Number, MemorySpaceType>;
+ Assert(dynamic_cast<const VectorType *>(&vv) != nullptr,
+ ExcVectorTypeNotCompatible());
+ const VectorType &v = dynamic_cast<const VectorType &>(vv);
+
+ AssertDimension(local_size(), v.local_size());
+
+ auto values = this->begin();
+ const auto values_other = v.begin();
+
+ for (unsigned int i = 0; i < partitioner->local_size(); i++)
+ values[i] -= values_other[i];
+
+ if (vector_is_ghosted)
+ update_ghost_values();
return *this;
}
const Number a,
const VectorSpaceVector<Number> &vv)
{
- Assert(false, ExcNotImplemented());
- (void)a;
- (void)vv;
+ // Downcast. Throws an exception if invalid.
+ using VectorType = Vector<Number, MemorySpaceType>;
+ Assert(dynamic_cast<const VectorType *>(&vv) != nullptr,
+ ExcVectorTypeNotCompatible());
+ const VectorType &v = dynamic_cast<const VectorType &>(vv);
+
+ AssertIsFinite(a);
+ AssertDimension(local_size(), v.local_size());
+
+ // nothing to do if a is zero
+ if (a == Number(0.))
+ return;
+
+ auto values = this->begin();
+ const auto values_other = v.begin();
+
+ for (unsigned int i = 0; i < partitioner->local_size(); i++)
+ values[i] += a * values_other[i];
}
Vector<Number, MemorySpaceType>::add(const Number a,
const VectorSpaceVector<Number> &vv)
{
- Assert(false, ExcNotImplemented());
- (void)a;
- (void)vv;
+ add_local(a, vv);
+
+ if (vector_is_ghosted)
+ update_ghost_values();
}
const Number a,
const VectorSpaceVector<Number> &vv)
{
- Assert(false, ExcNotImplemented());
- (void)x;
- (void)a;
- (void)vv;
+ // Downcast. Throws an exception if invalid.
+ using VectorType = Vector<Number, MemorySpaceType>;
+ Assert(dynamic_cast<const VectorType *>(&vv) != nullptr,
+ ExcVectorTypeNotCompatible());
+ const VectorType &v = dynamic_cast<const VectorType &>(vv);
+
+ AssertIsFinite(a);
+ AssertDimension(local_size(), v.local_size());
+
+ // nothing to do if a is zero
+ if (a == Number(0.))
+ return;
+
+ auto values = this->begin();
+ const auto values_other = v.begin();
+
+ for (unsigned int i = 0; i < partitioner->local_size(); i++)
+ values[i] = x * values[i] + a * values_other[i];
}
const Number a,
const VectorSpaceVector<Number> &vv)
{
- Assert(false, ExcNotImplemented());
- (void)x;
- (void)a;
- (void)vv;
+ sadd_local(x, a, vv);
+
+ if (vector_is_ghosted)
+ update_ghost_values();
}
typename Vector<Number, MemorySpaceType>::real_type
Vector<Number, MemorySpaceType>::linfty_norm_local() const
{
- Assert(false, ExcNotImplemented());
- return 0;
+ real_type max = 0.;
+
+ auto values = this->begin();
+
+ for (unsigned int i = 0; i < partitioner->local_size(); i++)
+ max = std::max(std::abs(values[i]), max);
+
+ return max;
}
inline typename Vector<Number, MemorySpaceType>::real_type
Vector<Number, MemorySpaceType>::linfty_norm() const
{
- Assert(false, ExcNotImplemented());
- return 0;
+ const real_type local_result = linfty_norm_local();
+ if (partitioner->n_mpi_processes() > 1)
+ return Utilities::MPI::max(local_result,
+ partitioner->get_mpi_communicator());
+ else
+ return local_result;
}
Vector<Number, MemorySpaceType>::partitioners_are_compatible(
const Utilities::MPI::Partitioner &part) const
{
- Assert(false, ExcNotImplemented());
return partitioner->is_compatible(part);
}
::dealii::MemorySpace::Host>(
const Vector<SCALAR, ::dealii::MemorySpace::Host> &,
VectorOperation::values);
- \}
- \}
- }
-
-for (S1 : REAL_SCALARS; S2 : REAL_SCALARS)
- {
- namespace LinearAlgebra
- \{
- namespace SharedMPI
- \{
template void
- Vector<S1, ::dealii::MemorySpace::Host>::reinit<S2>(
- const Vector<S2, ::dealii::MemorySpace::Host> &,
+ Vector<SCALAR, ::dealii::MemorySpace::Host>::reinit<SCALAR>(
+ const Vector<SCALAR, ::dealii::MemorySpace::Host> &,
const bool);
- template S1
- Vector<S1, ::dealii::MemorySpace::Host>::inner_product_local<S2>(
- const Vector<S2, ::dealii::MemorySpace::Host> &) const;
+ template SCALAR
+ Vector<SCALAR, ::dealii::MemorySpace::Host>::inner_product_local<
+ SCALAR>(const Vector<SCALAR, ::dealii::MemorySpace::Host> &) const;
template void
- Vector<S1, ::dealii::MemorySpace::Host>::copy_locally_owned_data_from<
- S2>(const Vector<S2, ::dealii::MemorySpace::Host> &);
+ Vector<SCALAR, ::dealii::MemorySpace::Host>::
+ copy_locally_owned_data_from<SCALAR>(
+ const Vector<SCALAR, ::dealii::MemorySpace::Host> &);
\}
\}
}
dealii::MatrixFree<dim, Number> matrix_free;
matrix_free.reinit(mapping, dof_handler, constraint, quad, additional_data);
- MPI_Comm comm = MPI_COMM_WORLD;
+ const MPI_Comm comm = MPI_COMM_WORLD;
const unsigned int rank = Utilities::MPI::this_mpi_process(comm);
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2019 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+
+// test the correctness of the matrix-free cell loop with additional
+// operation_before_loop and operation_after_loop lambdas
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/la_sm_vector.h>
+#include <deal.II/lac/sparse_matrix.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+#include "../matrix_free/matrix_vector_mf.h"
+
+
+
+template <int dim,
+ int fe_degree,
+ int n_q_points_1d,
+ typename Number,
+ typename VectorType>
+class Matrix
+{
+public:
+ Matrix(const MatrixFree<dim, Number> &data_in)
+ : data(data_in)
+ {}
+
+ void
+ vmult(VectorType &dst, const VectorType &src) const
+ {
+ const std::function<void(const MatrixFree<dim, Number> &,
+ VectorType &,
+ const VectorType &,
+ const std::pair<unsigned int, unsigned int> &)>
+ wrap = helmholtz_operator<dim, fe_degree, VectorType, n_q_points_1d>;
+ dst = 0;
+ data.cell_loop(wrap, dst, src);
+ for (auto i : data.get_constrained_dofs())
+ dst.local_element(i) = src.local_element(i);
+ }
+
+ void
+ vmult_with_update_basic(VectorType &dst,
+ VectorType &src,
+ VectorType &other_vector) const
+ {
+ src.add(1.5, other_vector);
+ other_vector.add(0.7, dst);
+ dst = 0;
+ vmult(dst, src);
+ dst.sadd(0.63, -1.3, other_vector);
+ }
+
+ void
+ vmult_with_update_merged(VectorType &dst,
+ VectorType &src,
+ VectorType &other_vector) const
+ {
+ const std::function<void(const MatrixFree<dim, Number> &,
+ VectorType &,
+ const VectorType &,
+ const std::pair<unsigned int, unsigned int> &)>
+ wrap = helmholtz_operator<dim, fe_degree, VectorType, n_q_points_1d>;
+ data.cell_loop(
+ wrap,
+ dst,
+ src,
+ // operation before cell operation
+ [&](const unsigned int start_range, const unsigned int end_range) {
+ for (unsigned int i = start_range; i < end_range; ++i)
+ {
+ src.local_element(i) += 1.5 * other_vector.local_element(i);
+ other_vector.local_element(i) += 0.7 * dst.local_element(i);
+ dst.local_element(i) = 0;
+ }
+ },
+ // operation after cell operation
+ [&](const unsigned int start_range, const unsigned int end_range) {
+ for (unsigned int i = start_range; i < end_range; ++i)
+ {
+ dst.local_element(i) =
+ 0.63 * dst.local_element(i) - 1.3 * other_vector.local_element(i);
+ }
+ });
+ }
+
+private:
+ const MatrixFree<dim, Number> &data;
+};
+
+
+
+template <int dim, int fe_degree, typename number>
+void
+test()
+{
+ const MPI_Comm comm = MPI_COMM_WORLD;
+
+ parallel::distributed::Triangulation<dim> tria(comm);
+ GridGenerator::hyper_cube(tria);
+ tria.refine_global(7 - dim);
+
+ FE_Q<dim> fe(fe_degree);
+ DoFHandler<dim> dof(tria);
+ dof.distribute_dofs(fe);
+ AffineConstraints<double> constraints;
+ constraints.close();
+
+ deallog << "Testing " << dof.get_fe().get_name() << std::endl;
+
+ MatrixFree<dim, number> mf_data;
+ {
+ const QGauss<1> quad(fe_degree + 1);
+ typename MatrixFree<dim, number>::AdditionalData data;
+ data.tasks_parallel_scheme = MatrixFree<dim, number>::AdditionalData::none;
+ mf_data.reinit(dof, constraints, quad, data);
+ }
+
+ const unsigned int group_size = 1;
+ const unsigned int rank = Utilities::MPI::this_mpi_process(comm);
+
+ MPI_Comm comm_sm;
+ MPI_Comm_split(comm, rank / group_size, rank, &comm_sm);
+
+ using VectorType = LinearAlgebra::SharedMPI::Vector<number>;
+
+ Matrix<dim, fe_degree, fe_degree + 1, number, VectorType> mf(mf_data);
+ VectorType vec1, vec2, vec3;
+ mf_data.initialize_dof_vector(vec1, comm_sm);
+ mf_data.initialize_dof_vector(vec2, comm_sm);
+ mf_data.initialize_dof_vector(vec3, comm_sm);
+
+ for (unsigned int i = 0; i < vec1.local_size(); ++i)
+ {
+ double entry = random_value<double>();
+ vec1.local_element(i) = entry;
+ entry = random_value<double>();
+ vec2.local_element(i) = entry;
+ entry = random_value<double>();
+ vec3.local_element(i) = entry;
+ }
+
+ VectorType ref1 = vec1;
+ VectorType ref2 = vec2;
+ VectorType ref3 = vec3;
+
+ mf.vmult_with_update_basic(ref3, ref2, ref1);
+ mf.vmult_with_update_basic(vec3, vec2, vec1);
+
+ vec3 -= ref3;
+ deallog << "Error in 1x merged operation: " << vec3.linfty_norm()
+ << std::endl;
+
+ ref3 = 0;
+
+ mf.vmult_with_update_basic(ref1, ref2, ref3);
+ mf.vmult_with_update_basic(vec1, vec2, vec3);
+
+ mf.vmult_with_update_basic(ref1, ref2, ref3);
+ mf.vmult_with_update_basic(vec1, vec2, vec3);
+
+ vec3 -= ref3;
+ deallog << "Error in 2x merged operation: " << vec3.linfty_norm()
+ << std::endl;
+}
+
+
+int
+main(int argc, char **argv)
+{
+ Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1);
+
+ mpi_initlog();
+ deallog << std::setprecision(3);
+
+ test<2, 3, double>();
+ test<2, 3, float>();
+ test<3, 3, double>();
+}
\ No newline at end of file
--- /dev/null
+
+DEAL::Testing FE_Q<2>(3)
+DEAL::Error in 1x merged operation: 0.00
+DEAL::Error in 2x merged operation: 0.00
+DEAL::Testing FE_Q<2>(3)
+DEAL::Error in 1x merged operation: 0.00
+DEAL::Error in 2x merged operation: 0.00
+DEAL::Testing FE_Q<3>(3)
+DEAL::Error in 1x merged operation: 0.00
+DEAL::Error in 2x merged operation: 0.00
--- /dev/null
+
+DEAL::Testing FE_Q<2>(3)
+DEAL::Error in 1x merged operation: 0.00
+DEAL::Error in 2x merged operation: 0.00
+DEAL::Testing FE_Q<2>(3)
+DEAL::Error in 1x merged operation: 0.00
+DEAL::Error in 2x merged operation: 0.00
+DEAL::Testing FE_Q<3>(3)
+DEAL::Error in 1x merged operation: 0.00
+DEAL::Error in 2x merged operation: 0.00
projects / dealii.git / commitdiff