--- /dev/null
+/*---------------------------- blocksparsematrix.h ---------------------*/
+/* $Id$ */
+/* Ralf Hartmann, University of Heidelberg */
+#ifndef __blocksparsematrix_H
+#define __blocksparsematrix_H
+/*---------------------------- blocksparsematrix.h ---------------------*/
+
+
+#include <lac/sparsematrix.h>
+#include <lac/fullmatrix.h>
+#include <vector.h>
+
+/**
+ * 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 <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
+/*---------------------------- blocksparsematrix.templates.h --------------------*/
+/* $Id$ */
+/* Ralf Hartmann, University of Heidelberg */
+/*---------------------------- blocksparsematrix.templates.h ---------------------*/
+
+
+
+#include <lac/blocksparsematrix.h>
+#include <lac/vector.h>
+
+
+template <typename number, typename inverse_type>
+BlockSparseMatrix<number, inverse_type>::BlockSparseMatrix ():
+ blocksize(0) {};
+
+
+template <typename number, typename inverse_type>
+BlockSparseMatrix<number, inverse_type>::~BlockSparseMatrix ()
+{
+ if (inverse.size()!=0)
+ inverse.erase(inverse.begin(), inverse.end());
+}
+
+
+template <typename number, typename inverse_type>
+void BlockSparseMatrix<number, inverse_type>::reinit ()
+{
+ 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)
+{
+ if (inverse.size()!=0)
+ inverse.erase(inverse.begin(), inverse.end());
+ blocksize=0;
+ SparseMatrix<number>::reinit (sparsity);
+}
+
+
+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;
+}
+
+template <typename number, typename inverse_type>
+void BlockSparseMatrix<number, inverse_type>::set_block_size(unsigned int bsize) {
+ blocksize=bsize;
+}
+
+
+template <typename number, typename inverse_type>
+unsigned int BlockSparseMatrix<number, inverse_type>::block_size() const {
+ return blocksize;
+}
+
+
+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
+{
+ 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<number2> 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<number> M_cell(blocksize);
+ for (unsigned int cell=0; cell<n_cells; ++cell)
+ {
+ for (row=cell*blocksize, row_cell=0; row_cell<blocksize; ++row_cell, ++row)
+ {
+ b_cell_row=src(row);
+ 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_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.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);
+
+ begin_diag_block+=blocksize;
+ }
+ }
+ else
+ for (unsigned int cell=0; cell<n_cells; ++cell)
+ {
+ for (row=cell*blocksize, row_cell=0; row_cell<blocksize; ++row_cell, ++row)
+ {
+ b_cell_row=src(row);
+ 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_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);
+
+ 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);
+ if (blocksize<=4)
+ inverse[cell].invert(M_cell);
+ else
+ {
+ M_cell.gauss_jordan();
+ inverse[cell]=M_cell;
+ }
+ }
+}
+
+
+
+/*---------------------------- blocksparsematrix.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>
+
+
+template <typename number, typename blocknumber>
+bool datatype_compatible ();
+
+
+// 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 sence 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 ---------------------------*/