Optimized constructor of Trilinos matrices on one processor. This makes assembly...

git-svn-id: https://svn.dealii.org/trunk@18322 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 5af9a75a7016b45052c8b847fb6bd58684396ada..580bcb552b52286b7a9bb23718b82c800af6f640 100755 (executable)
--- a/deal.II/lac/source/trilinos_sparse_matrix.cc
+++ b/deal.II/lac/source/trilinos_sparse_matrix.cc
@@ -93,7 +93,7 @@ namespace TrilinosWrappers
                  last_action (Zero),
                  compressed (true),
                  matrix (std::auto_ptr<Epetra_FECrsMatrix>
-                               (new Epetra_FECrsMatrix(View, row_map, 0)))
+                         (new Epetra_FECrsMatrix(View, row_map, col_map, 0)))
   {
     matrix->FillComplete();
   }
@@ -153,7 +153,7 @@ namespace TrilinosWrappers
                  last_action (Zero),
                  compressed (true),
                  matrix (std::auto_ptr<Epetra_FECrsMatrix>
-                   (new Epetra_FECrsMatrix(Copy, row_map, 
+                   (new Epetra_FECrsMatrix(Copy, row_map,
                      (int*)const_cast<unsigned int*>(&(n_entries_per_row[0])),
                                            false)))
   {}
@@ -173,9 +173,18 @@ namespace TrilinosWrappers
 #endif
                  last_action (Zero),
                  compressed (true),
+                                  // on one processor only, we know how the
+                                  // columns of the matrix will be
+                                  // distributed (everything on one
+                                  // processor), so we can hand in this
+                                  // information to the constructor. we
+                                  // can't do so in parallel, where the
+                                  // information from columns is only
+                                  // available when entries have been added
                  matrix (std::auto_ptr<Epetra_FECrsMatrix>
-                               (new Epetra_FECrsMatrix(Copy, row_map, 
-                                       int(n_max_entries_per_row), false)))
+                         (new Epetra_FECrsMatrix(Copy, row_map, col_map,
+                                                 int(n_max_entries_per_row), 
+                                                 false)))
   {}
 
 
@@ -194,9 +203,9 @@ namespace TrilinosWrappers
                  last_action (Zero),
                  compressed (true),
                  matrix (std::auto_ptr<Epetra_FECrsMatrix>
-                    (new Epetra_FECrsMatrix(Copy, row_map, 
-                       (int*)const_cast<unsigned int*>(&(n_entries_per_row[0])), 
-                                            false)))
+                         (new Epetra_FECrsMatrix(Copy, row_map, col_map, 
+                          (int*)const_cast<unsigned int*>(&(n_entries_per_row[0])), 
+                                                 false)))
   {}
 
 
@@ -232,8 +241,7 @@ namespace TrilinosWrappers
 
 
   SparseMatrix::~SparseMatrix ()
-  {
-  }
+  {}
 
 
 
@@ -331,9 +339,26 @@ namespace TrilinosWrappers
                                  // it were square. The call to
                                  // GlobalAssemble later will set the
                                  // correct values.
-    Epetra_CrsGraph graph (Copy, row_map, 
-                          &n_entries_per_row[input_row_map.MinMyGID()], true);
 
+                                  // for more than one processor, need to
+                                  // specify only row map first and let the
+                                  // matrix entries decide about the column
+                                  // map (which says which columns are
+                                  // present in the matrix, not to be
+                                  // confused with the col_map that tells
+                                  // how the domain dofs of the matrix will
+                                  // be distributed). for only one
+                                  // processor, we can directly assign the
+                                  // columns as well.
+    std::auto_ptr<Epetra_CrsGraph> graph;
+    if (row_map.Comm().NumProc() > 1)
+      graph = std::auto_ptr<Epetra_CrsGraph> 
+       (new Epetra_CrsGraph (Copy, row_map, 
+                             &n_entries_per_row[input_row_map.MinMyGID()], true));
+    else
+      graph = std::auto_ptr<Epetra_CrsGraph> 
+       (new Epetra_CrsGraph (Copy, row_map, col_map,
+                             &n_entries_per_row[input_row_map.MinMyGID()], true));
 
                                  // TODO: As of now, assume that the
                                  // sparsity pattern sits at the all
