From 94a0c5083e35707c3e20fd7e11c6fd0007068f33 Mon Sep 17 00:00:00 2001 From: kronbichler Date: Fri, 22 Aug 2008 16:17:11 +0000 Subject: [PATCH] Added a block sparse matrix class based on Trilinos (parallel) matrices. git-svn-id: https://svn.dealii.org/trunk@16653 0785d39b-7218-0410-832d-ea1e28bc413d --- .../lac/trilinos_block_sparse_matrix.h | 407 ++++++++++++++++++ .../source/trilinos_block_sparse_matrix.cc | 108 +++++ deal.II/lac/source/trilinos_sparse_matrix.cc | 18 + 3 files changed, 533 insertions(+) create mode 100644 deal.II/lac/include/lac/trilinos_block_sparse_matrix.h create mode 100644 deal.II/lac/source/trilinos_block_sparse_matrix.cc diff --git a/deal.II/lac/include/lac/trilinos_block_sparse_matrix.h b/deal.II/lac/include/lac/trilinos_block_sparse_matrix.h new file mode 100644 index 0000000000..c9a5a81552 --- /dev/null +++ b/deal.II/lac/include/lac/trilinos_block_sparse_matrix.h @@ -0,0 +1,407 @@ +//--------------------------------------------------------------------------- +// $Id: trilinos_block_sparse_matrix.h 14783 2007-06-18 14:52:01Z bangerth $ +// Version: $Name$ +// +// Copyright (C) 2004, 2005, 2006, 2007 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//--------------------------------------------------------------------------- +#ifndef __deal2__trilinos_block_sparse_matrix_h +#define __deal2__trilinos_block_sparse_matrix_h + + +#include +#include +#include +#include +#include +#include + +#include + + +#ifdef DEAL_II_USE_TRILINOS + +DEAL_II_NAMESPACE_OPEN + + + +namespace TrilinosWrappers +{ + +/*! @addtogroup TrilinosWrappers + *@{ + */ + +/** + * Blocked sparse matrix based on the TrilinosWrappers::SparseMatrix class. This + * class implements the functions that are specific to the Trilinos SparseMatrix + * base objects for a blocked sparse matrix, and leaves the actual work + * relaying most of the calls to the individual blocks to the functions + * implemented in the base class. See there also for a description of when + * this class is useful. + * + * In contrast to the deal.II-type SparseMatrix class, the Trilinos matrices do + * not have external objects for the sparsity patterns. Thus, one does not + * determine the size of the individual blocks of a block matrix of this type + * by attaching a block sparsity pattern, but by calling reinit() to set the + * number of blocks and then by setting the size of each block separately. In + * order to fix the data structures of the block matrix, it is then necessary + * to let it know that we have changed the sizes of the underlying + * matrices. For this, one has to call the collect_sizes() function, for much + * the same reason as is documented with the BlockSparsityPattern class. + * + * @ingroup Matrix1 + * @author Wolfgang Bangerth, 2004 + */ + class BlockSparseMatrix : public BlockMatrixBase + { + public: + /** + * Typedef the base class for simpler + * access to its own typedefs. + */ + typedef BlockMatrixBase BaseClass; + + /** + * Typedef the type of the underlying + * matrix. + */ + typedef BaseClass::BlockType BlockType; + + /** + * Import the typedefs from the base + * class. + */ + typedef BaseClass::value_type value_type; + typedef BaseClass::pointer pointer; + typedef BaseClass::const_pointer const_pointer; + typedef BaseClass::reference reference; + typedef BaseClass::const_reference const_reference; + typedef BaseClass::size_type size_type; + typedef BaseClass::iterator iterator; + typedef BaseClass::const_iterator const_iterator; + + /** + * Constructor; initializes the + * matrix to be empty, without + * any structure, i.e. the + * matrix is not usable at + * all. This constructor is + * therefore only useful for + * matrices which are members of + * a class. All other matrices + * should be created at a point + * in the data flow where all + * necessary information is + * available. + * + * You have to initialize the + * matrix before usage with + * reinit(BlockSparsityPattern). The + * number of blocks per row and + * column are then determined by + * that function. + */ + BlockSparseMatrix (); + + /** + * Destructor. + */ + ~BlockSparseMatrix (); + + /** + * Pseudo copy operator only copying + * empty objects. The sizes of the block + * matrices need to be the same. + */ + BlockSparseMatrix & + operator = (const BlockSparseMatrix &); + + /** + * 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 matrix=0, which + * sets all elements of the matrix to + * zero, but keep the sparsity pattern + * previously used. + */ + BlockSparseMatrix & + operator = (const double d); + + /** + * Resize the matrix, by setting + * the number of block rows and + * columns. This deletes all + * blocks and replaces them by + * unitialized ones, i.e. ones + * for which also the sizes are + * not yet set. You have to do + * that by calling the @p reinit + * functions of the blocks + * themselves. Do not forget to + * call collect_sizes() after + * that on this object. + * + * The reason that you have to + * set sizes of the blocks + * yourself is that the sizes may + * be varying, the maximum number + * of elements per row may be + * varying, etc. It is simpler + * not to reproduce the interface + * of the @p SparsityPattern + * class here but rather let the + * user call whatever function + * she desires. + */ + void reinit (const unsigned int n_block_rows, + const unsigned int n_block_columns); + + /** + * Resize the matrix, by using an + * array of Epetra maps to determine + * the distribution of the + * individual matrices. This function + * assumes that a quadratic block + * matrix is generated. + */ + //void reinit (const std::vector &input_maps); + + /** + * This function calls the compress() + * command of all matrices after + * the sparsity pattern has been + * generated or the assembly is + * completed. Note that all MPI + * processes need to call this + * command (in contrast to sparsity + * pattern generation, which is + * probably only called on one + * single process) before any + * can complete it. + */ + void compress (); + + /** + * This function collects the + * sizes of the sub-objects and + * stores them in internal + * arrays, in order to be able to + * relay global indices into the + * matrix to indices into the + * subobjects. You *must* call + * this function each time after + * you have changed the size of + * the sub-objects. + */ + void collect_sizes (); + + /** + * Matrix-vector multiplication: + * let $dst = M*src$ with $M$ + * being this matrix. + */ + void vmult (BlockVector &dst, + const BlockVector &src) const; + + /** + * Matrix-vector + * multiplication. Just like the + * previous function, but only + * applicable if the matrix has + * only one block column. + */ + void vmult (BlockVector &dst, + const Vector &src) const; + + /** + * Matrix-vector + * multiplication. Just like the + * previous function, but only + * applicable if the matrix has + * only one block row. + */ + void vmult (Vector &dst, + const BlockVector &src) const; + + /** + * Matrix-vector + * multiplication. Just like the + * previous function, but only + * applicable if the matrix has + * only one block. + */ + void vmult (Vector &dst, + const Vector &src) const; + + /** + * Matrix-vector multiplication: + * let $dst = M^T*src$ with $M$ + * being this matrix. This + * function does the same as + * vmult() but takes the + * transposed matrix. + */ + void Tvmult (BlockVector &dst, + const BlockVector &src) const; + + /** + * Matrix-vector + * multiplication. Just like the + * previous function, but only + * applicable if the matrix has + * only one block row. + */ + void Tvmult (BlockVector &dst, + const Vector &src) const; + + /** + * Matrix-vector + * multiplication. Just like the + * previous function, but only + * applicable if the matrix has + * only one block column. + */ + void Tvmult (Vector &dst, + const BlockVector &src) const; + + /** + * Matrix-vector + * multiplication. Just like the + * previous function, but only + * applicable if the matrix has + * only one block. + */ + void Tvmult (Vector &dst, + const Vector &src) const; + + /** + * Make the clear() function in the + * base class visible, though it is + * protected. + */ + using BlockMatrixBase::clear; + + /** @addtogroup Exceptions + * @{ + */ + + /** + * Exception + */ + DeclException4 (ExcIncompatibleRowNumbers, + int, int, int, int, + << "The blocks [" << arg1 << ',' << arg2 << "] and [" + << arg3 << ',' << arg4 << "] have differing row numbers."); + /** + * Exception + */ + DeclException4 (ExcIncompatibleColNumbers, + int, int, int, int, + << "The blocks [" << arg1 << ',' << arg2 << "] and [" + << arg3 << ',' << arg4 << "] have differing column numbers."); + ///@} + }; + + + +/*@}*/ + +// ------------- inline and template functions ----------------- + + inline + void + BlockSparseMatrix::vmult (BlockVector &dst, + const BlockVector &src) const + { + BaseClass::vmult_block_block (dst, src); + } + + + + inline + void + BlockSparseMatrix::vmult (BlockVector &dst, + const Vector &src) const + { + BaseClass::vmult_block_nonblock (dst, src); + } + + + + inline + void + BlockSparseMatrix::vmult (Vector &dst, + const BlockVector &src) const + { + BaseClass::vmult_nonblock_block (dst, src); + } + + + + inline + void + BlockSparseMatrix::vmult (Vector &dst, + const Vector &src) const + { + BaseClass::vmult_nonblock_nonblock (dst, src); + } + + + inline + void + BlockSparseMatrix::Tvmult (BlockVector &dst, + const BlockVector &src) const + { + BaseClass::Tvmult_block_block (dst, src); + } + + + + inline + void + BlockSparseMatrix::Tvmult (BlockVector &dst, + const Vector &src) const + { + BaseClass::Tvmult_block_nonblock (dst, src); + } + + + + inline + void + BlockSparseMatrix::Tvmult (Vector &dst, + const BlockVector &src) const + { + BaseClass::Tvmult_nonblock_block (dst, src); + } + + + + inline + void + BlockSparseMatrix::Tvmult (Vector &dst, + const Vector &src) const + { + BaseClass::Tvmult_nonblock_nonblock (dst, src); + } + +} + +DEAL_II_NAMESPACE_CLOSE + +#endif // DEAL_II_USE_TRILINOS + +#endif // __deal2__trilinos_block_sparse_matrix_h diff --git a/deal.II/lac/source/trilinos_block_sparse_matrix.cc b/deal.II/lac/source/trilinos_block_sparse_matrix.cc new file mode 100644 index 0000000000..314c8a158d --- /dev/null +++ b/deal.II/lac/source/trilinos_block_sparse_matrix.cc @@ -0,0 +1,108 @@ +//--------------------------------------------------------------------------- +// $Id: trilinos_block_sparse_matrix.cc 15631 2008-01-17 23:47:31Z bangerth $ +// Version: $Name$ +// +// Copyright (C) 2004, 2005, 2006, 2008 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//--------------------------------------------------------------------------- + +#include + + +#ifdef DEAL_II_USE_TRILINOS + +DEAL_II_NAMESPACE_OPEN + +namespace TrilinosWrappers +{ + BlockSparseMatrix::BlockSparseMatrix () + {} + + + + BlockSparseMatrix::~BlockSparseMatrix () + {} + + + + BlockSparseMatrix & + BlockSparseMatrix::operator = (const BlockSparseMatrix &m) + { + BaseClass::operator = (m); + + return *this; + } + + + + BlockSparseMatrix & + BlockSparseMatrix::operator = (const double d) + { + for (unsigned int r=0; rn_block_rows(); ++r) + for (unsigned int c=0; cn_block_cols(); ++c) + this->block(r,c) = d; + + return *this; + } + + + void + BlockSparseMatrix:: + reinit (const unsigned int n_block_rows, + const unsigned int n_block_columns) + { + // first delete previous content of + // the subobjects array + for (unsigned int r=0; rn_block_rows(); ++r) + for (unsigned int c=0; cn_block_cols(); ++c) + { + BlockType *p = this->sub_objects[r][c]; + this->sub_objects[r][c] = 0; + delete p; + } + + this->sub_objects.reinit (0,0); + + // then resize. set sizes of blocks to + // zero. user will later have to call + // collect_sizes for this + this->sub_objects.reinit (n_block_rows, + n_block_columns); + this->row_block_indices.reinit (n_block_rows, 0); + this->column_block_indices.reinit (n_block_columns, 0); + + // and reinitialize the blocks + for (unsigned int r=0; rn_block_rows(); ++r) + for (unsigned int c=0; cn_block_cols(); ++c) + { + BlockType *p = new BlockType(); + this->sub_objects[r][c] = p; + } + } + + + void + BlockSparseMatrix::compress() + { + for (unsigned int r=0; rn_block_rows(); ++r) + for (unsigned int c=0; cn_block_cols(); ++c) + this->block(r,c).compress(); + } + + void + BlockSparseMatrix::collect_sizes () + { + BaseClass::collect_sizes (); + } + +} + + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/deal.II/lac/source/trilinos_sparse_matrix.cc b/deal.II/lac/source/trilinos_sparse_matrix.cc index 6e3d0afa03..f8f8872550 100755 --- a/deal.II/lac/source/trilinos_sparse_matrix.cc +++ b/deal.II/lac/source/trilinos_sparse_matrix.cc @@ -194,6 +194,14 @@ namespace TrilinosWrappers Assert (matrix->NumGlobalRows() == (int)sparsity_pattern.n_rows(), ExcDimensionMismatch (matrix->NumGlobalRows(), sparsity_pattern.n_rows())); + + // Trilinos seems to have a bug for + // rectangular matrices at this point, + // so do not check for consistent + // column numbers here. +// Assert (matrix->NumGlobalCols() == (int)sparsity_pattern.n_cols(), +// ExcDimensionMismatch (matrix->NumGlobalCols(), +// sparsity_pattern.n_cols())); std::vector n_entries_per_row(n_rows); @@ -216,6 +224,7 @@ namespace TrilinosWrappers unsigned int n_rows = sparsity_pattern.n_rows(); matrix.reset(); row_map = input_map; + col_map = row_map; Assert (input_map.NumGlobalElements() == (int)sparsity_pattern.n_rows(), ExcDimensionMismatch (input_map.NumGlobalElements(), @@ -266,6 +275,15 @@ namespace TrilinosWrappers matrix = std::auto_ptr (new Epetra_FECrsMatrix(Copy, row_map, col_map, &n_entries_per_row[0], true)); + +// Assert (matrix->NumGlobalCols() == col_map.NumGlobalElements(), +// ExcInternalError()); +// Assert (matrix->NumGlobalRows() == row_map.NumGlobalElements(), +// ExcInternalError()); +// Assert (matrix->NumGlobalCols() == sparsity_pattern.n_cols(), +// ExcInternalError()); +// Assert (matrix->NumGlobalRows() == sparsity_pattern.n_rows(), +// ExcInternalError()); const unsigned int n_max_entries_per_row = *std::max_element ( &n_entries_per_row[0], &n_entries_per_row[n_rows-1]); -- 2.39.5