From: hartmann Date: Fri, 26 Feb 1999 15:00:13 +0000 (+0000) Subject: replace the old dBlockSMatrix by the new template BlockSparseMatrix with a second... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=4c06b81f175261bd25ef834f1234778903ec695e;p=dealii-svn.git replace the old dBlockSMatrix by the new template BlockSparseMatrix with a second template parameter denoting the datatype the inverse diagonal block matrices are stored with. git-svn-id: https://svn.dealii.org/trunk@917 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/lac/include/lac/blocksparsematrix.h b/deal.II/lac/include/lac/blocksparsematrix.h new file mode 100644 index 0000000000..5aa47cce88 --- /dev/null +++ b/deal.II/lac/include/lac/blocksparsematrix.h @@ -0,0 +1,171 @@ +/*---------------------------- blocksparsematrix.h ---------------------*/ +/* $Id$ */ +/* Ralf Hartmann, University of Heidelberg */ +#ifndef __blocksparsematrix_H +#define __blocksparsematrix_H +/*---------------------------- blocksparsematrix.h ---------------------*/ + + +#include +#include +#include + +/** + * Block sparse matrix. + * The block matrix assumes the matrix consisting of blocks on + * the diagonal. These diagonal blocks and the elements + * (of arbitray structure) below the + * diagonal blocks are used in the #precondition_BlockSOR#. + * + * This block matrix structure is given e.g. for the DG method + * for the transport equation and a downstream numbering. + * If (as for this DG method) the matrix is empty above the + * diagonal blocks then BlockSOR is a direct solver. + * + * This first implementation of the BlockMatrix assumes the + * matrix having blocks each of the same block size. Varying + * block sizes within the matrix must still be implemented if needed. + * + * The first template parameter denotes the type of the matrix, the second + * denotes the type of data in which the inverse diagonal block + * matrices are stored by #invert_diagblocks()#. + * + * @author Ralf Hartmann, 1999 + */ +template +class BlockSparseMatrix: public SparseMatrix +{ + public: + /** + * Constructor + */ + BlockSparseMatrix(); + + /** + * Destructor + */ + virtual ~BlockSparseMatrix(); + + /** + * Call #dSMatrix::reinit()# and + * delete the inverse diagonal block matrices if existent. + */ + virtual void reinit(); + + /** + * Call #dSMatrix::reinit + * (const dSMatrixStruct &sparsity)# and + * delete the inverse diagonal block matrices if existent. + */ + virtual void reinit (const SparseMatrixStruct &sparsity); + + /** + * Call #dSMatrix::clear# and + * delete the inverse diagonal block matrices if existent. + */ + virtual void clear (); + + /** + * Stores the inverse of + * the diagonal blocks matrices + * in #inverse#. This costs some + * additional memory - for DG + * methods about 1/3 (for double inverses) + * or 1/6 (for float inverses) of that + * used for the matrix - but it + * makes the preconditioning much faster. + * + * It is not allowed to call this function + * twice (will produce an error) before + * a call of #reinit(..)# + * because at the second time there already + * exist the inverse matrices. + */ + void invert_diagblocks(); + + /** + * Block SOR. Make sure that the right + * block size + * of the matrix is set by + * #set_block_size# + * before calling this function. + * + * BlockSOR will automatically use the + * inverse matrices if they exist, if not + * then BlockSOR will waste much time + * inverting the diagonal block + * matrices in each preconditioning step. + * + * For matrices which are + * empty above the diagonal blocks + * BlockSOR is a direct solver. + */ + template + void precondition_BlockSOR (Vector &dst, const Vector &src, + const number om = 1.) const; + + /** + * Set the right block size before calling + * #precondition_BlockSOR#. + * If block_size==1 BlockSOR is the + * same as SOR. + */ + void set_block_size (const unsigned int bsize); + + /** + * Gives back the size of the blocks. + */ + unsigned int block_size () const; + + /** + * Exception + */ + DeclException2 (ExcWrongBlockSize, + int, int, + << "The blocksize " << arg1 + << " and the size of the matrix " << arg2 + << " do not match."); + + /** + * Exception + */ + DeclException2 (ExcWrongNumberOfInverses, + int, int, + << "There are " << arg1 + << " inverse matrices but " << arg2 + << " cells."); + + /** + * Exception + */ + DeclException0 (ExcInverseMatricesAlreadyExist); + + /** + * Exception + */ + DeclException0 (ExcBlockSizeNotSet); + + private: + /** + * Size of the blocks. Each diagonal + * block is assumed to be of the + * same size. + */ + unsigned int blocksize; + + /** + * Storage of the inverse matrices of + * the diagonal blocks matrices as + * #FullMatrix# matrices. + * For BlockSOR as preconditioning + * using #inverse_type=float# saves memory + * in comparison with #inverse_type=double#. + */ + vector > inverse; +}; + + +/*---------------------------- blocksparsematrix.h ---------------------------*/ +/* end of #ifndef __blocksparsematrix_H */ +#endif +/*---------------------------- blocksparsematrix.h ---------------------------*/ diff --git a/deal.II/lac/include/lac/blocksparsematrix.templates.h b/deal.II/lac/include/lac/blocksparsematrix.templates.h new file mode 100644 index 0000000000..ae2835d2ce --- /dev/null +++ b/deal.II/lac/include/lac/blocksparsematrix.templates.h @@ -0,0 +1,183 @@ +/*---------------------------- blocksparsematrix.templates.h --------------------*/ +/* $Id$ */ +/* Ralf Hartmann, University of Heidelberg */ +/*---------------------------- blocksparsematrix.templates.h ---------------------*/ + + + +#include +#include + + +template +BlockSparseMatrix::BlockSparseMatrix (): + blocksize(0) {}; + + +template +BlockSparseMatrix::~BlockSparseMatrix () +{ + if (inverse.size()!=0) + inverse.erase(inverse.begin(), inverse.end()); +} + + +template +void BlockSparseMatrix::reinit () +{ + if (inverse.size()!=0) + inverse.erase(inverse.begin(), inverse.end()); + blocksize=0; + SparseMatrix::reinit (); +} + + +template +void BlockSparseMatrix::reinit (const SparseMatrixStruct &sparsity) +{ + if (inverse.size()!=0) + inverse.erase(inverse.begin(), inverse.end()); + blocksize=0; + SparseMatrix::reinit (sparsity); +} + + +template +void BlockSparseMatrix::clear () +{ + SparseMatrix::clear(); + if (inverse.size()!=0) + inverse.erase(inverse.begin(), inverse.end()); + blocksize=0; +} + +template +void BlockSparseMatrix::set_block_size(unsigned int bsize) { + blocksize=bsize; +} + + +template +unsigned int BlockSparseMatrix::block_size() const { + return blocksize; +} + + +template +template +void BlockSparseMatrix::precondition_BlockSOR (Vector &dst, const Vector &src, const number) const +{ + Assert (m() == n(), ExcMatrixNotSquare()); + Assert (blocksize!=0, ExcBlockSizeNotSet()); + Assert (m()%blocksize==0, ExcWrongBlockSize(blocksize, m())); + unsigned int n_cells=m()/blocksize; + Assert (inverse.size()==0 || inverse.size()==n_cells, + ExcWrongNumberOfInverses(inverse.size(), n_cells)); + + const SparseMatrixStruct &spars=get_sparsity_pattern(); + const unsigned int *rowstart = spars.get_rowstart_indices(); + const int *columns = spars.get_column_numbers(); + + Vector b_cell(blocksize), x_cell(blocksize); + + // cell_row, cell_column are the + // numbering of the blocks (cells). + // row_cell, column_cell are the local + // numbering of the unknowns in the + // blocks. + // row, column are the global numbering + // of the unkowns. + unsigned int row, column, row_cell, begin_diag_block=0; + number2 b_cell_row; + + if (inverse.size()==0) + { + FullMatrix M_cell(blocksize); + for (unsigned int cell=0; cell(columns[j])) + < begin_diag_block) + b_cell_row -= global_entry(j) * dst(column); + b_cell(row_cell)=b_cell_row; + for (unsigned int column_cell=0, column=cell*blocksize; + column_cell(columns[j])) < begin_diag_block) + { + b_cell_row -= global_entry(j) * dst(column); + } + b_cell(row_cell)=b_cell_row; + } + inverse[cell].vmult(x_cell, b_cell); + // distribute x_cell to dst + for (row=cell*blocksize, row_cell=0; row_cell +void BlockSparseMatrix::invert_diagblocks() +{ + Assert (m() == n(), ExcMatrixNotSquare()); + Assert (inverse.size()==0, ExcInverseMatricesAlreadyExist()); + + Assert (blocksize!=0, ExcBlockSizeNotSet()); + Assert (m()%blocksize==0, ExcWrongBlockSize(blocksize, m())); + + unsigned int n_cells=m()/blocksize; + + // cell_row, cell_column are the + // numbering of the blocks (cells). + // row_cell, column_cell are the local + // numbering of the unknowns in the + // blocks. + // row, column are the global numbering + // of the unkowns. + + inverse = vector > (n_cells, FullMatrix(blocksize)); + FullMatrix M_cell(blocksize); + + for (unsigned int cell=0, row=0; cell +#include +//#include + + +template +bool datatype_compatible (); + + +// explicit instantiations for "float" BlockSparseMatrix +template class BlockSparseMatrix; + +template void BlockSparseMatrix::precondition_BlockSOR ( + Vector &, const Vector &, float) const; +template void BlockSparseMatrix::precondition_BlockSOR ( + Vector &, const Vector &, float) const; + +// the instantiation for class BlockSparseMatrix is skipped +// because it does not make sence to have inverse block matrices with +// higher precision than the matrix itself + +// explicit instantiations for "double" BlockSparseMatrix +template class BlockSparseMatrix; + +template void BlockSparseMatrix::precondition_BlockSOR ( + Vector &, const Vector &, double) const; +template void BlockSparseMatrix::precondition_BlockSOR ( + Vector &, const Vector &, double) const; + +template class BlockSparseMatrix; + +template void BlockSparseMatrix::precondition_BlockSOR ( + Vector &, const Vector &, double) const; +template void BlockSparseMatrix::precondition_BlockSOR ( + Vector &, const Vector &, double) const; + + + + +/*---------------------------- blocksparsematrix.cc ---------------------------*/