]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Avoid using deprecated interfaces
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Wed, 29 Jul 2015 16:41:27 +0000 (18:41 +0200)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Thu, 30 Jul 2015 11:54:44 +0000 (13:54 +0200)
include/deal.II/lac/trilinos_sparse_matrix.h
source/lac/trilinos_sparse_matrix.cc
source/lac/trilinos_sparse_matrix.inst.in

index b9157880b04ed08f001420652333c12735b6f0aa..1096be7f0843fb69c744d16eb5e5c248c9c9ba06 100644 (file)
@@ -2540,37 +2540,12 @@ namespace TrilinosWrappers
                              const MPI_Comm      &communicator,
                              const bool           exchange_data)
   {
-    Epetra_Map map = parallel_partitioning.make_trilinos_map (communicator, false);
-    reinit (map, map, sparsity_pattern, exchange_data);
+    reinit (parallel_partitioning, parallel_partitioning,
+            sparsity_pattern, communicator, exchange_data);
   }
 
 
 
-  template <typename SparsityType>
-  inline
-  void SparseMatrix::reinit (const IndexSet      &row_parallel_partitioning,
-                             const IndexSet      &col_parallel_partitioning,
-                             const SparsityType  &sparsity_pattern,
-                             const MPI_Comm      &communicator,
-                             const bool           exchange_data)
-  {
-    Epetra_Map row_map =
-      row_parallel_partitioning.make_trilinos_map (communicator, false);
-    Epetra_Map col_map =
-      col_parallel_partitioning.make_trilinos_map (communicator, false);
-    reinit (row_map, col_map, sparsity_pattern, exchange_data);
-  }
-
-
-  // declare the existence of an explicit specialization
-  template <>
-  void
-  SparseMatrix::reinit (const Epetra_Map    &input_row_map,
-                        const Epetra_Map    &input_col_map,
-                        const DynamicSparsityPattern &sparsity_pattern,
-                        const bool           exchange_data);
-
-
   template <typename number>
   inline
   void SparseMatrix::reinit (const IndexSet      &parallel_partitioning,
@@ -2581,28 +2556,8 @@ namespace TrilinosWrappers
                              const ::dealii::SparsityPattern *use_this_sparsity)
   {
     Epetra_Map map = parallel_partitioning.make_trilinos_map (communicator, false);
-    reinit (map, map, sparse_matrix, drop_tolerance, copy_values,
-            use_this_sparsity);
-  }
-
-
-
-  template <typename number>
-  inline
-  void SparseMatrix::reinit (const IndexSet      &row_parallel_partitioning,
-                             const IndexSet      &col_parallel_partitioning,
-                             const ::dealii::SparseMatrix<number> &sparse_matrix,
-                             const MPI_Comm      &communicator,
-                             const double         drop_tolerance,
-                             const bool           copy_values,
-                             const ::dealii::SparsityPattern *use_this_sparsity)
-  {
-    Epetra_Map row_map =
-      row_parallel_partitioning.make_trilinos_map (communicator, false);
-    Epetra_Map col_map =
-      col_parallel_partitioning.make_trilinos_map (communicator, false);
-    reinit (row_map, col_map, sparse_matrix, drop_tolerance, copy_values,
-            use_this_sparsity);
+    reinit (parallel_partitioning, parallel_partitioning, sparse_matrix,
+            drop_tolerance, copy_values, use_this_sparsity);
   }
 
 
@@ -2625,24 +2580,6 @@ namespace TrilinosWrappers
 
 
 
-  inline
-  const Epetra_Map &
-  SparseMatrix::domain_partitioner () const
-  {
-    return matrix->DomainMap();
-  }
-
-
-
-  inline
-  const Epetra_Map &
-  SparseMatrix::range_partitioner () const
-  {
-    return matrix->RangeMap();
-  }
-
-
-
   inline
   IndexSet
   SparseMatrix::locally_owned_domain_indices () const
@@ -2661,24 +2598,6 @@ namespace TrilinosWrappers
 
 
 
-  inline
-  const Epetra_Map &
-  SparseMatrix::row_partitioner () const
-  {
-    return matrix->RowMap();
-  }
-
-
-
-  inline
-  const Epetra_Map &
-  SparseMatrix::col_partitioner () const
-  {
-    return matrix->ColMap();
-  }
-
-
-
   inline
   void
   SparseMatrix::prepare_add()
index e21bfeb8a2a3bfd56aee61060713083aed4f9205..f48945329ea3c4ceab66d0af8d215dcbddb1dfd1 100644 (file)
@@ -450,306 +450,360 @@ namespace TrilinosWrappers
 
 
 
