]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Relaxation parameter in PreconditionBlock
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 18 Sep 2000 19:32:52 +0000 (19:32 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 18 Sep 2000 19:32:52 +0000 (19:32 +0000)
git-svn-id: https://svn.dealii.org/trunk@3333 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/precondition_block.h
deal.II/lac/include/lac/precondition_block.templates.h
deal.II/lac/source/precondition_block.cc

index 49c3778fd810d8eb5ab83350f95332426ac151e7..509bab60faf6bb4240a2ce3f2bfbec99e57af5eb 100644 (file)
@@ -96,9 +96,14 @@ class PreconditionBlock
                                      * second step, the inverses of
                                      * the diagonal blocks may be
                                      * computed.
+                                     *
+                                     * Additionally, a relaxation
+                                     * parameter for derived classes
+                                     * may be provided.
                                      */
     void initialize (const SparseMatrix<number>& A,
-                    const unsigned int block_size);
+                    const unsigned int block_size,
+                    const double relaxation = 1.);
 
                                     /**
                                      * Deletes the inverse diagonal
@@ -192,6 +197,11 @@ class PreconditionBlock
                                      * Exception
                                      */
     DeclException0 (ExcMatrixNotSquare);
+
+                                    /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcDiagonalsNotStored);
    
   protected:
                                     /**
@@ -200,6 +210,12 @@ class PreconditionBlock
                                      */
     const FullMatrix<inverse_type>& inverse (unsigned int i) const;
     
+                                    /**
+                                     * Access to the diagonal
+                                     * blocks.
+                                     */
+    const FullMatrix<inverse_type>& diagonal (unsigned int i) const;
+    
                                     /**
                                      * Size of the blocks. Each
                                      * diagonal block is assumed to
@@ -220,6 +236,18 @@ class PreconditionBlock
                                      * derived classes.
                                      */
     SmartPointer<const SparseMatrix<number> > A;
+                                    /**
+                                     * Relaxation parameter to be
+                                     * used by derived classes.
+                                     */
+    double relaxation;
+
+
+                                    /**
+                                     * Flag for storing the diagonal
+                                     * blocks of the matrix.
+                                     */
+    bool store_diagonals;
 
   private:
                                     /**
@@ -227,14 +255,22 @@ class PreconditionBlock
                                      * matrices of the diagonal
                                      * blocks matrices as
                                      * @p{FullMatrix<inverse_type>}
-                                     * matrices.  For BlockSOR as
-                                     * preconditioning using
+                                     * matrices. Using
                                      * @p{inverse_type=float} saves
                                      * memory in comparison with
                                      * @p{inverse_type=double}.
                                      */
-    vector<FullMatrix<inverse_type> > _inverse;
+    vector<FullMatrix<inverse_type> > var_inverse;
 
+                                    /**
+                                     * Storage of the original diagonal blocks.
+                                     * These are only filled if @p{store_diagonals}
+                                     * is @p{true}.
+                                     *
+                                     * Used by the blocked SSOR method.
+                                     */
+    vector<FullMatrix<inverse_type> > var_diagonal;
+                                     
                                     /**
                                      * Flag for diagonal compression.
                                      * @see set_same_diagonal()
@@ -302,29 +338,19 @@ class PreconditionBlockJacobi : public Subscriptor,
 /**
  * Block SOR preconditioning.
  *
- * The diagonal blocks and the elements 
- * (of arbitray structure) below the diagonal blocks are used
- * in the @p{operator ()} function of this class.
+ * The functions @p{vmult} and @p{Tvmult} execute a (transposed)
+ * block-SOR step, based on the blocks in @ref{PreconditionBlock}. The
+ * elements outside the diagonal blocks may be distributed
+ * arbitrarily.
  *
  * See @ref{PreconditionBlock} for requirements on the matrix.
  * @author Ralf Hartmann, Guido Kanschat, 1999, 2000
  */
 template<typename number, typename inverse_type = number>
 class PreconditionBlockSOR : public  Subscriptor,
