]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Fix bug in PETScWrappers::MPI::BlockSparseMatrix by synchronizing the state (set...
authorTimo Heister <timo.heister@gmail.com>
Wed, 5 Aug 2009 10:25:07 +0000 (10:25 +0000)
committerTimo Heister <timo.heister@gmail.com>
Wed, 5 Aug 2009 10:25:07 +0000 (10:25 +0000)
git-svn-id: https://svn.dealii.org/trunk@19179 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.h
deal.II/lac/include/lac/block_matrix_base.h
deal.II/lac/include/lac/petsc_matrix_base.h
deal.II/lac/include/lac/petsc_parallel_sparse_matrix.h
deal.II/lac/include/lac/petsc_sparse_matrix.h
deal.II/lac/include/lac/sparse_matrix.h
deal.II/lac/include/lac/trilinos_sparse_matrix.h

index 1e5436f06a91822b22f096637307aeb7763b1146..9109f128dc98e069d9d83ac8fa6dae8726e26e89 100644 (file)
@@ -224,6 +224,18 @@ inconvenience this causes.
 <h3>lac</h3>
 
 <ol>
+  <li>
+  <p>
+  Fixed: Crash or strange behaviour (wrong matrix entries written) in
+  PETScWrappers::MPI::BlockSparseMatrix when adding or setting elements
+  through any of the set() and add() routines. This happened when different
+  CPUs access different blocks at the start of assembly or when switching
+  between adding and setting.
+  <br>
+  (Timo Heister 2009/08/05)
+  </p>
+  </li>
+  
   <li> <p>New: The relaxation preconditioners PreconditionJacobi, PreconditionSOR and
   PreconditionSSOR, as well as their blocked versions PreconditionBlockJacobi,
   PreconditionBlockSOR and PreconditionBlockSSOR now have functions <code>step</code>
@@ -231,6 +243,7 @@ inconvenience this causes.
   <br>
   (GK 2009/08/04)
   </p>
+  </li>
 
   <li>
   <p>
index 663803a7ae29d5f8c4de91a8a59303ec3d21eb52..5af050a1b380daaf971eaf25b615226844903636 100644 (file)
@@ -1244,6 +1244,29 @@ class BlockMatrixBase : public Subscriptor
     void Tvmult_nonblock_nonblock (VectorType       &dst,
                                   const VectorType &src) const;
     
+
+  protected:
+
+                                    /**
+                                     * Some matrix types, in particular PETSc,
+                                     * need to synchronize set and add
+                                     * operations. This has to be done for all
+                                     * matrices in the BlockMatrix.
+                                     * This routine prepares adding of elements
+                                     * by notifying all blocks. Called by all
+                                     * internal routines before adding
+                                     * elements.
+                                     */ 
+    void prepare_add_operation();
+
+                                    /**
+                                     * Notifies all blocks to let them prepare
+                                     * for setting elements, see
+                                     * prepare_add_operation().
+                                     */
+    void prepare_set_operation();
+
+
   private:
                                     /**
                                      * Temporary vector for counting the
@@ -1868,6 +1891,7 @@ BlockMatrixBase<MatrixType>::set (const unsigned int i,
                                   const unsigned int j,
                                   const value_type value)
 {
+  prepare_set_operation();
 
   Assert (numbers::is_finite(value), 
           ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
@@ -1953,6 +1977,8 @@ BlockMatrixBase<MatrixType>::set (const unsigned int  row,
                                  const number       *values,
                                  const bool          elide_zero_values)
 {
+  prepare_set_operation();
+       
                                   // Resize scratch arrays
   if (column_indices.size() < this->n_block_cols())
     {
@@ -2053,6 +2079,8 @@ BlockMatrixBase<MatrixType>::add (const unsigned int i,
   Assert (numbers::is_finite(value), 
           ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
 
+  prepare_add_operation();
+
                                    // save some cycles for zero additions, but
                                    // only if it is safe for the matrix we are
                                    // working with
@@ -2144,6 +2172,8 @@ BlockMatrixBase<MatrixType>::add (const unsigned int   row,
                                  const bool           elide_zero_values,
                                  const bool           col_indices_are_sorted)
 {
+  prepare_add_operation();
+       
                                   // TODO: Look over this to find out
                                   // whether we can do that more
                                   // efficiently.
@@ -2824,6 +2854,28 @@ BlockMatrixBase<MatrixType>::collect_sizes ()
   this->column_block_indices.reinit (col_sizes);
 }
 
+
+
+template <class MatrixType>
+void
+BlockMatrixBase<MatrixType>::prepare_add_operation ()
+{
+  for (unsigned int row=0; row<n_block_rows(); ++row)
+    for (unsigned int col=0; col<n_block_cols(); ++col)
+      block(row, col).prepare_add();
+
+}
+
+template <class MatrixType>
+void
+BlockMatrixBase<MatrixType>::prepare_set_operation ()
+{
+  for (unsigned int row=0; row<n_block_rows(); ++row)
+    for (unsigned int col=0; col<n_block_cols(); ++col)
+      block(row, col).prepare_set();
+
+}
+
 #endif // DOXYGEN
 
 
index 0a5589fb2c4bf4f49dff5370e7dec65790d940ac..786773dbeff4ec76466587f281049dc22926d048 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2004, 2005, 2006, 2007, 2008 by the deal.II authors
+//    Copyright (C) 2004, 2005, 2006, 2007, 2008, 2009 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -1157,6 +1157,49 @@ namespace PETScWrappers
                                         */
       LastAction::Values last_action;
 
+                                      /**
+                                       * Flush buffers on all CPUs when
+                                       * switching between inserting and
+                                       * adding to elements, no-op otherwise.
+                                       * Should be called from all internal
+                                       * functions accessing matrix elements.
+                                       */
+      void prepare_action(const LastAction::Values new_action);
+
+                                      /**
+                                       * For some matrix storage
+                                        * formats, in particular for the
+                                        * PETSc distributed blockmatrices,
+                                        * set and add operations on
+                                        * individual elements can not be
+                                        * freely mixed. Rather, one has
+                                        * to synchronize operations when
+                                        * one wants to switch from
+                                        * setting elements to adding to
+                                        * elements.
+                                        * BlockMatrixBase automatically
+                                        * synchronizes the access by
+                                        * calling this helper function
+                                        * for each block.
+                                        * This function ensures that the
+                                        * matrix is in a state that
+                                        * allows adding elements; if it
+                                        * previously already was in this
+                                        * state, the function does
+                                        * nothing.
+                                        */
+      void prepare_add();
+                                      /**
+                                        * Same as prepare_add() but
+                                        * prepare the matrix for setting
+                                        * elements if the representation
+                                        * of elements in this class
+                                        * requires such an operation.
+                                        */
+      void prepare_set();
+      
+
+      
     private:
                                       /**
                                        * An internal array of integer
@@ -1176,6 +1219,13 @@ namespace PETScWrappers
                                        */
       std::vector<PetscScalar> column_values;
 
+
+                                      /**
+                                       *  To allow calling protected prepare_add()
+                                       *  and prepare_set().
+                                       */ 
+      template <class> friend class BlockMatrixBase;
+
   };
 
 
