<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>
<br>
(GK 2009/08/04)
</p>
+ </li>
<li>
<p>
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
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)"));
const number *values,
const bool elide_zero_values)
{
+ prepare_set_operation();
+
// Resize scratch arrays
if (column_indices.size() < this->n_block_cols())
{
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
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.
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
// $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
*/
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
*/
std::vector<PetscScalar> column_values;
+
+ /**
+ * To allow calling protected prepare_add()
+ * and prepare_set().
+ */
+ template <class> friend class BlockMatrixBase;
+
};
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;
// 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;
}
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;
(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
}
// $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
{
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>;
};
// $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
#ifdef DEAL_II_USE_PETSC
DEAL_II_NAMESPACE_OPEN
+ // forward declaration
+template <typename MatrixType>
+class BlockMatrixBase;
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>;
};
}
* 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
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;
};
/*@}*/
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
DEAL_II_NAMESPACE_OPEN
-// forward declaration
+ // forward declarations
+template <typename MatrixType>
+class BlockMatrixBase;
+
template <typename number> class SparseMatrix;
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
*/
std::auto_ptr<Epetra_FECrsMatrix> matrix;
+ /**
+ * To allow calling protected prepare_add()
+ * and prepare_set().
+ */
+ friend class BlockMatrixBase<SparseMatrix>;
};
return matrix->ColMap();
}
+
+
+ inline
+ void
+ SparseMatrix::prepare_add()
+ {
+ //nothing to do here
+ }
+
+
+
+ inline
+ void
+ SparseMatrix::prepare_set()
+ {
+ //nothing to do here
+ }
+
+
+
#endif // DOXYGEN
}