@@ -341,8 +366,8 @@ namespace TrilinosWrappers
                                  // each processor set its rows. Since
                                  // this is wasteful, a better solution
                                  // should be found in the future.
-    Assert (graph.NumGlobalRows() == (int)sparsity_pattern.n_rows(),
-           ExcDimensionMismatch (graph.NumGlobalRows(),
+    Assert (graph->NumGlobalRows() == (int)sparsity_pattern.n_rows(),
+           ExcDimensionMismatch (graph->NumGlobalRows(),
                                  sparsity_pattern.n_rows()));
 
                                  // Trilinos has a bug for rectangular
@@ -368,16 +393,16 @@ namespace TrilinosWrappers
          for (int col=0; col < row_length; ++col)
            row_indices[col] = sparsity_pattern.column_number (row, col);
 
-         graph.InsertGlobalIndices (row, row_length, &row_indices[0]);
+         graph->InsertGlobalIndices (row, row_length, &row_indices[0]);
        }
 
                                  // Now, fill the graph (sort indices, make
                                  // memory contiguous, etc).
-    graph.FillComplete(col_map, row_map);
+    graph->FillComplete(col_map, row_map);
 
                                  // And now finally generate the matrix.
     matrix =  std::auto_ptr<Epetra_FECrsMatrix>
-                       (new Epetra_FECrsMatrix(Copy, graph, false));
+                       (new Epetra_FECrsMatrix(Copy, *graph, false));
 
     last_action = Zero;
 
@@ -443,9 +468,15 @@ namespace TrilinosWrappers
                                  // it were square. The call to
                                  // GlobalAssemble later will set the
                                  // correct values.
-    Epetra_CrsGraph graph (Copy, row_map, 
-                          &n_entries_per_row[input_row_map.MinMyGID()], true);
-
+    std::auto_ptr<Epetra_CrsGraph> graph;
+    if (row_map.Comm().NumProc() > 1)
+      graph = std::auto_ptr<Epetra_CrsGraph> 
+       (new Epetra_CrsGraph (Copy, row_map, 
+                             &n_entries_per_row[input_row_map.MinMyGID()], true));
+    else
+      graph = std::auto_ptr<Epetra_CrsGraph> 
+       (new Epetra_CrsGraph (Copy, row_map, col_map,
+                             &n_entries_per_row[input_row_map.MinMyGID()], true));
 
                                  // TODO: As of now, assume that the
                                  // sparsity pattern sits at the all
@@ -453,8 +484,8 @@ namespace TrilinosWrappers
                                  // each processor set its rows. Since
                                  // this is wasteful, a better solution
                                  // should be found in the future.
-    Assert (graph.NumGlobalRows() == (int)sparsity_pattern.n_rows(),
-           ExcDimensionMismatch (graph.NumGlobalRows(),
+    Assert (graph->NumGlobalRows() == (int)sparsity_pattern.n_rows(),
+           ExcDimensionMismatch (graph->NumGlobalRows(),
                                  sparsity_pattern.n_rows()));
 
                                  // Trilinos has a bug for rectangular
@@ -485,16 +516,16 @@ namespace TrilinosWrappers
               ++col_num, ++col)
            row_indices[col] = *col_num;
 
-         graph.InsertGlobalIndices (row, row_length, &row_indices[0]);
+         graph->InsertGlobalIndices (row, row_length, &row_indices[0]);
        }
 
                                  // Now, fill the graph (sort indices, make
                                  // memory contiguous, etc).
-    graph.FillComplete(col_map, row_map);
+    graph->FillComplete(col_map, row_map);
 
                                  // And now finally generate the matrix.
     matrix =  std::auto_ptr<Epetra_FECrsMatrix>
-                       (new Epetra_FECrsMatrix(Copy, graph, false));
+                       (new Epetra_FECrsMatrix(Copy, *graph, false));
     
     last_action = Zero;
 

diff --git a/deal.II/lac/source/trilinos_sparsity_pattern.cc b/deal.II/lac/source/trilinos_sparsity_pattern.cc
index 7b005997cf69c4cefa354e7c882d31219f66e1de..a4de59a8f6923b3d9a13ab05acf5b2a90eb190fc 100755 (executable)
--- a/deal.II/lac/source/trilinos_sparsity_pattern.cc
+++ b/deal.II/lac/source/trilinos_sparsity_pattern.cc
@@ -90,7 +90,7 @@ namespace TrilinosWrappers
                  col_map (row_map),
                  compressed (true),
                  graph (std::auto_ptr<Epetra_FECrsGraph>
-                        (new Epetra_FECrsGraph(View, row_map, 0))),
+                        (new Epetra_FECrsGraph(View, row_map, col_map, 0))),
                  cached_row_indices (1),
                  n_cached_elements (0),
                  row_in_cache (0)
