]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Extended the rudimentary Trilinos wrappers a little bit. Rectangular matrices do...
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 21 Aug 2008 18:03:55 +0000 (18:03 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 21 Aug 2008 18:03:55 +0000 (18:03 +0000)
git-svn-id: https://svn.dealii.org/trunk@16629 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/trilinos_sparse_matrix.h
deal.II/lac/include/lac/trilinos_vector.h
deal.II/lac/source/trilinos_sparse_matrix.cc
deal.II/lac/source/trilinos_vector.cc

index dee47baf3885bc03d26d8553fa6549e7485f4237..becabaa03335490d05d5b87c78b7377140b366b9 100755 (executable)
@@ -337,29 +337,52 @@ namespace TrilinosWrappers
                    const std::vector<unsigned int> &n_entries_per_row);
 
                                        /**
-                                        * Destructor. Made virtual so that one
-                                        * can use pointers to this class.
+                                        * This function is similar to the
+                                       * one above, but it now takes two
+                                       * different Epetra maps for rows
+                                       * and columns. This interface is
+                                       * meant to be used for generating
+                                       * rectangular matrices, where one
+                                       * map specifies the parallel 
+                                       * distribution of rows and the other
+                                       * the one of the columns. This is
+                                       * in contrast to the first 
+                                       * constructor, where the same map
+                                       * is used for both the number of 
+                                       * rows and the number of columns.
+                                       * The number
+                                       * of columns per row is specified
+                                       * by the maximum number of entries.
                                         */
-      virtual ~SparseMatrix ();
+      SparseMatrix (const Epetra_Map   &InputRowMap,
+                   const Epetra_Map   &InputColMap,
+                   const unsigned int  n_max_entries_per_row);
+
                                        /**
-                                        * This operator assigns a scalar to a
-                                        * matrix. Since this does usually not
-                                        * make much sense (should we set all
-                                        * matrix entries to this value? Only
-                                        * the nonzero entries of the sparsity
-                                        * pattern?), this operation is only
-                                        * allowed if the actual value to be
-                                        * assigned is zero. This operator only
-                                        * exists to allow for the obvious
-                                        * notation <tt>matrix=0</tt>, which
-                                        * sets all elements of the matrix to
-                                        * zero, but keeps the sparsity pattern
-                                        * previously used.
+                                        * This function is similar to the
+                                       * one above, but it now takes two
+                                       * different Epetra maps for rows
+                                       * and columns. This interface is
+                                       * meant to be used for generating
+                                       * rectangular matrices, where one
+                                       * map specifies the parallel 
+                                       * distribution of rows and the other
+                                       * the one of the columns. The 
+                                       * vector n_entries_per_row specifies
+                                       * the number of entries in each 
+                                       * row of the newly generated matrix.
                                         */
-      SparseMatrix &
-      operator = (const double d);
+      SparseMatrix (const Epetra_Map                &InputRowMap,
+                   const Epetra_Map                &InputColMap,
+                   const std::vector<unsigned int> &n_entries_per_row);
 
                                        /**
+                                        * Destructor. Made virtual so that one
+                                        * can use pointers to this class.
+                                        */
+      virtual ~SparseMatrix ();
+                                       
+                                      /**
                                         * This function initializes the
                                        * Trilinos matrix by attaching all
                                        * the elements to the sparsity
@@ -394,6 +417,19 @@ namespace TrilinosWrappers
                                        * e.g. when the grid is refined.
                                         */
       void reinit (const Epetra_Map                &input_map,
+                  const CompressedSparsityPattern &sparsity_pattern);
+
+                                      /**
+                                        * This function is similar to the
+                                       * other initialization function above,
+                                       * but now also reassigns the matrix 
+                                       * rows and columns according to 
+                                       * two user-supplied Epetra maps.
+                                       * To be used e.g. for rectangular 
+                                       * matrices after remeshing.
+                                        */
+      void reinit (const Epetra_Map                &input_row_map,
+                  const Epetra_Map                &input_col_map,
                   const CompressedSparsityPattern &sparsity_pattern);
 
                                        /**
@@ -404,6 +440,24 @@ namespace TrilinosWrappers
                                         */
       void clear ();
 
+                                       /**
+                                        * This operator assigns a scalar to a
+                                        * matrix. Since this does usually not
+                                        * make much sense (should we set all
+                                        * matrix entries to this value? Only
+                                        * the nonzero entries of the sparsity
+                                        * pattern?), this operation is only
+                                        * allowed if the actual value to be
+                                        * assigned is zero. This operator only
+                                        * exists to allow for the obvious
+                                        * notation <tt>matrix=0</tt>, which
+                                        * sets all elements of the matrix to
+                                        * zero, but keeps the sparsity pattern
+                                        * previously used.
+                                        */
+      SparseMatrix &
+      operator = (const double d);
+
                                        /**
                                         * Set the element (<i>i,j</i>)
                                         * to @p value.
@@ -929,7 +983,7 @@ namespace TrilinosWrappers
                                        /**
                                         * Exception
                                         */
