From: Timo Heister Date: Wed, 5 Aug 2009 10:25:07 +0000 (+0000) Subject: Fix bug in PETScWrappers::MPI::BlockSparseMatrix by synchronizing the state (set... X-Git-Tag: v8.0.0~7387 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=57e5aec191a54acc7005d292b05cfdee1eff60f3;p=dealii.git Fix bug in PETScWrappers::MPI::BlockSparseMatrix by synchronizing the state (set/add) between all blocks. git-svn-id: https://svn.dealii.org/trunk@19179 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 1e5436f06a..9109f128dc 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -224,6 +224,18 @@ inconvenience this causes.

lac

    +
  1. +

    + 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. +
    + (Timo Heister 2009/08/05) +

    +
    2. +

  2. New: The relaxation preconditioners PreconditionJacobi, PreconditionSOR and PreconditionSSOR, as well as their blocked versions PreconditionBlockJacobi, PreconditionBlockSOR and PreconditionBlockSSOR now have functions step @@ -231,6 +243,7 @@ inconvenience this causes.
    (GK 2009/08/04)

    +

  3. diff --git a/deal.II/lac/include/lac/block_matrix_base.h b/deal.II/lac/include/lac/block_matrix_base.h index 663803a7ae..5af050a1b3 100644 --- a/deal.II/lac/include/lac/block_matrix_base.h +++ b/deal.II/lac/include/lac/block_matrix_base.h @@ -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::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::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::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::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::collect_sizes () this->column_block_indices.reinit (col_sizes); } + + +template +void +BlockMatrixBase::prepare_add_operation () +{ + for (unsigned int row=0; row +void +BlockMatrixBase::prepare_set_operation () +{ + for (unsigned int row=0; row column_values; + + /** + * To allow calling protected prepare_add() + * and prepare_set(). + */ + template 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(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 } diff --git a/deal.II/lac/include/lac/petsc_parallel_sparse_matrix.h b/deal.II/lac/include/lac/petsc_parallel_sparse_matrix.h index 6d81587c7f..a60e9c9bad 100644 --- a/deal.II/lac/include/lac/petsc_parallel_sparse_matrix.h +++ b/deal.II/lac/include/lac/petsc_parallel_sparse_matrix.h @@ -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 @@ -26,6 +26,10 @@ DEAL_II_NAMESPACE_OPEN + // forward declaration +template +class BlockMatrixBase; + namespace PETScWrappers { @@ -423,6 +427,12 @@ namespace PETScWrappers const std::vector &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; }; diff --git a/deal.II/lac/include/lac/petsc_sparse_matrix.h b/deal.II/lac/include/lac/petsc_sparse_matrix.h index 54a3978001..37d3b79e17 100644 --- a/deal.II/lac/include/lac/petsc_sparse_matrix.h +++ b/deal.II/lac/include/lac/petsc_sparse_matrix.h @@ -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 +class BlockMatrixBase; namespace PETScWrappers @@ -318,6 +321,12 @@ namespace PETScWrappers template void do_reinit (const SparsityType &sparsity_pattern, const bool preset_nonzero_locations); + + /** + * To allow calling protected prepare_add() + * and prepare_set(). + */ + friend class BlockMatrixBase; }; } diff --git a/deal.II/lac/include/lac/sparse_matrix.h b/deal.II/lac/include/lac/sparse_matrix.h index f80ac00f6d..1bb763cf7b 100644 --- a/deal.II/lac/include/lac/sparse_matrix.h +++ b/deal.II/lac/include/lac/sparse_matrix.h @@ -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 friend class SparseMatrix; template friend class SparseLUDecomposition; template friend class SparseILU; + + /** + * To allow it calling private + * prepare_add() and + * prepare_set(). + */ + template friend class BlockMatrixBase; }; /*@}*/ @@ -2862,6 +2904,28 @@ SparseMatrix::end (const unsigned int r) return end(); } + + +template +inline +void +SparseMatrix:: +prepare_add() +{ + //nothing to do here +} + + + +template +inline +void +SparseMatrix:: +prepare_set() +{ + //nothing to do here +} + #endif // DOXYGEN diff --git a/deal.II/lac/include/lac/trilinos_sparse_matrix.h b/deal.II/lac/include/lac/trilinos_sparse_matrix.h index 6f96a1378e..9a220c1072 100755 --- a/deal.II/lac/include/lac/trilinos_sparse_matrix.h +++ b/deal.II/lac/include/lac/trilinos_sparse_matrix.h @@ -40,7 +40,10 @@ DEAL_II_NAMESPACE_OPEN -// forward declaration + // forward declarations +template +class BlockMatrixBase; + template 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 matrix; + /** + * To allow calling protected prepare_add() + * and prepare_set(). + */ + friend class BlockMatrixBase; }; @@ -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 }