@@ -1417,18 +1467,7 @@ namespace PETScWrappers
                   const PetscScalar  *values,
                   const bool          elide_zero_values)
   {
-
-    if (last_action != LastAction::insert)
-      {
-        int ierr;
-        ierr = MatAssemblyBegin(matrix,MAT_FLUSH_ASSEMBLY);
-        AssertThrow (ierr == 0, ExcPETScError(ierr));
-
-        ierr = MatAssemblyEnd(matrix,MAT_FLUSH_ASSEMBLY);
-        AssertThrow (ierr == 0, ExcPETScError(ierr));
-
-       last_action = LastAction::insert;
-      }
+    prepare_action(LastAction::insert);
     
     const signed int petsc_i = row;
     int * col_index_ptr;
@@ -1504,17 +1543,7 @@ namespace PETScWrappers
                                  // zero. However, these actions are done
                                  // in case we pass on to the other
                                  // function.
-       if (last_action != LastAction::add)
-         {
-           int ierr;
-           ierr = MatAssemblyBegin(matrix,MAT_FLUSH_ASSEMBLY);
-           AssertThrow (ierr == 0, ExcPETScError(ierr));
-
-           ierr = MatAssemblyEnd(matrix,MAT_FLUSH_ASSEMBLY);
-           AssertThrow (ierr == 0, ExcPETScError(ierr));
-
-           last_action = LastAction::add;
-         }
+       prepare_action(LastAction::add);
 
        return;
       }