-                            private PreconditionBlock<number,inverse_type>
+                            protected PreconditionBlock<number,inverse_type>
 {
   public:
-                                    /**
-                                     * Constructor. Takes the damping
-                                     * Parameter as 
-                                     */
-    PreconditionBlockSOR (const number omega=1.);
-    
-                                    /**
-                                     * Destructor.
-                                     */
-    ~PreconditionBlockSOR();
-
                                     /**
                                      * Make initialization function
                                      * publicly available.
@@ -346,11 +372,6 @@ class PreconditionBlockSOR : public  Subscriptor,
                                      */
     PreconditionBlock<number,inverse_type>::invert_diagblocks;
 
-                                    /**
-                                     * Set the relaxation parameter.
-                                     */
-    void set_omega(number omega);
-    
                                     /**
                                      * Execute block SOR
                                      * preconditioning.
@@ -375,12 +396,68 @@ class PreconditionBlockSOR : public  Subscriptor,
                                      */
     template <typename number2>
     void Tvmult (Vector<number2>&, const Vector<number2>&) const;
+};
 
-  private:
+
+/**
+ * Block SSOR preconditioning.
+ *
+ * The functions @p{vmult} and @p{Tvmult} execute a block-SSOR step,
+ * based on the implementation in @ref{PreconditionBlockSOR}.  This
+ * class requires storage of the diagonal blocks and their inverses.
+ *
+ * See @ref{PreconditionBlock} for requirements on the matrix.
+ * @author Ralf Hartmann, Guido Kanschat, 1999, 2000
+ */
+template<typename number, typename inverse_type = number>
+class PreconditionBlockSSOR : private PreconditionBlockSOR<number,inverse_type>
+{
+  public:
                                     /**
-                                     * Relaxation parameter.
+                                     * Constructor.
+                                     */
+    PreconditionBlockSSOR ();
+                                    /**
+                                     * Make initialization function
+                                     * publicly available.
+                                     */
+    PreconditionBlock<number,inverse_type>::initialize;
+    
+                                    /**
+                                     * Make function of base class public again.
+                                     */
+    PreconditionBlock<number,inverse_type>::clear;
+
+                                    /**
+                                     * Make function of base class public again.
+                                     */
+    PreconditionBlock<number,inverse_type>::set_same_diagonal;
+
+                                    /**
+                                     * Make function of base class public again.
                                      */
-    number omega;
+    PreconditionBlock<number,inverse_type>::invert_diagblocks;
+
+                                    /**
+                                     * Execute block SSOR
+                                     * preconditioning.
+                                     *
+                                     * 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.
+                                     */
+    template <typename number2>
+    void vmult (Vector<number2>&, const Vector<number2>&) const;
+
+                                    /**
+                                     * Same as @ref{vmult}
+                                     */
+    template <typename number2>
+    void Tvmult (Vector<number2>&, const Vector<number2>&) const;
 };
 
 
index 6d53a57a2a031e483909ec8692154d63a92cf7b7..42aafca092f14bbb31bd76e2cd4c2842601f8539 100644 (file)
@@ -26,6 +26,7 @@ template <typename number, typename inverse_type>
 PreconditionBlock<number,inverse_type>::PreconditionBlock ():
                blocksize(0),
                A(0),
+               store_diagonals(false),
                same_diagonal(false)
 {};
 
