]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
PreconditionBlock... operates on SparseMatrixEZ
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 18 Sep 2002 14:58:47 +0000 (14:58 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 18 Sep 2002 14:58:47 +0000 (14:58 +0000)
git-svn-id: https://svn.dealii.org/trunk@6451 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/2002/c-3-4.html
deal.II/lac/include/lac/precondition_block.h
deal.II/lac/include/lac/precondition_block.templates.h
deal.II/lac/include/lac/sparse_matrix_ez.h
deal.II/lac/source/precondition_block.cc
tests/lac/sparse_matrices.cc

index c85c2133730ce9169ab1ec560bb8338f335e526e..51bbb5f41481cb666addb632f550055cbda76d4b 100644 (file)
@@ -234,6 +234,17 @@ contributor's names are abbreviated by WB (Wolfgang Bangerth), GK
 <h3>lac</h3>
 
 <ol>
+  <li> <p> Changed: The classes <code
+  class="class">PreconditionBlock</code>, <code
+  class="class">PreconditionBlockJacobi</code>, <code
+  class="class">PreconditionBlockSOR</code>, and <code
+  class="class">PreconditionBlockSSOR</code> have changed their
+  template signature. The first template argument is now the matrix
+  type, not just a number type.
+  <br>
+  (GK 2002/09/18)
+  </p>
+
   <li> <p> New: class <code class="class">BlockVector</code> has a
   function <code class="member">collect_sizes()</code>, very much as
   <code class="class">BlockSparsityPattern</code>. This allows
index c00ae7f56429492b1f1ebc46ca70663b0a2ab292..b251d65737d20b8a1bf604925f6fc46e2c7baba7 100644 (file)
@@ -75,9 +75,15 @@ template <typename number> class Vector;
  *
  * @author Ralf Hartmann, Guido Kanschat, 1999, 2000
  */
-template<typename number, typename inverse_type = number>
+template<class MATRIX, typename inverse_type = typename MATRIX::value_type>
 class PreconditionBlock : public virtual Subscriptor
 {
+  private:
+                                    /**
+                                     * Define number type of matrix.
+                                     */
+    typedef typename MATRIX::value_type number;
+    
   public:
                                     /**
                                      * Constructor.
@@ -102,7 +108,7 @@ class PreconditionBlock : public virtual Subscriptor
                                      * parameter for derived classes
                                      * may be provided.
                                      */
-    void initialize (const SparseMatrix<number>& A,
+    void initialize (const MATRIX& A,
                     const unsigned int block_size,
                     const double relaxation = 1.);
 
@@ -243,7 +249,7 @@ class PreconditionBlock : public virtual Subscriptor
                                      * @p{vmult} function of the
                                      * derived classes.
                                      */
-    SmartPointer<const SparseMatrix<number> > A;
+    SmartPointer<const MATRIX> A;
                                     /**
                                      * Relaxation parameter to be
                                      * used by derived classes.
@@ -293,31 +299,37 @@ class PreconditionBlock : public virtual Subscriptor
  * See @ref{PreconditionBlock} for requirements on the matrix.
  * @author Ralf Hartmann, Guido Kanschat, 1999, 2000
  */
-template<typename number, typename inverse_type = number>
+template<class MATRIX, typename inverse_type = typename MATRIX::value_type>
 class PreconditionBlockJacobi : public virtual Subscriptor,
-                               private PreconditionBlock<number,inverse_type>
+                               private PreconditionBlock<MATRIX, inverse_type>
 {
+  private:
+                                    /**
+                                     * Define number type of matrix.
+                                     */
+    typedef typename MATRIX::value_type number;
+    
   public:
                                     /**
                                      * Make initialization function
                                      * publicly available.
                                      */
-    PreconditionBlock<number,inverse_type>::initialize;
+    PreconditionBlock<MATRIX, inverse_type>::initialize;
     
                                     /**
                                      * Make function of base class public again.
                                      */
-    PreconditionBlock<number,inverse_type>::clear;
+    PreconditionBlock<MATRIX, inverse_type>::clear;
 
                                     /**
                                      * Make function of base class public again.
                                      */
-    PreconditionBlock<number,inverse_type>::set_same_diagonal;
+    PreconditionBlock<MATRIX, inverse_type>::set_same_diagonal;
 
                                     /**
                                      * Make function public.
                                      */
-    PreconditionBlock<number,inverse_type>::invert_diagblocks;
+    PreconditionBlock<MATRIX, inverse_type>::invert_diagblocks;
     
                                     /**
                                      * Execute block Jacobi
@@ -384,31 +396,37 @@ class PreconditionBlockJacobi : public virtual Subscriptor,
  * See @ref{PreconditionBlock} for requirements on the matrix.
  * @author Ralf Hartmann, Guido Kanschat, 1999, 2000
  */
-template<typename number, typename inverse_type = number>
+template<class MATRIX, typename inverse_type = typename MATRIX::value_type>
 class PreconditionBlockSOR : public virtual Subscriptor,
-                            protected PreconditionBlock<number,inverse_type>
+                            protected PreconditionBlock<MATRIX, inverse_type>
 {
+  private:
+                                    /**
+                                     * Define number type of matrix.
+                                     */
+    typedef typename MATRIX::value_type number;
+    
   public:
                                     /**
                                      * Make initialization function
                                      * publicly available.
                                      */
-    PreconditionBlock<number,inverse_type>::initialize;
+    PreconditionBlock<MATRIX, inverse_type>::initialize;
     
                                     /**
                                      * Make function of base class public again.
                                      */
-    PreconditionBlock<number,inverse_type>::clear;
+    PreconditionBlock<MATRIX, inverse_type>::clear;
 
                                     /**
                                      * Make function of base class public again.
                                      */
-    PreconditionBlock<number,inverse_type>::set_same_diagonal;
+    PreconditionBlock<MATRIX, inverse_type>::set_same_diagonal;
 
                                     /**
                                      * Make function of base class public again.
                                      */
-    PreconditionBlock<number,inverse_type>::invert_diagblocks;
+    PreconditionBlock<MATRIX, inverse_type>::invert_diagblocks;
 
                                     /**
                                      * Execute block SOR
@@ -447,10 +465,16 @@ class PreconditionBlockSOR : public virtual Subscriptor,
  * See @ref{PreconditionBlock} for requirements on the matrix.
  * @author Ralf Hartmann, Guido Kanschat, 1999, 2000
  */
-template<typename number, typename inverse_type = number>
+template<class MATRIX, typename inverse_type = typename MATRIX::value_type>
 class PreconditionBlockSSOR : public virtual Subscriptor,
-                             private PreconditionBlockSOR<number,inverse_type>
+                             private PreconditionBlockSOR<MATRIX, inverse_type>
 {
+  private:
+                                    /**
+                                     * Define number type of matrix.
+                                     */
+    typedef typename MATRIX::value_type number;
+    
   public:
                                     /**
                                      * Constructor.
@@ -460,22 +484,22 @@ class PreconditionBlockSSOR : public virtual Subscriptor,
                                      * Make initialization function
                                      * publicly available.
                                      */
-    PreconditionBlockSOR<number,inverse_type>::initialize;
+    PreconditionBlockSOR<MATRIX,inverse_type>::initialize;
     
                                     /**
                                      * Make function of base class public again.
                                      */
-    PreconditionBlockSOR<number,inverse_type>::clear;
+    PreconditionBlockSOR<MATRIX,inverse_type>::clear;
 
                                     /**
                                      * Make function of base class public again.
                                      */
-    PreconditionBlockSOR<number,inverse_type>::set_same_diagonal;
+    PreconditionBlockSOR<MATRIX,inverse_type>::set_same_diagonal;
 
                                     /**
                                      * Make function of base class public again.
                                      */
-    PreconditionBlockSOR<number,inverse_type>::invert_diagblocks;
+    PreconditionBlockSOR<MATRIX,inverse_type>::invert_diagblocks;
 
                                     /**
                                      * Execute block SSOR
index 8ea6a1ea4268ccbc09bb25e5826dfe63298ea259..0def0191bea6846b70d3cbd57a1c124ceb6937d7 100644 (file)
@@ -24,8 +24,8 @@
 #include <lac/sparse_matrix.h>
 
 
-template <typename number, typename inverse_type>
-PreconditionBlock<number,inverse_type>::PreconditionBlock ():
+template <class MATRIX, typename inverse_type>
+PreconditionBlock<MATRIX,inverse_type>::PreconditionBlock ():
                blocksize(0),
                A(0),
                store_diagonals(false),
@@ -33,8 +33,8 @@ PreconditionBlock<number,inverse_type>::PreconditionBlock ():
 {};
 
 
-template <typename number, typename inverse_type>
-PreconditionBlock<number,inverse_type>::~PreconditionBlock ()
+template <class MATRIX, typename inverse_type>
+PreconditionBlock<MATRIX,inverse_type>::~PreconditionBlock ()
 {
   if (var_inverse.size()!=0)
     var_inverse.erase(var_inverse.begin(), var_inverse.end());
@@ -43,8 +43,8 @@ PreconditionBlock<number,inverse_type>::~PreconditionBlock ()
 }
 
 
-template <typename number, typename inverse_type>
-void PreconditionBlock<number,inverse_type>::clear ()
+template <class MATRIX, typename inverse_type>
+void PreconditionBlock<MATRIX,inverse_type>::clear ()
 {
   if (var_inverse.size()!=0)
     var_inverse.erase(var_inverse.begin(), var_inverse.end());
@@ -56,8 +56,8 @@ void PreconditionBlock<number,inverse_type>::clear ()
 }
 
 
-template <typename number, typename inverse_type>
-void PreconditionBlock<number,inverse_type>::initialize (const SparseMatrix<number> &M,
+template <class MATRIX, typename inverse_type>
+void PreconditionBlock<MATRIX,inverse_type>::initialize (const MATRIX &M,
                                                         const unsigned int bsize,
                                                         const double relax)
 {
@@ -71,9 +71,9 @@ void PreconditionBlock<number,inverse_type>::initialize (const SparseMatrix<numb
 }
 
 
-template <typename number, typename inverse_type>
+template <class MATRIX, typename inverse_type>
 const FullMatrix<inverse_type>&
-PreconditionBlock<number,inverse_type>::inverse(unsigned int i) const
+PreconditionBlock<MATRIX,inverse_type>::inverse(unsigned int i) const
 {
   if (same_diagonal)
     return var_inverse[0];
@@ -83,9 +83,9 @@ PreconditionBlock<number,inverse_type>::inverse(unsigned int i) const
 }
 
 
-template <typename number, typename inverse_type>
+template <class MATRIX, typename inverse_type>
 const FullMatrix<inverse_type>&
-PreconditionBlock<number,inverse_type>::diagonal(unsigned int i) const
+PreconditionBlock<MATRIX,inverse_type>::diagonal(unsigned int i) const
 {
   Assert(store_diagonals, ExcDiagonalsNotStored());
   
@@ -97,37 +97,37 @@ PreconditionBlock<number,inverse_type>::diagonal(unsigned int i) const
 }
 
 
-template <typename number, typename inverse_type>
-unsigned int PreconditionBlock<number,inverse_type>::block_size() const
+template <class MATRIX, typename inverse_type>
+unsigned int PreconditionBlock<MATRIX,inverse_type>::block_size() const
 {
   return blocksize;
 }
 
 
-template <typename number, typename inverse_type>
+template <class MATRIX, typename inverse_type>
 void
-PreconditionBlock<number,inverse_type>::set_same_diagonal()
+PreconditionBlock<MATRIX,inverse_type>::set_same_diagonal()
 {
   Assert(var_inverse.size()==0, ExcInverseMatricesAlreadyExist());
   same_diagonal = true;
 }
 
 
-template <typename number, typename inverse_type>
+template <class MATRIX, typename inverse_type>
 bool
-PreconditionBlock<number,inverse_type>::inverses_ready() const
+PreconditionBlock<MATRIX,inverse_type>::inverses_ready() const
 {
   return (var_inverse.size() != 0);
 }
 
 
-template <typename number, typename inverse_type>
-void PreconditionBlock<number,inverse_type>::invert_diagblocks()
+template <class MATRIX, typename inverse_type>
+void PreconditionBlock<MATRIX,inverse_type>::invert_diagblocks()
 {
   Assert (A!=0, ExcNotInitialized());
   Assert (blocksize!=0, ExcNotInitialized());
 
-  const SparseMatrix<number> &M=*A;
+  const MATRIX &M=*A;
   Assert (var_inverse.size()==0, ExcInverseMatricesAlreadyExist());
 
   const unsigned int n_cells = M.m()/blocksize;
@@ -147,11 +147,12 @@ void PreconditionBlock<number,inverse_type>::invert_diagblocks()
 
       for (unsigned int row_cell=0; row_cell<blocksize; ++row_cell)
        {
-         SparseMatrix<number>::const_iterator entry = M.begin(row_cell);
-         const SparseMatrix<number>::const_iterator row_end = M.end(row_cell);
-         while(entry != row_end && entry->column() < blocksize)
+         typename MATRIX::const_iterator entry = M.begin(row_cell);
+         const typename MATRIX::const_iterator row_end = M.end(row_cell);
+         while(entry != row_end)
            {
-             M_cell(row_cell, entry->column()) = entry->value();
+             if (entry->column() < blocksize)
+               M_cell(row_cell, entry->column()) = entry->value();
              ++entry;
            }
        }
@@ -202,18 +203,18 @@ void PreconditionBlock<number,inverse_type>::invert_diagblocks()
          for (unsigned int row_cell=0; row_cell<blocksize; ++row_cell)
            {
              const unsigned int row = row_cell + cell_start;
-             SparseMatrix<number>::const_iterator entry = M.begin(row);
-             const SparseMatrix<number>::const_iterator row_end = M.end(row);
+             typename MATRIX::const_iterator entry = M.begin(row);
+             const typename MATRIX::const_iterator row_end = M.end(row);
 
-             while (entry != row_end && entry->column()<cell_start)
-               ++entry;
-             while(entry != row_end)
+             for (;entry != row_end; ++entry)
                {
+                 if (entry->column()<cell_start)
+                   continue;
+                 
                  const unsigned int column_cell = entry->column()-cell_start;
                  if (column_cell >= blocksize)
-                   break;
+                   continue;
                  M_cell(row_cell, column_cell) = entry->value();
-                 ++entry;
                }
            }
 
@@ -226,9 +227,9 @@ void PreconditionBlock<number,inverse_type>::invert_diagblocks()
 
 
 
-template <typename number, typename inverse_type>
+template <class MATRIX, typename inverse_type>
 unsigned int
-PreconditionBlock<number,inverse_type>::memory_consumption () const
+PreconditionBlock<MATRIX,inverse_type>::memory_consumption () const
 {
   unsigned int mem = sizeof(*this);
   for (unsigned int i=0; i<var_inverse.size(); ++i)
@@ -244,9 +245,9 @@ PreconditionBlock<number,inverse_type>::memory_consumption () const
 /*--------------------- PreconditionBlockJacobi -----------------------*/
 
 
-template <typename number, typename inverse_type>
+template <class MATRIX, typename inverse_type>
 template <typename number2>
-void PreconditionBlockJacobi<number,inverse_type>
+void PreconditionBlockJacobi<MATRIX,inverse_type>
 ::do_vmult (Vector<number2>       &dst,
            const Vector<number2> &src,
            bool adding) const
@@ -267,10 +268,10 @@ void PreconditionBlockJacobi<number,inverse_type>
                                   // parse error for gcc for the
                                   // exceptions that do not take
                                   // args...
-  typedef PreconditionBlock<number,inverse_type> BaseClass;
+  typedef PreconditionBlock<MATRIX,inverse_type> BaseClass;
   Assert(this->A!=0, ExcNotInitialized());
   
-  const SparseMatrix<number> &M=*this->A;
+  const MATRIX &M=*this->A;
   const unsigned int n_cells=M.m()/this->blocksize;
 
   Vector<number2> b_cell(this->blocksize), x_cell(this->blocksize);
@@ -337,9 +338,9 @@ void PreconditionBlockJacobi<number,inverse_type>
 }
 
 
-template <typename number, typename inverse_type>
+template <class MATRIX, typename inverse_type>
 template <typename number2>
-void PreconditionBlockJacobi<number,inverse_type>
+void PreconditionBlockJacobi<MATRIX,inverse_type>
 ::vmult (Vector<number2>       &dst,
         const Vector<number2> &src) const
 {
@@ -347,9 +348,9 @@ void PreconditionBlockJacobi<number,inverse_type>
 }
 
 
-template <typename number, typename inverse_type>
+template <class MATRIX, typename inverse_type>
 template <typename number2>
-void PreconditionBlockJacobi<number,inverse_type>
+void PreconditionBlockJacobi<MATRIX,inverse_type>
 ::Tvmult (Vector<number2>       &dst,
          const Vector<number2> &src) const
 {
@@ -357,9 +358,9 @@ void PreconditionBlockJacobi<number,inverse_type>
 }
 
 
-template <typename number, typename inverse_type>
+template <class MATRIX, typename inverse_type>
 template <typename number2>
-void PreconditionBlockJacobi<number,inverse_type>
+void PreconditionBlockJacobi<MATRIX,inverse_type>
 ::vmult_add (Vector<number2>       &dst,
             const Vector<number2> &src) const
 {
@@ -367,9 +368,9 @@ void PreconditionBlockJacobi<number,inverse_type>
 }
 
 
-template <typename number, typename inverse_type>
+template <class MATRIX, typename inverse_type>
 template <typename number2>
-void PreconditionBlockJacobi<number,inverse_type>
+void PreconditionBlockJacobi<MATRIX,inverse_type>
 ::Tvmult_add (Vector<number2>       &dst,
              const Vector<number2> &src) const
 {
@@ -379,9 +380,9 @@ void PreconditionBlockJacobi<number,inverse_type>
 
 /*--------------------- PreconditionBlockSOR -----------------------*/
 
-template <typename number, typename inverse_type>
+template <class MATRIX, typename inverse_type>
 template <typename number2>
-void PreconditionBlockSOR<number,inverse_type>::vmult (Vector<number2>       &dst,
+void PreconditionBlockSOR<MATRIX,inverse_type>::vmult (Vector<number2>       &dst,
                                                       const Vector<number2> &src) const
 {
                                   // introduce the following typedef
@@ -400,17 +401,13 @@ void PreconditionBlockSOR<number,inverse_type>::vmult (Vector<number2>       &ds
                                   // parse error for gcc for the
                                   // exceptions that do not take
                                   // args...
-  typedef PreconditionBlock<number,inverse_type> BaseClass;
+  typedef PreconditionBlock<MATRIX,inverse_type> BaseClass;
 
   Assert (this->A!=0, ExcNotInitialized());
   
-  const SparseMatrix<number> &M=*this->A;
+  const MATRIX &M=*this->A;
   const unsigned int n_cells=M.m()/this->blocksize;
 
-  const SparsityPattern &spars    = M.get_sparsity_pattern();
-  const unsigned int    *rowstart = spars.get_rowstart_indices();
-  const unsigned int    *columns  = spars.get_column_numbers();
-
   Vector<number2> b_cell(this->blocksize), x_cell(this->blocksize);
 
                                       // cell_row, cell_column are the
@@ -420,7 +417,7 @@ void PreconditionBlockSOR<number,inverse_type>::vmult (Vector<number2>       &ds
                                       // blocks.
                                       // row, column are the global numbering
                                       // of the unkowns.
-  unsigned int row, column, row_cell, begin_diag_block=0;
+  unsigned int row, column, row_cell, block_start=0;
   number2 b_cell_row;
 
   if (!this->inverses_ready())
@@ -428,44 +425,55 @@ void PreconditionBlockSOR<number,inverse_type>::vmult (Vector<number2>       &ds
       FullMatrix<number> M_cell(this->blocksize);
       for (unsigned int cell=0; cell<n_cells; ++cell)
        {
-         for (row=cell*this->blocksize, row_cell=0;
-              row_cell<this->blocksize;
+         for (row = block_start, row_cell = 0;
+              row_cell < this->blocksize;
               ++row_cell, ++row)
            {
+             const typename MATRIX::const_iterator row_end = M.end(row);
+             typename MATRIX::const_iterator entry = M.begin(row);
+             
              b_cell_row=src(row);
-             for (unsigned int j=rowstart[row]; j<rowstart[row+1]; ++j)
-               if ((column=columns[j]) < begin_diag_block)
-                   b_cell_row -= M.global_entry(j) * dst(column);
+             for (; entry != row_end; ++entry)
+               {
+                 const unsigned int column = entry->column();
+                 if (column < block_start)
+                   b_cell_row -= entry->value() * dst(column);
+                 else if (column < block_start + this->blocksize)
+                   {
+                     const unsigned int column_cell = column - block_start;
+                     M_cell(row_cell, column_cell) = entry->value();
+                   }
+               }
              b_cell(row_cell)=b_cell_row;
-             for (unsigned int column_cell=0, column=cell*this->blocksize;
-                  column_cell<this->blocksize;
-                  ++column_cell, ++column)
-                 M_cell(row_cell,column_cell)=M(row,column);
            }
          M_cell.householder(b_cell);
          M_cell.backward(x_cell,b_cell);
                                           // distribute x_cell to dst
-         for (row=cell*this->blocksize, row_cell=0;
+         for (row=block_start, row_cell=0;
               row_cell<this->blocksize;
               ++row_cell, ++row)
            dst(row)=this->relaxation*x_cell(row_cell);
          
-         begin_diag_block+=this->blocksize;
+         block_start+=this->blocksize;
        }
     }
   else
     for (unsigned int cell=0; cell<n_cells; ++cell)
       {
-       for (row=cell*this->blocksize, row_cell=0;
+       for (row = block_start, row_cell = 0;
             row_cell<this->blocksize;
             ++row_cell, ++row)
          {
+           const typename MATRIX::const_iterator row_end = M.end(row);
+           typename MATRIX::const_iterator entry = M.begin(row);
+           
            b_cell_row=src(row);
-           for (unsigned int j=rowstart[row]; j<rowstart[row+1]; ++j)
-             if ((column=columns[j]) < begin_diag_block)
-               {
-                 b_cell_row -= M.global_entry(j) * dst(column);
-               }
+           for (; entry != row_end; ++entry)
+             {
+               const unsigned int column = entry->column();
+               if (column < block_start)
+                 b_cell_row -= entry->value() * dst(column);
+             }               
            b_cell(row_cell)=b_cell_row;
          }
        this->inverse(cell).vmult(x_cell, b_cell);
@@ -475,14 +483,14 @@ void PreconditionBlockSOR<number,inverse_type>::vmult (Vector<number2>       &ds
             ++row_cell, ++row)
          dst(row)=this->relaxation*x_cell(row_cell);
        
-       begin_diag_block+=this->blocksize;
+       block_start+=this->blocksize;
       }
 }
 
 
-template <typename number, typename inverse_type>
+template <class MATRIX, typename inverse_type>
 template <typename number2>
-void PreconditionBlockSOR<number,inverse_type>::Tvmult (Vector<number2>       &dst,
+void PreconditionBlockSOR<MATRIX,inverse_type>::Tvmult (Vector<number2>       &dst,
                                                       const Vector<number2> &src) const
 {
                                   // introduce the following typedef
@@ -501,17 +509,13 @@ void PreconditionBlockSOR<number,inverse_type>::Tvmult (Vector<number2>       &d
                                   // parse error for gcc for the
                                   // exceptions that do not take
                                   // args...
-  typedef PreconditionBlock<number,inverse_type> BaseClass;
+  typedef PreconditionBlock<MATRIX,inverse_type> BaseClass;
 
   Assert (this->A!=0, ExcNotInitialized());
   
-  const SparseMatrix<number> &M=*this->A;
+  const MATRIX &M=*this->A;
   const unsigned int n_cells=M.m()/this->blocksize;
 
-  const SparsityPattern &spars    = M.get_sparsity_pattern();
-  const unsigned int    *rowstart = spars.get_rowstart_indices();
-  const unsigned int    *columns  = spars.get_column_numbers();
-
   Vector<number2> b_cell(this->blocksize), x_cell(this->blocksize);
 
                                       // cell_row, cell_column are the
@@ -521,8 +525,8 @@ void PreconditionBlockSOR<number,inverse_type>::Tvmult (Vector<number2>       &d
                                       // blocks.
                                       // row, column are the global numbering
                                       // of the unkowns.
-  unsigned int row, column, row_cell;
-  unsigned int end_diag_block=this->blocksize *n_cells;
+  unsigned int row, column, row_cell, block_start = 0;
+  unsigned int block_end=this->blocksize *n_cells;
   number2 b_cell_row;
 
   if (!this->inverses_ready())
@@ -532,23 +536,26 @@ void PreconditionBlockSOR<number,inverse_type>::Tvmult (Vector<number2>       &d
        {
          unsigned int cell = (unsigned int) icell;
                                           // Collect upper triangle
-         for (row=cell*this->blocksize, row_cell=0;
+         for (row = block_start, row_cell = 0;
               row_cell<this->blocksize;
               ++row_cell, ++row)
            {
+           const typename MATRIX::const_iterator row_end = M.end(row);
+           typename MATRIX::const_iterator entry = M.begin(row);
+           
              b_cell_row=src(row);
-             for (unsigned int j=rowstart[row]; j<rowstart[row+1]; ++j)
+             for (; entry != row_end; ++entry)
                {
-                 column=columns[j];
-                 if (column >= end_diag_block)
-                   b_cell_row -= M.global_entry(j) * dst(column);
+                 const unsigned int column = entry->column();
+                 if (column >= block_end)
+                   b_cell_row -= entry->value() * dst(column);
+                 else if (column >= block_start)
+                   {
+                     const unsigned int column_cell = column - block_start;
+                     M_cell(row_cell, column_cell) = entry->value();
+                   }
                }
-             
              b_cell(row_cell)=b_cell_row;
-                                              // Collect diagonal block
-             for (unsigned int column_cell=0, column=cell*this->blocksize;
-                  column_cell<this->blocksize; ++column_cell, ++column)
-                 M_cell(row_cell,column_cell)=M(row,column);
            }
          M_cell.householder(b_cell);
          M_cell.backward(x_cell,b_cell);
@@ -558,7 +565,8 @@ void PreconditionBlockSOR<number,inverse_type>::Tvmult (Vector<number2>       &d
               ++row_cell, ++row)
            dst(row)=this->relaxation*x_cell(row_cell);
          
-         end_diag_block-=this->blocksize;
+         block_end   -= this->blocksize;
+         block_start += this->blocksize;
        }
     }
   else
@@ -569,13 +577,16 @@ void PreconditionBlockSOR<number,inverse_type>::Tvmult (Vector<number2>       &d
             row_cell<this->blocksize;
             ++row_cell, ++row)
          {
+           const typename MATRIX::const_iterator row_end = M.end(row);
+           typename MATRIX::const_iterator entry = M.begin(row);
+           
            b_cell_row=src(row);
-           for (unsigned int j=rowstart[row]; j<rowstart[row+1]; ++j)
+           for (; entry != row_end; ++entry)
              {
-               column=columns[j];
-               if (column >= end_diag_block)
-                 b_cell_row -= M.global_entry(j) * dst(column);
-             }
+               const unsigned int column = entry->column();
+               if (column >= block_end)
+                 b_cell_row -= entry->value() * dst(column);
+             }               
            b_cell(row_cell)=b_cell_row;
          }
        this->inverse(cell).vmult(x_cell, b_cell);
@@ -585,29 +596,30 @@ void PreconditionBlockSOR<number,inverse_type>::Tvmult (Vector<number2>       &d
             ++row_cell, ++row)
          dst(row)=this->relaxation*x_cell(row_cell);
        
-       end_diag_block-=this->blocksize;
+       block_end   -= this->blocksize;
+       block_start += this->blocksize;
       }
 }
 
 //----------------------------------------------------------------------//
 
 
-template <typename number, typename inverse_type>
-PreconditionBlockSSOR<number,inverse_type>::PreconditionBlockSSOR ()
+template <class MATRIX, typename inverse_type>
+PreconditionBlockSSOR<MATRIX,inverse_type>::PreconditionBlockSSOR ()
 {
   this->store_diagonals = 1;
 }
 
 
-template <typename number, typename inverse_type>
+template <class MATRIX, typename inverse_type>
 template <typename number2>
-void PreconditionBlockSSOR<number,inverse_type>::vmult (Vector<number2>       &dst,
+void PreconditionBlockSSOR<MATRIX,inverse_type>::vmult (Vector<number2>       &dst,
                                                        const Vector<number2> &src) const
 {
   Vector<number2> help;
   help.reinit(dst);
   
-  PreconditionBlockSOR<number,inverse_type>::vmult(help, src);
+  PreconditionBlockSOR<MATRIX,inverse_type>::vmult(help, src);
 
   Vector<inverse_type> cell_src(this->blocksize);
   Vector<inverse_type> cell_dst(this->blocksize);
@@ -628,18 +640,18 @@ void PreconditionBlockSSOR<number,inverse_type>::vmult (Vector<number2>       &d
        help(row+row_cell)=this->relaxation*cell_dst(row_cell);
     }
   
-  PreconditionBlockSOR<number,inverse_type>::Tvmult(dst, help);
+  PreconditionBlockSOR<MATRIX,inverse_type>::Tvmult(dst, help);
 }
 
-template <typename number, typename inverse_type>
+template <class MATRIX, typename inverse_type>
 template <typename number2>
-void PreconditionBlockSSOR<number,inverse_type>::Tvmult (Vector<number2>       &dst,
+void PreconditionBlockSSOR<MATRIX,inverse_type>::Tvmult (Vector<number2>       &dst,
                                                        const Vector<number2> &src) const
 {
   Vector<number2> help;
   help.reinit(dst);
   
-  PreconditionBlockSOR<number,inverse_type>::Tvmult(help, src);
+  PreconditionBlockSOR<MATRIX,inverse_type>::Tvmult(help, src);
 
   Vector<inverse_type> cell_src(this->blocksize);
   Vector<inverse_type> cell_dst(this->blocksize);
@@ -660,7 +672,7 @@ void PreconditionBlockSSOR<number,inverse_type>::Tvmult (Vector<number2>       &
        help(row+row_cell)=this->relaxation*cell_dst(row_cell);
     }
   
-  PreconditionBlockSOR<number,inverse_type>::vmult(dst, help);
+  PreconditionBlockSOR<MATRIX,inverse_type>::vmult(dst, help);
 }
 
 #endif
index 20053ebdc4dbacdf5acf970892d707c851e46186..1e9f9cd3d910c7e9490d71221ca3daed5f387c87 100644 (file)
@@ -1292,7 +1292,7 @@ inline
 number SparseMatrixEZ<number>::operator() (const unsigned int i,
                                           const unsigned int j) const
 {
-  Entry* entry = locate(i,j);
+  const Entry* entry = locate(i,j);
   if (entry)
     return entry->value;
   Assert(false, ExcInvalidEntry(i,j));
index 3c5a66587199739c826dd757a3292610ca25ee19..744547422dbf7139fab3b9dc36be133feef0035a 100644 (file)
 
 
 // explicit instantiations for "float" PreconditionBlock
-template class PreconditionBlock<float, float>;
+template class PreconditionBlock<SparseMatrix<float>, float>;
 
-// the instantiation for class PreconditionBlock<float, double> is skipped
+// the instantiation for class PreconditionBlock<SparseMatrix<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" PreconditionBlock
-template class PreconditionBlock<double, float>;
+template class PreconditionBlock<SparseMatrix<double>, float>;
 
-template class PreconditionBlock<double, double>;
+template class PreconditionBlock<SparseMatrix<double>, double>;
 
 
 /*--------------------- PreconditionBlockJacobi -----------------------*/
 
 
 // explicit instantiations for "float" PreconditionBlock
-template class PreconditionBlockJacobi<float, float>;
+template class PreconditionBlockJacobi<SparseMatrix<float>, float>;
 
-template void PreconditionBlockJacobi<float, float>::vmult<float>
+template void PreconditionBlockJacobi<SparseMatrix<float>, float>::vmult<float>
 (Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockJacobi<float, float>::vmult<double>
+template void PreconditionBlockJacobi<SparseMatrix<float>, float>::vmult<double>
 (Vector<double> &, const Vector<double> &) const;
-template void PreconditionBlockJacobi<float, float>::Tvmult<float>
+template void PreconditionBlockJacobi<SparseMatrix<float>, float>::Tvmult<float>
 (Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockJacobi<float, float>::Tvmult<double>
+template void PreconditionBlockJacobi<SparseMatrix<float>, float>::Tvmult<double>
 (Vector<double> &, const Vector<double> &) const;
-template void PreconditionBlockJacobi<float, float>::vmult_add<float>
+template void PreconditionBlockJacobi<SparseMatrix<float>, float>::vmult_add<float>
 (Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockJacobi<float, float>::vmult_add<double>
+template void PreconditionBlockJacobi<SparseMatrix<float>, float>::vmult_add<double>
 (Vector<double> &, const Vector<double> &) const;
-template void PreconditionBlockJacobi<float, float>::Tvmult_add<float>
+template void PreconditionBlockJacobi<SparseMatrix<float>, float>::Tvmult_add<float>
 (Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockJacobi<float, float>::Tvmult_add<double>
+template void PreconditionBlockJacobi<SparseMatrix<float>, float>::Tvmult_add<double>
 (Vector<double> &, const Vector<double> &) const;
 
-template class PreconditionBlockJacobi<double, float>;
+template class PreconditionBlockJacobi<SparseMatrix<double>, float>;
 
-template void PreconditionBlockJacobi<double, float>::vmult<float>
+template void PreconditionBlockJacobi<SparseMatrix<double>, float>::vmult<float>
 (Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockJacobi<double, float>::vmult<double>
+template void PreconditionBlockJacobi<SparseMatrix<double>, float>::vmult<double>
 (Vector<double> &, const Vector<double> &) const;
-template void PreconditionBlockJacobi<double, float>::Tvmult<float>
+template void PreconditionBlockJacobi<SparseMatrix<double>, float>::Tvmult<float>
 (Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockJacobi<double, float>::Tvmult<double>
+template void PreconditionBlockJacobi<SparseMatrix<double>, float>::Tvmult<double>
 (Vector<double> &, const Vector<double> &) const;
-template void PreconditionBlockJacobi<double, float>::vmult_add<float>
+template void PreconditionBlockJacobi<SparseMatrix<double>, float>::vmult_add<float>
 (Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockJacobi<double, float>::vmult_add<double>
+template void PreconditionBlockJacobi<SparseMatrix<double>, float>::vmult_add<double>
 (Vector<double> &, const Vector<double> &) const;
-template void PreconditionBlockJacobi<double, float>::Tvmult_add<float>
+template void PreconditionBlockJacobi<SparseMatrix<double>, float>::Tvmult_add<float>
 (Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockJacobi<double, float>::Tvmult_add<double>
+template void PreconditionBlockJacobi<SparseMatrix<double>, float>::Tvmult_add<double>
 (Vector<double> &, const Vector<double> &) const;
 
-template class PreconditionBlockJacobi<double, double>;
+template class PreconditionBlockJacobi<SparseMatrix<double>, double>;
 
-template void PreconditionBlockJacobi<double, double>::vmult<float>
+template void PreconditionBlockJacobi<SparseMatrix<double>, double>::vmult<float>
 (Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockJacobi<double, double>::vmult<double>
+template void PreconditionBlockJacobi<SparseMatrix<double>, double>::vmult<double>
 (Vector<double> &, const Vector<double> &) const;
-template void PreconditionBlockJacobi<double, double>::Tvmult<float>
+template void PreconditionBlockJacobi<SparseMatrix<double>, double>::Tvmult<float>
 (Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockJacobi<double, double>::Tvmult<double>
+template void PreconditionBlockJacobi<SparseMatrix<double>, double>::Tvmult<double>
 (Vector<double> &, const Vector<double> &) const;
-template void PreconditionBlockJacobi<double, double>::vmult_add<float>
+template void PreconditionBlockJacobi<SparseMatrix<double>, double>::vmult_add<float>
 (Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockJacobi<double, double>::vmult_add<double>
+template void PreconditionBlockJacobi<SparseMatrix<double>, double>::vmult_add<double>
 (Vector<double> &, const Vector<double> &) const;
-template void PreconditionBlockJacobi<double, double>::Tvmult_add<float>
+template void PreconditionBlockJacobi<SparseMatrix<double>, double>::Tvmult_add<float>
 (Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockJacobi<double, double>::Tvmult_add<double>
+template void PreconditionBlockJacobi<SparseMatrix<double>, double>::Tvmult_add<double>
 (Vector<double> &, const Vector<double> &) const;
 
 /*--------------------- PreconditionBlockGaussSeidel -----------------------*/
 
 
 // explicit instantiations for "float" PreconditionBlock
-template class PreconditionBlockSOR<float, float>;
+template class PreconditionBlockSOR<SparseMatrix<float>, float>;
 
-template void PreconditionBlockSOR<float, float>::vmult<float> (
+template void PreconditionBlockSOR<SparseMatrix<float>, float>::vmult<float> (
   Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockSOR<float, float>::vmult<double> (
+template void PreconditionBlockSOR<SparseMatrix<float>, float>::vmult<double> (
   Vector<double> &, const Vector<double> &) const;
-template void PreconditionBlockSOR<float, float>::Tvmult<float> (
+template void PreconditionBlockSOR<SparseMatrix<float>, float>::Tvmult<float> (
   Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockSOR<float, float>::Tvmult<double> (
+template void PreconditionBlockSOR<SparseMatrix<float>, float>::Tvmult<double> (
   Vector<double> &, const Vector<double> &) const;
 
 
-// the instantiation for class PreconditionBlockSOR<float, double> is skipped
+// the instantiation for class PreconditionBlockSOR<SparseMatrix<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" PreconditionBlockSOR
-template class PreconditionBlockSOR<double, float>;
+template class PreconditionBlockSOR<SparseMatrix<double>, float>;
 
 
-template void PreconditionBlockSOR<double, float>::vmult<float> (
+template void PreconditionBlockSOR<SparseMatrix<double>, float>::vmult<float> (
   Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockSOR<double, float>::vmult<double> (
+template void PreconditionBlockSOR<SparseMatrix<double>, float>::vmult<double> (
   Vector<double> &, const Vector<double> &) const;
-template void PreconditionBlockSOR<double, float>::Tvmult<float> (
+template void PreconditionBlockSOR<SparseMatrix<double>, float>::Tvmult<float> (
   Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockSOR<double, float>::Tvmult<double> (
+template void PreconditionBlockSOR<SparseMatrix<double>, float>::Tvmult<double> (
   Vector<double> &, const Vector<double> &) const;
 
-template class PreconditionBlockSOR<double, double>;
+template class PreconditionBlockSOR<SparseMatrix<double>, double>;
 
-template void PreconditionBlockSOR<double, double>::vmult<float> (
+template void PreconditionBlockSOR<SparseMatrix<double>, double>::vmult<float> (
   Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockSOR<double, double>::vmult<double> (
+template void PreconditionBlockSOR<SparseMatrix<double>, double>::vmult<double> (
   Vector<double> &, const Vector<double> &) const;
-template void PreconditionBlockSOR<double, double>::Tvmult<float> (
+template void PreconditionBlockSOR<SparseMatrix<double>, double>::Tvmult<float> (
   Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockSOR<double, double>::Tvmult<double> (
+template void PreconditionBlockSOR<SparseMatrix<double>, double>::Tvmult<double> (
   Vector<double> &, const Vector<double> &) const;
 
 
@@ -140,43 +140,43 @@ template void PreconditionBlockSOR<double, double>::Tvmult<double> (
 
 
 // explicit instantiations for "float" PreconditionBlock
-template class PreconditionBlockSSOR<float, float>;
+template class PreconditionBlockSSOR<SparseMatrix<float>, float>;
 
-template void PreconditionBlockSSOR<float, float>::vmult<float> (
+template void PreconditionBlockSSOR<SparseMatrix<float>, float>::vmult<float> (
   Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockSSOR<float, float>::vmult<double> (
+template void PreconditionBlockSSOR<SparseMatrix<float>, float>::vmult<double> (
   Vector<double> &, const Vector<double> &) const;
-template void PreconditionBlockSSOR<float, float>::Tvmult<float> (
+template void PreconditionBlockSSOR<SparseMatrix<float>, float>::Tvmult<float> (
   Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockSSOR<float, float>::Tvmult<double> (
+template void PreconditionBlockSSOR<SparseMatrix<float>, float>::Tvmult<double> (
   Vector<double> &, const Vector<double> &) const;
 
 
-// the instantiation for class PreconditionBlockSSOR<float, double> is skipped
+// the instantiation for class PreconditionBlockSSOR<SparseMatrix<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 class PreconditionBlockSSOR<SparseMatrix<double>, float>;
 
 
-template void PreconditionBlockSSOR<double, float>::vmult<float> (
+template void PreconditionBlockSSOR<SparseMatrix<double>, float>::vmult<float> (
   Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockSSOR<double, float>::vmult<double> (
+template void PreconditionBlockSSOR<SparseMatrix<double>, float>::vmult<double> (
   Vector<double> &, const Vector<double> &) const;
-template void PreconditionBlockSSOR<double, float>::Tvmult<float> (
+template void PreconditionBlockSSOR<SparseMatrix<double>, float>::Tvmult<float> (
   Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockSSOR<double, float>::Tvmult<double> (
+template void PreconditionBlockSSOR<SparseMatrix<double>, float>::Tvmult<double> (
   Vector<double> &, const Vector<double> &) const;
 
-template class PreconditionBlockSSOR<double, double>;
+template class PreconditionBlockSSOR<SparseMatrix<double>, double>;
 
-template void PreconditionBlockSSOR<double, double>::vmult<float> (
+template void PreconditionBlockSSOR<SparseMatrix<double>, double>::vmult<float> (
   Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockSSOR<double, double>::vmult<double> (
+template void PreconditionBlockSSOR<SparseMatrix<double>, double>::vmult<double> (
   Vector<double> &, const Vector<double> &) const;
-template void PreconditionBlockSSOR<double, double>::Tvmult<float> (
+template void PreconditionBlockSSOR<SparseMatrix<double>, double>::Tvmult<float> (
   Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockSSOR<double, double>::Tvmult<double> (
+template void PreconditionBlockSSOR<SparseMatrix<double>, double>::Tvmult<double> (
   Vector<double> &, const Vector<double> &) const;
index 241aadbd412aa0e2f4b57a2a9417c13880cec10e..ccfab17f9397066a9de234e5b54227075b2cfff2 100644 (file)
 #include <lac/vector.h>
 #include <lac/solver_richardson.h>
 #include <lac/precondition.h>
+#include <lac/precondition_block.h>
 
 #include <lac/sparse_matrix_ez.templates.h>
+#include <lac/precondition_block.templates.h>
 
 #include <fstream>
 #include <cstdio>
@@ -54,6 +56,15 @@ check_vmult_quadratic(std::vector<double>& residuals,
   sor.initialize(A, 1.2);
   PreconditionSSOR<MATRIX> ssor;
   ssor.initialize(A, 1.2);
+  PreconditionBlockJacobi<MATRIX, float> block_jacobi;
+  block_jacobi.initialize(A, (unsigned int) (sqrt(A.n()+.3)), 1.2);
+  block_jacobi.invert_diagblocks();
+  PreconditionBlockSSOR<MATRIX, float> block_ssor;
+  block_ssor.initialize(A, (unsigned int) (sqrt(A.n()+.3)), 1.2);
+  block_ssor.invert_diagblocks();
+  PreconditionBlockSOR<MATRIX, float> block_sor;
+  block_sor.initialize(A, (unsigned int) (sqrt(A.n()+.3)), 1.2);
+  block_sor.invert_diagblocks();
   
   u = 0.;
   f = 1.;
@@ -62,6 +73,10 @@ check_vmult_quadratic(std::vector<double>& residuals,
   PREC_CHECK(prich, solve, jacobi);
   PREC_CHECK(prich, solve, ssor);
   PREC_CHECK(prich, solve, sor);
+  u = 0.;
+  PREC_CHECK(prich, solve, block_jacobi);
+  PREC_CHECK(prich, solve, block_ssor);
+  PREC_CHECK(prich, solve, block_sor);
   
   u = 0.;
   deallog << "Transpose" << std::endl;
@@ -69,6 +84,10 @@ check_vmult_quadratic(std::vector<double>& residuals,
   PREC_CHECK(prich, Tsolve, jacobi);
   PREC_CHECK(prich, Tsolve, ssor);
   PREC_CHECK(prich, Tsolve, sor);
+  u = 0.;
+  PREC_CHECK(prich, Tsolve, block_jacobi);
+  PREC_CHECK(prich, Tsolve, block_ssor);
+  PREC_CHECK(prich, Tsolve, block_sor);
   deallog.pop();
 }
 
@@ -110,7 +129,7 @@ int main()
   deallog.depth_console(1000);
   deallog.log_execution_time(true);
   deallog.log_time_differences(true);
-  const unsigned int size = 500;
+  const unsigned int size = 50;
   const unsigned int row_length = 9;
 #endif
   
@@ -156,7 +175,7 @@ int main()
   testproblem.nine_point(E);
   check_vmult_quadratic(E_res, E, "9-SparseMatrixEZ<double>");
 
-  for (unsigned int i=0;i<A_res.size();++i)
+  for (unsigned int i=0;i<E_res.size();++i)
     if (std::fabs(A_res[i] - E_res[i]) > 1.e-14)
       deallog << "SparseMatrix and SparseMatrixEZ differ!!!"
              << std::endl;

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.