-  template <typename SparsityType>
-  void
-  SparseMatrix::reinit (const SparsityType &sparsity_pattern)
+  namespace
   {
-    const Epetra_Map rows (static_cast<TrilinosWrappers::types::int_type>(sparsity_pattern.n_rows()),
-                           0,
-                           Utilities::Trilinos::comm_self());
-    const Epetra_Map columns (static_cast<TrilinosWrappers::types::int_type>(sparsity_pattern.n_cols()),
-                              0,
-                              Utilities::Trilinos::comm_self());
+    typedef SparseMatrix::size_type size_type;
 
-    reinit (rows, columns, sparsity_pattern);
-  }
+    template <typename SparsityType>
+    void
+    reinit_matrix (const Epetra_Map                          &input_row_map,
+                   const Epetra_Map                          &input_col_map,
+                   const SparsityType                        &sparsity_pattern,
+                   const bool                                 exchange_data,
+                   std_cxx11::shared_ptr<Epetra_Map>         &column_space_map,
+                   std_cxx11::shared_ptr<Epetra_FECrsMatrix> &matrix,
+                   std_cxx11::shared_ptr<Epetra_CrsMatrix>   &nonlocal_matrix,
+                   std_cxx11::shared_ptr<Epetra_Export>      &nonlocal_matrix_exporter)
+    {
+      // release memory before reallocation
+      matrix.reset();
+      nonlocal_matrix.reset();
+      nonlocal_matrix_exporter.reset();
 
+      if (input_row_map.Comm().MyPID() == 0)
+        {
+          AssertDimension (sparsity_pattern.n_rows(),
+                           static_cast<size_type>(n_global_elements(input_row_map)));
+          AssertDimension (sparsity_pattern.n_cols(),
+                           static_cast<size_type>(n_global_elements(input_col_map)));
+        }
 
+      column_space_map.reset (new Epetra_Map (input_col_map));
 
-  template <typename SparsityType>
-  void
-  SparseMatrix::reinit (const Epetra_Map    &input_map,
-                        const SparsityType  &sparsity_pattern,
-                        const bool           exchange_data)
-  {
-    reinit (input_map, input_map, sparsity_pattern, exchange_data);
-  }
+      // if we want to exchange data, build a usual Trilinos sparsity pattern
+      // and let that handle the exchange. otherwise, manually create a
+      // CrsGraph, which consumes considerably less memory because it can set
+      // correct number of indices right from the start
+      if (exchange_data)
+        {
+          SparsityPattern trilinos_sparsity;
+          trilinos_sparsity.reinit (input_row_map, input_col_map,
+                                    sparsity_pattern, exchange_data);
+          matrix.reset (new Epetra_FECrsMatrix
+                        (Copy, trilinos_sparsity.trilinos_sparsity_pattern(), false));
 
+          return;
+        }
 
+      const size_type first_row = min_my_gid(input_row_map),
+                      last_row = max_my_gid(input_row_map)+1;
+      std::vector<int> n_entries_per_row(last_row-first_row);
+
+      for (size_type row=first_row; row<last_row; ++row)
+        n_entries_per_row[row-first_row] = sparsity_pattern.row_length(row);
+
+      // 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. Moreover, it consumes less memory, since the
+      // internal reordering is done on ints only, and we can leave the
+      // doubles aside.
+
+      // 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. Compare this with bug # 4123 in the Sandia Bugzilla.
+      std_cxx11::shared_ptr<Epetra_CrsGraph> graph;
+      if (input_row_map.Comm().NumProc() > 1)
+        graph.reset (new Epetra_CrsGraph (Copy, input_row_map,
+                                          &n_entries_per_row[0], true));
+      else
+        graph.reset (new Epetra_CrsGraph (Copy, input_row_map, input_col_map,
+                                          &n_entries_per_row[0], true));
 
-  template <typename SparsityType>
-  void
-  SparseMatrix::reinit (const Epetra_Map    &input_row_map,
-                        const Epetra_Map    &input_col_map,
-                        const SparsityType  &sparsity_pattern,
-                        const bool           exchange_data)
-  {
-    // release memory before reallocation
-    matrix.reset();
-    nonlocal_matrix.reset();
-    nonlocal_matrix_exporter.reset();
+      // This functions assumes that the sparsity pattern sits on all
+      // processors (completely). The parallel version uses an Epetra graph
+      // that is already distributed.
 
-    // if we want to exchange data, build a usual Trilinos sparsity pattern
-    // and let that handle the exchange. otherwise, manually create a
-    // CrsGraph, which consumes considerably less memory because it can set
-    // correct number of indices right from the start
-    if (exchange_data)
-      {
-        SparsityPattern trilinos_sparsity;
-        trilinos_sparsity.reinit (input_row_map, input_col_map,
-                                  sparsity_pattern, exchange_data);
-        reinit (trilinos_sparsity);
+      // now insert the indices
+      std::vector<TrilinosWrappers::types::int_type>   row_indices;
 
-        return;
-      }
+      for (size_type row=first_row; row<last_row; ++row)
+        {
+          const int row_length = sparsity_pattern.row_length(row);
+          if (row_length == 0)
+            continue;
 
-    Assert (exchange_data == false, ExcNotImplemented());
-    if (input_row_map.Comm().MyPID() == 0)
-      {
-        AssertDimension (sparsity_pattern.n_rows(),
-                         static_cast<size_type>(n_global_elements(input_row_map)));
-        AssertDimension (sparsity_pattern.n_cols(),
-                         static_cast<size_type>(n_global_elements(input_col_map)));
-      }
+          row_indices.resize (row_length, -1);
+          {
+            typename SparsityType::iterator p = sparsity_pattern.begin(row);
+            for (size_type col=0; p != sparsity_pattern.end(row); ++p, ++col)
+              row_indices[col] = p->column();
+          }
+          graph->Epetra_CrsGraph::InsertGlobalIndices (row, row_length,
+                                                       &row_indices[0]);
+        }
 
-    column_space_map.reset (new Epetra_Map (input_col_map));
-
-    const size_type first_row = min_my_gid(input_row_map),
-                    last_row = max_my_gid(input_row_map)+1;
-    std::vector<int> n_entries_per_row(last_row-first_row);
-
-    for (size_type row=first_row; row<last_row; ++row)
-      n_entries_per_row[row-first_row] = sparsity_pattern.row_length(row);
-
-    // 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. Moreover, it consumes less memory, since the internal
-    // reordering is done on ints only, and we can leave the doubles aside.
-
-    // 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. Compare
-    // this with bug # 4123 in the Sandia Bugzilla.
-    std_cxx11::shared_ptr<Epetra_CrsGraph> graph;
-    if (input_row_map.Comm().NumProc() > 1)
-      graph.reset (new Epetra_CrsGraph (Copy, input_row_map,
-                                        &n_entries_per_row[0], true));
-    else
-      graph.reset (new Epetra_CrsGraph (Copy, input_row_map, input_col_map,
-                                        &n_entries_per_row[0], true));
+      // Eventually, optimize the graph structure (sort indices, make memory
+      // contiguous, etc). note that the documentation of the function indeed
+      // states that we first need to provide the column (domain) map and then
+      // the row (range) map
+      graph->FillComplete(input_col_map, input_row_map);
+      graph->OptimizeStorage();
 
-    // This functions assumes that the sparsity pattern sits on all processors
-    // (completely). The parallel version uses an Epetra graph that is already
-    // distributed.
+      // check whether we got the number of columns right.
+      AssertDimension (sparsity_pattern.n_cols(),
+                       static_cast<size_type>(n_global_cols(*graph)));
+      (void)n_global_cols;
 
-    // now insert the indices
-    std::vector<TrilinosWrappers::types::int_type>   row_indices;
+      // And now finally generate the matrix.
+      matrix.reset (new Epetra_FECrsMatrix(Copy, *graph, false));
+    }
 
-    for (size_type row=first_row; row<last_row; ++row)
-      {
-        const int row_length = sparsity_pattern.row_length(row);
-        if (row_length == 0)
-          continue;
 
-        row_indices.resize (row_length, -1);
+
+    // specialization for DynamicSparsityPattern which can provide us with
+    // more information about the non-locally owned rows
+    template <>
+    void
+    reinit_matrix (const Epetra_Map                          &input_row_map,
+                   const Epetra_Map                          &input_col_map,
+                   const DynamicSparsityPattern              &sparsity_pattern,
+                   const bool                                 exchange_data,
+                   std_cxx11::shared_ptr<Epetra_Map>         &column_space_map,
+                   std_cxx11::shared_ptr<Epetra_FECrsMatrix> &matrix,
+                   std_cxx11::shared_ptr<Epetra_CrsMatrix>   &nonlocal_matrix,
+                   std_cxx11::shared_ptr<Epetra_Export>      &nonlocal_matrix_exporter)
+    {
+      matrix.reset();
+      nonlocal_matrix.reset();
+      nonlocal_matrix_exporter.reset();
+
+      AssertDimension (sparsity_pattern.n_rows(),
+                       static_cast<size_type>(n_global_elements(input_row_map)));
+      AssertDimension (sparsity_pattern.n_cols(),
+                       static_cast<size_type>(n_global_elements(input_col_map)));
+
+      column_space_map.reset (new Epetra_Map (input_col_map));
+
+      IndexSet relevant_rows (sparsity_pattern.row_index_set());
+      // serial case
+      if (relevant_rows.size() == 0)
         {
-          typename SparsityType::iterator p = sparsity_pattern.begin(row);
-          for (size_type col=0; p != sparsity_pattern.end(row); ++p, ++col)
-            row_indices[col] = p->column();
+          relevant_rows.set_size(n_global_elements(input_row_map));
+          relevant_rows.add_range(0, n_global_elements(input_row_map));
         }
-        graph->Epetra_CrsGraph::InsertGlobalIndices (row, row_length,
-                                                     &row_indices[0]);
+      relevant_rows.compress();
+      Assert(relevant_rows.n_elements() >= static_cast<unsigned int>(input_row_map.NumMyElements()),
+             ExcMessage("Locally relevant rows of sparsity pattern must contain "
+                        "all locally owned rows"));
+
+      // check whether the relevant rows correspond to exactly the same map as
+      // the owned rows. In that case, do not create the nonlocal graph and
+      // fill the columns by demand
+      bool have_ghost_rows = false;
+      {
+        std::vector<dealii::types::global_dof_index> indices;
+        relevant_rows.fill_index_vector(indices);
+        Epetra_Map relevant_map (TrilinosWrappers::types::int_type(-1),
+                                 TrilinosWrappers::types::int_type(relevant_rows.n_elements()),
+                                 (indices.empty() ? 0 :
+                                  reinterpret_cast<TrilinosWrappers::types::int_type *>(&indices[0])),
+                                 0, input_row_map.Comm());
+        if (relevant_map.SameAs(input_row_map))
+          have_ghost_rows = false;
+        else
+          have_ghost_rows = true;
       }
 
-    // Eventually, optimize the graph structure (sort indices, make memory
-    // contiguous, etc). note that the documentation of the function indeed
-    // states that we first need to provide the column (domain) map and then
-    // the row (range) map
-    graph->FillComplete(input_col_map, input_row_map);
-    graph->OptimizeStorage();
+      const unsigned int n_rows = relevant_rows.n_elements();
+      std::vector<TrilinosWrappers::types::int_type> ghost_rows;
+      std::vector<int> n_entries_per_row(input_row_map.NumMyElements());
+      std::vector<int> n_entries_per_ghost_row;
+      for (unsigned int i=0, own=0; i<n_rows; ++i)
+        {
+          const TrilinosWrappers::types::int_type global_row =
+            relevant_rows.nth_index_in_set(i);
+          if (input_row_map.MyGID(global_row))
+            n_entries_per_row[own++] = sparsity_pattern.row_length(global_row);
+          else if (sparsity_pattern.row_length(global_row) > 0)
+            {
+              ghost_rows.push_back(global_row);
+              n_entries_per_ghost_row.push_back(sparsity_pattern.row_length(global_row));
+            }
+        }
 
-    // check whether we got the number of columns right.
-    AssertDimension (sparsity_pattern.n_cols(),static_cast<size_type>(
-                       n_global_cols(*graph)));
-    (void)n_global_cols;
+      // make sure all processors create an off-processor matrix with at least
+      // one entry
+      if (have_ghost_rows == true && ghost_rows.empty() == true)
+        {
+          ghost_rows.push_back(0);
+          n_entries_per_ghost_row.push_back(1);
+        }
 
-    // And now finally generate the matrix.
-    matrix.reset (new Epetra_FECrsMatrix(Copy, *graph, false));
-    last_action = Zero;
+      Epetra_Map off_processor_map(-1, ghost_rows.size(),
+                                   (ghost_rows.size()>0)?(&ghost_rows[0]):NULL,
+                                   0, input_row_map.Comm());
 
-    // In the end, the matrix needs to be compressed in order to be really
-    // ready.
-    compress(VectorOperation::insert);
-  }
+      std_cxx11::shared_ptr<Epetra_CrsGraph> graph, nonlocal_graph;
+      if (input_row_map.Comm().NumProc() > 1)
+        {
+          graph.reset (new Epetra_CrsGraph (Copy, input_row_map,
+                                            (n_entries_per_row.size()>0)?(&n_entries_per_row[0]):NULL,
+                                            exchange_data ? false : true));
+          if (have_ghost_rows == true)
+            nonlocal_graph.reset (new Epetra_CrsGraph (Copy, off_processor_map,
+                                                       &n_entries_per_ghost_row[0],
+                                                       true));
+        }
+      else
+        graph.reset (new Epetra_CrsGraph (Copy, input_row_map, input_col_map,
+                                          (n_entries_per_row.size()>0)?(&n_entries_per_row[0]):NULL,
+                                          true));
 
+      // now insert the indices, select between the right matrix
+      std::vector<TrilinosWrappers::types::int_type> row_indices;
 
+      for (unsigned int i=0; i<n_rows; ++i)
+        {
+          const TrilinosWrappers::types::int_type global_row =
+            relevant_rows.nth_index_in_set(i);
+          const int row_length = sparsity_pattern.row_length(global_row);
+          if (row_length == 0)
+            continue;
 
-  // specialization for CompressedSimpleSparsityPattern which can provide us
-  // with more information about the non-locally owned rows
-  template <>
-  void
-  SparseMatrix::reinit (const Epetra_Map    &input_row_map,
-                        const Epetra_Map    &input_col_map,
-                        const DynamicSparsityPattern &sparsity_pattern,
-                        const bool           exchange_data)
-  {
-    matrix.reset();
-    nonlocal_matrix.reset();
-    nonlocal_matrix_exporter.reset();
+          row_indices.resize (row_length, -1);
+          {
+            dealii::DynamicSparsityPattern::iterator p = sparsity_pattern.begin(global_row);
+            for (size_type col=0; p != sparsity_pattern.end(global_row); ++p, ++col)
+              row_indices[col] = p->column();
+          }
 
-    AssertDimension (sparsity_pattern.n_rows(),
-                     static_cast<size_type>(n_global_elements(input_row_map)));
-    AssertDimension (sparsity_pattern.n_cols(),
-                     static_cast<size_type>(n_global_elements(input_col_map)));
+          if (input_row_map.MyGID(global_row))
+            graph->InsertGlobalIndices (global_row, row_length, &row_indices[0]);
+          else
+            {
+              Assert(nonlocal_graph.get() != 0, ExcInternalError());
+              nonlocal_graph->InsertGlobalIndices (global_row, row_length,
+                                                   &row_indices[0]);
+            }
+        }
 
-    column_space_map.reset (new Epetra_Map (input_col_map));
+      // finalize nonlocal graph and create nonlocal matrix
+      if (nonlocal_graph.get() != 0)
+        {
+          if (nonlocal_graph->IndicesAreGlobal() == false &&
+              nonlocal_graph->RowMap().NumMyElements() > 0)
+            {
+              // insert dummy element
+              TrilinosWrappers::types::int_type row =
+                nonlocal_graph->RowMap().MyGID(TrilinosWrappers::types::int_type(0));
+              nonlocal_graph->InsertGlobalIndices(row, 1, &row);
+            }
+          Assert(nonlocal_graph->IndicesAreGlobal() == true,
+                 ExcInternalError());
+          nonlocal_graph->FillComplete(input_col_map, input_row_map);
+          nonlocal_graph->OptimizeStorage();
 
-    IndexSet relevant_rows (sparsity_pattern.row_index_set());
-    // serial case
-    if (relevant_rows.size() == 0)
-      {
-        relevant_rows.set_size(n_global_elements(input_row_map));
-        relevant_rows.add_range(0, n_global_elements(input_row_map));
-      }
-    relevant_rows.compress();
-    Assert(relevant_rows.n_elements() >= static_cast<unsigned int>(input_row_map.NumMyElements()),
-           ExcMessage("Locally relevant rows of sparsity pattern must contain "
-                      "all locally owned rows"));
-
-    // check whether the relevant rows correspond to exactly the same map as
-    // the owned rows. In that case, do not create the nonlocal graph and fill
-    // the columns by demand
-    bool have_ghost_rows = false;
-    {
-      std::vector<dealii::types::global_dof_index> indices;
-      relevant_rows.fill_index_vector(indices);
-      Epetra_Map relevant_map (TrilinosWrappers::types::int_type(-1),
-                               TrilinosWrappers::types::int_type(relevant_rows.n_elements()),
-                               (indices.empty() ? 0 :
-                                reinterpret_cast<TrilinosWrappers::types::int_type *>(&indices[0])),
-                               0, input_row_map.Comm());
-      if (relevant_map.SameAs(input_row_map))
-        have_ghost_rows = false;
-      else
-        have_ghost_rows = true;
+          // insert data from nonlocal graph into the final sparsity pattern
+          if (exchange_data)
+            {
+              Epetra_Export exporter(nonlocal_graph->RowMap(), input_row_map);
+              int ierr = graph->Export(*nonlocal_graph, exporter, Add);
+              (void)ierr;
+              Assert (ierr==0, ExcTrilinosError(ierr));
+            }
+
+          nonlocal_matrix.reset (new Epetra_CrsMatrix(Copy, *nonlocal_graph));
+        }
+
+      graph->FillComplete(input_col_map, input_row_map);
+      graph->OptimizeStorage();
+
+      AssertDimension (sparsity_pattern.n_cols(),static_cast<size_type>(
+                         n_global_cols(*graph)));
+
+      matrix.reset (new Epetra_FECrsMatrix(Copy, *graph, false));
     }
+  }
 
-    const unsigned int n_rows = relevant_rows.n_elements();
-    std::vector<TrilinosWrappers::types::int_type> ghost_rows;
-    std::vector<int> n_entries_per_row(input_row_map.NumMyElements());
-    std::vector<int> n_entries_per_ghost_row;
-    for (unsigned int i=0, own=0; i<n_rows; ++i)
-      {
-        const TrilinosWrappers::types::int_type global_row =
-          relevant_rows.nth_index_in_set(i);
-        if (input_row_map.MyGID(global_row))
-          n_entries_per_row[own++] = sparsity_pattern.row_length(global_row);
-        else if (sparsity_pattern.row_length(global_row) > 0)
-          {
-            ghost_rows.push_back(global_row);
-            n_entries_per_ghost_row.push_back(sparsity_pattern.row_length(global_row));
-          }
-      }
 
-    // make sure all processors create an off-processor matrix with at least
-    // one entry
-    if (have_ghost_rows == true && ghost_rows.empty() == true)
-      {
-        ghost_rows.push_back(0);
-        n_entries_per_ghost_row.push_back(1);
-      }
 
-    Epetra_Map off_processor_map(-1, ghost_rows.size(),
-                                 (ghost_rows.size()>0)?(&ghost_rows[0]):NULL,
-                                 0, input_row_map.Comm());
+  template <typename SparsityType>
+  void
+  SparseMatrix::reinit (const SparsityType &sparsity_pattern)
+  {
+    const Epetra_Map rows (static_cast<TrilinosWrappers::types::int_type>(sparsity_pattern.n_rows()),
+                           0,
+                           Utilities::Trilinos::comm_self());
+    const Epetra_Map columns (static_cast<TrilinosWrappers::types::int_type>(sparsity_pattern.n_cols()),
+                              0,
+                              Utilities::Trilinos::comm_self());
 
-    std_cxx11::shared_ptr<Epetra_CrsGraph> graph, nonlocal_graph;
-    if (input_row_map.Comm().NumProc() > 1)
-      {
-        graph.reset (new Epetra_CrsGraph (Copy, input_row_map,
-                                          (n_entries_per_row.size()>0)?(&n_entries_per_row[0]):NULL,
-                                          exchange_data ? false : true));
-        if (have_ghost_rows == true)
-          nonlocal_graph.reset (new Epetra_CrsGraph (Copy, off_processor_map,
-                                                     &n_entries_per_ghost_row[0],
-                                                     true));
-      }
-    else
-      graph.reset (new Epetra_CrsGraph (Copy, input_row_map, input_col_map,
-                                        (n_entries_per_row.size()>0)?(&n_entries_per_row[0]):NULL,
-                                        true));
+    reinit_matrix (rows, columns, sparsity_pattern, false,
+                   column_space_map, matrix, nonlocal_matrix,
+                   nonlocal_matrix_exporter);
+  }
 
-    // now insert the indices, select between the right matrix
-    std::vector<TrilinosWrappers::types::int_type> row_indices;
 
-    for (unsigned int i=0; i<n_rows; ++i)
-      {
-        const TrilinosWrappers::types::int_type global_row =
-          relevant_rows.nth_index_in_set(i);
-        const int row_length = sparsity_pattern.row_length(global_row);
-        if (row_length == 0)
-          continue;
 
-        row_indices.resize (row_length, -1);
-        {
-          dealii::DynamicSparsityPattern::iterator p = sparsity_pattern.begin(global_row);
-          for (size_type col=0; p != sparsity_pattern.end(global_row); ++p, ++col)
-            row_indices[col] = p->column();
-        }
+  template <typename SparsityType>
+  void
+  SparseMatrix::reinit (const Epetra_Map    &input_map,
+                        const SparsityType  &sparsity_pattern,
+                        const bool           exchange_data)
+  {
+    reinit_matrix (input_map, input_map, sparsity_pattern, exchange_data,
+                   column_space_map, matrix, nonlocal_matrix,
+                   nonlocal_matrix_exporter);
+  }
 
-        if (input_row_map.MyGID(global_row))
-          graph->InsertGlobalIndices (global_row, row_length, &row_indices[0]);
-        else
-          {
-            Assert(nonlocal_graph.get() != 0, ExcInternalError());
-            nonlocal_graph->InsertGlobalIndices (global_row, row_length,
-                                                 &row_indices[0]);
-          }
-      }
 
-    // finalize nonlocal graph and create nonlocal matrix
-    if (nonlocal_graph.get() != 0)
-      {
-        if (nonlocal_graph->IndicesAreGlobal() == false &&
-            nonlocal_graph->RowMap().NumMyElements() > 0)
-          {
-            // insert dummy element
-            TrilinosWrappers::types::int_type row =
-              nonlocal_graph->RowMap().MyGID(TrilinosWrappers::types::int_type(0));
-            nonlocal_graph->InsertGlobalIndices(row, 1, &row);
-          }
-        Assert(nonlocal_graph->IndicesAreGlobal() == true,
-               ExcInternalError());
-        nonlocal_graph->FillComplete(input_col_map, input_row_map);
-        nonlocal_graph->OptimizeStorage();
 
-        // insert data from nonlocal graph into the final sparsity pattern
-        if (exchange_data)
-          {
-            Epetra_Export exporter(nonlocal_graph->RowMap(), input_row_map);
-            int ierr = graph->Export(*nonlocal_graph, exporter, Add);
-            (void)ierr;
-            Assert (ierr==0, ExcTrilinosError(ierr));
-          }
 
-        nonlocal_matrix.reset (new Epetra_CrsMatrix(Copy, *nonlocal_graph));
-      }
 
-    graph->FillComplete(input_col_map, input_row_map);
-    graph->OptimizeStorage();
 
-    AssertDimension (sparsity_pattern.n_cols(),static_cast<size_type>(
-                       n_global_cols(*graph)));
+  template <typename SparsityType>
+  inline
+  void SparseMatrix::reinit (const IndexSet      &row_parallel_partitioning,
+                             const IndexSet      &col_parallel_partitioning,
+                             const SparsityType  &sparsity_pattern,
+                             const MPI_Comm      &communicator,
+                             const bool           exchange_data)
+  {
+    Epetra_Map row_map =
+      row_parallel_partitioning.make_trilinos_map (communicator, false);
+    Epetra_Map col_map =
+      col_parallel_partitioning.make_trilinos_map (communicator, false);
+    reinit_matrix (row_map, col_map, sparsity_pattern, exchange_data,
+                   column_space_map, matrix, nonlocal_matrix,
+                   nonlocal_matrix_exporter);
 
-    matrix.reset (new Epetra_FECrsMatrix(Copy, *graph, false));
+    // In the end, the matrix needs to be compressed in order to be really
+    // ready.
     last_action = Zero;
+    compress(VectorOperation::insert);
+  }
+
+
+
+  template <typename SparsityType>
+  inline
+  void SparseMatrix::reinit (const Epetra_Map    &row_map,
+                             const Epetra_Map    &col_map,
+                             const SparsityType  &sparsity_pattern,
+                             const bool           exchange_data)
+  {
+    reinit_matrix (row_map, col_map, sparsity_pattern, exchange_data,
+                   column_space_map, matrix, nonlocal_matrix,
+                   nonlocal_matrix_exporter);
 
     // In the end, the matrix needs to be compressed in order to be really
     // ready.
+    last_action = Zero;
     compress(VectorOperation::insert);
   }
 
@@ -804,71 +858,33 @@ namespace TrilinosWrappers
 
 
   template <typename number>
-  void
-  SparseMatrix::reinit (const ::dealii::SparseMatrix<number> &dealii_sparse_matrix,
-                        const double                          drop_tolerance,
-                        const bool                            copy_values,
-                        const ::dealii::SparsityPattern      *use_this_sparsity)
-  {
-    const Epetra_Map rows (static_cast<TrilinosWrappers::types::int_type>(dealii_sparse_matrix.m()),
-                           0,
-                           Utilities::Trilinos::comm_self());
-    const Epetra_Map columns (static_cast<TrilinosWrappers::types::int_type>(dealii_sparse_matrix.n()),
-                              0,
-                              Utilities::Trilinos::comm_self());
-    reinit (rows, columns, dealii_sparse_matrix, drop_tolerance,
-            copy_values, use_this_sparsity);
-  }
-
-
-
-  template <typename number>
-  void
-  SparseMatrix::reinit (const Epetra_Map                     &input_map,
-                        const ::dealii::SparseMatrix<number> &dealii_sparse_matrix,
-                        const double                          drop_tolerance,
-                        const bool                            copy_values,
-                        const ::dealii::SparsityPattern      *use_this_sparsity)
-  {
-    reinit (input_map, input_map, dealii_sparse_matrix, drop_tolerance,
-            copy_values, use_this_sparsity);
-  }
-
-
-
-  template <typename number>
-  void
-  SparseMatrix::reinit (const Epetra_Map                     &input_row_map,
-                        const Epetra_Map                     &input_col_map,
-                        const ::dealii::SparseMatrix<number> &dealii_sparse_matrix,
-                        const double                          drop_tolerance,
-                        const bool                            copy_values,
-                        const ::dealii::SparsityPattern      *use_this_sparsity)
+  inline
+  void SparseMatrix::reinit (const IndexSet      &row_parallel_partitioning,
+                             const IndexSet      &col_parallel_partitioning,
+                             const ::dealii::SparseMatrix<number> &dealii_sparse_matrix,
+                             const MPI_Comm      &communicator,
+                             const double         drop_tolerance,
+                             const bool           copy_values,
+                             const ::dealii::SparsityPattern *use_this_sparsity)
   {
     if (copy_values == false)
       {
         // in case we do not copy values, just
         // call the other function.
         if (use_this_sparsity == 0)
-          reinit (input_row_map, input_col_map,
-                  dealii_sparse_matrix.get_sparsity_pattern());
+          reinit (row_parallel_partitioning, col_parallel_partitioning,
+                  dealii_sparse_matrix.get_sparsity_pattern(),
+                  communicator, false);
         else
-          reinit (input_row_map, input_col_map,
-                  *use_this_sparsity);
+          reinit (row_parallel_partitioning, col_parallel_partitioning,
+                  *use_this_sparsity, communicator, false);
         return;
       }
 
     const size_type n_rows = dealii_sparse_matrix.m();
 
-    Assert (static_cast<size_type>(n_global_elements(input_row_map)) == n_rows,
-            ExcDimensionMismatch (n_global_elements(input_row_map),
-                                  n_rows));
-    Assert (n_global_elements(input_row_map) == (TrilinosWrappers::types::int_type)n_rows,
-            ExcDimensionMismatch (n_global_elements(input_row_map),
-                                  n_rows));
-    Assert (n_global_elements(input_col_map) == (TrilinosWrappers::types::int_type)dealii_sparse_matrix.n(),
-            ExcDimensionMismatch (n_global_elements(input_col_map),
-                                  dealii_sparse_matrix.n()));
+    AssertDimension (row_parallel_partitioning.size(), n_rows);
+    AssertDimension (col_parallel_partitioning.size(), dealii_sparse_matrix.n());
 
     const ::dealii::SparsityPattern &sparsity_pattern =
       (use_this_sparsity!=0)? *use_this_sparsity :
@@ -878,23 +894,22 @@ namespace TrilinosWrappers
         m() != n_rows ||
         n_nonzero_elements() != sparsity_pattern.n_nonzero_elements())
       {
-        SparsityPattern trilinos_sparsity;
-        trilinos_sparsity.reinit (input_row_map, input_col_map, sparsity_pattern);
-        reinit (trilinos_sparsity);
+        reinit (row_parallel_partitioning, col_parallel_partitioning,
+                sparsity_pattern, communicator, false);
       }
 
-    // fill the values. the same as above: go through all rows of the matrix,
-    // and then all columns. since the sparsity patterns of the input matrix
-    // and the specified sparsity pattern might be different, need to go
-    // through the row for both these sparsity structures simultaneously in
-    // order to really set the correct values.
+    // fill the values. the same as above: go through all rows of the
+    // matrix, and then all columns. since the sparsity patterns of the
+    // input matrix and the specified sparsity pattern might be different,
+    // need to go through the row for both these sparsity structures
+    // simultaneously in order to really set the correct values.
     size_type maximum_row_length = matrix->MaxNumEntries();
     std::vector<size_type> row_indices (maximum_row_length);
     std::vector<TrilinosScalar> values (maximum_row_length);
 
     for (size_type row=0; row<n_rows; ++row)
       // see if the row is locally stored on this processor
-      if (input_row_map.MyGID(static_cast<TrilinosWrappers::types::int_type>(row)) == true)
+      if (row_parallel_partitioning.is_element(row) == true)
         {
           ::dealii::SparsityPattern::iterator select_index =
             sparsity_pattern.begin(row);
@@ -937,12 +952,56 @@ namespace TrilinosWrappers
           set (row, col, reinterpret_cast<size_type *>(&row_indices[0]),
                &values[0], false);
         }
-
     compress(VectorOperation::insert);
   }
 
 
 
+  template <typename number>
+  void
+  SparseMatrix::reinit (const ::dealii::SparseMatrix<number> &dealii_sparse_matrix,
+                        const double                          drop_tolerance,
+                        const bool                            copy_values,
+                        const ::dealii::SparsityPattern      *use_this_sparsity)
+  {
+    reinit (complete_index_set(dealii_sparse_matrix.m()),
+            complete_index_set(dealii_sparse_matrix.n()),
+            dealii_sparse_matrix, MPI_COMM_SELF, drop_tolerance,
+            copy_values, use_this_sparsity);
+  }
+
+
+
+  template <typename number>
+  void
+  SparseMatrix::reinit (const Epetra_Map                     &input_map,
+                        const ::dealii::SparseMatrix<number> &dealii_sparse_matrix,
+                        const double                          drop_tolerance,
+                        const bool                            copy_values,
+                        const ::dealii::SparsityPattern      *use_this_sparsity)
+  {
+    reinit (IndexSet(input_map), IndexSet(input_map), dealii_sparse_matrix,
+            MPI_COMM_SELF, drop_tolerance, copy_values, use_this_sparsity);
+  }
+
+
+
+  template <typename number>
+  void
+  SparseMatrix::reinit (const Epetra_Map                     &input_row_map,
+                        const Epetra_Map                     &input_col_map,
+                        const ::dealii::SparseMatrix<number> &dealii_sparse_matrix,
+                        const double                          drop_tolerance,
+                        const bool                            copy_values,
+                        const ::dealii::SparsityPattern      *use_this_sparsity)
+  {
+    reinit (IndexSet(input_row_map), IndexSet(input_col_map),
+            dealii_sparse_matrix, MPI_COMM_SELF,
+            drop_tolerance, copy_values, use_this_sparsity);
+  }
+
+
+
   void
   SparseMatrix::reinit (const Epetra_CrsMatrix &input_matrix,
                         const bool              copy_values)
@@ -972,6 +1031,7 @@ namespace TrilinosWrappers
                      my_nonzeros*sizeof (TrilinosScalar));
       }
 
