From: kronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Tue, 25 Nov 2008 09:59:11 +0000 (+0000)
Subject: When initializing the Trilinos sparse matrix, use an Epetra_CrsGraph as in intermedia... 
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=09537d22aa55d8d1ef1810a7607240d3c09744c5;p=dealii-svn.git

When initializing the Trilinos sparse matrix, use an Epetra_CrsGraph as in intermediate object. That saves some computing time in setup.

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

diff --git a/deal.II/lac/source/trilinos_sparse_matrix.cc b/deal.II/lac/source/trilinos_sparse_matrix.cc
index d45ebca374..d8fed0969d 100755
--- a/deal.II/lac/source/trilinos_sparse_matrix.cc
+++ b/deal.II/lac/source/trilinos_sparse_matrix.cc
@@ -256,7 +256,7 @@ namespace TrilinosWrappers
   {
     matrix.reset();
 
-    unsigned int n_rows = sparsity_pattern.n_rows();
+    const unsigned int n_rows = sparsity_pattern.n_rows();
 
     if (input_row_map.Comm().MyPID() == 0)
       {
@@ -276,21 +276,28 @@ namespace TrilinosWrappers
     for (unsigned int row=0; row<n_rows; ++row)
       n_entries_per_row[row] = sparsity_pattern.row_length(row);
 
-				  // TODO: There seems to be problem
-				  // in Epetra when a matrix is
-				  // initialized with both row and
-				  // column map. Maybe find something
-				  // more out about this... It could
-				  // be related to the bug #4123. For
-				  // the moment, just ignore the
-				  // column information and generate
-				  // the matrix as if it were
-				  // square. The call to
+				  // The deal.II notation of a Sparsity
+				  // pattern corresponds to the Epetra
+				  // concept of a Graph. Hence, we generate
+				  // a graph by copying the sparsity pattern
+				  // into it, and then build up the matrix
+				  // from the graph. This is considerable
+				  // faster than directly filling elements
+				  // into the matrix.
+
+				  // TODO: There seems to be problem in
+				  // Epetra when a Graph/matrix is
+				  // initialized with both row and column
+				  // map. Maybe find something more out
+				  // about this... It could be related to
+				  // the bug #4123. For the moment, just
+				  // ignore the column information and
+				  // generate the graph to the matrix as if
+				  // it were square. The call to
 				  // GlobalAssemble later will set the
 				  // correct values.
-    matrix = std::auto_ptr<Epetra_FECrsMatrix>
-	        (new Epetra_FECrsMatrix(Copy, row_map,
-					&n_entries_per_row[0], false));
+    Epetra_CrsGraph graph (Copy, row_map, 
+			   &n_entries_per_row[input_row_map.MinMyGID()], true);
 
 
 				  // TODO: As of now, assume that the
@@ -299,10 +306,10 @@ namespace TrilinosWrappers
 				  // each processor set its rows. Since
 				  // this is wasteful, a better solution
 				  // should be found in the future.
-    Assert (matrix->NumGlobalRows() == (int)sparsity_pattern.n_rows(),
-	    ExcDimensionMismatch (matrix->NumGlobalRows(),
-				  sparsity_pattern.n_rows()));
-    
+    Assert (graph.NumGlobalRows() == (int)sparsity_pattern.n_rows(),
+    	    ExcDimensionMismatch (graph.NumGlobalRows(),
+    				  sparsity_pattern.n_rows()));
+
 				  // Trilinos has a bug for rectangular
 				  // matrices at this point, so do not
 				  // check for consistent column numbers
@@ -311,26 +318,32 @@ namespace TrilinosWrappers
 				  // this bug is filed in the Sandia
 				  // bugzilla under #4123 and should be
 				  // fixed for version 9.0
-//        Assert (matrix->NumGlobalCols() == (int)sparsity_pattern.n_cols(),
-//		ExcDimensionMismatch (matrix->NumGlobalCols(),
+//        Assert (graph.NumGlobalCols() == (int)sparsity_pattern.n_cols(),
+//		ExcDimensionMismatch (graph.NumGlobalCols(),
 //				      sparsity_pattern.n_cols()));
 
-    std::vector<TrilinosScalar> values;
-    std::vector<unsigned int>   row_indices;
+    std::vector<int>   row_indices;
     
     for (unsigned int row=0; row<n_rows; ++row)
       if (row_map.MyGID(row))
 	{
 	  const int row_length = sparsity_pattern.row_length(row);
 	  row_indices.resize (row_length, -1);
-	  values.resize (row_length, 0.);
 
 	  for (int col=0; col < row_length; ++col)
 	    row_indices[col] = sparsity_pattern.column_number (row, col);
-	  
-	  set (row, row_length, &row_indices[0], &values[0], false);
+
+	  graph.InsertGlobalIndices (row, row_length, &row_indices[0]);
 	}
-    
+
+				  // Now, fill the graph (sort indices, make
+				  // memory contiguous, etc).
+    graph.FillComplete(col_map, row_map);
+
+				  // And now finally generate the matrix.
+    matrix =  std::auto_ptr<Epetra_FECrsMatrix>
+                       (new Epetra_FECrsMatrix(Copy, graph, false));
+
     last_action = Zero;
 
 				  // In the end, the matrix needs to
@@ -354,7 +367,7 @@ namespace TrilinosWrappers
   {
     matrix.reset();
 
-    unsigned int n_rows = sparsity_pattern.n_rows();
+    const unsigned int n_rows = sparsity_pattern.n_rows();
 
     if (input_row_map.Comm().MyPID() == 0)
       {
@@ -374,21 +387,29 @@ namespace TrilinosWrappers
     for (unsigned int row=0; row<n_rows; ++row)
       n_entries_per_row[row] = sparsity_pattern.row_length(row);
 
-				  // TODO: There seems to be problem
-				  // in Epetra when a matrix is
-				  // initialized with both row and
-				  // column map. Maybe find something
-				  // more out about this... It could
-				  // be related to the bug #4123. For
-				  // the moment, just ignore the
-				  // column information and generate
-				  // the matrix as if it were
-				  // square. The call to
+
+				  // The deal.II notation of a Sparsity
+				  // pattern corresponds to the Epetra
+				  // concept of a Graph. Hence, we generate
+				  // a graph by copying the sparsity pattern
+				  // into it, and then build up the matrix
+				  // from the graph. This is considerable
+				  // faster than directly filling elements
+				  // into the matrix.
+
+				  // TODO: There seems to be problem in
+				  // Epetra when a Graph/matrix is
+				  // initialized with both row and column
+				  // map. Maybe find something more out
+				  // about this... It could be related to
+				  // the bug #4123. For the moment, just
+				  // ignore the column information and
+				  // generate the graph to the matrix as if
+				  // it were square. The call to
 				  // GlobalAssemble later will set the
 				  // correct values.
-    matrix = std::auto_ptr<Epetra_FECrsMatrix>
-	        (new Epetra_FECrsMatrix(Copy, row_map,
-					&n_entries_per_row[0], false));
+    Epetra_CrsGraph graph (Copy, row_map, 
+			   &n_entries_per_row[input_row_map.MinMyGID()], true);
 
 
 				  // TODO: As of now, assume that the
@@ -397,10 +418,10 @@ namespace TrilinosWrappers
 				  // each processor set its rows. Since
 				  // this is wasteful, a better solution
 				  // should be found in the future.
-    Assert (matrix->NumGlobalRows() == (int)sparsity_pattern.n_rows(),
-	    ExcDimensionMismatch (matrix->NumGlobalRows(),
-				  sparsity_pattern.n_rows()));
-    
+    Assert (graph.NumGlobalRows() == (int)sparsity_pattern.n_rows(),
+    	    ExcDimensionMismatch (graph.NumGlobalRows(),
+    				  sparsity_pattern.n_rows()));
+
 				  // Trilinos has a bug for rectangular
 				  // matrices at this point, so do not
 				  // check for consistent column numbers
@@ -409,19 +430,17 @@ namespace TrilinosWrappers
 				  // this bug is filed in the Sandia
 				  // bugzilla under #4123 and should be
 				  // fixed for version 9.0
-//        Assert (matrix->NumGlobalCols() == (int)sparsity_pattern.n_cols(),
-//		ExcDimensionMismatch (matrix->NumGlobalCols(),
+//        Assert (graph.NumGlobalCols() == (int)sparsity_pattern.n_cols(),
+//		ExcDimensionMismatch (graph.NumGlobalCols(),
 //				      sparsity_pattern.n_cols()));
 
-    std::vector<TrilinosScalar> values;
-    std::vector<unsigned int>   row_indices;
+    std::vector<int>   row_indices;
     
     for (unsigned int row=0; row<n_rows; ++row)
       if (row_map.MyGID(row))
 	{
 	  const int row_length = sparsity_pattern.row_length(row);
 	  row_indices.resize (row_length, -1);
-	  values.resize (row_length, 0.);
 
 	  CompressedSetSparsityPattern::row_iterator col_num = 
 	    sparsity_pattern.row_begin (row);
@@ -431,8 +450,16 @@ namespace TrilinosWrappers
 	       ++col_num, ++col)
 	    row_indices[col] = *col_num;
 
-	  set (row, row_length, &row_indices[0], &values[0], false);
+	  graph.InsertGlobalIndices (row, row_length, &row_indices[0]);
 	}
+
+				  // Now, fill the graph (sort indices, make
+				  // memory contiguous, etc).
+    graph.FillComplete(col_map, row_map);
+
+				  // And now finally generate the matrix.
+    matrix =  std::auto_ptr<Epetra_FECrsMatrix>
+                       (new Epetra_FECrsMatrix(Copy, graph, false));
     
     last_action = Zero;
 
diff --git a/deal.II/lac/source/trilinos_vector_base.cc b/deal.II/lac/source/trilinos_vector_base.cc
index 37ebf97a76..6147ec10ed 100644
--- a/deal.II/lac/source/trilinos_vector_base.cc
+++ b/deal.II/lac/source/trilinos_vector_base.cc
@@ -315,7 +315,7 @@ namespace TrilinosWrappers
   TrilinosScalar
   VectorBase::operator * (const VectorBase &vec) const
   {
-    Assert (size() == vec.size(),
+    Assert (local_range() == vec.local_range(),
 	    ExcDimensionMismatch(size(), vec.size()));
 
     TrilinosScalar result;
@@ -766,9 +766,31 @@ namespace TrilinosWrappers
 	    ExcMessage("The given value is not finite but "
 		       "either infinite or Not A Number (NaN)"));
 
-    *vector = *v.vector;
+				   // If we don't have the same map, copy.
+    if (local_range() != v.local_range())
+      {
+	vector.reset();
+	last_action = Zero;
+	*vector = *v.vector;
+	*this *= a;
+      }
+    else
+      {
+				   // Otherwise, just update
+	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));
+
+	ierr = vector->Update(a, *v.vector, 0.0);
+	AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+	last_action = Zero;
+      }
 
-    *this *= a;
   }
 
 
@@ -790,9 +812,30 @@ namespace TrilinosWrappers
 	    ExcMessage("The given value is not finite but "
 		       "either infinite or Not A Number (NaN)"));
 
-    *vector = *v.vector;
+				   // If we don't have the same map, copy.
+    if (local_range() != v.local_range())
+      {
+	vector.reset();
+	last_action = Zero;
+	*vector = *v.vector;
+	sadd (a, b, w);
+      }
+    else
+      {
+				   // Otherwise, just update
+	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));
 
-    sadd (a, b, w);
+	ierr = vector->Update(a, *v.vector, b, *w.vector, 0.0);
+	AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+	last_action = Zero;
+      }
   }