@@ -33,16 +34,20 @@ PreconditionBlock<number,inverse_type>::PreconditionBlock ():
 template <typename number, typename inverse_type>
 PreconditionBlock<number,inverse_type>::~PreconditionBlock ()
 {
-  if (_inverse.size()!=0)
-    _inverse.erase(_inverse.begin(), _inverse.end());
+  if (var_inverse.size()!=0)
+    var_inverse.erase(var_inverse.begin(), var_inverse.end());
+  if (var_diagonal.size()!=0)
+    var_diagonal.erase(var_diagonal.begin(), var_diagonal.end());
 }
 
 
 template <typename number, typename inverse_type>
 void PreconditionBlock<number,inverse_type>::clear ()
 {
-  if (_inverse.size()!=0)
-    _inverse.erase(_inverse.begin(), _inverse.end());
+  if (var_inverse.size()!=0)
+    var_inverse.erase(var_inverse.begin(), var_inverse.end());
+  if (var_diagonal.size()!=0)
+    var_diagonal.erase(var_diagonal.begin(), var_diagonal.end());
   blocksize     = 0;
   same_diagonal = false;
   A = 0;
@@ -51,7 +56,8 @@ void PreconditionBlock<number,inverse_type>::clear ()
 
 template <typename number, typename inverse_type>
 void PreconditionBlock<number,inverse_type>::initialize (const SparseMatrix<number> &M,
-                                                        unsigned int bsize)
+                                                        unsigned int bsize,
+                                                        const double relax)
 {
   clear();
   Assert (M.m() == M.n(), ExcMatrixNotSquare());
@@ -59,6 +65,7 @@ void PreconditionBlock<number,inverse_type>::initialize (const SparseMatrix<numb
   Assert (bsize>0, ExcIndexRange(bsize, 1, M.m()));
   Assert (A->m()%bsize==0, ExcWrongBlockSize(bsize, A->m()));
   blocksize=bsize;
+  relaxation = relax;
 }
 
 
@@ -67,10 +74,24 @@ const FullMatrix<inverse_type>&
 PreconditionBlock<number,inverse_type>::inverse(unsigned int i) const
 {
   if (same_diagonal)
-    return _inverse[0];
+    return var_inverse[0];
   
-  Assert (i < _inverse.size(), ExcIndexRange(i,0,_inverse.size()));
-  return _inverse[i];
+  Assert (i < var_inverse.size(), ExcIndexRange(i,0,var_inverse.size()));
+  return var_inverse[i];
+}
+
+
+template <typename number, typename inverse_type>
+const FullMatrix<inverse_type>&
+PreconditionBlock<number,inverse_type>::diagonal(unsigned int i) const
+{
+  Assert(store_diagonals, ExcDiagonalsNotStored());
+  
+  if (same_diagonal)
+    return var_diagonal[0];
+  
+  Assert (i < var_diagonal.size(), ExcIndexRange(i,0,var_diagonal.size()));
+  return var_diagonal[i];
 }
 
 
@@ -85,7 +106,7 @@ template <typename number, typename inverse_type>
 void
 PreconditionBlock<number,inverse_type>::set_same_diagonal()
 {
-  Assert(_inverse.size()==0, ExcInverseMatricesAlreadyExist());
+  Assert(var_inverse.size()==0, ExcInverseMatricesAlreadyExist());
   same_diagonal = true;
 }
 
@@ -94,7 +115,7 @@ template <typename number, typename inverse_type>
 bool
 PreconditionBlock<number,inverse_type>::inverses_ready() const
 {
-  return (_inverse.size() != 0);
+  return (var_inverse.size() != 0);
 }
 
 
@@ -105,7 +126,7 @@ void PreconditionBlock<number,inverse_type>::invert_diagblocks()
   Assert (blocksize!=0, ExcNotInitialized());
 
   const SparseMatrix<number> &M=*A;
-  Assert (_inverse.size()==0, ExcInverseMatricesAlreadyExist());
+  Assert (var_inverse.size()==0, ExcInverseMatricesAlreadyExist());
 
   const unsigned int n_cells = M.m()/blocksize;
 
@@ -117,13 +138,17 @@ void PreconditionBlock<number,inverse_type>::invert_diagblocks()
                                       // Invert only the first block
                                       // This is a copy of the code in the
                                       // 'else' part, stripped of the outer loop
-      _inverse.resize(1);
+      if (store_diagonals)
+       var_diagonal.resize(1);
+      var_inverse.resize(1);
 
       for (unsigned int row_cell=0; row_cell<blocksize; ++row_cell)
        for (unsigned int column_cell=0; column_cell<blocksize; ++column_cell)
          M_cell(row_cell,column_cell)=M.el(row_cell,column_cell);
 
-      _inverse[0].invert(M_cell);
+      if (store_diagonals)
+       var_diagonal[0] = M_cell;
+      var_inverse[0].invert(M_cell);
     }
   else
     {
@@ -135,17 +160,24 @@ void PreconditionBlock<number,inverse_type>::invert_diagblocks()
                                       // row, column are the global numbering
                                       // of the unkowns.
 
-                                      // set the @p{_inverse} array to the right
+                                      // set the @p{var_inverse} array to the right
                                       // size. we could do it like this:
-                                      // _inverse = vector<>(n_cells,FullMatrix<>())
+                                      // var_inverse = vector<>(n_cells,FullMatrix<>())
                                       // but this would involve copying many
                                       // FullMatrix objects.
                                       //
                                       // the following is a neat trick which
                                       // avoids copying
+      if (store_diagonals)
+       {
+         vector<FullMatrix<inverse_type> > tmp(n_cells,
+                                               FullMatrix<inverse_type>(blocksize));
+         var_diagonal.swap (tmp);
+       }
+      
       vector<FullMatrix<inverse_type> > tmp(n_cells,
                                            FullMatrix<inverse_type>(blocksize));
-      _inverse.swap (tmp);
+      var_inverse.swap (tmp);
 
       M_cell.clear ();
       
@@ -156,7 +188,9 @@ void PreconditionBlock<number,inverse_type>::invert_diagblocks()
                 column_cell<blocksize; ++column_cell, ++column)
              M_cell(row_cell,column_cell)=M.el(row,column);
 