+    last_action = Zero;
     compress(VectorOperation::insert);
   }
 
@@ -2313,13 +2373,47 @@ namespace TrilinosWrappers
             matrix->NumMyNonzeros() + sizeof(int)*local_size() + static_memory);
   }
 
+
+
+  const Epetra_Map &
+  SparseMatrix::domain_partitioner () const
+  {
+    return matrix->DomainMap();
+  }
+
+
+
+  const Epetra_Map &
+  SparseMatrix::range_partitioner () const
+  {
+    return matrix->RangeMap();
+  }
+
+
+
+  const Epetra_Map &
+  SparseMatrix::row_partitioner () const
+  {
+    return matrix->RowMap();
+  }
+
+
+
+  const Epetra_Map &
+  SparseMatrix::col_partitioner () const
+  {
+    return matrix->ColMap();
+  }
+
+
+
   MPI_Comm SparseMatrix::get_mpi_communicator () const
   {
 
 #ifdef DEAL_II_WITH_MPI
 
     const Epetra_MpiComm *mpi_comm
-      = dynamic_cast<const Epetra_MpiComm *>(&range_partitioner().Comm());
+      = dynamic_cast<const Epetra_MpiComm *>(&matrix->RangeMap().Comm());
     return mpi_comm->Comm();
 #else
 
@@ -2353,13 +2447,28 @@ namespace TrilinosWrappers
   SparseMatrix::reinit (const Epetra_Map &,
                         const DynamicSparsityPattern &,
                         const bool);
-
-
   template void
   SparseMatrix::reinit (const Epetra_Map &,
                         const Epetra_Map &,
                         const dealii::SparsityPattern &,
                         const bool);
+  template void
+  SparseMatrix::reinit (const Epetra_Map &,
+                        const Epetra_Map &,
+                        const DynamicSparsityPattern &,
+                        const bool);
+  template void
+  SparseMatrix::reinit (const IndexSet &,
+                        const IndexSet &,
+                        const dealii::SparsityPattern &,
+                        const MPI_Comm &,
+                        const bool);
+  template void
+  SparseMatrix::reinit (const IndexSet &,
+                        const IndexSet &,
+                        const DynamicSparsityPattern &,
+                        const MPI_Comm &,
+                        const bool);
 
   template void
   SparseMatrix::vmult (VectorBase &,
index 7e52ae22aaed6850e2fc0d4d273c39c2f940c316..8583f7047c140e79e816a402e9300ad25da9bec1 100644 (file)
@@ -37,5 +37,13 @@ for (S : REAL_SCALARS)
                             const double,
                             const bool,
                             const dealii::SparsityPattern *);
+      template void
+      SparseMatrix::reinit (const IndexSet &,
+                            const IndexSet &,
+                            const dealii::SparseMatrix<S> &,
+                            const MPI_Comm &,
+                            const double,
+                            const bool,
+                            const dealii::SparsityPattern *);
     \}
   }

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.