Corrected an error in add() and set function of Trilinos sparse matrix. Some more...

git-svn-id: https://svn.dealii.org/trunk@17694 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/lac/include/lac/trilinos_block_vector.h b/deal.II/lac/include/lac/trilinos_block_vector.h
index ce144e0079453ff4dc72d9c1a3e8f7b946918b9a..13e27772f52f3e106ccec83bfa4ad590b0bd3ab0 100644 (file)
--- a/deal.II/lac/include/lac/trilinos_block_vector.h
+++ b/deal.II/lac/include/lac/trilinos_block_vector.h
@@ -277,8 +277,8 @@ namespace TrilinosWrappers
                                        * for example done during the
                                        * solution of linear systems.
                                        */
-       void do_data_exchange (const TrilinosWrappers::BlockSparseMatrix &m,
-                              const BlockVector                         &v);
+       void import_nonlocal_data_for_fe (const TrilinosWrappers::BlockSparseMatrix &m,
+                                         const BlockVector                         &v);
 
                                         /**
                                          * Compresses all the components

diff --git a/deal.II/lac/include/lac/trilinos_sparse_matrix.h b/deal.II/lac/include/lac/trilinos_sparse_matrix.h
index ff52dd92488918f068862130b514b6f403bed56e..d4a7b18ce6f1779c826d2e0bbfef47fb18149ac9 100755 (executable)
--- a/deal.II/lac/include/lac/trilinos_sparse_matrix.h
+++ b/deal.II/lac/include/lac/trilinos_sparse_matrix.h
@@ -2013,29 +2013,42 @@ namespace TrilinosWrappers
                                   // present in the input data.
     for (unsigned int i=0; i<n_rows; ++i)
       {
+       const int row = row_indices[i];
+
                                   // If the calling matrix owns the row to
-                                  // which we want to add values, we
+                                  // which we want to insert values, we
                                   // can directly call the Epetra_CrsMatrix
                                   // input function, which is much faster
-                                  // than the Epetra_FECrsMatrix function.
-       if (row_map.MyGID(i) == true)
+                                  // than the Epetra_FECrsMatrix
+                                   // function. We distinguish between two 
+                                  // cases: the first one is when the matrix
+                                  // is not filled (i.e., it is possible to 
+                                  // add new elements to the sparsity pattern), 
+                                  // and the second one is when the pattern is
+                                  // already fixed. In the former case, we 
+                                  // add the possibility to insert new values,
+                                  // and in the second we just replace
+                                  // data.
+       if (row_map.MyGID(row) == true)
          {
            if (matrix->Filled() == false)
              {
                ierr = matrix->Epetra_CrsMatrix::InsertGlobalValues(
-                                  (int)row_indices[i], (int)n_cols,
+                                  row, (int)n_cols,
                                   const_cast<TrilinosScalar*>(&values[i*n_cols]),
                                   (int*)&col_indices[0]);
 
-                                       // When adding up elements, we do
+                                       // When inserting elements, we do
                                        // not want to create exceptions in
-                                       // the case when adding elements.
+                                       // the case when inserting non-local
+                                       // data (since that's what we want 
+                                       // to do right now).
                if (ierr > 0)
                  ierr = 0;
              }
            else
              ierr = matrix->Epetra_CrsMatrix::ReplaceGlobalValues(
-                                  (int)row_indices[i], (int)n_cols,
+                                  row, (int)n_cols,
                                   const_cast<TrilinosScalar*>(&values[i*n_cols]),
                                   (int*)&col_indices[0]);
          }
@@ -2043,7 +2056,7 @@ namespace TrilinosWrappers
          {
                                   // When we're at off-processor data, we
                                   // have to stick with the standard
-                                  // SumIntoGlobalValues
+                                  // Insert/ReplaceGlobalValues
                                   // function. Nevertheless, the way we
                                   // call it is the fastest one (any other
                                   // will lead to repeated allocation and
@@ -2057,7 +2070,7 @@ namespace TrilinosWrappers
 
            if (matrix->Filled() == false)
              {
-               ierr = matrix->InsertGlobalValues (1, (int*)&i, 
+               ierr = matrix->InsertGlobalValues (1, &row, 
                                            (int)n_cols, (int*)&col_indices[0],
                                            &value_ptr, 
                                            Epetra_FECrsMatrix::ROW_MAJOR);
@@ -2065,7 +2078,7 @@ namespace TrilinosWrappers
                  ierr = 0;
              }
            else
-             ierr = matrix->ReplaceGlobalValues (1, (int*)&i, 
+             ierr = matrix->ReplaceGlobalValues (1, &row, 
                                            (int)n_cols, (int*)&col_indices[0],
                                            &value_ptr, 
                                            Epetra_FECrsMatrix::ROW_MAJOR);
@@ -2174,15 +2187,16 @@ namespace TrilinosWrappers
                                   // present in the input data.
     for (unsigned int i=0; i<n_rows; ++i)
       {
+       const int row = row_indices[i];
                                   // If the calling matrix owns the row to
                                   // which we want to add values, we
                                   // can directly call the Epetra_CrsMatrix
                                   // input function, which is much faster
                                   // than the Epetra_FECrsMatrix function.
-       if (row_map.MyGID(i) == true)
+       if (row_map.MyGID(row_indices[i]) == true)
          {
            ierr = matrix->Epetra_CrsMatrix::SumIntoGlobalValues(
-                                  (int)row_indices[i], (int)n_cols,
+                                  row, (int)n_cols,
                                   const_cast<TrilinosScalar*>(&values[i*n_cols]),
                                   (int*)&col_indices[0]);
          }
@@ -2202,7 +2216,7 @@ namespace TrilinosWrappers
       
            const TrilinosScalar* value_ptr = &values[i*n_cols];
 
-           ierr = matrix->SumIntoGlobalValues (1, (int*)&i, 
+           ierr = matrix->SumIntoGlobalValues (1, &row, 
                                                (int)n_cols, (int*)&col_indices[0],
                                                &value_ptr, 
                                                Epetra_FECrsMatrix::ROW_MAJOR);

diff --git a/deal.II/lac/include/lac/trilinos_vector.h b/deal.II/lac/include/lac/trilinos_vector.h
index d142b9f479e8c222ccad0897452d0dc8ed8766a2..e1240798809ce321669112fa82056a3b57e3b3e6 100755 (executable)
--- a/deal.II/lac/include/lac/trilinos_vector.h
+++ b/deal.II/lac/include/lac/trilinos_vector.h
@@ -368,8 +368,9 @@ namespace TrilinosWrappers
                                        * for example done during the
                                        * solution of linear systems.
                                        */