@@ -1585,18 +1614,7 @@ namespace PETScWrappers
                   const bool          elide_zero_values,
                   const bool          /*col_indices_are_sorted*/)
   {
-    if (last_action != LastAction::add)
-      {
-        int ierr;
-        ierr = MatAssemblyBegin(matrix,MAT_FLUSH_ASSEMBLY);
-        AssertThrow (ierr == 0, ExcPETScError(ierr));
-
-        ierr = MatAssemblyEnd(matrix,MAT_FLUSH_ASSEMBLY);
-        AssertThrow (ierr == 0, ExcPETScError(ierr));
-
-        last_action = LastAction::add;
-      }
-    
+    prepare_action(LastAction::add);
     
     const signed int petsc_i = row;
     int * col_index_ptr;
@@ -1723,6 +1741,44 @@ namespace PETScWrappers
             (index < static_cast<unsigned int>(end)));
   }
 
+  
+
+  inline
+  void
+  MatrixBase::prepare_action(const LastAction::Values new_action)
+  {
+                                    // flush PETSc buffers when switching
+                                    // actions, otherwise just return.
+    if (last_action == new_action) return;
+
+    int ierr;
+    ierr = MatAssemblyBegin(matrix, MAT_FLUSH_ASSEMBLY);
+    AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+    ierr = MatAssemblyEnd(matrix, MAT_FLUSH_ASSEMBLY);
+    AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+    last_action = new_action;
+  }
+
+
+  
+  inline
+  void
+  MatrixBase::prepare_add()
+  {
+    prepare_action(LastAction::add);
+  }
+
+
+  
+  inline
+  void
+  MatrixBase::prepare_set()
+  {
+    prepare_action(LastAction::insert);
+  }
+  
 #endif // DOXYGEN      
 }
 
index 6d81587c7f1887748543f85eba6e241b67237ac6..a60e9c9bad949c85b126a45f68b1c2f5e6ccab9c 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2004, 2005, 2006, 2007 by the deal.II authors
+//    Copyright (C) 2004, 2005, 2006, 2007, 2009 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
 DEAL_II_NAMESPACE_OPEN
 
 
+                                 // forward declaration
+template <typename MatrixType>
+class BlockMatrixBase;
+
 
 namespace PETScWrappers
 {
@@ -423,6 +427,12 @@ namespace PETScWrappers
                         const std::vector<unsigned int> &local_columns_per_process,
                         const unsigned int               this_process,
                         const bool                       preset_nonzero_locations);
+
+                                      /**
+                                       *  To allow calling protected prepare_add()
+                                       *  and prepare_set().
+                                       */  
+       friend class BlockMatrixBase<SparseMatrix>;
     };
 
 
index 54a397800196b907de9230704c1be51dee217a59..37d3b79e179a30020a7d838ca963a28eb1fc577e 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2004, 2005, 2006, 2007 by the deal.II authors
+//    Copyright (C) 2004, 2005, 2006, 2007, 2009 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -23,6 +23,9 @@
 #ifdef DEAL_II_USE_PETSC
 
 DEAL_II_NAMESPACE_OPEN
+                                 // forward declaration
+template <typename MatrixType>
+class BlockMatrixBase;
 
 
 namespace PETScWrappers
@@ -318,6 +321,12 @@ namespace PETScWrappers
       template <typename SparsityType>
       void do_reinit (const SparsityType &sparsity_pattern,
                       const bool          preset_nonzero_locations);
+                      
+                                      /**
+                                       *  To allow calling protected prepare_add()
+                                       *  and prepare_set().
+                                       */                       
+         friend class BlockMatrixBase<SparseMatrix>;
   };
 }
 
index f80ac00f6d93e13fd6e5989df255a3783b842a54..1bb763cf7bdd09192c188ca12069730c2a42f4d1 100644 (file)
@@ -2008,7 +2008,42 @@ class SparseMatrix : public virtual Subscriptor
                                       * Exception
                                       */
     DeclException0 (ExcSourceEqualsDestination);
