From 4c0d03db12b59dbd7bc905576284872caf712f41 Mon Sep 17 00:00:00 2001 From: hartmann Date: Wed, 10 Feb 1999 18:04:00 +0000 Subject: [PATCH] Block matrix for using of block SOR as direct solver or preconditioning, optional storing of the inverses of the diagonal matrices git-svn-id: https://svn.dealii.org/trunk@777 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/include/lac/dblocksmatrix.h | 154 +++++++++++++++++++++ deal.II/lac/source/dblocksmatrix.cc | 175 ++++++++++++++++++++++++ 2 files changed, 329 insertions(+) create mode 100644 deal.II/lac/include/lac/dblocksmatrix.h create mode 100644 deal.II/lac/source/dblocksmatrix.cc diff --git a/deal.II/lac/include/lac/dblocksmatrix.h b/deal.II/lac/include/lac/dblocksmatrix.h new file mode 100644 index 0000000000..5501acc8d1 --- /dev/null +++ b/deal.II/lac/include/lac/dblocksmatrix.h @@ -0,0 +1,154 @@ +/*---------------------------- dblocksmatrix.h ---------------------------*/ +/* $Id$ */ +#ifndef __dblocksmatrix_H +#define __dblocksmatrix_H +/*---------------------------- dblocksmatrix.h ---------------------------*/ + + +#include +#include +#include + +/** + * Double precision block sparse matrix. + * The block matrix assumes the matrix consisting of blocks on + * the diagonal. These diagonal blocks and the elements 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 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. + * @author Ralf Hartmann, 1999 + */ +class dBlockSMatrix: public dSMatrix +{ + public: + /** + * Constructor + */ + dBlockSMatrix(); + + /** + * Destructor + */ + virtual ~dBlockSMatrix(); + + /** + * Call #dSMatrix::reinit()# and + * delete the inverse matrices if existent. + */ + + virtual void reinit(); + + /** + * Call #dSMatrix::reinit + * (const dSMatrixStruct &sparsity)# and + * delete the inverse matrices if existent. + */ + virtual void reinit (const dSMatrixStruct &sparsity); + + /** + * Call #dSMatrix::clear# and + * delete the inverse matrices if existent. + */ + virtual void clear (); + + /** + * Stores the inverse matrices of + * the diagonal blocks matrices + * in #inverse#. This costs some + * additional memory (for DG + * methods about 1/3 of that used for + * the matrix) but it + * makes the preconditioning much faster. + */ + 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. + */ + void precondition_BlockSOR (dVector &dst, const dVector &src) 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."); + + DeclException2 (ExcWrongInverses, + int, int, + << "There are " << arg1 + << " inverse matrices but " << arg2 + << " cells."); + + /** + * Exception + */ + DeclException0 (ExcInverseMatricesDoNotExist); + + /** + * Exception + */ + DeclException0 (ExcInverseMatricesAlreadyExist); + + /** + * Exception + */ + DeclException0 (ExcBlockSizeNotSet); + + /** + * Exception + */ + DeclException0 (ExcInternalError); + + private: + /** + * size of the blocks. + */ + unsigned int blocksize; + + /** + * stores the inverse matrices of + * the diagonal blocks matrices + */ + vector inverse; +}; + + + +/*---------------------------- dblocksmatrix.h ---------------------------*/ +/* end of #ifndef __dblocksmatrix_H */ +#endif +/*---------------------------- dblocksmatrix.h ---------------------------*/ diff --git a/deal.II/lac/source/dblocksmatrix.cc b/deal.II/lac/source/dblocksmatrix.cc new file mode 100644 index 0000000000..4e119332c6 --- /dev/null +++ b/deal.II/lac/source/dblocksmatrix.cc @@ -0,0 +1,175 @@ +/*---------------------------- dblocksmatrix.cc ---------------------------*/ +/* $Id$ */ +/*---------------------------- dblocksmatrix.cc ---------------------------*/ + +#include +#include + + +dBlockSMatrix::dBlockSMatrix (): + blocksize(0) {}; + +dBlockSMatrix::~dBlockSMatrix () +{ + if (inverse.size()!=0) + inverse.erase(inverse.begin(), inverse.end()); +} + + +void dBlockSMatrix::reinit () +{ + if (inverse.size()!=0) + inverse.erase(inverse.begin(), inverse.end()); + blocksize=0; + dSMatrix::reinit (); +} + + +void dBlockSMatrix::reinit (const dSMatrixStruct &sparsity) +{ + if (inverse.size()!=0) + inverse.erase(inverse.begin(), inverse.end()); + blocksize=0; + dSMatrix::reinit (sparsity); +} + + +void dBlockSMatrix::clear () +{ + dSMatrix::clear(); + if (inverse.size()!=0) + inverse.erase(inverse.begin(), inverse.end()); + blocksize=0; +} + + +void dBlockSMatrix::set_block_size(unsigned int bsize) { + blocksize=bsize; +} + + + +unsigned int dBlockSMatrix::block_size() const { + return blocksize; +} + + + +void dBlockSMatrix::precondition_BlockSOR (dVector &dst, const dVector &src) 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, + ExcWrongInverses(inverse.size(), n_cells)); + + const dSMatrixStruct &spars=get_sparsity_pattern(); + const unsigned int *rowstart = spars.get_rowstart_indices(); + const int *columns = spars.get_column_numbers(); + + dVector 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; + double b_cell_row; + + if (inverse.size()==0) + { + dFMatrix 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 4 + if (blocksize<=4) + inverse[cell].invert(M_cell); + else + { + M_cell.gauss_jordan(); + inverse[cell]=M_cell; + } + } +} + + + +/*---------------------------- dblocksmatrix.cc ---------------------------*/ -- 2.39.5