-      DeclException4 (ExcAccessToNonlocalElement,
+      DeclException4 (ExcAccessToNonLocalElement,
                      int, int, int, int,
                      << "You tried to access element (" << arg1
                      << "/" << arg2 << ")"
@@ -938,11 +992,33 @@ namespace TrilinosWrappers
                      << " are stored locally and can be accessed.");
 
                                        /**
-                                        * The Epetra Trilinos mapping that
+                                        * Exception
+                                        */
+      DeclException2 (ExcAccessToNonPresentElement,
+                     int, int,
+                     << "You tried to access element (" << arg1
+                     << "/" << arg2 << ")"
+                     << " of a sparse matrix, but it appears to not "
+                     << " exist in the Trilinos sparsity pattern.");
+
+                                       /**
+                                       * Epetra Trilinos mapping of the 
+                                       * matrix rows that
+                                       * assigns parts of the matrix to
+                                       * the individual processes.
+                                       * TODO: is it possible to only use
+                                       * a pointer?
+                                       */
+      Epetra_Map row_map;
+
+                                       /**
+                                        * Pointer to the user-supplied 
+                                       * Epetra Trilinos mapping of the 
+                                       * matrix columns that
                                        * assigns parts of the matrix to
                                        * the individual processes.
                                        */
-      Epetra_Map map;
+      Epetra_Map col_map;
 
                                        /**
                                         * A sparse matrix object in
@@ -1040,7 +1116,7 @@ namespace TrilinosWrappers
 
     inline
     const_iterator::
-    const_iterator(const SparseMatrix   *matrix,
+    const_iterator(const SparseMatrix *matrix,
                    const unsigned int  row,
                    const unsigned int  index)
                     :
index 7998b13a1345ed2a45f76665a0c07529e7625053..39a68a3f5ae291b3330968297491b56e4f4add22 100755 (executable)
@@ -162,7 +162,7 @@ namespace TrilinosWrappers
                                          /**
                                           * Exception
                                           */
-        DeclException3 (ExcAccessToNonlocalElement,
+        DeclException3 (ExcAccessToNonLocalElement,
                         int, int, int,
                         << "You tried to access element " << arg1
                         << " of a distributed vector, but only elements "
@@ -824,12 +824,6 @@ namespace TrilinosWrappers
                                         * the communicator and data distribution
                                        * object common to all
                                         * Trilinos objects used by deal.II.
-                                       * TODO: we probably only need a pointer
-                                       * to the map, since the information
-                                       * is provided from outside and there
-                                       * is no need to copy the map. 
-                                       * Especially not when we have many
-                                       * vectors based on the same map.
                                         */
       Epetra_Map map;
 
index 5ce2e21e09aabc94612757263053bd302b38d2c4..51d0d58d601881f7f0f2256dc4fdadc15eddf612 100755 (executable)
@@ -80,22 +80,23 @@ namespace TrilinosWrappers
   SparseMatrix::SparseMatrix ()
                   :
 #ifdef DEAL_II_COMPILER_SUPPORTS_MPI
-                 map (0,0,Epetra_MpiComm(MPI_COMM_WORLD)),
+                 row_map (0, 0, Epetra_MpiComm(MPI_COMM_WORLD)),
 #else
-                 map (0,0,Epetra_SerialComm()),
+                 row_map (0, 0, Epetra_SerialComm()),
 #endif
+                 col_map (row_map),
                  matrix (std::auto_ptr<Epetra_FECrsMatrix>
-                                (new Epetra_FECrsMatrix(Copy, map, 0))),
+                                (new Epetra_FECrsMatrix(Copy, row_map, 0))),
                   last_action (Insert)
   {}
 
-
   SparseMatrix::SparseMatrix (const Epetra_Map  &InputMap,
                              const unsigned int n_max_entries_per_row)
                   :
-                 map (InputMap),
+                 row_map (InputMap),
+                 col_map (row_map),
                  matrix (std::auto_ptr<Epetra_FECrsMatrix>
-                                (new Epetra_FECrsMatrix(Copy, map, 
+                                (new Epetra_FECrsMatrix(Copy, row_map, 
                                        int(n_max_entries_per_row), false))),
                   last_action (Insert)
   {}
@@ -103,15 +104,42 @@ namespace TrilinosWrappers
   SparseMatrix::SparseMatrix (const Epetra_Map  &InputMap,
                              const std::vector<unsigned int> &n_entries_per_row)
                   :
-                 map (InputMap),
+                 row_map (InputMap),
+                 col_map (row_map),
+                 matrix (std::auto_ptr<Epetra_FECrsMatrix>
+                   (new Epetra_FECrsMatrix(Copy, row_map, 
+                     (int*)const_cast<unsigned int*>(&(n_entries_per_row[0])),
+                     false))),
+                  last_action (Insert)
+  {}
+
+  SparseMatrix::SparseMatrix (const Epetra_Map  &InputRowMap,
+                             const Epetra_Map  &InputColMap,
+                             const unsigned int n_max_entries_per_row)
+                  :
+                 row_map (InputRowMap),
+                 col_map (InputColMap),
+                 matrix (std::auto_ptr<Epetra_FECrsMatrix>
+                                (new Epetra_FECrsMatrix(Copy, row_map, col_map, 
+                                       int(n_max_entries_per_row), false))),
+                  last_action (Insert)
+  {}
+
+  SparseMatrix::SparseMatrix (const Epetra_Map  &InputRowMap,
+                             const Epetra_Map  &InputColMap,
+                             const std::vector<unsigned int> &n_entries_per_row)
+                  :
+                 row_map (InputRowMap),
+                 col_map (InputColMap),
                  matrix (std::auto_ptr<Epetra_FECrsMatrix>
-                   (new Epetra_FECrsMatrix(Copy, map, 
+                   (new Epetra_FECrsMatrix(Copy, row_map, col_map, 
                      (int*)const_cast<unsigned int*>(&(n_entries_per_row[0])),
-                     true))),
+                     false))),
                   last_action (Insert)
   {}
 
 
+
   SparseMatrix::~SparseMatrix ()
   {
   }
@@ -124,14 +152,14 @@ namespace TrilinosWrappers
   {
 
     unsigned int n_rows = sparsity_pattern.n_rows();
-    
+
     Assert (matrix->NumGlobalRows() == (int)sparsity_pattern.n_rows(),
            ExcDimensionMismatch (matrix->NumGlobalRows(),
                                  sparsity_pattern.n_rows()));
-           
+
     std::vector<double> values(n_max_entries_per_row, 0.);
     std::vector<int>    row_indices(n_max_entries_per_row);
-    
+
     for (unsigned int row=0; row<n_rows; ++row)
       {
        const int row_length = sparsity_pattern.row_length(row);
@@ -145,7 +173,7 @@ namespace TrilinosWrappers
                                    &values[0], &row_indices[0]);
       }
 
-                                 // In the end, the matrix is to
+                                 // In the end, the matrix needs to
                                  // be compressed in order to be
                                  // really ready. However, that is
                                  // a collective operation, so it
@@ -156,8 +184,8 @@ namespace TrilinosWrappers
                                  // types are all serial.
   }
 
-  
-  
+
+
   void
   SparseMatrix::reinit (const CompressedSparsityPattern &sparsity_pattern)
   {
@@ -177,9 +205,9 @@ namespace TrilinosWrappers
 
     reinit (sparsity_pattern, n_max_entries_per_row);
   }
-  
-  
-  
+
+
+
   void
   SparseMatrix::reinit (const Epetra_Map                &input_map,
                        const CompressedSparsityPattern &sparsity_pattern)
@@ -187,10 +215,48 @@ namespace TrilinosWrappers
 
     unsigned int n_rows = sparsity_pattern.n_rows();
     matrix.reset();
-    map = input_map;
+    row_map = input_map;
+    
     Assert (input_map.NumGlobalElements() == (int)sparsity_pattern.n_rows(),
            ExcDimensionMismatch (input_map.NumGlobalElements(),
                                  sparsity_pattern.n_rows()));
+    Assert (input_map.NumGlobalElements() == (int)sparsity_pattern.n_cols(),
+           ExcDimensionMismatch (input_map.NumGlobalElements(),
+                                 sparsity_pattern.n_cols()));
+    
+    std::vector<int> n_entries_per_row(n_rows);
+    
+    for (unsigned int row=0; row<n_rows; ++row)
+      n_entries_per_row[(int)row] = sparsity_pattern.row_length(row);
+    
+    matrix = std::auto_ptr<Epetra_FECrsMatrix> (new Epetra_FECrsMatrix (
+                               Copy, row_map, &n_entries_per_row[0], true));
+    
+    const unsigned int n_max_entries_per_row = *std::max_element (
+                   &n_entries_per_row[0], &n_entries_per_row[n_rows-1]);
+
+    reinit (sparsity_pattern, n_max_entries_per_row);
+  }
+
+
+
+  void
+  SparseMatrix::reinit (const Epetra_Map                &input_row_map,
+                       const Epetra_Map                &input_col_map,
+                       const CompressedSparsityPattern &sparsity_pattern)
+  {
+
+    unsigned int n_rows = sparsity_pattern.n_rows();
+    matrix.reset();
+    row_map = input_row_map;
+    col_map = input_col_map;
+    
+    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()));
     
     std::vector<int> n_entries_per_row(n_rows);
     
