--- /dev/null
+//---------------------------------------------------------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 1998, 1999, 2000, 2001, 2002 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__block_sparse_matrix_ez_h
+#define __deal2__block_sparse_matrix_ez_h
+
+
+#include <base/config.h>
+#include <base/exceptions.h>
+#include <base/subscriptor.h>
+#include <base/smartpointer.h>
+#include <lac/block_indices.h>
+
+template <typename Number> class BlockVector;
+
+/**
+ * A block matrix consisting of blocks of type @p{SparseMatrixEZ}.
+ *
+ * Like the other Block-objects, this matrix can be used like a
+ * @p{SparseMatrixEZ}, when it comes to access to entries. Then, there
+ * are functions for the multiplication with @p{BlockVector} and
+ * access to the individual blocks.
+ *
+ * @author Guido Kanschat, 2002
+ */
+template<typename Number>
+class BlockSparseMatrixEZ : public Subscriptor
+{
+ public:
+ /**
+ * Default constructor. The
+ * result is an empty object with
+ * zero dimensions.
+ */
+ BlockSparseMatrixEZ ();
+
+ /**
+ * Constructor setting up an
+ * object with given unmber of
+ * block rows and columns. The
+ * blocks themselves still have
+ * zero dimension.
+ */
+ BlockSparseMatrixEZ (const unsigned int block_rows,
+ const unsigned int block_cols);
+
+ /**
+ * Copy constructor. This is
+ * needed for some container
+ * classes. It creates an object
+ * of the same number of block
+ * rows and columns. Since it
+ * calls the copy constructor of
+ * @ref{SparseMatrixEZ}, the
+ * block s must be empty.
+ */
+ BlockSparseMatrixEZ (const BlockSparseMatrixEZ<Number>&);
+
+ /**
+ * Copy operator. Like the copy
+ * constructor, this may be
+ * called for objects with empty
+ * blocks only.
+ */
+ BlockSparseMatrixEZ operator = (const BlockSparseMatrixEZ<Number>&);
+
+ /**
+ * Initialize to given block
+ * numbers. After this
+ * operation, the matrix will
+ * have the block dimensions
+ * provided. Each block will have
+ * zero dimensions and must be
+ * initialized
+ * subsequently. After setting
+ * the sizes of the blocks,
+ * @ref{collect_sizes} must be
+ * called to update internal data
+ * structures.
+ */
+ void reinit (const unsigned int n_block_rows,
+ const unsigned int n_block_cols);
+ /**
+ * 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 ();
+
+ /**
+ * Access the block with the
+ * given coordinates.
+ */
+ SparseMatrixEZ<Number>&
+ block (const unsigned int row,
+ const unsigned int column);
+
+
+ /**
+ * Access the block with the
+ * given coordinates. Version for
+ * constant objects.
+ */
+ const SparseMatrixEZ<Number>&
+ block (const unsigned int row,
+ const unsigned int column) const;
+
+ /**
+ * Return the number of blocks in a
+ * column.
+ */
+ unsigned int n_block_rows () const;
+
+ /**
+ * Return the number of blocks in a
+ * row.
+ */
+ unsigned int n_block_cols () const;
+
+ /**
+ * Return whether the object is
+ * empty. It is empty if no
+ * memory is allocated, which is
+ * the same as that both
+ * dimensions are zero. This
+ * function is just the
+ * concatenation of the
+ * respective call to all
+ * sub-matrices.
+ */
+ bool empty () const;
+
+ /**
+ * Return number of rows of this
+ * matrix, which equals the
+ * dimension of the image
+ * space. It is the sum of rows
+ * of the rows of sub-matrices.
+ */
+ unsigned int n_rows () const;
+
+ /**
+ * Return number of columns of
+ * this matrix, which equals the
+ * dimension of the range
+ * space. It is the sum of
+ * columns of the columns of
+ * sub-matrices.
+ */
+ unsigned int n_cols () const;
+
+ /**
+ * Return the dimension of the
+ * image space. To remember: the
+ * matrix is of dimension
+ * $m \times n$.
+ */
+ unsigned int m () const;
+
+ /**
+ * Return the dimension of the
+ * range space. To remember: the
+ * matrix is of dimension
+ * $m \times n$.
+ */
+ unsigned int n () const;
+
+ /**
+ * Set the element @p{(i,j)} to
+ * @p{value}. Throws an error if
+ * the entry does not
+ * exist. Still, it is allowed to
+ * store zero values in
+ * non-existent fields.
+ */
+ void set (const unsigned int i,
+ const unsigned int j,
+ const number value);
+
+ /**
+ * Add @p{value} to the element
+ * @p{(i,j)}. Throws an error if
+ * the entry does not
+ * exist. Still, it is allowed to
+ * store zero values in
+ * non-existent fields.
+ */
+ void add (const unsigned int i, const unsigned int j,
+ const number value);
+
+
+ /**
+ * Matrix-vector multiplication:
+ * let $dst = M*src$ with $M$
+ * being this matrix.
+ */
+ template <typename somenumber>
+ void vmult (BlockVector<somenumber> &dst,
+ const BlockVector<somenumber> &src) const;
+
+ /**
+ * Matrix-vector multiplication:
+ * let $dst = M^T*src$ with $M$
+ * being this matrix. This
+ * function does the same as
+ * @p{vmult} but takes the
+ * transposed matrix.
+ */
+ template <typename somenumber>
+ void Tvmult (BlockVector<somenumber> &dst,
+ const BlockVector<somenumber> &src) const;
+
+ /**
+ * Adding Matrix-vector
+ * multiplication. Add $M*src$ on
+ * $dst$ with $M$ being this
+ * matrix.
+ */
+ template <typename somenumber>
+ void vmult_add (BlockVector<somenumber> &dst,
+ const BlockVector<somenumber> &src) const;
+
+ /**
+ * Adding Matrix-vector
+ * multiplication. Add $M^T*src$
+ * to $dst$ with $M$ being this
+ * matrix. This function does the
+ * same as @p{vmult_add} but takes
+ * the transposed matrix.
+ */
+ template <typename somenumber>
+ void Tvmult_add (BlockVector<somenumber> &dst,
+ const BlockVector<somenumber> &src) const;
+
+
+ private:
+ /**
+ * Object storing and managing
+ * the transformation of row
+ * indices to indices of the
+ * sub-objects.
+ */
+ BlockIndices row_indices;
+
+ /**
+ * Object storing and managing
+ * the transformation of column
+ * indices to indices of the
+ * sub-objects.
+ */
+ BlockIndices column_indices;
+
+ /**
+ * The actual matrices
+ */
+ Table<2, SparseMatrixEZ<Number> > blocks;
+};
+
+
+/*----------------------------------------------------------------------*/
+
+
+template <typename number>
+inline
+unsigned int
+BlockSparseMatrixEZ<number>::n_block_rows () const
+{
+ return row_indices.size()
+}
+
+
+
+template <typename number>
+inline
+unsigned int
+BlockSparseMatrixEZ<number>::n_rows () const
+{
+ return row_indices.total_size()
+}
+
+
+
+template <typename number>
+inline
+unsigned int
+BlockSparseMatrixEZ<number>::n_block_cols () const
+{
+ return column_indices.size()
+}
+
+
+
+template <typename number>
+inline
+unsigned int
+BlockSparseMatrixEZ<number>::n_cols () const
+{
+ return column_indices.total_size()
+}
+
+
+
+template <typename number>
+inline
+SparseMatrixEZ<number> &
+BlockSparseMatrixEZ<number>::block (const unsigned int row,
+ const unsigned int column)
+{
+ Assert (row<n_rows(), ExcIndexRange (row, 0, n_rows()));
+ Assert (column<n_cols(), ExcIndexRange (column, 0, n_cols()));
+
+ return blocks[row][column];
+};
+
+
+
+template <typename number>
+inline
+const SparseMatrixEZ<number> &
+BlockSparseMatrixEZ<number>::block (const unsigned int row,
+ const unsigned int column) const
+{
+ Assert (row<n_rows(), ExcIndexRange (row, 0, n_rows()));
+ Assert (column<n_cols(), ExcIndexRange (column, 0, n_cols()));
+
+ return blocks[row][column];
+};
+
+
+
+template <typename number>
+inline
+unsigned int
+BlockSparseMatrixEZ<number>::m () const
+{
+ return n_rows();
+};
+
+
+
+template <typename number>
+inline
+unsigned int
+BlockSparseMatrixEZ<number>::n () const
+{
+ return n_cols();
+};
+
+
+
+template <typename number>
+inline
+void
+BlockSparseMatrixEZ<number>::set (const unsigned int i,
+ const unsigned int j,
+ const number value)
+{
+ const std::pair<unsigned int,unsigned int>
+ row_index = row_indices.global_to_local (i),
+ col_index = column_indices.global_to_local (j);
+ block(row_index.first,col_index.first).set (row_index.second,
+ col_index.second,
+ value);
+};
+
+
+
+template <typename number>
+inline
+void
+BlockSparseMatrixEZ<number>::add (const unsigned int i,
+ const unsigned int j,
+ const number value)
+{
+ const std::pair<unsigned int,unsigned int>
+ row_index = row_indices.global_to_local (i),
+ col_index = column_indices.global_to_local (j);
+ block(row_index.first,col_index.first).add (row_index.second,
+ col_index.second,
+ value);
+};
+
+
+template <typename number>
+template <typename somenumber>
+void
+BlockSparseMatrixEZ<number>::vmult (BlockVector<somenumber> &dst,
+ const BlockVector<somenumber> &src) const
+{
+ Assert (dst.n_blocks() == n_block_rows(),
+ ExcDimensionMismatch(dst.n_blocks(), n_block_rows()));
+ Assert (src.n_blocks() == n_block_cols(),
+ ExcDimensionMismatch(src.n_blocks(), n_block_cols()));
+
+ dst = 0.;
+
+ for (unsigned int row=0; row<n_block_rows(); ++row)
+ {
+ for (unsigned int col=0; col<n_block_cols(); ++col)
+ block(row,col).vmult_add (dst.block(row),
+ src.block(col));
+ };
+};
+
+
+
+template <typename number>
+template <typename somenumber>
+void
+BlockSparseMatrix<number>::vmult_add (BlockVector<somenumber> &dst,
+ const BlockVector<somenumber> &src) const
+{
+ Assert (dst.n_blocks() == n_block_rows(),
+ ExcDimensionMismatch(dst.n_blocks(), n_block_rows()));
+ Assert (src.n_blocks() == n_block_cols(),
+ ExcDimensionMismatch(src.n_blocks(), n_block_cols()));
+
+ for (unsigned int row=0; row<n_block_rows(); ++row)
+ {
+ for (unsigned int col=0; col<n_block_cols(); ++col)
+ block(row,col).vmult_add (dst.block(row),
+ src.block(col));
+ };
+};
+
+
+
+
+template <typename number>
+template <typename somenumber>
+void
+BlockSparseMatrix<number>::Tvmult (BlockVector<somenumber> &dst,
+ const BlockVector<somenumber>& src) const
+{
+ Assert (dst.n_blocks() == n_block_cols(),
+ ExcDimensionMismatch(dst.n_blocks(), n_block_cols()));
+ Assert (src.n_blocks() == n_block_rows(),
+ ExcDimensionMismatch(src.n_blocks(), n_block_rows()));
+
+ dst = 0.;
+
+ for (unsigned int row=0; row<n_block_rows(); ++row)
+ {
+ for (unsigned int col=0; col<n_block_cols(); ++col)
+ block(row,col).Tvmult_add (dst.block(col),
+ src.block(row));
+ };
+};
+
+
+
+template <typename number>
+template <typename somenumber>
+void
+BlockSparseMatrix<number>::Tvmult_add (BlockVector<somenumber> &dst,
+ const BlockVector<somenumber> &src) const
+{
+ Assert (dst.n_blocks() == n_block_cols(),
+ ExcDimensionMismatch(dst.n_blocks(), n_block_cols()));
+ Assert (src.n_blocks() == n_block_rows(),
+ ExcDimensionMismatch(src.n_blocks(), n_block_rows()));
+
+ for (unsigned int row=0; row<n_block_rows(); ++row)
+ {
+ for (unsigned int col=0; col<n_block_cols(); ++col)
+ block(row,col).Tvmult_add (dst.block(col),
+ src.block(row));
+ };
+};
--- /dev/null
+//-------------------------------------------------------------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 1998, 1999, 2000, 2001, 2002 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__block_sparse_matrix_ez_templates_h
+#define __deal2__block_sparse_matrix_ez_templates_h
+
+
+#include <base/config.h>
+#include <base/memory_consumption.h>
+#include <lac/block_sparse_matrix_ez.h>
+
+
+
+template <typename number>
+BlockSparseMatrixEZ<number>::BlockSparseMatrixEZ () :
+ row_indices (0),
+ column_indices (0)
+{};
+
+
+
+template <typename number>
+BlockSparseMatrixEZ<number>::
+BlockSparseMatrixEZ (unsigned int rows, unsigned int cols) :
+ row_indices (rows),
+ column_indices (cols)
+{};
+
+
+
+// template <typename number>
+// BlockSparseMatrixEZ<number>::~BlockSparseMatrixEZ ()
+// {
+// // delete previous content of
+// // the subobjects array
+// for (unsigned int r=0; r<rows; ++r)
+// for (unsigned int c=0; c<columns; ++c)
+// {
+// SparseMatrixEZ<number> *p = sub_objects[r][c];
+// sub_objects[r][c] = 0;
+// delete p;
+// };
+// };
+
+
+
+template <typename number>
+BlockSparseMatrixEZ<number> &
+BlockSparseMatrixEZ<number>::
+operator = (const BlockSparseMatrixEZ<number> &m)
+{
+ Assert (n_block_rows() == m.n_block_rows(), ExcIncompatibleObjects());
+ Assert (n_block_cols() == m.n_block_cols(), ExcIncompatibleObjects());
+ // this operator does not do
+ // anything except than checking
+ // whether the base objects want to
+ // do something
+ for (unsigned int r=0; r<n_block_rows(); ++r)
+ for (unsigned int c=0; c<n_block_cols(); ++c)
+ block(r,c) = m.block(r,c);
+ return *this;
+};
+
+
+
+template <typename number>
+BlockSparseMatrixEZ<number>::BlockSparseMatrixEZ (
+ const BlockSparseMatrixEZ<number> &m)
+ :
+ row_indices(m.row_indices),
+ column_indices(m.column_indices),
+ blocks(m.blocks)
+{};
+
+
+
+template <typename number>
+void
+BlockSparseMatrixEZ<number>::reinit (unsigned int rows,
+ unsigned int cols)
+{
+ row_indices.reinit(rows, 0);
+ column_indices.reinit(cols, 0);
+ blocks.reinit(rows, cols);
+};
+
+
+
+template <typename number>
+void
+BlockSparseMatrixEZ<number>::clear ()
+{
+ row_indices.reinit(0, 0);
+ column_indices.reinit(0, 0);
+ blocks.reinit(0, 0);
+};
+
+
+
+template <typename number>
+bool
+BlockSparseMatrixEZ<number>::empty () const
+{
+ for (unsigned int r=0; r<n_block_rows(); ++r)
+ for (unsigned int c=0; c<n_block_cols(); ++c)
+ if (block(r,c).empty () == false)
+ return false;
+ return true;
+};
+
+
+
+template <typename number>
+void
+BlockSparseMatrixEZ<number>::collect_sizes ()
+{
+ const unsigned int rows = n_block_rows();
+ const unsigned int columns = n_block_cols();
+ std::vector<unsigned int> row_sizes (rows);
+ std::vector<unsigned int> col_sizes (columns);
+
+ // first find out the row sizes
+ // from the first block column
+ for (unsigned int r=0; r<rows; ++r)
+ row_sizes[r] = blocks[r][0].n_rows();
+ // then check that the following
+ // block columns have the same
+ // sizes
+ for (unsigned int c=1; c<columns; ++c)
+ for (unsigned int r=0; r<rows; ++r)
+ Assert (row_sizes[r] == blocks[r][c].n_rows(),
+ ExcIncompatibleRowNumbers (r,0,r,c));
+
+ // finally initialize the row
+ // indices with this array
+ row_indices.reinit (row_sizes);
+
+
+ // then do the same with the columns
+ for (unsigned int c=0; c<columns; ++c)
+ col_sizes[c] = blocks[0][c].n_cols();
+ for (unsigned int r=1; r<rows; ++r)
+ for (unsigned int c=0; c<columns; ++c)
+ Assert (col_sizes[c] == blocks[r][c].n_cols(),
+ ExcIncompatibleRowNumbers (0,c,r,c));
+
+ // finally initialize the row
+ // indices with this array
+ column_indices.reinit (col_sizes);
+};
+
+
+
+
+// template <typename number>
+// unsigned int
+// BlockSparseMatrixEZ<number>::n_nonzero_elements () const
+// {
+// return sparsity_pattern->n_nonzero_elements ();
+// };
+
+
+
+// template <typename number>
+// unsigned int
+// BlockSparseMatrixEZ<number>::n_actually_nonzero_elements () const
+// {
+// unsigned int count = 0;
+// for (unsigned int i=0; i<rows; ++i)
+// for (unsigned int j=0; j<columns; ++j)
+// count += sub_objects[i][j]->n_actually_nonzero_elements ();
+// return count;
+// };
+
+
+
+// template <typename number>
+// unsigned int
+// BlockSparseMatrixEZ<number>::memory_consumption () const
+// {
+// unsigned int mem = sizeof(*this);
+// mem += MemoryConsumption::memory_consumption (sub_objects);
+// for (unsigned int r=0; r<rows; ++r)
+// for (unsigned int c=0; c<columns; ++c)
+// mem += MemoryConsumption::memory_consumption(*sub_objects[r][c]);
+// return mem;
+// };
+
+
+
+#endif // ifdef block_sparse_matrix_templates_h