From: hartmann Date: Mon, 3 May 1999 15:31:22 +0000 (+0000) Subject: Replace the BlockSparseMatrix by the PreconditionBlock X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=b9dcc138bb3eb9a62206106df0fe25db9c30832c;p=dealii-svn.git Replace the BlockSparseMatrix by the PreconditionBlock git-svn-id: https://svn.dealii.org/trunk@1250 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/lac/include/lac/blocksparsematrix.h b/deal.II/lac/include/lac/blocksparsematrix.h deleted file mode 100644 index 5d414208ee..0000000000 --- a/deal.II/lac/include/lac/blocksparsematrix.h +++ /dev/null @@ -1,176 +0,0 @@ -/*---------------------------- 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# function. - * - * 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 number representation in - * the sparse matrix, the second denotes the type of number representation 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/precondition_block.h b/deal.II/lac/include/lac/precondition_block.h new file mode 100644 index 0000000000..16c6ec8514 --- /dev/null +++ b/deal.II/lac/include/lac/precondition_block.h @@ -0,0 +1,232 @@ +/*---------------------------- precondition_block.h ---------------------------*/ +/* $Id$ */ +#ifndef __precondition_block_H +#define __precondition_block_H +/*---------------------------- precondition_block.h ---------------------------*/ + +#include +#include +#include +#include + +#include + + +/** + * Base class for #PreconditionBlockJacobi#, #PreconditionBlockSOR#, ... + * This class assumes the #SparseMatrix# consisting of invertible blocks + * of #blocksize# on the diagonal and provides the inversion of the diagonal blocks + * of the matrix. NOT only diagonal block matrices are allowed but all + * matrices of arbitrary structure with the minimal property of having got + * invertible blocks on the diagonal! + * + * This block matrix structure is given e.g. for the DG method + * for the transport equation. For a downstream numbering the matrices + * even have got a lower block diagonal matrix structure, i.e. the matrices + * are empty above the diagonal blocks. + * + * For all matrices that are empty above and below the diagonal + * blocks (i.e. for all block diagonal matrices) the #BlockJacobi# preconditioner + * is a direct solver. For all matrices that are empty only above the diagonal blocks + * (e.g. the matrices one get by the DG method with downstream numbering) the + * #BlockSOR# is a direct solver. + * + * This first implementation of the #PreconditionBlock# 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 number representation in + * the sparse matrix, the second denotes the type of number representation in + * which the inverse diagonal block + * matrices are stored by #invert_diagblocks(const SparseMatrix &A)#. + */ +template +class PreconditionBlock: public Subscriptor +{ + public: + /** + * Constructor. Takes the matrix as + * argument. + */ + PreconditionBlock(); + + /** + * Destructor. + */ + virtual ~PreconditionBlock(); + + /** + * Deletes the inverse diagonal block + * matrices if existent, sets the + * blocksize to 0, hence leaves the + * class in the state that it had + * directly after + * calling the constructor. + */ + virtual void clear(); + + /** + * Takes the matrix that should be used + * for the preconditioning. + */ + void use_matrix(const SparseMatrix &M); + + /** + * Set the right block size before calling + * #invert_diagblocks# or calling the + * #operator ()# function + * of #PreconditionBlockSOR# + * or #PreconditionBlockJacobi#. + * If #block_size==1# BlockSOR or + * BlockJacobi are equal to the + * standard SOR or Jacobi + * preconditioner. + */ + void set_block_size (const unsigned int bsize); + + /** + * 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(); + + /** + * 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); + + /** + * Exception + */ + DeclException0 (ExcMatrixNotSquare); + + /** + * Exception + */ + DeclException0 (ExcNoMatrixGivenToUse); + + protected: + /** + * 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; + + /** + * Pointer to the matrix. Make sure that + * the matrix exists as long as this class + * needs it, i.e. until calling #invert_diagblocks#, + * or (if the inverse matrices should not be + * stored) until the last call of the + * preconditoining #operator()# function of the + * derived classes. + */ + SmartPointer > A; +}; + + + +/** + * Block SOR preconditioning. + * + * The diagonal blocks and the elements + * (of arbitray structure) below the diagonal blocks are used + * in the #operator ()# function of this class. + */ +template +class PreconditionBlockSOR : public PreconditionBlock +{ + public: + /** + * Constructor. Takes the damping + * Parameter as + */ + PreconditionBlockSOR(number omega=1.); + + /** + * Destructor. + */ + virtual ~PreconditionBlockSOR(); + + /** + * Execute block SOR preconditioning. + * Make sure that the right + * block size + * of the matrix is set by + * #set_block_size# + * before calling this function. + * + * This function 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 operator() (Vector&, const Vector&) const; + + private: + number omega; +}; + + + +/*------------------------ precondition_block.h ---------------------------*/ +/* end of #ifndef __precondition_block_H */ +#endif +/*------------------------ precondition_block.h ---------------------------*/ diff --git a/deal.II/lac/include/lac/blocksparsematrix.templates.h b/deal.II/lac/include/lac/precondition_block.templates.h similarity index 61% rename from deal.II/lac/include/lac/blocksparsematrix.templates.h rename to deal.II/lac/include/lac/precondition_block.templates.h index 7967fb9fd8..a15886ff34 100644 --- a/deal.II/lac/include/lac/blocksparsematrix.templates.h +++ b/deal.II/lac/include/lac/precondition_block.templates.h @@ -1,81 +1,131 @@ -/*---------------------------- blocksparsematrix.templates.h --------------------*/ +/*----------------- precondition_block.templates.h ----------------------*/ /* $Id$ */ -/* Ralf Hartmann, University of Heidelberg */ -/*---------------------------- blocksparsematrix.templates.h ---------------------*/ +/* Ralf Hartmann, University of Heidelberg */ - -#include +#include #include #include +#include template -BlockSparseMatrix::BlockSparseMatrix (): - blocksize(0) {}; +PreconditionBlock::PreconditionBlock (): + blocksize(0), + A(0) {}; template -BlockSparseMatrix::~BlockSparseMatrix () +PreconditionBlock::~PreconditionBlock () { if (inverse.size()!=0) inverse.erase(inverse.begin(), inverse.end()); } - template -void BlockSparseMatrix::reinit () +void PreconditionBlock::clear () { if (inverse.size()!=0) inverse.erase(inverse.begin(), inverse.end()); blocksize=0; - SparseMatrix::reinit (); } - template -void BlockSparseMatrix::reinit (const SparseMatrixStruct &sparsity) +void PreconditionBlock::use_matrix( + const SparseMatrix &M) { - if (inverse.size()!=0) - inverse.erase(inverse.begin(), inverse.end()); - blocksize=0; - SparseMatrix::reinit (sparsity); + A=0; + A=&M; } template -void BlockSparseMatrix::clear () -{ - SparseMatrix::clear(); - if (inverse.size()!=0) - inverse.erase(inverse.begin(), inverse.end()); - blocksize=0; +void PreconditionBlock::set_block_size(unsigned int bsize) { + blocksize=bsize; } + template -void BlockSparseMatrix::set_block_size(unsigned int bsize) { - blocksize=bsize; +unsigned int PreconditionBlock::block_size() const { + return blocksize; } template -unsigned int BlockSparseMatrix::block_size() const { - return blocksize; +void PreconditionBlock::invert_diagblocks() +{ + Assert (A!=0, ExcNoMatrixGivenToUse()); + const SparseMatrix &M=*A; + Assert (M.m() == M.n(), ExcMatrixNotSquare()); + Assert (inverse.size()==0, ExcInverseMatricesAlreadyExist()); + + Assert (blocksize!=0, ExcBlockSizeNotSet()); + Assert (M.m()%blocksize==0, ExcWrongBlockSize(blocksize, M.m())); + + unsigned int n_cells=M.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 +PreconditionBlockSOR::PreconditionBlockSOR(const number omega): + omega(omega) {} + + +template +PreconditionBlockSOR::~PreconditionBlockSOR(){} + + template template -void BlockSparseMatrix::precondition_BlockSOR (Vector &dst, const Vector &src, const number) const +void PreconditionBlockSOR::operator() (Vector &dst, const Vector &src) const { - Assert (m() == n(), ExcMatrixNotSquare()); + Assert(A!=0, ExcNoMatrixGivenToUse()); + const SparseMatrix &M=*A; + Assert (M.m() == M.n(), ExcMatrixNotSquare()); Assert (blocksize!=0, ExcBlockSizeNotSet()); - Assert (m()%blocksize==0, ExcWrongBlockSize(blocksize, m())); - unsigned int n_cells=m()/blocksize; + Assert (M.m()%blocksize==0, ExcWrongBlockSize(blocksize, M.m())); + unsigned int n_cells=M.m()/blocksize; Assert (inverse.size()==0 || inverse.size()==n_cells, ExcWrongNumberOfInverses(inverse.size(), n_cells)); - const SparseMatrixStruct &spars=get_sparsity_pattern(); + const SparseMatrixStruct &spars=M.get_sparsity_pattern(); const unsigned int *rowstart = spars.get_rowstart_indices(); const int *columns = spars.get_column_numbers(); @@ -102,20 +152,20 @@ void BlockSparseMatrix::precondition_BlockSOR (Vector(columns[j])) < begin_diag_block) - b_cell_row -= global_entry(j) * dst(column); + b_cell_row -= M.global_entry(j) * dst(column); b_cell(row_cell)=b_cell_row; for (unsigned int column_cell=0, column=cell*blocksize; column_cell::precondition_BlockSOR (Vector(columns[j])) < begin_diag_block) { - b_cell_row -= global_entry(j) * dst(column); + b_cell_row -= M.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::precondition_BlockSOR (Vector -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 - - - -// 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 sense 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 ---------------------------*/ diff --git a/deal.II/lac/source/precondition_block.cc b/deal.II/lac/source/precondition_block.cc new file mode 100644 index 0000000000..528dd90c02 --- /dev/null +++ b/deal.II/lac/source/precondition_block.cc @@ -0,0 +1,65 @@ +/*---------------------------- precondition_block.cc ---------------------------*/ +/* $Id$ */ +/* Ralf Hartmann, University of Heidelberg */ + + + +#include + + + + +// explicit instantiations for "float" PreconditionBlock +template class PreconditionBlock; + +// the instantiation for class PreconditionBlock is skipped +// because it does not make sense to have inverse block matrices with +// higher precision than the matrix itself + + +// explicit instantiations for "double" PreconditionBlock +template class PreconditionBlock; + +template class PreconditionBlock; + + +/*--------------------- PreconditionBlockGaussSeidel -----------------------*/ + + +// explicit instantiations for "float" PreconditionBlock +template class PreconditionBlockSOR; + +template void PreconditionBlockSOR::operator() ( + Vector &, const Vector &) const; +template void PreconditionBlockSOR::operator() ( + Vector &, const Vector &) const; + + +// the instantiation for class PreconditionBlockSOR is skipped +// because it does not make sense to have inverse block matrices with +// higher precision than the matrix itself + + + +// explicit instantiations for "double" PreconditionBlockSOR +template class PreconditionBlockSOR; + + +template void PreconditionBlockSOR::operator() ( + Vector &, const Vector &) const; +template void PreconditionBlockSOR::operator() ( + Vector &, const Vector &) const; + +template class PreconditionBlockSOR; + +template void PreconditionBlockSOR::operator() ( + Vector &, const Vector &) const; +template void PreconditionBlockSOR::operator() ( + Vector &, const Vector &) const; + + + + + + +/*---------------------------- precondition_block.cc ---------------------------*/