--- /dev/null
+/*---------------------------- dblocksmatrix.h ---------------------------*/
+/* $Id$ */
+#ifndef __dblocksmatrix_H
+#define __dblocksmatrix_H
+/*---------------------------- dblocksmatrix.h ---------------------------*/
+
+
+#include <lac/dsmatrix.h>
+#include <lac/dfmatrix.h>
+#include <vector.h>
+
+/**
+ * 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<dFMatrix> inverse;
+};
+
+
+
+/*---------------------------- dblocksmatrix.h ---------------------------*/
+/* end of #ifndef __dblocksmatrix_H */
+#endif
+/*---------------------------- dblocksmatrix.h ---------------------------*/
--- /dev/null
+/*---------------------------- dblocksmatrix.cc ---------------------------*/
+/* $Id$ */
+/*---------------------------- dblocksmatrix.cc ---------------------------*/
+
+#include <lac/dblocksmatrix.h>
+#include <lac/dvector.h>
+
+
+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<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;
+ }
+}
+
+
+void dBlockSMatrix::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;
+
+ inverse.insert(inverse.begin(), n_cells, dFMatrix(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.
+ dFMatrix 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);
+
+ // perhaps #dFMatrix::invert# should
+ // be change such that it calls
+ // #gauss_jordan()# automatically
+ // if blocksize > 4
+ if (blocksize<=4)
+ inverse[cell].invert(M_cell);
+ else
+ {
+ M_cell.gauss_jordan();
+ inverse[cell]=M_cell;
+ }
+ }
+}
+
+
+
+/*---------------------------- dblocksmatrix.cc ---------------------------*/