-       void do_data_exchange (const dealii::TrilinosWrappers::SparseMatrix &matrix,
-                              const Vector                                 &vector);
+       void import_nonlocal_data_for_fe 
+         (const dealii::TrilinosWrappers::SparseMatrix &matrix,
+          const Vector                                 &vector);
 
       private:
                                        /**

diff --git a/deal.II/lac/include/lac/trilinos_vector_base.h b/deal.II/lac/include/lac/trilinos_vector_base.h
index d36cc7ff8d462c8f7589cb757a0154a7583555c1..1d12d5a8f4785506ee1970725b8a66af5fc30cb4 100644 (file)
--- a/deal.II/lac/include/lac/trilinos_vector_base.h
+++ b/deal.II/lac/include/lac/trilinos_vector_base.h
@@ -1185,6 +1185,16 @@ namespace TrilinosWrappers
 
 
 
+  inline
+  std::pair<unsigned int, unsigned int>
+  VectorBase::local_range () const
+  {
+    int begin, end;
+    begin = vector->Map().MinMyGID();
+    end = vector->Map().MaxMyGID()+1;
+    return std::make_pair (begin, end);
+  }
+
 
 #endif // DOXYGEN
 

diff --git a/deal.II/lac/source/trilinos_block_vector.cc b/deal.II/lac/source/trilinos_block_vector.cc
index 369b99ecc941586981d1428d01ee363de861d1d4..a2d7bae017c4ba0a57589f8bab6552472d5c1477 100644 (file)
--- a/deal.II/lac/source/trilinos_block_vector.cc
+++ b/deal.II/lac/source/trilinos_block_vector.cc
@@ -139,8 +139,9 @@ namespace TrilinosWrappers
 
 
     void
-    BlockVector::do_data_exchange (const TrilinosWrappers::BlockSparseMatrix &m,
-                                  const BlockVector                         &v)
+    BlockVector::import_nonlocal_data_for_fe 
+      (const TrilinosWrappers::BlockSparseMatrix &m,
+       const BlockVector                         &v)
     {
       Assert (m.n_block_rows() == v.n_blocks(),
              ExcDimensionMismatch(m.n_block_rows(),v.n_blocks()));
@@ -154,7 +155,7 @@ namespace TrilinosWrappers
        }
 
       for (unsigned int i=0; i<this->n_blocks(); ++i)
-       components[i].do_data_exchange(m.block(i,i), v.block(i));
+       components[i].import_nonlocal_data_for_fe(m.block(i,i), v.block(i));
 
       collect_sizes();
     }

diff --git a/deal.II/lac/source/trilinos_sparse_matrix.cc b/deal.II/lac/source/trilinos_sparse_matrix.cc
index 674f7b3c97f8fde964c04e881dff98a31609e46f..18570f3f9db739bd6446c3bc11f7f511f0834652 100755 (executable)
--- a/deal.II/lac/source/trilinos_sparse_matrix.cc
+++ b/deal.II/lac/source/trilinos_sparse_matrix.cc
@@ -599,6 +599,7 @@ namespace TrilinosWrappers
   }
 
 
+
   void
   SparseMatrix::clear ()
   {

diff --git a/deal.II/lac/source/trilinos_vector.cc b/deal.II/lac/source/trilinos_vector.cc
index e61f8d68dee0af0ad36962ad3f70e63b67856a1d..ef8ac334046e37e7bb0bec8627b421b6efbfc91f 100755 (executable)
--- a/deal.II/lac/source/trilinos_vector.cc
+++ b/deal.II/lac/source/trilinos_vector.cc
@@ -124,7 +124,7 @@ namespace TrilinosWrappers
                                        // generate the vector.
       if (allow_different_maps == false)
         {
-         if (map.SameAs(v.vector->Map()) == false)
+         if (local_range() != v.local_range())
            {
              vector.reset();
              map = Epetra_Map(v.vector->Map().NumGlobalElements(),
@@ -138,6 +138,10 @@ namespace TrilinosWrappers
          else if (fast == false)
            {
              int ierr;
+             Assert (vector->Map().SameAs(v.vector->Map()) == true,
+                     ExcMessage ("The Epetra maps in the assignment operator ="
+                                 " do not match, even though the local_range "
+                                 " seems to be the same. Check vector setup!"));
              ierr = vector->GlobalAssemble (last_action);      
              AssertThrow (ierr == 0, ExcTrilinosError(ierr));
 
@@ -181,8 +185,12 @@ namespace TrilinosWrappers
     Vector &
     Vector::operator = (const Vector &v)
     {
-      if (vector->Map().SameAs(v.vector->Map()) == true)
+      if (local_range() == v.local_range())
        {
+         Assert (vector->Map().SameAs(v.vector->Map()) == true,
+                 ExcMessage ("The Epetra maps in the assignment operator ="
+                             " do not match, even though the local_range "
+                             " seems to be the same. Check vector setup!"));
          vector->GlobalAssemble(last_action);
          const int ierr = vector->Update(1.0, *v.vector, 0.0);
          AssertThrow (ierr == 0, ExcTrilinosError(ierr));
@@ -224,8 +232,8 @@ namespace TrilinosWrappers
 
 
     void
-    Vector::do_data_exchange (const TrilinosWrappers::SparseMatrix &m,
-                             const Vector                         &v)
+    Vector::import_nonlocal_data_for_fe (const TrilinosWrappers::SparseMatrix &m,
+                                        const Vector                         &v)
     {
       Assert (m.matrix->Filled() == true,
              ExcMessage ("Matrix is not compressed. "
@@ -382,7 +390,7 @@ namespace TrilinosWrappers
                                        // generate the vector.
     if (allow_different_maps == false)
       {
-       if (map.SameAs(v.vector->Map()) == false)
+       if (local_range() != v.local_range())
          {
            vector.reset();
            map = Epetra_LocalMap (v.vector->GlobalLength(),
@@ -394,6 +402,11 @@ namespace TrilinosWrappers
        else
          {
            int ierr;
+           Assert (vector->Map().SameAs(v.vector->Map()) == true,
+                   ExcMessage ("The Epetra maps in the assignment operator ="
+                               " do not match, even though the local_range "
+                               " seems to be the same. Check vector setup!"));
+
            ierr = vector->GlobalAssemble(last_action);
            AssertThrow (ierr == 0, ExcTrilinosError(ierr));
 
@@ -452,7 +465,7 @@ namespace TrilinosWrappers
   Vector &
   Vector::operator = (const Vector &v)
   {
-    if (vector->Map().SameAs(v.vector->Map()) == false)
+    if (size() != v.size())
       {
        map = Epetra_LocalMap (v.vector->Map().NumGlobalElements(), 
                               v.vector->Map().IndexBase(),

diff --git a/deal.II/lac/source/trilinos_vector_base.cc b/deal.II/lac/source/trilinos_vector_base.cc
index 3fa69f8bc611f5413c2ebf3c81e8cffbbb18a14a..37ebf97a76bc0e905f5fa31263485bdd5257cbea 100644 (file)
--- a/deal.II/lac/source/trilinos_vector_base.cc
+++ b/deal.II/lac/source/trilinos_vector_base.cc
@@ -110,7 +110,7 @@ namespace TrilinosWrappers
     Assert (&*vector != 0,
            ExcMessage("Vector has not been constructed properly."));
 
-    if (fast == false || vector->Map().SameAs(v.vector->Map()) == false)
+    if (fast == false || local_range() != v.local_range())
       vector = std::auto_ptr<Epetra_FEVector>(new Epetra_FEVector(*v.vector));
   }
 
@@ -155,7 +155,7 @@ namespace TrilinosWrappers
     Assert (&*vector != 0,
            ExcMessage("Vector is not constructed properly."));
 
-    if (vector->Map().SameAs(v.vector->Map()) == false)
+    if (local_range() != v.local_range())
       {
        vector.reset();
        last_action = Zero;
@@ -163,6 +163,10 @@ namespace TrilinosWrappers
       }
     else
       {
+       Assert (vector->Map().SameAs(v.vector->Map()) == true,
+               ExcMessage ("The Epetra maps in the assignment operator ="
+                           " do not match, even though the local_range "
+                           " seems to be the same. Check vector setup!"));
        int ierr;
        ierr = vector->GlobalAssemble(last_action);
        AssertThrow (ierr == 0, ExcTrilinosError(ierr));
@@ -254,17 +258,6 @@ namespace TrilinosWrappers
 
 
 
-  std::pair<unsigned int, unsigned int>
-  VectorBase::local_range () const
-  {
-    int begin, end;
-    begin = vector->Map().MinMyGID();
-    end = vector->Map().MaxMyGID()+1;
-    return std::make_pair (begin, end);
-  }
-
-
-
   TrilinosScalar
   VectorBase::el (const unsigned int index) const
   {