]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
replace the old dBlockSMatrix by the new template BlockSparseMatrix with a second...
authorhartmann <hartmann@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 26 Feb 1999 15:00:13 +0000 (15:00 +0000)
committerhartmann <hartmann@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 26 Feb 1999 15:00:13 +0000 (15:00 +0000)
git-svn-id: https://svn.dealii.org/trunk@917 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/blocksparsematrix.h [new file with mode: 0644]
deal.II/lac/include/lac/blocksparsematrix.templates.h [new file with mode: 0644]
deal.II/lac/source/blocksparsematrix.cc [new file with mode: 0644]

diff --git a/deal.II/lac/include/lac/blocksparsematrix.h b/deal.II/lac/include/lac/blocksparsematrix.h
new file mode 100644 (file)
index 0000000..5aa47cc
--- /dev/null
@@ -0,0 +1,171 @@
+/*----------------------------   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     ---------------------------*/
diff --git a/deal.II/lac/include/lac/blocksparsematrix.templates.h b/deal.II/lac/include/lac/blocksparsematrix.templates.h
new file mode 100644 (file)
index 0000000..ae2835d
--- /dev/null
@@ -0,0 +1,183 @@
+/*----------------------------   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     ---------------------*/
diff --git a/deal.II/lac/source/blocksparsematrix.cc b/deal.II/lac/source/blocksparsematrix.cc
new file mode 100644 (file)
index 0000000..d80607c
--- /dev/null
@@ -0,0 +1,47 @@
+/*----------------------------   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     ---------------------------*/

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.