From: wolf Date: Fri, 5 May 2000 15:36:13 +0000 (+0000) Subject: Bring this into the tree, albeit untested, for possible work at home. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=464ccb87ace4f4667982cbab6631e92a9c9ca1e1;p=dealii-svn.git Bring this into the tree, albeit untested, for possible work at home. git-svn-id: https://svn.dealii.org/trunk@2804 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/lac/include/lac/block_sparse_matrix.h b/deal.II/lac/include/lac/block_sparse_matrix.h new file mode 100644 index 0000000000..201f37add1 --- /dev/null +++ b/deal.II/lac/include/lac/block_sparse_matrix.h @@ -0,0 +1,447 @@ +//---------------------------- block_sparse_matrix.h --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 1998, 1999, 2000 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. +// +//---------------------------- block_sparse_matrix.h --------------------------- +#ifndef __deal2__block_sparse_matrix_h +#define __deal2__block_sparse_matrix_h + + +#include +#include + + + +template +class BlockSparseMatrix : public Subscriptor +{ + public: + /** + * Type of matrix entries. In analogy to + * the STL container classes. + */ + typedef number value_type; + + /** + * 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)#. + */ + BlockSparseMatrix (); + + /** + * Constructor. Takes the given + * matrix sparsity structure to + * represent the sparsity pattern + * of this matrix. You can change + * the sparsity pattern later on + * by calling the #reinit# + * function. + * + * This constructor initializes + * all sub-matrices with the + * sub-sparsity pattern within + * the argument. + * + * You have to make sure that the + * lifetime of the sparsity + * structure is at least as long + * as that of this matrix or as + * long as #reinit# is not called + * with a new sparsity structure. + */ + BlockSparseMatrix (const BlockSparsityPattern &sparsity); + + + + /** + * Pseudo operator only copying + * empty objects. + */ + BlockSparseMatrix & + operator = (const BlockSparseMatrix &); + + + /** + * Reinitialize the object but + * keep to the sparsity pattern + * previously used. This may be + * necessary if you #reinit#'d + * the sparsity structure and + * want to update the size of the + * matrix. It only calls #reinit# + * on the sub-matrices. + */ + virtual void reinit (); + + /** + * Reinitialize the sparse matrix + * with the given sparsity + * pattern. The latter tells the + * matrix how many nonzero + * elements there need to be + * reserved. + * + * Basically, this function only + * calls #reinit# of the + * sub-matrices with the block + * sparsity patterns of the + * parameter. + */ + virtual void reinit (const BlockSparsityPattern &sparsity); + + /** + * Release all memory and return + * to a state just like after + * having called the default + * constructor. It also forgets + * the sparsity pattern it was + * previously tied to. + * + * This calls #clear# on all + * sub-matrices. + */ + virtual void clear (); + + /** + * Return whether the object is + * empty. It is empty if either + * both dimensions are zero or no + * #SparsityPattern# is + * associated. + */ + bool empty () 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; + + /** + * Return the number of nonzero + * elements of this + * matrix. Actually, it returns + * the number of entries in the + * sparsity pattern; if any of + * the entries should happen to + * be zero, it is counted anyway. + */ + unsigned int n_nonzero_elements () const; + + /** + * Set the element #(i,j)# to #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 #value# to the element + * #(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); + + /** + * Copy the given matrix to this + * one. The operation throws an + * error if the sparsity patterns + * of the two involved matrices + * do not point to the same + * object, since in this case the + * copy operation is + * cheaper. Since this operation + * is notheless not for free, we + * do not make it available + * through #operator =#, since + * this may lead to unwanted + * usage, e.g. in copy arguments + * to functions, which should + * really be arguments by + * reference. + * + * The source matrix may be a + * matrix of arbitrary type, as + * long as its data type is + * convertible to the data type + * of this matrix. + * + * The function returns a + * reference to #this#. + */ + template + BlockSparseMatrix & + copy_from (const BlockSparseMatrix &source); + + /** + * Add #matrix# scaled by + * #factor# to this matrix. The + * function throws an error if + * the sparsity patterns of the + * two involved matrices do not + * point to the same object, + * since in this case the + * operation is cheaper. + * + * The source matrix may be a matrix + * of arbitrary type, as long as its + * data type is convertible to the + * data type of this matrix. + */ + template + void add_scaled (const number factor, + const BlockSparseMatrix &matrix); + + /** + * Return the value of the entry (i,j). + * This may be an expensive operation + * and you should always take care + * where to call this function. + * In order to avoid abuse, this function + * throws an exception if the wanted + * element does not exist in the matrix. + */ + number operator () (const unsigned int i, const unsigned int j) const; + + private: + /** + * Pointer to the block sparsity + * pattern used for this + * matrix. In order to guarantee + * that it is not deleted while + * still in use, we subscribe to + * it using the #SmartPointer# + * class. + */ + SmartPointer > sparsity_pattern; + + /** + * Array of sub-matrices. + */ + SparseMatrix sub_objects[rows][columns]; +}; + + + + +/* ------------------------- Template functions ---------------------- */ + + +template +BlockSparseMatrix::BlockSparseMatrix () : + sparsity_pattern (0) +{}; + + + +template +BlockSparseMatrix:: +BlockSparseMatrix (const BlockSparsityPattern &sparsity) +{ + reinit (sparsity); +}; + + + +template +BlockSparseMatrix & +BlockSparseMatrix:: +operator = (const BlockSparseMatrix &m) +{ + *sparsity_pattern = *m.sparsity_pattern; + for (unsigned int r=0; r +void +BlockSparseMatrix::reinit () +{ + for (unsigned int r=0; r +void +BlockSparseMatrix::reinit (const BlockSparsityPattern &sparsity) +{ + sparsity_pattern = &sparsity; + + for (unsigned int r=0; r +void +BlockSparseMatrix::clear () +{ + sparsity_pattern = 0; + for (unsigned int r=0; r +bool +BlockSparseMatrix::empty () const +{ + for (unsigned int r=0; r +unsigned int +BlockSparseMatrix::m () const +{ + return sparsity_pattern->n_rows(); +}; + + + +template +unsigned int +BlockSparseMatrix::n () const +{ + return sparsity_pattern->n_cols(); +}; + + + +template +unsigned int +BlockSparseMatrix::n_nonzero_elements () const +{ + return sparsity_pattern->n_nonzero_elements (); +}; + + + +template +void +BlockSparseMatrix::set (const unsigned int i, + const unsigned int j, + const number value) +{ + const pair + row_index = sparsity_pattern.row_indices.global_to_local (i), + col_index = sparsity_pattern.column_indices.global_to_local (j); + sub_objects[row_index.first][col_index.first].set (row_index.second, + col_index.second, + value); +}; + + + +template +void +BlockSparseMatrix::add (const unsigned int i, + const unsigned int j, + const number value) +{ + const pair + row_index = sparsity_pattern.row_indices.global_to_local (i), + col_index = sparsity_pattern.column_indices.global_to_local (j); + sub_objects[row_index.first][col_index.first].add (row_index.second, + col_index.second, + value); +}; + + + +template +template +BlockSparseMatrix & +BlockSparseMatrix:: +copy_from (const BlockSparseMatrix &source) +{ + for (unsigned int r=0; r +template +void +BlockSparseMatrix:: +add_scaled (const number factor, + const BlockSparseMatrix &matrix) +{ + for (unsigned int r=0; r +number +BlockSparseMatrix::operator () (const unsigned int i, + const unsigned int j) const +{ + const pair + row_index = sparsity_pattern.row_indices.global_to_local (i), + col_index = sparsity_pattern.column_indices.global_to_local (j); + return sub_objects[row_index.first][col_index.first] (row_index.second, + col_index.second); +}; + + + + +#endif // __deal2__block_sparse_matrix_h