@@ -198,7 +264,8 @@ namespace TrilinosWrappers
       n_entries_per_row[(int)row] = sparsity_pattern.row_length(row);
     
     matrix = std::auto_ptr<Epetra_FECrsMatrix>
-             (new Epetra_FECrsMatrix(Copy, map, &n_entries_per_row[0], true));
+             (new Epetra_FECrsMatrix(Copy, row_map, col_map,
+                                     &n_entries_per_row[0], true));
     
     const unsigned int n_max_entries_per_row = *std::max_element (
                    &n_entries_per_row[0], &n_entries_per_row[n_rows-1]);
@@ -216,13 +283,13 @@ namespace TrilinosWrappers
                                     // generate an empty matrix.
     matrix.reset();
 #ifdef DEAL_II_COMPILER_SUPPORTS_MPI
-    map = Epetra_Map (0,0,Epetra_MpiComm(MPI_COMM_WORLD)),
+    row_map = Epetra_Map (0,0,Epetra_MpiComm(MPI_COMM_WORLD)),
 #else
-    map = Epetra_Map (0,0,Epetra_SerialComm()),
+    row_map = Epetra_Map (0,0,Epetra_SerialComm()),
 #endif
 
     matrix = std::auto_ptr<Epetra_FECrsMatrix> 
