+++ /dev/null
-/*---------------------------- blocksparsematrix.h ---------------------*/
-/* $Id$ */
-/* Ralf Hartmann, University of Heidelberg */
-#ifndef __blocksparsematrix_H
-#define __blocksparsematrix_H
-/*---------------------------- blocksparsematrix.h ---------------------*/
-
-
-#include <lac/sparsematrix.h>
-#include <lac/forward-declarations.h>
-
-#include <vector>
-
-
-
-
-/**
- * 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 <typename number, typename inverse_type>
-class BlockSparseMatrix: public SparseMatrix<number>
-{
- 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 <typename number2>
- void precondition_BlockSOR (Vector<number2> &dst, const Vector<number2> &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<inverse_type># matrices.
- * For BlockSOR as preconditioning
- * using #inverse_type=float# saves memory
- * in comparison with #inverse_type=double#.
- */
- vector<FullMatrix<inverse_type> > inverse;
-};
-
-
-/*---------------------------- blocksparsematrix.h ---------------------------*/
-/* end of #ifndef __blocksparsematrix_H */
-#endif
-/*---------------------------- blocksparsematrix.h ---------------------------*/
--- /dev/null
+/*---------------------------- precondition_block.h ---------------------------*/
+/* $Id$ */
+#ifndef __precondition_block_H
+#define __precondition_block_H
+/*---------------------------- precondition_block.h ---------------------------*/
+
+#include <lac/forward-declarations.h>
+#include <base/exceptions.h>
+#include <base/subscriptor.h>
+#include <base/smartpointer.h>
+
+#include <vector>
+
+
+/**
+ * Base class for #PreconditionBlockJacobi#, #PreconditionBlockSOR#, ...
+ * This class assumes the #SparseMatrix<number># 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<number> &A)#.
+ */
+template<typename number, typename inverse_type>
+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<number> &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<inverse_type># matrices.
+ * For BlockSOR as preconditioning
+ * using #inverse_type=float# saves memory
+ * in comparison with #inverse_type=double#.
+ */
+ vector<FullMatrix<inverse_type> > 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<const SparseMatrix<number> > 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<typename number, typename inverse_type>
+class PreconditionBlockSOR : public PreconditionBlock<number, inverse_type>
+{
+ 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 <typename number2>
+ void operator() (Vector<number2>&, const Vector<number2>&) const;
+
+ private:
+ number omega;
+};
+
+
+
+/*------------------------ precondition_block.h ---------------------------*/
+/* end of #ifndef __precondition_block_H */
+#endif
+/*------------------------ precondition_block.h ---------------------------*/
-/*---------------------------- blocksparsematrix.templates.h --------------------*/
+/*----------------- precondition_block.templates.h ----------------------*/
/* $Id$ */
-/* Ralf Hartmann, University of Heidelberg */
-/*---------------------------- blocksparsematrix.templates.h ---------------------*/
+/* Ralf Hartmann, University of Heidelberg */
-
-#include <lac/blocksparsematrix.h>
+#include <lac/precondition_block.h>
#include <lac/vector.h>
#include <lac/fullmatrix.h>
+#include <lac/sparsematrix.h>
template <typename number, typename inverse_type>
-BlockSparseMatrix<number, inverse_type>::BlockSparseMatrix ():
- blocksize(0) {};
+PreconditionBlock<number, inverse_type>::PreconditionBlock ():
+ blocksize(0),
+ A(0) {};
template <typename number, typename inverse_type>
-BlockSparseMatrix<number, inverse_type>::~BlockSparseMatrix ()
+PreconditionBlock<number, inverse_type>::~PreconditionBlock ()
{
if (inverse.size()!=0)
inverse.erase(inverse.begin(), inverse.end());
}
-
template <typename number, typename inverse_type>
-void BlockSparseMatrix<number, inverse_type>::reinit ()
+void PreconditionBlock<number, inverse_type>::clear ()
{
if (inverse.size()!=0)
inverse.erase(inverse.begin(), inverse.end());
blocksize=0;
- SparseMatrix<number>::reinit ();
}
-
template <typename number, typename inverse_type>
-void BlockSparseMatrix<number, inverse_type>::reinit (const SparseMatrixStruct &sparsity)
+void PreconditionBlock<number, inverse_type>::use_matrix(
+ const SparseMatrix<number> &M)
{
- if (inverse.size()!=0)
- inverse.erase(inverse.begin(), inverse.end());
- blocksize=0;
- SparseMatrix<number>::reinit (sparsity);
+ A=0;
+ A=&M;
}
template <typename number, typename inverse_type>
-void BlockSparseMatrix<number, inverse_type>::clear ()
-{
- SparseMatrix<number>::clear();
- if (inverse.size()!=0)
- inverse.erase(inverse.begin(), inverse.end());
- blocksize=0;
+void PreconditionBlock<number, inverse_type>::set_block_size(unsigned int bsize) {
+ blocksize=bsize;
}
+
template <typename number, typename inverse_type>
-void BlockSparseMatrix<number, inverse_type>::set_block_size(unsigned int bsize) {
- blocksize=bsize;
+unsigned int PreconditionBlock<number, inverse_type>::block_size() const {
+ return blocksize;
}
template <typename number, typename inverse_type>
-unsigned int BlockSparseMatrix<number, inverse_type>::block_size() const {
- return blocksize;
+void PreconditionBlock<number, inverse_type>::invert_diagblocks()
+{
+ Assert (A!=0, ExcNoMatrixGivenToUse());
+ const SparseMatrix<number> &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<FullMatrix<inverse_type> > (
+ n_cells, FullMatrix<inverse_type>(blocksize));
+ FullMatrix<inverse_type> M_cell(blocksize);
+
+ for (unsigned int cell=0, row=0; cell<n_cells; ++cell)
+ {
+ for (unsigned int row_cell=0; row_cell<blocksize; ++row_cell, ++row)
+ for (unsigned int column_cell=0, column=cell*blocksize;
+ column_cell<blocksize; ++column_cell, ++column)
+ M_cell(row_cell,column_cell)=M(row,column);
+// try
+// {
+ if (blocksize <=4)
+// {
+ inverse[cell].invert(M_cell);
+// }
+// catch (ExcNotImplemented &)
+ else
+ {
+ M_cell.gauss_jordan();
+ inverse[cell]=M_cell;
+ }
+ }
+ // release the lock on A
+ A=0;
}
+
+/*--------------------- PreconditionBlockSOR -----------------------*/
+
+template<typename number, typename inverse_type>
+PreconditionBlockSOR<number, inverse_type>::PreconditionBlockSOR(const number omega):
+ omega(omega) {}
+
+
+template<typename number, typename inverse_type>
+PreconditionBlockSOR<number, inverse_type>::~PreconditionBlockSOR(){}
+
+
template <typename number, typename inverse_type>
template <typename number2>
-void BlockSparseMatrix<number, inverse_type>::precondition_BlockSOR (Vector<number2> &dst, const Vector<number2> &src, const number) const
+void PreconditionBlockSOR<number, inverse_type>::operator() (Vector<number2> &dst, const Vector<number2> &src) const
{
- Assert (m() == n(), ExcMatrixNotSquare());
+ Assert(A!=0, ExcNoMatrixGivenToUse());
+ const SparseMatrix<number> &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();
for (unsigned int j=rowstart[row]; j<rowstart[row+1]; ++j)
if ((column=static_cast<unsigned int>(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<blocksize; ++column_cell, ++column)
- M_cell(row_cell,column_cell)=(*this)(row,column);
+ M_cell(row_cell,column_cell)=M(row,column);
}
M_cell.householder(b_cell);
M_cell.backward(x_cell,b_cell);
// distribute x_cell to dst
for (row=cell*blocksize, row_cell=0; row_cell<blocksize; ++row_cell, ++row)
- dst(row)=x_cell(row_cell);
+ dst(row)=omega*x_cell(row_cell);
begin_diag_block+=blocksize;
- }
+ }
}
else
for (unsigned int cell=0; cell<n_cells; ++cell)
for (unsigned int j=rowstart[row]; j<rowstart[row+1]; ++j)
if ((column=static_cast<unsigned int>(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<blocksize; ++row_cell, ++row)
- dst(row)=x_cell(row_cell);
+ dst(row)=omega*x_cell(row_cell);
begin_diag_block+=blocksize;
}
-template <typename number, typename inverse_type>
-void BlockSparseMatrix<number, inverse_type>::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<FullMatrix<inverse_type> > (n_cells, FullMatrix<inverse_type>(blocksize));
- FullMatrix<inverse_type> M_cell(blocksize);
-
- for (unsigned int cell=0, row=0; cell<n_cells; ++cell)
- {
- for (unsigned int row_cell=0; row_cell<blocksize; ++row_cell, ++row)
- for (unsigned int column_cell=0, column=cell*blocksize;
- column_cell<blocksize; ++column_cell, ++column)
- M_cell(row_cell,column_cell)=(*this)(row,column);
- try
- {
- inverse[cell].invert(M_cell);
- }
- catch (ExcNotImplemented &)
- {
- M_cell.gauss_jordan();
- inverse[cell]=M_cell;
- }
- }
-}
-
-/*---------------------------- blocksparsematrix.templates.h ---------------------*/
+/*---------------------- precondition_block.templates.h ----------------*/
+++ /dev/null
-/*---------------------------- blocksparsematrix.cc ---------------------------*/
-/* $Id$ */
-/* Ralf Hartmann, University of Heidelberg */
-/*---------------------------- blocksparsematrix.cc ---------------------------*/
-
-
-
-#include <lac/blocksparsematrix.h>
-#include <lac/blocksparsematrix.templates.h>
-//#include <lac/vector.h>
-
-
-
-// explicit instantiations for "float" BlockSparseMatrix
-template class BlockSparseMatrix<float, float>;
-
-template void BlockSparseMatrix<float, float>::precondition_BlockSOR (
- Vector<float> &, const Vector<float> &, float) const;
-template void BlockSparseMatrix<float, float>::precondition_BlockSOR (
- Vector<double> &, const Vector<double> &, float) const;
-
-// the instantiation for class BlockSparseMatrix<float, double> 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<double, float>;
-
-template void BlockSparseMatrix<double, float>::precondition_BlockSOR (
- Vector<float> &, const Vector<float> &, double) const;
-template void BlockSparseMatrix<double, float>::precondition_BlockSOR (
- Vector<double> &, const Vector<double> &, double) const;
-
-template class BlockSparseMatrix<double, double>;
-
-template void BlockSparseMatrix<double, double>::precondition_BlockSOR (
- Vector<float> &, const Vector<float> &, double) const;
-template void BlockSparseMatrix<double, double>::precondition_BlockSOR (
- Vector<double> &, const Vector<double> &, double) const;
-
-
-
-
-/*---------------------------- blocksparsematrix.cc ---------------------------*/
--- /dev/null
+/*---------------------------- precondition_block.cc ---------------------------*/
+/* $Id$ */
+/* Ralf Hartmann, University of Heidelberg */
+
+
+
+#include <lac/precondition_block.templates.h>
+
+
+
+
+// explicit instantiations for "float" PreconditionBlock
+template class PreconditionBlock<float, float>;
+
+// the instantiation for class PreconditionBlock<float, double> 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<double, float>;
+
+template class PreconditionBlock<double, double>;
+
+
+/*--------------------- PreconditionBlockGaussSeidel -----------------------*/
+
+
+// explicit instantiations for "float" PreconditionBlock
+template class PreconditionBlockSOR<float, float>;
+
+template void PreconditionBlockSOR<float, float>::operator() (
+ Vector<float> &, const Vector<float> &) const;
+template void PreconditionBlockSOR<float, float>::operator() (
+ Vector<double> &, const Vector<double> &) const;
+
+
+// the instantiation for class PreconditionBlockSOR<float, double> 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<double, float>;
+
+
+template void PreconditionBlockSOR<double, float>::operator() (
+ Vector<float> &, const Vector<float> &) const;
+template void PreconditionBlockSOR<double, float>::operator() (
+ Vector<double> &, const Vector<double> &) const;
+
+template class PreconditionBlockSOR<double, double>;
+
+template void PreconditionBlockSOR<double, double>::operator() (
+ Vector<float> &, const Vector<float> &) const;
+template void PreconditionBlockSOR<double, double>::operator() (
+ Vector<double> &, const Vector<double> &) const;
+
+
+
+
+
+
+/*---------------------------- precondition_block.cc ---------------------------*/