]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Block matrix for using of block SOR as direct solver or preconditioning, optional...
authorhartmann <hartmann@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 10 Feb 1999 18:04:00 +0000 (18:04 +0000)
committerhartmann <hartmann@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 10 Feb 1999 18:04:00 +0000 (18:04 +0000)
git-svn-id: https://svn.dealii.org/trunk@777 0785d39b-7218-0410-832d-ea1e28bc413d

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

diff --git a/deal.II/lac/include/lac/dblocksmatrix.h b/deal.II/lac/include/lac/dblocksmatrix.h
new file mode 100644 (file)
index 0000000..5501acc
--- /dev/null
@@ -0,0 +1,154 @@
+/*----------------------------   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     ---------------------------*/
diff --git a/deal.II/lac/source/dblocksmatrix.cc b/deal.II/lac/source/dblocksmatrix.cc
new file mode 100644 (file)
index 0000000..4e11933
--- /dev/null
@@ -0,0 +1,175 @@
+/*----------------------------   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     ---------------------------*/

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.