-             (new Epetra_FECrsMatrix(Copy, map, 0));
+             (new Epetra_FECrsMatrix(Copy, row_map, 0));
   }
 
 
@@ -232,7 +299,11 @@ namespace TrilinosWrappers
   {
                                      // flush buffers
     int ierr;
-    ierr = matrix->GlobalAssemble (true);
+    if (row_map.SameAs(col_map))
+      ierr = matrix->GlobalAssemble ();
+    else
+      ierr = matrix->GlobalAssemble (row_map, col_map);
+    
     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
 
     ierr = matrix->OptimizeStorage ();
@@ -269,7 +340,11 @@ namespace TrilinosWrappers
     if (last_action == Add)
       {
         int ierr;
-        ierr = matrix->GlobalAssemble(false);
+       if (row_map.SameAs(col_map))
+         ierr = matrix->GlobalAssemble (false);
+       else
+         ierr = matrix->GlobalAssemble(row_map, col_map, false);
+       
         AssertThrow (ierr == 0, ExcTrilinosError(ierr));
        last_action = Insert;
       }
@@ -281,6 +356,7 @@ namespace TrilinosWrappers
                                                  const_cast<double*>(&value), 
                                                  &trilinos_j);
 
+    AssertThrow (ierr <= 0, ExcAccessToNonPresentElement(i,j));
     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
   }
 
@@ -299,8 +375,12 @@ namespace TrilinosWrappers
     if (last_action == Insert)
       {
         int ierr;
-        ierr = matrix->GlobalAssemble(false);
-        AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+       if (row_map.SameAs(col_map))
+         ierr = matrix->GlobalAssemble (false);
+       else
+         ierr = matrix->GlobalAssemble(row_map, col_map, false);
+        
+       AssertThrow (ierr == 0, ExcTrilinosError(ierr));
 
         last_action = Add;
       }
@@ -322,6 +402,7 @@ namespace TrilinosWrappers
                                                  const_cast<double*>(&value), 
                                                  &trilinos_j);
 
+    AssertThrow (ierr <= 0, ExcAccessToNonPresentElement(i,j));
     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
   }
 
@@ -341,7 +422,7 @@ namespace TrilinosWrappers
                                       // continue.
     if ((trilinos_i == -1 ) || (trilinos_j == -1))
       {
-        Assert (false, ExcAccessToNonlocalElement(i, j, local_range().first,
+        Assert (false, ExcAccessToNonLocalElement(i, j, local_range().first,
                                                  local_range().second));
       }
     else
index 18f8afc3a895e729688edcf87b62b787fc430aa0..8634dbb26990ab8cc7cb08d70a762766c17a2460 100755 (executable)
@@ -37,7 +37,7 @@ namespace TrilinosWrappers
 
       AssertThrow ((static_cast<signed int>(index) >= vector.map.MinMyGID()) &&
                   (static_cast<signed int>(index) <= vector.map.MaxMyGID()),
-                  ExcAccessToNonlocalElement (index, vector.map.MinMyGID(),
+                  ExcAccessToNonLocalElement (index, vector.map.MinMyGID(),
                                               vector.map.MaxMyGID()-1));
       return *(*(vector.vector))[index];
     }

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.