]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Make tests work again.
authorBruno Turcksin <bruno.turcksin@gmail.com>
Mon, 15 Feb 2016 16:18:45 +0000 (11:18 -0500)
committerBruno Turcksin <bruno.turcksin@gmail.com>
Wed, 23 Mar 2016 13:33:42 +0000 (09:33 -0400)
include/deal.II/lac/la_vector.h
include/deal.II/lac/la_vector.templates.h
include/deal.II/lac/read_write_vector.h
include/deal.II/lac/read_write_vector.templates.h
include/deal.II/lac/trilinos_epetra_vector.h
include/deal.II/lac/vector_space_vector.h
source/lac/trilinos_epetra_communication_pattern.cc
source/lac/trilinos_epetra_vector.cc
tests/trilinos/epetra_vector_01.cc
tests/trilinos/epetra_vector_02.cc
tests/trilinos/readwritevector.cc

index 2ce6a28c82584c5f612cd45ae734efce49bcc1a0..fbdb0a75c73cea2c37dfa74269188ea8b40416cf 100644 (file)
@@ -132,6 +132,13 @@ namespace LinearAlgebra
      */
     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.
      */
index bc014e9daeb23e75d25b2cf904fe233755615d9a..a6ea7c6b52d0f7d060c7c87d8e12e53dd2625852 100644 (file)
@@ -114,6 +114,16 @@ namespace LinearAlgebra
 
 
 
+  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)
   {
index cb61c2c9a705a3055453a2c251955497a9210394..197ead3db0fd48fa16219f2dcaced65c5c4e1c1c 100644 (file)
@@ -233,7 +233,7 @@ namespace LinearAlgebra
      */
     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
@@ -247,7 +247,7 @@ namespace LinearAlgebra
      */
     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
@@ -259,7 +259,7 @@ namespace LinearAlgebra
      */
     void import(const EpetraWrappers::Vector &epetra_vec,
                 VectorOperation::values operation,
-                const CommunicationPatternBase *communication_pattern = NULL);
+                const CommunicationPatternBase *communication_pattern = nullptr);
 #endif
 
     /**
@@ -446,7 +446,7 @@ namespace LinearAlgebra
   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.
      */
index e3a7ea8a56d7d30c8a62a2f991c996aba053e69a..26c0c360deeb2aab511b9bcad4d606405adc32e0 100644 (file)
@@ -191,10 +191,11 @@ namespace LinearAlgebra
         // 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));
@@ -323,8 +324,8 @@ namespace LinearAlgebra
     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;
   }
index 67fd2efe714bf08ec6302e68735daa7743479d9f..b46e0530095fa233020641d0ae634cc308aafb51 100644 (file)
@@ -202,6 +202,8 @@ namespace LinearAlgebra
        * 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,
index 9f1f57ddc509ad0ba3c36f8612ee0b70d5ac5a93..577da9b9d212e0c7bb311e7688159552c3ec40a6 100644 (file)
@@ -70,7 +70,7 @@ namespace LinearAlgebra
     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
@@ -79,7 +79,7 @@ namespace LinearAlgebra
      */
     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.
index ef93fc9db1793f1035c8a01f0ea65e5847047515..992d2d9c6841abfd18e4b72e6220b6de29764beb 100644 (file)
@@ -43,8 +43,10 @@ namespace LinearAlgebra
     {
       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
index b9ef42e3aeb1f056487424f44084a2e3b0ee9c6f..b2b4cab8c6a66f8cc228ed0acd7056d8cccac02f 100644 (file)
@@ -33,7 +33,7 @@ namespace LinearAlgebra
   {
     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)
     {}
 
@@ -110,7 +110,9 @@ namespace LinearAlgebra
           // 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());
@@ -414,8 +416,30 @@ namespace LinearAlgebra
 
     ::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;
     }
 
 
index 3059910224c34fb9afb480e443aa3e63af8325c6..1192ef53826af80552c3f256d0af23d599578354 100644 (file)
@@ -86,16 +86,18 @@ void test()
         }
     }
 
-  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
     {
@@ -104,7 +106,7 @@ void test()
                     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)
@@ -118,7 +120,7 @@ void test()
                     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)
@@ -134,7 +136,7 @@ void test()
 
 
   a *= 2;
-  read_write_3 = a;
+  read_write_3.import(a, VectorOperation::insert);
   if (rank==0)
     {
       for (unsigned int i=0; i<5; ++i)
@@ -149,7 +151,7 @@ void test()
     }
 
   c /= 2.;
-  read_write_3 = c;
+  read_write_3.import(c, VectorOperation::insert);
   if (rank==0)
     {
       for (unsigned int i=0; i<5; ++i)
@@ -164,7 +166,7 @@ void test()
     }
 
   b += a;
-  read_write_3 = b;
+  read_write_3.import(b, VectorOperation::insert);
   if (rank==0)
     {
       for (unsigned int i=0; i<5; ++i)
@@ -179,7 +181,7 @@ void test()
     }
 
   b -= c;
-  read_write_3 = b;
+  read_write_3.import(b, VectorOperation::insert);
   if (rank==0)
     {
       for (unsigned int i=0; i<5; ++i)
@@ -193,8 +195,8 @@ void test()
                     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 *."));
 }
@@ -209,6 +211,8 @@ int main(int argc, char **argv)
 
   Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1);
 
+  test();
+
   deallog << "OK" <<std::endl;
 
   return 0;
index 687038a39397dc98fc3ac61545f0a8702b165049..437e11284ca75d8d3e6bec9dc0afda5521140302 100644 (file)
@@ -71,12 +71,12 @@ void test()
         }
     }
 
-  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)
@@ -91,7 +91,7 @@ void test()
     }
 
   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)
@@ -108,41 +108,41 @@ void test()
 
   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)
@@ -158,7 +158,7 @@ void test()
 
 
   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)
@@ -173,15 +173,16 @@ void test()
     }
 
 
-  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"));
 }
 
@@ -195,6 +196,8 @@ int main(int argc, char **argv)
 
   Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1);
 
+  test();
+
   deallog << "OK" <<std::endl;
 
   return 0;
index 1b55c2d77f37d78fb755dc5626115c78abac03ed..069320ca24bfbd0c22c7913fa0e695fec182c7c6 100644 (file)
@@ -54,7 +54,7 @@ void test()
     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);
@@ -78,6 +78,33 @@ void test()
         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;
 }
 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.