When initializing the Trilinos sparse matrix, use an Epetra_CrsGraph as in intermedia...

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 d45ebca374019c2b4a9b561369006b498aeeff56..d8fed0969dc3aaad986f97e19799bc417979e03c 100755 (executable)
--- 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 37ebf97a76bc0e905f5fa31263485bdd5257cbea..6147ec10edcc10071991001ce766d93c761ad52a 100644 (file)
--- 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;
+      }
   }