@@ -104,24 +104,51 @@ namespace TrilinosWrappers
                   row_map (InputMap),
                  col_map (row_map),
                  compressed (false),
-                 graph (std::auto_ptr<Epetra_FECrsGraph>
-                                 (new Epetra_FECrsGraph(Copy, row_map, 
-                                       int(n_entries_per_row), false))),
+                                  // for more than one processor, need to
+                                  // specify only row map first and let the
+                                  // matrix entries decide about the column
+                                  // map (which says which columns are
+                                  // present in the matrix, not to be
+                                  // confused with the col_map that tells
+                                  // how the domain dofs of the matrix will
+                                  // be distributed). for only one
+                                  // processor, we can directly assign the
+                                  // columns as well.
+                 graph (row_map.Comm().NumProc() > 1 ?
+                        (std::auto_ptr<Epetra_FECrsGraph>
+                                  (new Epetra_FECrsGraph(Copy, row_map, 
+                                                         (int)n_entries_per_row, 
+                                                         false)))
+                        :
+                        (std::auto_ptr<Epetra_FECrsGraph>
+                                  (new Epetra_FECrsGraph(Copy, row_map, col_map, 
+                                                         (int)n_entries_per_row, 
+                                                         false)))
+                        ),
                  cached_row_indices (1),
                  n_cached_elements (0),
                  row_in_cache (0)
   {}
-
+ 
   SparsityPattern::SparsityPattern (const Epetra_Map                &InputMap,
                                    const std::vector<unsigned int> &n_entries_per_row)
                  :
                   row_map (InputMap),
                  col_map (row_map),
                  compressed (false),
-                 graph (std::auto_ptr<Epetra_FECrsGraph>
-                   (new Epetra_FECrsGraph(Copy, row_map, 
-                     (int*)const_cast<unsigned int*>(&(n_entries_per_row[0])),
-                                           false))),
+                 graph (row_map.Comm().NumProc() > 1 ?
+                        (std::auto_ptr<Epetra_FECrsGraph>
+                         (new Epetra_FECrsGraph(Copy, row_map, 
+                                                (int*)const_cast<unsigned int*>
+                                                (&(n_entries_per_row[row_map.MinMyGID()])),
+                                                false)))
+                        :
+                        (std::auto_ptr<Epetra_FECrsGraph>
+                         (new Epetra_FECrsGraph(Copy, row_map, col_map, 
+                                                (int*)const_cast<unsigned int*>
+                                                (&(n_entries_per_row[row_map.MinMyGID()])),
+                                                false)))
+                        ),
                  cached_row_indices (1),
                  n_cached_elements (0),
                  row_in_cache (0)
@@ -134,9 +161,17 @@ namespace TrilinosWrappers
                   row_map (InputRowMap),
                   col_map (InputColMap),
                  compressed (false),
