*/
virtual Number operator* (const VectorSpaceVector<Number> &V) const override;
+ /**
+ * This function is not implemented and will throw an exception.
+ */
+ virtual void import(const ReadWriteVector<Number> &V,
+ VectorOperation::values operation,
+ const CommunicationPatternBase *communication_pattern = NULL) override;
+
/**
* Add @p a to all components. Note that @p a is a scalar not a vector.
*/
+ template <typename Number>
+ void Vector<Number>::import(const ReadWriteVector<Number> &,
+ VectorOperation::values ,
+ const CommunicationPatternBase *)
+ {
+ AssertThrow(false, ExcMessage("This function is not implemented."));
+ }
+
+
+
template <typename Number>
void Vector<Number>::add(const Number a, const VectorSpaceVector<Number> &V)
{
*/
void import(const PETScWrappers::MPI::Vector &petsc_vec,
VectorOperation::values operation,
- const CommunicationPatternBase *communication_pattern = NULL);
+ const CommunicationPatternBase *communication_pattern = nullptr);
#endif
#ifdef DEAL_II_WITH_TRILINOS
*/
void import(const TrilinosWrappers::MPI::Vector &trilinos_vec,
VectorOperation::values operation,
- const CommunicationPatternBase *communication_pattern = NULL);
+ const CommunicationPatternBase *communication_pattern = nullptr);
/**
* Imports all the elements present in the vector's IndexSet from the input
*/
void import(const EpetraWrappers::Vector &epetra_vec,
VectorOperation::values operation,
- const CommunicationPatternBase *communication_pattern = NULL);
+ const CommunicationPatternBase *communication_pattern = nullptr);
#endif
/**
protected:
#ifdef DEAL_II_WITH_TRILINOS
/**
- * Imports all the elements present in the vector's IndexSet from the input
+ * Import all the elements present in the vector's IndexSet from the input
* vector @p multivector. This is an helper function and it should not be
* used directly.
*/
// The first time import is called, we create a communication pattern.
// Check if the communication pattern already exists and if it can be
// reused.
- if (source_elements == source_stored_elements)
+ if ((source_elements.size() == source_stored_elements.size()) &&
+ (source_elements == source_stored_elements))
{
- epetra_comm_pattern.reset(
- dynamic_cast<const EpetraWrappers::CommunicationPattern *> (comm_pattern.get()));
+ epetra_comm_pattern =
+ std::dynamic_pointer_cast<const EpetraWrappers::CommunicationPattern> (comm_pattern);
if (epetra_comm_pattern == nullptr)
epetra_comm_pattern = std::make_shared<const EpetraWrappers::CommunicationPattern>(
create_epetra_comm_pattern(source_elements, mpi_comm));
source_stored_elements = source_index_set;
EpetraWrappers::CommunicationPattern epetra_comm_pattern(
source_stored_elements, stored_elements, mpi_comm);
- comm_pattern = std::make_shared<EpetraWrappers::CommunicationPattern>(
- epetra_comm_pattern);
+ comm_pattern.reset(new EpetraWrappers::CommunicationPattern(
+ source_stored_elements, stored_elements, mpi_comm));
return epetra_comm_pattern;
}
* calling separate methods means to load the calling vector @p this
* twice. Since most vector operations are memory transfer limited, this
* reduces the time by 25\% (or 50\% if @p W equals @p this).
+ *
+ * The vectors need to have the same layout.
*/
virtual double add_and_dot(const double a,
const VectorSpaceVector<double> &V,
virtual VectorSpaceVector<Number> &operator-= (const VectorSpaceVector<Number> &V) = 0;
/**
- * Imports all the elements present in the vector's IndexSet from the input
+ * Import all the elements present in the vector's IndexSet from the input
* vector @p V. VectorOperation::values @p operation is used to decide if
* the elements in @p V should be added to the current vector or replace the
* current elements. The last parameter can be used if the same
*/
virtual void import(const ReadWriteVector<Number> &V,
VectorOperation::values operation,
- const CommunicationPatternBase *communication_pattern = NULL);
+ const CommunicationPatternBase *communication_pattern = NULL) = 0;
/**
* Return the scalar product of two vectors.
{
comm = std::make_shared<const MPI_Comm>(communicator);
- Epetra_Map vector_space_vector_map = vector_space_vector_index_set.make_trilinos_map(*comm);
- Epetra_Map read_write_vector_map = read_write_vector_index_set.make_trilinos_map(*comm);
+ Epetra_Map vector_space_vector_map = vector_space_vector_index_set.make_trilinos_map(*comm,
+ false);
+ Epetra_Map read_write_vector_map = read_write_vector_index_set.make_trilinos_map(*comm,
+ true);
// Target map is read_write_vector_map
// Source map is vector_space_vector_map. This map must have uniquely
{
Vector::Vector()
:
- vector(new Epetra_FEVector(Epetra_Map(0,1,0,Epetra_MpiComm(MPI_COMM_SELF)))),
+ vector(new Epetra_FEVector(Epetra_Map(0,0,0,Utilities::Trilinos::comm_self()))),
epetra_comm_pattern(nullptr)
{}
// The first time import is called, a communication pattern is created.
// Check if the communication pattern already exists and if it can be
// reused.
- if (source_stored_elements == V.get_stored_elements())
+ if ((source_stored_elements.size() != V.get_stored_elements().size()) ||
+ ((source_stored_elements.size() == V.get_stored_elements().size()) &&
+ (source_stored_elements != V.get_stored_elements())))
{
create_epetra_comm_pattern(V.get_stored_elements(),
dynamic_cast<const Epetra_MpiComm &>(vector->Comm()).Comm());
::dealii::IndexSet Vector::locally_owned_elements() const
{
- const Epetra_Map *map = dynamic_cast<const Epetra_Map *>(&(vector->Map()));
- return IndexSet(*map);
+ IndexSet is (size());
+
+ // easy case: local range is contiguous
+ if (vector->Map().LinearMap())
+ {
+#ifndef DEAL_II_WITH_64BIT_INDICES
+ is.add_range(vector->Map().MinMyGID(), vector->Map().MaxMyGID()+1);
+#else
+ is.add_range(vector->Map().MinMyGID64(), vector->Map().MaxMyGID64()+1);
+#endif
+ }
+ else if (vector->Map().NumMyElements() > 0)
+ {
+ const size_type n_indices = vector->Map().NumMyElements();
+#ifndef DEAL_II_WITH_64BIT_INDICES
+ unsigned int *vector_indices = (unsigned int *)vector->Map().MyGlobalElements();
+#else
+ size_type *vector_indices = (size_type *)vector->Map().MyGlobalElements64();
+#endif
+ is.add_indices(vector_indices, vector_indices+n_indices);
+ }
+ is.compress();
+
+ return is;
}
}
}
- a = read_write_2;
- b = read_write_1;
- c = read_write_2;
+ a.import(read_write_2, VectorOperation::insert);
+ b.import(read_write_1, VectorOperation::insert);
+ c.import(read_write_2, VectorOperation::insert);
- read_write_3 = a;
+ read_write_3.import(a, VectorOperation::insert);
if (rank==0)
{
for (unsigned int i=0; i<5; ++i)
- AssertThrow(read_write_2[i]==read_write_3[i],
- ExcMessage("Vector a has been modified."));
+ {
+ AssertThrow(read_write_2[i]==read_write_3[i],
+ ExcMessage("Vector a has been modified."));
+ }
}
else
{
ExcMessage("Vector a has been modified."));
}
- read_write_3 = b;
+ read_write_3.import(b, VectorOperation::insert);
if (rank==0)
{
for (unsigned int i=0; i<5; ++i)
ExcMessage("Vector b has been modified."));
}
- read_write_3 = c;
+ read_write_3.import(c, VectorOperation::insert);
if (rank==0)
{
for (unsigned int i=0; i<5; ++i)
a *= 2;
- read_write_3 = a;
+ read_write_3.import(a, VectorOperation::insert);
if (rank==0)
{
for (unsigned int i=0; i<5; ++i)
}
c /= 2.;
- read_write_3 = c;
+ read_write_3.import(c, VectorOperation::insert);
if (rank==0)
{
for (unsigned int i=0; i<5; ++i)
}
b += a;
- read_write_3 = b;
+ read_write_3.import(b, VectorOperation::insert);
if (rank==0)
{
for (unsigned int i=0; i<5; ++i)
}
b -= c;
- read_write_3 = b;
+ read_write_3.import(b, VectorOperation::insert);
if (rank==0)
{
for (unsigned int i=0; i<5; ++i)
ExcMessage("Problem in operator -=."));
}
- b = read_write_1;
- c = read_write_1;
+ b.import(read_write_1, VectorOperation::insert);
+ c.import(read_write_1, VectorOperation::insert);
const double val = b*c;
AssertThrow(val==285., ExcMessage("Problem in operator *."));
}
Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1);
+ test();
+
deallog << "OK" <<std::endl;
return 0;
}
}
- a = read_write_1;
- b = read_write_2;
- c = read_write_2;
+ a.import(read_write_1, VectorOperation::insert);
+ b.import(read_write_2, VectorOperation::insert);
+ c.import(read_write_2, VectorOperation::insert);
a.add(1.);
- read_write_3 = a;
+ read_write_3.import(a, VectorOperation::insert);
if (rank==0)
{
for (unsigned int i=0; i<5; ++i)
}
a.add(2.,b);
- read_write_3 = a;
+ read_write_3.import(a, VectorOperation::insert);
if (rank==0)
{
for (unsigned int i=0; i<5; ++i)
LinearAlgebra::EpetraWrappers::Vector d(a);
a.add(2.,b,3.,d);
- read_write_3 = a;
+ read_write_3.import(a, VectorOperation::insert);
if (rank==0)
{
for (unsigned int i=0; i<5; ++i)
- AssertThrow(3.+read_write_1[i]+8.*read_write_2[i]==read_write_3[i],
+ AssertThrow(4.+4.*read_write_1[i]+10.*read_write_2[i]==read_write_3[i],
ExcMessage("Problem in add(scalar,Vector,scalar,Vector)."));
}
else
{
for (unsigned int i=5; i<10; ++i)
- AssertThrow(3.+read_write_1[i]+8.*read_write_2[i]==read_write_3[i],
+ AssertThrow(4.+4.*read_write_1[i]+10.*read_write_2[i]==read_write_3[i],
ExcMessage("Problem in add(scalar,Vector,scalar,Vector)."));
}
- a = read_write_1;
+ a.import(read_write_1, VectorOperation::insert);
a.sadd(3.,2.,c);
- read_write_3 = a;
+ read_write_3.import(a, VectorOperation::insert);
if (rank==0)
{
for (unsigned int i=0; i<5; ++i)
- AssertThrow(3.+read_write_1[i]+2.*read_write_2[i]==read_write_3[i],
+ AssertThrow(3.*read_write_1[i]+2.*read_write_2[i]==read_write_3[i],
ExcMessage("Problem in sadd(scalar,scalar,Vector)."));
}
else
{
for (unsigned int i=5; i<10; ++i)
- AssertThrow(3.+read_write_1[i]+2.*read_write_2[i]==read_write_3[i],
+ AssertThrow(3.*read_write_1[i]+2.*read_write_2[i]==read_write_3[i],
ExcMessage("Problem in sadd(scalar,scalar,Vector)."));
}
- a = read_write_1;
+ a.import(read_write_1, VectorOperation::insert);
a.scale(b);
- read_write_3 = a;
+ read_write_3.import(a, VectorOperation::insert);
if (rank==0)
{
for (unsigned int i=0; i<5; ++i)
a.equ(2.,c);
- read_write_3 = a;
+ read_write_3.import(a, VectorOperation::insert);
if (rank==0)
{
for (unsigned int i=0; i<5; ++i)
}
- AssertThrow(b.l1_norm()==45., ExcMessage("Problem in l1_norm."));
+ AssertThrow(b.l1_norm()==95., ExcMessage("Problem in l1_norm."));
const double eps=1e-6;
- AssertThrow(std::fabs(b.l2_norm()-16.881943016)<eps,
+ AssertThrow(std::fabs(b.l2_norm()-31.3847096)<eps,
ExcMessage("Problem in l2_norm"));
- AssertThrow(b.linfty_norm()==9., ExcMessage("Problem in linfty_norm."));
+ AssertThrow(b.linfty_norm()==14., ExcMessage("Problem in linfty_norm."));
- const double val = a.add_and_dot(2.,c,b);
+ a.import(read_write_1, VectorOperation::insert);
+ const double val = a.add_and_dot(2.,a,b);
AssertThrow(val==1530., ExcMessage("Problem in add_and_dot"));
}
Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1);
+ test();
+
deallog << "OK" <<std::endl;
return 0;
readwrite_is.add_range(2,6);
readwrite_is.compress();
LinearAlgebra::ReadWriteVector<double> readwrite(readwrite_is);
- readwrite = tril_vector;
+ readwrite.import(tril_vector,VectorOperation::insert);
if (rank==0)
{
std::vector<double> comp(4);
AssertThrow(readwrite.local_element(i)==comp[i],
ExcMessage("Element not copied correctly"));
}
+
+ readwrite.import(tril_vector,VectorOperation::add);
+
+ if (rank==0)
+ {
+ std::vector<double> comp(4);
+ comp[0] = 0.;
+ comp[1] = 2.;
+ comp[2] = 12.;
+ comp[3] = 14.;
+ for (unsigned int i=0; i<4; ++i)
+ AssertThrow(readwrite.local_element(i)==comp[i],
+ ExcMessage("Element not copied correctly"));
+ }
+ MPI_Barrier(MPI_COMM_WORLD);
+ if (rank==1)
+ {
+ std::vector<double> comp(4);
+ comp[0] = 4.;
+ comp[1] = 6.;
+ comp[2] = 8.;
+ comp[3] = 10.;
+ for (unsigned int i=0; i<4; ++i)
+ AssertThrow(readwrite.local_element(i)==comp[i],
+ ExcMessage("Element not copied correctly"));
+ }
+
deallog << "OK" <<std::endl;
}
projects / dealii.git / commitdiff