]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Trilinos sparse matrices can now also be reinitialized from the compressed sparsity...
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 22 Oct 2008 10:15:10 +0000 (10:15 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 22 Oct 2008 10:15:10 +0000 (10:15 +0000)
git-svn-id: https://svn.dealii.org/trunk@17303 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/trilinos_block_sparse_matrix.h
deal.II/lac/include/lac/trilinos_precondition_block.h
deal.II/lac/include/lac/trilinos_sparse_matrix.h
deal.II/lac/source/trilinos_block_sparse_matrix.cc
deal.II/lac/source/trilinos_sparse_matrix.cc

index 661f6b6d98b75e902ef96818dad7ecb43a149215..90133e1bdc68174051d80304cfa6bbdaa2e77fc9 100644 (file)
@@ -18,7 +18,6 @@
 #include <base/table.h>
 #include <lac/block_matrix_base.h>
 #include <lac/block_sparse_matrix.h>
-#include <lac/block_sparsity_pattern.h>
 #include <lac/trilinos_sparse_matrix.h>
 #include <lac/trilinos_block_vector.h>
 #include <lac/exceptions.h>
 
 DEAL_II_NAMESPACE_OPEN
 
+                                   // forward declarations
+class BlockSparsityPattern;
+class BlockCompressedSparsityPattern;
+class BlockCompressedSetSparsityPattern;
+class BlockCompressedSimpleSparsityPattern;
 
 
 namespace TrilinosWrappers
@@ -181,8 +185,9 @@ namespace TrilinosWrappers
                                        * assumes that a quadratic block
                                        * matrix is generated.
                                         */
+      template <typename BlockSparsityType>
       void reinit (const std::vector<Epetra_Map> &input_maps,
-                  const BlockSparsityPattern    &block_sparsity_pattern);
+                  const BlockSparsityType       &block_sparsity_pattern);
 
                                        /**
                                         * Resize the matrix and initialize it
@@ -191,8 +196,9 @@ namespace TrilinosWrappers
                                         * result is a block matrix for which
                                         * all elements are stored locally.
                                         */
-      void reinit (const BlockSparsityPattern &block_sparsity_pattern);
-      
+      template <typename BlockSparsityType>
+      void reinit (const BlockSparsityType &block_sparsity_pattern);
+
                                        /**
                                         * This function initializes the
                                        * Trilinos matrix using the deal.II
@@ -208,7 +214,25 @@ namespace TrilinosWrappers
                   const ::dealii::BlockSparseMatrix<double> &deal_ii_sparse_matrix,
                   const double                               drop_tolerance=1e-13);
 
-                                      /** 
+                                       /**
+                                        * This function initializes
+                                       * the Trilinos matrix using
+                                       * the deal.II sparse matrix
+                                       * and the entries stored
+                                       * therein. It uses a threshold
+                                       * to copy only elements whose
+                                       * modulus is larger than the
+                                       * threshold (so zeros in the
+                                       * deal.II matrix can be
+                                       * filtered away). Since no
+                                       * Epetra_Map is given, all the
+                                       * elements will be locally
+                                       * stored.
+                                        */
+      void reinit (const ::dealii::BlockSparseMatrix<double> &deal_ii_sparse_matrix,
+                  const double                               drop_tolerance=1e-13);
+
+                                      /**
                                        * This function calls the compress()
                                        * command of all matrices after 
                                        * the assembly is 
index ae4909388aef1a0bdc9c6adbea34bb143175308e..7615f9c5942aa490899380b10b7e1e2b4cceb51b 100644 (file)
@@ -286,8 +286,8 @@ namespace TrilinosWrappers
                                        * block, so that the number of
                                        * outer iterations is
                                        * reduced. See the discussion
-                                       * in the @step_22 step-22
-                                       * tutorial program.
+                                       * in the @ref step_22
+                                       * "step-22" tutorial program.
                                        */
        unsigned int Av_do_solve;
          
index 9ecf706bdfccd26a2bdfb8a90cc743196c85e032..2bf080ecf02ea87e0919a2ac75a7c481a10d8f2e 100755 (executable)
@@ -16,7 +16,6 @@
 
 #include <base/config.h>
 #include <base/subscriptor.h>