-                 graph (std::auto_ptr<Epetra_FECrsGraph>
-                                 (new Epetra_FECrsGraph(Copy, row_map, 
-                                       int(n_entries_per_row), false))),
+                 graph (row_map.Comm().NumProc() > 1 ?
+                        (std::auto_ptr<Epetra_FECrsGraph>
+                                  (new Epetra_FECrsGraph(Copy, row_map, 
+                                                         (int)n_entries_per_row, 
+                                                         false)))
+                        :
+                        (std::auto_ptr<Epetra_FECrsGraph>
+                                  (new Epetra_FECrsGraph(Copy, row_map, col_map, 
+                                                         (int)n_entries_per_row, 
+                                                         false)))
+                        ),
                  cached_row_indices (1),
                  n_cached_elements (0),
                  row_in_cache (0)
@@ -149,10 +184,19 @@ namespace TrilinosWrappers
                   row_map (InputRowMap),
                   col_map (InputColMap),
                  compressed (false),
-                 graph (std::auto_ptr<Epetra_FECrsGraph>
-                   (new Epetra_FECrsGraph(Copy, row_map, 
-                     (int*)const_cast<unsigned int*>(&(n_entries_per_row[0])),
-                                           false))),
+                 graph (row_map.Comm().NumProc() > 1 ?
+                        (std::auto_ptr<Epetra_FECrsGraph>
+                         (new Epetra_FECrsGraph(Copy, row_map, 
+                                                (int*)const_cast<unsigned int*>
+                                                (&(n_entries_per_row[row_map.MinMyGID()])),
+                                                false)))
+                        :
+                        (std::auto_ptr<Epetra_FECrsGraph>
+                         (new Epetra_FECrsGraph(Copy, row_map, col_map, 
+                                                (int*)const_cast<unsigned int*>
+                                                (&(n_entries_per_row[row_map.MinMyGID()])),
+                                                false)))
+                        ),
                  cached_row_indices (1),
                  n_cached_elements (0),
                  row_in_cache (0)
@@ -171,8 +215,8 @@ namespace TrilinosWrappers
 #endif
                  compressed (false),
                  graph (std::auto_ptr<Epetra_FECrsGraph>
-                               (new Epetra_FECrsGraph(Copy, row_map, 
-                                       int(n_entries_per_row), false))),
+                        (new Epetra_FECrsGraph(Copy, row_map, col_map,
+                                               int(n_entries_per_row), false))),
                  cached_row_indices (1),
                  n_cached_elements (0),
                  row_in_cache (0)
@@ -191,9 +235,9 @@ namespace TrilinosWrappers
 #endif
                  compressed (false),
                  graph (std::auto_ptr<Epetra_FECrsGraph>
-                    (new Epetra_FECrsGraph(Copy, row_map, 
-                       (int*)const_cast<unsigned int*>(&(n_entries_per_row[0])), 
-                                            false))),
+                        (new Epetra_FECrsGraph(Copy, row_map, col_map, 
+                         (int*)const_cast<unsigned int*>(&(n_entries_per_row[0])), 
+                                               false))),
                  cached_row_indices (1),
                  n_cached_elements (0),
                  row_in_cache (0)
@@ -201,8 +245,8 @@ namespace TrilinosWrappers
 
                                   // Copy function is currently not working
                                   // because the Trilinos Epetra_FECrsGraph