-         _inverse[cell].invert(M_cell);
+         if (store_diagonals)
+           var_diagonal[cell] = M_cell;
+         var_inverse[cell].invert(M_cell);
        }
     }
 }
@@ -239,6 +273,7 @@ void PreconditionBlockJacobi<number,inverse_type>
        
        begin_diag_block+=blocksize;
       }
+  dst.scale(relaxation);
 }
 
 
@@ -254,27 +289,6 @@ void PreconditionBlockJacobi<number,inverse_type>
 
 /*--------------------- 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>
-void
-PreconditionBlockSOR<number,inverse_type>::set_omega(number om)
-{
-  omega = om;
-}
-
-
-
 template <typename number, typename inverse_type>
 template <typename number2>
 void PreconditionBlockSOR<number,inverse_type>::vmult (Vector<number2>       &dst,
@@ -339,7 +353,7 @@ void PreconditionBlockSOR<number,inverse_type>::vmult (Vector<number2>       &ds
          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)=omega*x_cell(row_cell);
+           dst(row)=relaxation*x_cell(row_cell);
          
          begin_diag_block+=blocksize;
        }
@@ -360,7 +374,7 @@ void PreconditionBlockSOR<number,inverse_type>::vmult (Vector<number2>       &ds
        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)=omega*x_cell(row_cell);
+         dst(row)=relaxation*x_cell(row_cell);
        
        begin_diag_block+=blocksize;
       }
@@ -438,7 +452,7 @@ void PreconditionBlockSOR<number,inverse_type>::Tvmult (Vector<number2>       &d
          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)=omega*x_cell(row_cell);
+           dst(row)=relaxation*x_cell(row_cell);
          
          end_diag_block-=blocksize;
        }
@@ -461,11 +475,84 @@ void PreconditionBlockSOR<number,inverse_type>::Tvmult (Vector<number2>       &d
        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)=omega*x_cell(row_cell);
+         dst(row)=relaxation*x_cell(row_cell);
        
        end_diag_block-=blocksize;
       }
 }
 
+//----------------------------------------------------------------------//
+
+
+template <typename number, typename inverse_type>
+PreconditionBlockSSOR<number,inverse_type>::PreconditionBlockSSOR ()
+{
+  store_diagonals = 1;
+}
+
+
+template <typename number, typename inverse_type>
+template <typename number2>
+void PreconditionBlockSSOR<number,inverse_type>::vmult (Vector<number2>       &dst,
+                                                       const Vector<number2> &src) const
+{
+  Vector<number2> help;
+  help.reinit(dst);
+  
+  PreconditionBlockSOR<number,inverse_type>::vmult(help, src);
+
+  Vector<inverse_type> cell_src(blocksize);
+  Vector<inverse_type> cell_dst(blocksize);
+  
+  const unsigned int n_cells = A->m()/blocksize;
+
+                                  // Multiply with diagonal blocks
+  for (unsigned int cell=0; cell<n_cells; ++cell)
+    {
+      unsigned int row = cell*blocksize;
+      
+      for (unsigned int row_cell=0; row_cell<blocksize; ++row_cell)
+       cell_src(row_cell)=help(row+row_cell);
+
+      diagonal(cell).vmult(cell_dst, cell_src);
+
+      for (unsigned int row_cell=0; row_cell<blocksize; ++row_cell)
+       help(row+row_cell)=relaxation*cell_dst(row_cell);
+    }
+  
+  PreconditionBlockSOR<number,inverse_type>::Tvmult(dst, help);
+}
+
+template <typename number, typename inverse_type>
+template <typename number2>
+void PreconditionBlockSSOR<number,inverse_type>::Tvmult (Vector<number2>       &dst,
+                                                       const Vector<number2> &src) const
+{
+  Vector<number2> help;
+  help.reinit(dst);
+  
+  PreconditionBlockSOR<number,inverse_type>::Tvmult(help, src);
+
+  Vector<inverse_type> cell_src(blocksize);
+  Vector<inverse_type> cell_dst(blocksize);
+  
+  const unsigned int n_cells = A->m()/blocksize;
+
+                                  // Multiply with diagonal blocks
+  for (unsigned int cell=0; cell<n_cells; ++cell)
+    {
+      unsigned int row = cell*blocksize;
+      
+      for (unsigned int row_cell=0; row_cell<blocksize; ++row_cell)
+       cell_src(row_cell)=help(row+row_cell);
+
+      diagonal(cell).vmult(cell_dst, cell_src);
+
+      for (unsigned int row_cell=0; row_cell<blocksize; ++row_cell)
+       help(row+row_cell)=relaxation*cell_dst(row_cell);
+    }
+  
+  PreconditionBlockSOR<number,inverse_type>::vmult(dst, help);
+}
 
 #endif
index eed558300dd1707009f1ec94309d187025bad751..0427e542e8a7f00decdd9674c2c664f2527b1e6d 100644 (file)
@@ -112,3 +112,49 @@ template void PreconditionBlockSOR<double, double>::Tvmult (
   Vector<double> &, const Vector<double> &) const;
 
 
+/*--------------------- PreconditionBlockSSOR -----------------------*/
+
+
+// explicit instantiations for "float" PreconditionBlock
+template class PreconditionBlockSSOR<float, float>;
+
+template void PreconditionBlockSSOR<float, float>::vmult (
+  Vector<float> &, const Vector<float> &) const;
+template void PreconditionBlockSSOR<float, float>::vmult (
+  Vector<double> &, const Vector<double> &) const;
+template void PreconditionBlockSSOR<float, float>::Tvmult (
+  Vector<float> &, const Vector<float> &) const;
+template void PreconditionBlockSSOR<float, float>::Tvmult (
+  Vector<double> &, const Vector<double> &) const;
+
+
+// the instantiation for class PreconditionBlockSSOR<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" PreconditionBlockSSOR
+template class PreconditionBlockSSOR<double, float>;
+
+
+template void PreconditionBlockSSOR<double, float>::vmult (
+  Vector<float> &, const Vector<float> &) const;
+template void PreconditionBlockSSOR<double, float>::vmult (
+  Vector<double> &, const Vector<double> &) const;
+template void PreconditionBlockSSOR<double, float>::Tvmult (
+  Vector<float> &, const Vector<float> &) const;
+template void PreconditionBlockSSOR<double, float>::Tvmult (
+  Vector<double> &, const Vector<double> &) const;
+
+template class PreconditionBlockSSOR<double, double>;
+
+template void PreconditionBlockSSOR<double, double>::vmult (
+  Vector<float> &, const Vector<float> &) const;
+template void PreconditionBlockSSOR<double, double>::vmult (
+  Vector<double> &, const Vector<double> &) const;
+template void PreconditionBlockSSOR<double, double>::Tvmult (
+  Vector<float> &, const Vector<float> &) const;
+template void PreconditionBlockSSOR<double, double>::Tvmult (
+  Vector<double> &, const Vector<double> &) const;
+
+

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.