From: Martin Kronbichler <kronbichler@lnm.mw.tum.de>
Date: Mon, 24 Nov 2008 12:02:40 +0000 (+0000)
Subject: Corrected an error in add() and set function of Trilinos sparse matrix. Some more... 
X-Git-Tag: v8.0.0~8350
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=923475b3e7d28e0bf53992bfda19c702598292f2;p=dealii.git

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

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 ce144e0079..13e27772f5 100644
--- 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 ff52dd9248..d4a7b18ce6 100755
--- 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 d142b9f479..e124079880 100755
--- 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 d36cc7ff8d..1d12d5a8f4 100644
--- 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 369b99ecc9..a2d7bae017 100644
--- 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 674f7b3c97..18570f3f9d 100755
--- 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 e61f8d68de..ef8ac33404 100755
--- 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 3fa69f8bc6..37ebf97a76 100644
--- 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
   {