-                                  // does not implement a reinit function
-                                  // from another graph.
+                                  // does not implement a constructor from
+                                  // another graph.
   /*
   SparsityPattern::SparsityPattern (const SparsityPattern &InputSP)
                  :
@@ -261,8 +305,22 @@ namespace TrilinosWrappers
     row_map = input_row_map;
     col_map = input_col_map;
 
-    graph = std::auto_ptr<Epetra_FECrsGraph>
-      (new Epetra_FECrsGraph(Copy, row_map, n_entries_per_row, false));
+                                  // for more than one processor, need to
+                                  // specify only row map first and let the
+                                  // matrix entries decide about the column
+                                  // map (which says which columns are
+                                  // present in the matrix, not to be
+                                  // confused with the col_map that tells
+                                  // how the domain dofs of the matrix will
+                                  // be distributed). for only one
+                                  // processor, we can directly assign the
+                                  // columns as well.
+    if (row_map.Comm().NumProc() > 1)
+      graph = std::auto_ptr<Epetra_FECrsGraph>
+       (new Epetra_FECrsGraph(Copy, row_map, n_entries_per_row, false));
+    else
+      graph = std::auto_ptr<Epetra_FECrsGraph>
+       (new Epetra_FECrsGraph(Copy, row_map, col_map, n_entries_per_row, false));
 
     n_cached_elements = 0;
     row_in_cache = 0;
@@ -312,10 +370,16 @@ namespace TrilinosWrappers
     row_map = input_row_map;
     col_map = input_col_map;
 
-    graph = std::auto_ptr<Epetra_FECrsGraph>
-      (new Epetra_FECrsGraph(Copy, row_map, 
-                            n_entries_per_row[input_row_map.MinMyGID()], 
-                            false));
+    if (row_map.Comm().NumProc() > 1)
+      graph = std::auto_ptr<Epetra_FECrsGraph>
+       (new Epetra_FECrsGraph(Copy, row_map, 
+                              n_entries_per_row[input_row_map.MinMyGID()], 
+                              false));
+    else
+      graph = std::auto_ptr<Epetra_FECrsGraph>
+       (new Epetra_FECrsGraph(Copy, row_map, col_map,
+                              n_entries_per_row[input_row_map.MinMyGID()], 
+                              false));      
 
     n_cached_elements = 0;
     row_in_cache = 0;
@@ -360,10 +424,16 @@ namespace TrilinosWrappers
     for (unsigned int row=0; row<n_rows; ++row)
       n_entries_per_row[row] = sp.row_length(row);
 
-    graph = std::auto_ptr<Epetra_FECrsGraph>
-      (new Epetra_FECrsGraph(Copy, row_map, 
-                            n_entries_per_row[input_row_map.MinMyGID()], 
-                            false));
+    if (row_map.Comm().NumProc() > 1)
+      graph = std::auto_ptr<Epetra_FECrsGraph>
+       (new Epetra_FECrsGraph(Copy, row_map, 
+                              n_entries_per_row[input_row_map.MinMyGID()], 
+                              false));
+    else
+      graph = std::auto_ptr<Epetra_FECrsGraph>
+       (new Epetra_FECrsGraph(Copy, row_map, col_map,
+                              n_entries_per_row[input_row_map.MinMyGID()], 
+                              false));      
 
     Assert (graph->NumGlobalRows() == (int)sp.n_rows(),
            ExcDimensionMismatch (graph->NumGlobalRows(),
@@ -417,10 +487,16 @@ namespace TrilinosWrappers
     for (unsigned int row=0; row<n_rows; ++row)
       n_entries_per_row[row] = sp.row_length(row);
 
-    graph = std::auto_ptr<Epetra_FECrsGraph>
-      (new Epetra_FECrsGraph(Copy, row_map, 
-                            n_entries_per_row[input_row_map.MinMyGID()], 
-                            false));
+    if (row_map.Comm().NumProc() > 1)
+      graph = std::auto_ptr<Epetra_FECrsGraph>
+       (new Epetra_FECrsGraph(Copy, row_map, 
+                              n_entries_per_row[input_row_map.MinMyGID()], 
+                              false));
+    else
+      graph = std::auto_ptr<Epetra_FECrsGraph>
+       (new Epetra_FECrsGraph(Copy, row_map, col_map,
+                              n_entries_per_row[input_row_map.MinMyGID()], 
+                              false));      
 
     Assert (graph->NumGlobalRows() == (int)sp.n_rows(),
            ExcDimensionMismatch (graph->NumGlobalRows(),