-                                    //@}    
+                                    //@}
+
+  protected:
+                                     /**
+                                      * For some matrix storage
+                                      * formats, in particular for the
+                                      * PETSc distributed blockmatrices,
+                                      * set and add operations on
+                                      * individual elements can not be
+                                      * freely mixed. Rather, one has
+                                      * to synchronize operations when
+                                      * one wants to switch from
+                                      * setting elements to adding to
+                                      * elements.
+                                      * BlockMatrixBase automatically
+                                      * synchronizes the access by
+                                      * calling this helper function
+                                      * for each block.
+                                      * This function ensures that the
+                                      * matrix is in a state that
+                                      * allows adding elements; if it
+                                      * previously already was in this
+                                      * state, the function does
+                                      * nothing.
+                                      */
+    void prepare_add();
+
+                                     /**
+                                      * Same as prepare_add() but
+                                      * prepare the matrix for setting
+                                      * elements if the representation
+                                      * of elements in this class
+                                      * requires such an operation.
+                                      */
+    void prepare_set();
+    
   private:
                                     /**
                                      * Pointer to the sparsity
@@ -2053,6 +2088,13 @@ class SparseMatrix : public virtual Subscriptor
     template <typename somenumber> friend class SparseMatrix;
     template <typename somenumber> friend class SparseLUDecomposition;
     template <typename> friend class SparseILU;
+    
+                                      /**
+                                       * To allow it calling private
+                                        * prepare_add() and
+                                        * prepare_set().
+                                       */  
+    template <typename> friend class BlockMatrixBase;
 };
 
 /*@}*/
@@ -2862,6 +2904,28 @@ SparseMatrix<number>::end (const unsigned int r)
   return end();
 }
 
+
+
+template <typename number>
+inline
+void
+SparseMatrix<number>::
+prepare_add()
+{
+                                   //nothing to do here
+}
+
+
+
+template <typename number>
+inline
+void
+SparseMatrix<number>::
+prepare_set()
+{
+                                   //nothing to do here
+}
+
 #endif // DOXYGEN
 
 
index 6f96a1378e016398e5e8d22d5f8c27f711f209f2..9a220c1072f7caa50e549072f973d79b40ee5f0e 100755 (executable)
 
 DEAL_II_NAMESPACE_OPEN
 
-// forward declaration
+                                // forward declarations
+template <typename MatrixType>
+class BlockMatrixBase;
+
 template <typename number> class SparseMatrix;
 
 namespace TrilinosWrappers
@@ -1862,9 +1865,50 @@ namespace TrilinosWrappers
                      << "/" << arg2 << ")"
                      << " of a sparse matrix, but it appears to not"
                      << " exist in the Trilinos sparsity pattern.");
-                                    //@}    
+                                    //@}
+                                
+
+
+    protected:    
+
+                                      /**
+                                      * For some matrix storage
+                                      * formats, in particular for the
+                                      * PETSc distributed blockmatrices,
+                                      * set and add operations on
+                                      * individual elements can not be
+                                      * freely mixed. Rather, one has
+                                      * to synchronize operations when
+                                      * one wants to switch from
+                                      * setting elements to adding to
+                                      * elements.
+                                      * BlockMatrixBase automatically
+                                      * synchronizes the access by
+                                      * calling this helper function
+                                      * for each block.
+                                      * This function ensures that the
+                                      * matrix is in a state that
+                                      * allows adding elements; if it
+                                      * previously already was in this
+                                      * state, the function does
+                                      * nothing.
+                                      */
+      void prepare_add();
+
+                                     /**
+                                      * Same as prepare_add() but
+                                      * prepare the matrix for setting
+                                      * elements if the representation
+                                      * of elements in this class
+                                      * requires such an operation.
+                                      */
+      void prepare_set();
+
+
+      
     private:
-                                       /**
+      
+                                      /**
                                        * Epetra Trilinos
                                        * mapping of the matrix rows
                                        * that assigns parts of the
@@ -1959,6 +2003,11 @@ namespace TrilinosWrappers
                                         */
       std::auto_ptr<Epetra_FECrsMatrix> matrix;
 
+                                      /**
+                                       *  To allow calling protected prepare_add()
+                                       *  and prepare_set().
+                                       */ 
+         friend class BlockMatrixBase<SparseMatrix>;
   };
 
 
@@ -2960,6 +3009,26 @@ namespace TrilinosWrappers
     return matrix->ColMap();
   }
 
+  
+
+  inline
+  void
+  SparseMatrix::prepare_add()
+  {
+                                   //nothing to do here
+  }
+
+
+  
+  inline
+  void
+  SparseMatrix::prepare_set()
+  {
+                                   //nothing to do here
+  }
+
+  
+
 #endif // DOXYGEN      
 }
 

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.