-#include <lac/sparsity_pattern.h>
 #include <lac/sparse_matrix.h>
 #include <lac/exceptions.h>
 #include <lac/trilinos_vector_base.h>
 
 DEAL_II_NAMESPACE_OPEN
 
+                                   // forward declarations
+class SparsityPattern;
+class CompressedSparsityPattern;
+class CompressedSetSparsityPattern;
+class CompressedSimpleSparsityPattern;
 
 namespace TrilinosWrappers
 {
@@ -466,28 +470,33 @@ namespace TrilinosWrappers
                                        * Epetra matrix know the
                                        * position of nonzero entries
                                        * according to the sparsity
-                                       * pattern. Note that, when
-                                       * using this function, the
-                                       * matrix must already be
-                                       * initialized with a suitable
-                                       * Epetra_Map that describes
-                                       * the distribution of the
-                                       * matrix among the MPI
-                                       * processes. Otherwise, an
-                                       * error will be thrown. In a
-                                       * parallel run, it is currently
-                                       * necessary that each processor
-                                       * holds the sparsity_pattern
-                                       * structure because each
-                                       * processor sets its rows.
+                                       * pattern. This function is
+                                       * meant for use in serial
+                                       * programs, where there is no
+                                       * need to specify how the
+                                       * matrix is going to be
+                                       * distributed among the
+                                       * processors. This function
+                                       * works in parallel, too, but
+                                       * it is recommended to
+                                       * manually specify the
+                                       * parallel partioning of the
+                                       * matrix using an
+                                       * Epetra_Map. When run in
+                                       * parallel, it is currently
+                                       * necessary that each
+                                       * processor holds the
+                                       * sparsity_pattern structure
+                                       * because each processor sets
+                                       * its rows.
                                        *
-                                       * This is a
-                                       * collective operation that
-                                       * needs to be called on all
-                                       * processors in order to avoid a
-                                       * dead lock.
+                                       * This is a collective
+                                       * operation that needs to be
+                                       * called on all processors in
+                                       * order to avoid a dead lock.
                                         */
-      void reinit (const SparsityPattern &sparsity_pattern);
+      template<typename SparsityType>
+      void reinit (const SparsityType &sparsity_pattern);
 
                                       /**
                                         * This function is initializes
@@ -526,8 +535,9 @@ namespace TrilinosWrappers
                                        * processors in order to avoid a
                                        * dead lock.
                                         */
-      void reinit (const Epetra_Map       &input_map,
-                  const SparsityPattern  &sparsity_pattern);
+      template<typename SparsityType>
+      void reinit (const Epetra_Map    &input_map,
+                  const SparsityType  &sparsity_pattern);
 
                                       /**
                                         * This function is similar to
@@ -545,9 +555,10 @@ namespace TrilinosWrappers
                                        * processors in order to avoid a
                                        * dead lock.
                                         */
-      void reinit (const Epetra_Map       &input_row_map,
-                  const Epetra_Map       &input_col_map,
-                  const SparsityPattern  &sparsity_pattern);
+      template<typename SparsityType>
+      void reinit (const Epetra_Map    &input_row_map,
+                  const Epetra_Map    &input_col_map,
+                  const SparsityType  &sparsity_pattern);
 
                                       /**
                                        * This function copies the
@@ -581,6 +592,34 @@ namespace TrilinosWrappers
                                        * processors in order to avoid a
                                        * dead lock.
                                         */
+      void reinit (const ::dealii::SparseMatrix<double> &dealii_sparse_matrix,
+                  const double                          drop_tolerance=1e-13);
+
+                                      /**
+                                        * This function initializes
+                                       * the Trilinos matrix using
+                                       * the deal.II sparse matrix
+                                       * and the entries stored
+                                       * therein. It uses a threshold
+                                       * to copy only elements with
+                                       * modulus larger than the
+                                       * threshold (so zeros in the
+                                       * deal.II matrix can be
+                                       * filtered away). In contrast
+                                       * to the other reinit function
+                                       * with deal.II sparse matrix
+                                       * argument, this function
+                                       * takes a parallel
+                                       * partitioning specified by
+                                       * the user instead of
+                                       * internally generating one.
+                                       *
+                                       * This is a
+                                       * collective operation that
+                                       * needs to be called on all
+                                       * processors in order to avoid a
+                                       * dead lock.
+                                        */
       void reinit (const Epetra_Map                     &input_map,
                   const ::dealii::SparseMatrix<double> &dealii_sparse_matrix,
                   const double                          drop_tolerance=1e-13);
index 0995b26e7bd9dbc95b262023e26e0fb41078ac98..6c8272208634c608c15fbf12d366f84d2def4b24 100644 (file)
@@ -13,6 +13,7 @@
 
 #include <lac/trilinos_block_sparse_matrix.h>
 
+#include <lac/block_sparsity_pattern.h>
 
 #ifdef DEAL_II_USE_TRILINOS
 
@@ -88,10 +89,11 @@ namespace TrilinosWrappers
 
 
 
+  template <typename BlockSparsityType>  
   void
   BlockSparseMatrix::
   reinit (const std::vector<Epetra_Map> &input_maps,
-         const BlockSparsityPattern    &block_sparsity_pattern)
+         const BlockSparsityType       &block_sparsity_pattern)
   {
     Assert (input_maps.size() == block_sparsity_pattern.n_block_rows(),
            ExcDimensionMismatch (input_maps.size(),
@@ -129,9 +131,10 @@ namespace TrilinosWrappers
 
 
 
+  template <typename BlockSparsityType>
   void
   BlockSparseMatrix::
-  reinit (const BlockSparsityPattern    &block_sparsity_pattern)
+  reinit (const BlockSparsityType &block_sparsity_pattern)
   {
     Assert (block_sparsity_pattern.n_block_rows() ==
            block_sparsity_pattern.n_block_cols(),
@@ -158,7 +161,7 @@ namespace TrilinosWrappers
 
     reinit (input_maps, block_sparsity_pattern);
   }
-  
+
 
 
   void
@@ -194,6 +197,39 @@ namespace TrilinosWrappers
 
 
 
+  void
+  BlockSparseMatrix::
+  reinit (const ::dealii::BlockSparseMatrix<double> &dealii_block_sparse_matrix,
+         const double                               drop_tolerance)
+  {
+    Assert (dealii_block_sparse_matrix.n_block_rows() ==
+           dealii_block_sparse_matrix.n_block_cols(),
+           ExcDimensionMismatch (dealii_block_sparse_matrix.n_block_rows(),
+                                 dealii_block_sparse_matrix.n_block_cols()));
+    Assert (dealii_block_sparse_matrix.m() ==
+           dealii_block_sparse_matrix.n(),
+           ExcDimensionMismatch (dealii_block_sparse_matrix.m(),
+                                 dealii_block_sparse_matrix.n()));
+    
+                                    // produce a dummy local map and pass it
+                                    // off to the other function
+#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
+    Epetra_MpiComm    trilinos_communicator (MPI_COMM_WORLD);
+#else
+    Epetra_SerialComm trilinos_communicator;
+#endif
+
+    std::vector<Epetra_Map> input_maps;
+    for (unsigned int i=0; i<dealii_block_sparse_matrix.n_block_rows(); ++i)
+      input_maps.push_back (Epetra_Map(dealii_block_sparse_matrix.block(i,0).m(),
+                                      0,
+                                      trilinos_communicator));
+
+    reinit (input_maps, dealii_block_sparse_matrix, drop_tolerance);
+  }
+
+
+
   void
   BlockSparseMatrix::compress()
   {
@@ -326,6 +362,35 @@ namespace TrilinosWrappers
     return dst.l2_norm();
   }
 
+
+
+
+
+  // -------------------- explicit instantiations -----------------------
+  //
+  template void
+  BlockSparseMatrix::reinit (const BlockSparsityPattern &);
+  template void
+  BlockSparseMatrix::reinit (const BlockCompressedSparsityPattern &);
+  template void
+  BlockSparseMatrix::reinit (const BlockCompressedSetSparsityPattern &);
+  template void
+  BlockSparseMatrix::reinit (const BlockCompressedSimpleSparsityPattern &);
+
+
+  template void
+  BlockSparseMatrix::reinit (const std::vector<Epetra_Map> &,
+                            const BlockSparsityPattern    &);
+  template void
+  BlockSparseMatrix::reinit (const std::vector<Epetra_Map> &,
+                            const BlockCompressedSparsityPattern &);
+  template void
+  BlockSparseMatrix::reinit (const std::vector<Epetra_Map> &,
+                            const BlockCompressedSetSparsityPattern &);
+  template void
+  BlockSparseMatrix::reinit (const std::vector<Epetra_Map> &,
+                            const BlockCompressedSimpleSparsityPattern &);
+
 }
 
 
index 4953d8e3dbabdeac9388991b5d7e9694147f247f..36a3844dfd622724105600379dba20476264697b 100755 (executable)
 
 #include <lac/trilinos_sparse_matrix.h>
 
+#include <lac/sparsity_pattern.h>
+#include <lac/compressed_sparsity_pattern.h>
+#include <lac/compressed_set_sparsity_pattern.h>
+#include <lac/compressed_simple_sparsity_pattern.h>
+
 #ifdef DEAL_II_USE_TRILINOS
 
 DEAL_II_NAMESPACE_OPEN
@@ -172,8 +177,9 @@ namespace TrilinosWrappers
   
 
 
+  template <typename SparsityType>
   void
-  SparseMatrix::reinit (const SparsityPattern &sparsity_pattern)
+  SparseMatrix::reinit (const SparsityType &sparsity_pattern)
   {
 #ifdef DEAL_II_COMPILER_SUPPORTS_MPI
     Epetra_MpiComm    trilinos_communicator (MPI_COMM_WORLD);
@@ -193,19 +199,21 @@ namespace TrilinosWrappers
 
 
 
+  template <typename SparsityType>
   void
-  SparseMatrix::reinit (const Epetra_Map       &input_map,
-                       const SparsityPattern  &sparsity_pattern)
+  SparseMatrix::reinit (const Epetra_Map    &input_map,
+                       const SparsityType  &sparsity_pattern)
   {
     reinit (input_map, input_map, sparsity_pattern);
   }
 
 
 
+  template <typename SparsityType>
   void
-  SparseMatrix::reinit (const Epetra_Map       &input_row_map,
-                       const Epetra_Map       &input_col_map,
-                       const SparsityPattern  &sparsity_pattern)
+  SparseMatrix::reinit (const Epetra_Map    &input_row_map,
+                       const Epetra_Map    &input_col_map,
+                       const SparsityType  &sparsity_pattern)
   {
     matrix.reset();
 
@@ -295,6 +303,110 @@ namespace TrilinosWrappers
 
 
 
+                                  // The CompressedSetSparsityPattern
+                                  // class stores the columns
+                                  // differently, so we need to
+                                  // explicitly rewrite this function
+                                  // in that case.
+  template<>
+  void
+  SparseMatrix::reinit (const Epetra_Map                    &input_row_map,
+                       const Epetra_Map                    &input_col_map,
+                       const CompressedSetSparsityPattern  &sparsity_pattern)
+  {
+    matrix.reset();
+
+    unsigned int n_rows = sparsity_pattern.n_rows();
+
+    if (input_row_map.Comm().MyPID() == 0)
+      {
+       Assert (input_row_map.NumGlobalElements() == (int)sparsity_pattern.n_rows(),
+               ExcDimensionMismatch (input_row_map.NumGlobalElements(),
+                                     sparsity_pattern.n_rows()));
+       Assert (input_col_map.NumGlobalElements() == (int)sparsity_pattern.n_cols(),
+               ExcDimensionMismatch (input_col_map.NumGlobalElements(),
+                                     sparsity_pattern.n_cols()));
+      }
+
+    row_map = input_row_map;
+    col_map = input_col_map;
+
+    std::vector<int> n_entries_per_row(n_rows);
+
+    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
+                                 // 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));
+
+
+                                 // TODO: As of now, assume that the
+                                 // sparsity pattern sits at the all
+                                 // processors (completely), and let
+                                 // 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()));
+    
+                                 // Trilinos has a bug for rectangular
+                                 // matrices at this point, so do not
+                                 // check for consistent column numbers
+                                 // here.
+                                 //
+                                 // 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(),
+//                                   sparsity_pattern.n_cols()));
+
+    std::vector<double> values;
+    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);
+
+         for (unsigned int col = 0; 
+              col_num != sparsity_pattern.row_end (row); 
+              ++col_num, ++col)
+           row_indices[col] = *col_num;
+
+         matrix->InsertGlobalValues (row, row_length,
+                                     &values[0], &row_indices[0]);
+       }
+    
+    last_action = Zero;
+
+                                 // In the end, the matrix needs to
+                                 // be compressed in order to be
+                                 // really ready.
+    compress();
+  }
+
+
+
   void
   SparseMatrix::reinit (const SparseMatrix &sparse_matrix)
   {
@@ -309,6 +421,27 @@ namespace TrilinosWrappers
 
 
 
+  void
+  SparseMatrix::reinit (const ::dealii::SparseMatrix<double> &dealii_sparse_matrix,
+                       const double                          drop_tolerance)
+  {
+#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
+    Epetra_MpiComm    trilinos_communicator (MPI_COMM_WORLD);
+#else
+    Epetra_SerialComm trilinos_communicator;
+#endif
+
+    const Epetra_Map rows (dealii_sparse_matrix.m(),
+                          0,
+                          trilinos_communicator);
+    const Epetra_Map columns (dealii_sparse_matrix.n(),
+                             0,
+                             trilinos_communicator);
+    reinit (rows, columns, dealii_sparse_matrix, drop_tolerance);
+  }
+
+
+
   void
   SparseMatrix::reinit (const Epetra_Map                     &input_map,
                        const ::dealii::SparseMatrix<double> &dealii_sparse_matrix,
@@ -1104,7 +1237,8 @@ namespace TrilinosWrappers
                                  // As of now, no particularly neat
                                  // ouput is generated in case of
                                  // multiple processors.
-  void SparseMatrix::print (std::ostream &out) const
+  void
+  SparseMatrix::print (std::ostream &out) const
   {
     double * values;
     int * indices;
@@ -1121,6 +1255,48 @@ namespace TrilinosWrappers
     AssertThrow (out, ExcIO());
   }
 
+
+
+
+  // explicit instantiations
+  //
+  template void
+  SparseMatrix::reinit (const SparsityPattern &);
+  template void
+  SparseMatrix::reinit (const CompressedSparsityPattern &);
+  template void
+  SparseMatrix::reinit (const CompressedSetSparsityPattern &);
+  template void
+  SparseMatrix::reinit (const CompressedSimpleSparsityPattern &);
+
+
+  template void
+  SparseMatrix::reinit (const Epetra_Map &,
+                       const SparsityPattern &);
+  template void
+  SparseMatrix::reinit (const Epetra_Map &,
+                       const CompressedSparsityPattern &);
+  template void
+  SparseMatrix::reinit (const Epetra_Map &,
+                       const CompressedSetSparsityPattern &);
+  template void
+  SparseMatrix::reinit (const Epetra_Map &,
+                       const CompressedSimpleSparsityPattern &);
+
+
+  template void
+  SparseMatrix::reinit (const Epetra_Map &,
+                       const Epetra_Map &,
+                       const SparsityPattern &);
+  template void
+  SparseMatrix::reinit (const Epetra_Map &,
+                       const Epetra_Map &,
+                       const CompressedSparsityPattern &);
+  template void
+  SparseMatrix::reinit (const Epetra_Map &,
+                       const Epetra_Map &,
+                       const CompressedSimpleSparsityPattern &);
+
 }
 
 DEAL_II_NAMESPACE_CLOSE

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.