]> https://gitweb.dealii.org/ - dealii.git/commitdiff
introduce a base class for PreconditionBlock which only stores matrices
authorGuido Kanschat <dr.guido.kanschat@gmail.com>
Tue, 20 Jul 2010 05:54:46 +0000 (05:54 +0000)
committerGuido Kanschat <dr.guido.kanschat@gmail.com>
Tue, 20 Jul 2010 05:54:46 +0000 (05:54 +0000)
git-svn-id: https://svn.dealii.org/trunk@21542 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/include/lac/precondition_block_base.h [new file with mode: 0644]

index 4ac9f842871943a4c804c0f83ebbc3f09ae405a5..100be2e1d824473ccf5301a5cef0e43788977ac1 100644 (file)
 #include <base/exceptions.h>
 #include <base/subscriptor.h>
 #include <base/smartpointer.h>
+#include <lac/precondition_block_base.h>
 
 #include <vector>
 
 DEAL_II_NAMESPACE_OPEN
 
-template <typename number> class FullMatrix;
-template <typename number> class Vector;
-
 template<class MATRIX, typename inverse_type>
 class PreconditionBlockJacobi;
 
@@ -35,15 +33,14 @@ class PreconditionBlockJacobi;
 
 
 /**
- * Base class for @p PreconditionBlockJacobi, @p PreconditionBlockSOR,
- * ...  This class assumes the <tt>MATRIX</tt> consisting of
- * invertible blocks of @p blocksize on the diagonal and provides the
- * inversion of the diagonal blocks of the matrix. NOT only block
- * diagonal matrices are allowed but all matrices of arbitrary
- * structure with the minimal property of having invertible blocks on
- * the diagonal! Still the matrix must have access to single matrix
- * entries. therefore, BlockMatrixArray is not a possible matrix
- * class.
+ * Base class for actual block preconditioners. This class assumes the
+ * <tt>MATRIX</tt> consisting of invertible blocks of @p blocksize on
+ * the diagonal and provides the inversion of the diagonal blocks of
+ * the matrix. NOT only block diagonal matrices are allowed but all
+ * matrices of arbitrary structure with the minimal property of having
+ * invertible blocks on the diagonal! Still the matrix must have
+ * access to single matrix entries. therefore, BlockMatrixArray is not
+ * a possible matrix class.
  *
  * This block matrix structure is given e.g. for the DG method for the
  * transport equation. For a downstream numbering the matrices even
@@ -72,10 +69,13 @@ class PreconditionBlockJacobi;
  * choice.
  *
  * @see @ref GlossBlockLA "Block (linear algebra)"
- * @author Ralf Hartmann, Guido Kanschat, 1999, 2000
+ * @author Ralf Hartmann, Guido Kanschat
+ * @date 1999, 2000, 2010
  */
 template<class MATRIX, typename inverse_type = typename MATRIX::value_type>
-class PreconditionBlock : public virtual Subscriptor
+class PreconditionBlock
+  : public virtual Subscriptor,
+    protected PreconditionBlockBase<inverse_type>
 {
   private:
                                     /**
@@ -132,7 +132,7 @@ class PreconditionBlock : public virtual Subscriptor
                                     /**
                                      * Constructor.
                                      */
-    PreconditionBlock();
+    PreconditionBlock(bool store_diagonals = false);
 
                                     /**
                                      * Destructor.
@@ -225,27 +225,6 @@ class PreconditionBlock : public virtual Subscriptor
     value_type el(unsigned int i,
                  unsigned int j) const;
 
-                                    /**
-                                     * Use only the inverse of the
-                                     * first diagonal block to save
-                                     * memory and computation time.
-                                     *
-                                     * Possible applications:
-                                     * computing on a cartesian grid,
-                                     * all diagonal blocks are the
-                                     * same or all diagonal blocks
-                                     * are at least similar and
-                                     * inversion of one of them still
-                                     * yields a preconditioner.
-                                     */
-    void set_same_diagonal ();
-
-                                    /**
-                                     * Does the matrix use only one
-                                     * diagonal block?
-                                     */
-    bool same_diagonal () const;
-
                                     /**
                                      * Stores the inverse of the
                                      * diagonal blocks in
@@ -342,12 +321,6 @@ class PreconditionBlock : public virtual Subscriptor
                                      */
     unsigned int memory_consumption () const;
 
-                                    /**
-                                     * Determine, whether inverses
-                                     * have been computed.
-                                     */
-    bool inverses_ready () const;
-
                                     /** @addtogroup Exceptions
                                      * @{ */
 
@@ -360,39 +333,13 @@ class PreconditionBlock : public virtual Subscriptor
                    << " 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 (ExcDiagonalsNotStored);
-
                                     //@}
 
-                                    /**
-                                     * Access to the inverse diagonal
-                                     * blocks.
-                                     */
-    const FullMatrix<inverse_type>& inverse (unsigned int i) const;
-
-                                    /**
-                                     * Access to the diagonal
-                                     * blocks.
-                                     */
-    const FullMatrix<inverse_type>& diagonal (unsigned int i) const;
-
   protected:
                                     /**
                                      * Size of the blocks. Each
@@ -419,14 +366,7 @@ class PreconditionBlock : public virtual Subscriptor
                                      * used by derived classes.
                                      */
     double relaxation;
-
-
-                                    /**
-                                     * Flag for storing the diagonal
-                                     * blocks of the matrix.
-                                     */
-    bool store_diagonals;
-
+    
                                     /**
                                      * Number of blocks.
                                      */
@@ -446,30 +386,6 @@ class PreconditionBlock : public virtual Subscriptor
                                      * Flag for diagonal compression.
                                      * @ref set_same_diagonal()
                                      */
-  private:
-
-                                    /**
-                                     * Storage of the inverse
-                                     * matrices of the diagonal
-                                     * blocks matrices as
-                                     * <tt>FullMatrix<inverse_type></tt>
-                                     * matrices. Using
-                                     * <tt>inverse_type=float</tt> saves
-                                     * memory in comparison with
-                                     * <tt>inverse_type=double</tt>.
-                                     */
-    std::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.
-                                     */
-    std::vector<FullMatrix<inverse_type> > var_diagonal;
-
-    bool var_same_diagonal;
 };
 
 
@@ -808,6 +724,11 @@ class PreconditionBlockSOR : public virtual Subscriptor,
                             protected PreconditionBlock<MATRIX, inverse_type>
 {
   public:
+                                    /**
+                                     * Default constructor.
+                                     */
+    PreconditionBlockSOR();
+    
                                     /**
                                      * Define number type of matrix.
                                      */
@@ -944,6 +865,12 @@ class PreconditionBlockSOR : public virtual Subscriptor,
     void Tstep (Vector<number2>& dst, const Vector<number2>& rhs) const;
 
   protected:
+                                    /**
+                                     * Constructor to be used by
+                                     * PreconditionBlockSSOR.
+                                     */
+    PreconditionBlockSOR(bool store);
+    
                                     /**
                                      * Implementation of the forward
                                      * substitution loop called by
@@ -1138,7 +1065,7 @@ PreconditionBlock<MATRIX, inverse_type>::el (
   const unsigned int bs = blocksize;
   const unsigned int nb = i/bs;
 
-  const FullMatrix<inverse_type>& B = inverse(nb);
+  const FullMatrix<inverse_type>& B = this->inverse(nb);
 
   const unsigned int ib = i % bs;
   const unsigned int jb = j % bs;
index e86b7e04c51110b378b975beec7914a4524caf84..4ad3844a75fbf395c5eae7422c3f8a913a33bcd9 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009 by the deal.II authors
+//    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009, 2010 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -27,33 +27,23 @@ DEAL_II_NAMESPACE_OPEN
 
 
 template <class MATRIX, typename inverse_type>
-PreconditionBlock<MATRIX,inverse_type>::PreconditionBlock ():
-               blocksize(0),
-               A(0, typeid(*this).name()),
-               store_diagonals(false),
-               var_same_diagonal(false)
+PreconditionBlock<MATRIX,inverse_type>::PreconditionBlock (bool store)
+               : PreconditionBlockBase<inverse_type>(store),
+                 blocksize(0),
+                 A(0, typeid(*this).name())
 {}
 
 
 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());
-  if (var_diagonal.size()!=0)
-    var_diagonal.erase(var_diagonal.begin(), var_diagonal.end());
-}
+{}
 
 
 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());
-  if (var_diagonal.size()!=0)
-    var_diagonal.erase(var_diagonal.begin(), var_diagonal.end());
+  PreconditionBlockBase<inverse_type>::clear();
   blocksize     = 0;
-  var_same_diagonal = false;
   A = 0;
 }
 
@@ -72,13 +62,14 @@ void PreconditionBlock<MATRIX,inverse_type>::initialize (
   Assert (A->m()%bsize==0, ExcWrongBlockSize(bsize, A->m()));
   blocksize=bsize;
   relaxation = parameters.relaxation;
-  var_same_diagonal = parameters.same_diagonal;
   nblocks = A->m()/bsize;
+  this->reinit(nblocks, blocksize, parameters.same_diagonal);
   
   if (parameters.invert_diagonal)
     invert_diagblocks();
 }
 
+
 template <class MATRIX, typename inverse_type>
 void PreconditionBlock<MATRIX,inverse_type>::initialize (
   const MATRIX &M,
@@ -96,8 +87,8 @@ void PreconditionBlock<MATRIX,inverse_type>::initialize (
   Assert (A->m()%bsize==0, ExcWrongBlockSize(bsize, A->m()));
   blocksize=bsize;
   relaxation = parameters.relaxation;
-  var_same_diagonal = parameters.same_diagonal;
   nblocks = A->m()/bsize;
+  this->reinit(nblocks, blocksize, parameters.same_diagonal);
   
   if (parameters.invert_diagonal)
     invert_permuted_diagblocks(permutation, inverse_permutation);
@@ -113,21 +104,14 @@ void PreconditionBlock<MATRIX,inverse_type>::invert_permuted_diagblocks(
   Assert (blocksize!=0, ExcNotInitialized());
 
   const MATRIX &M=*A;
-  Assert (var_inverse.size()==0, ExcInverseMatricesAlreadyExist());
+  Assert (this->inverses_ready()==0, ExcInverseMatricesAlreadyExist());
 
   FullMatrix<inverse_type> M_cell(blocksize);
 
-  if (same_diagonal())
+  if (this->same_diagonal())
     {
       deallog << "PreconditionBlock uses only one diagonal block" << std::endl;
-                                      // Invert only the first block
-                                      // This is a copy of the code in the
-                                      // 'else' part, stripped of the outer loop
-      if (store_diagonals)
-       var_diagonal.resize(1);
-      var_inverse.resize(1);
-      var_inverse[0].reinit(blocksize, blocksize);
-
+      
       for (unsigned int row_cell=0; row_cell<blocksize; ++row_cell)
        {
          typename MATRIX::const_iterator entry = M.begin(row_cell);
@@ -139,9 +123,9 @@ void PreconditionBlock<MATRIX,inverse_type>::invert_permuted_diagblocks(
              ++entry;
            }
        }
-      if (store_diagonals)
-       var_diagonal[0] = M_cell;
-      var_inverse[0].invert(M_cell);
+      if (this->store_diagonals())
+       this->diagonal(0) = M_cell;
+      this->inverse(0).invert(M_cell);
     }
   else
     {
@@ -152,32 +136,6 @@ void PreconditionBlock<MATRIX,inverse_type>::invert_permuted_diagblocks(
                                       // blocks.
                                       // row, column are the global numbering
                                       // of the unkowns.
-
-                                      // set the var_inverse array to the right
-                                      // size. we could do it like this:
-                                      // var_inverse = vector<>(nblocks,FullMatrix<>())
-                                      // but this would involve copying many
-                                      // FullMatrix objects.
-                                      //
-                                      // the following is a neat trick which
-                                      // avoids copying
-      if (store_diagonals)
-       {
-         std::vector<FullMatrix<inverse_type> >
-           tmp(nblocks, FullMatrix<inverse_type>(blocksize));
-         var_diagonal.swap (tmp);
-       }
-
-      if (true)
-       {
-         std::vector<FullMatrix<inverse_type> >
-           tmp(nblocks, FullMatrix<inverse_type>(blocksize));
-         var_inverse.swap (tmp);
-                                          // make sure the tmp object
-                                          // goes out of scope as
-                                          // soon as possible
-       };
-
       M_cell = 0;
       
       for (unsigned int cell=0; cell<nblocks; ++cell)
@@ -204,12 +162,13 @@ void PreconditionBlock<MATRIX,inverse_type>::invert_permuted_diagblocks(
                  M_cell(row_cell, column_cell) = entry->value();
                }
            }
-
-         if (store_diagonals)
-           var_diagonal[cell] = M_cell;
-         var_inverse[cell].invert(M_cell);
+         
+         if (this->store_diagonals())
+           this->diagonal(cell) = M_cell;
+         this->inverse(cell).invert(M_cell);
        }
     }
+  this->inverses_computed(true);
 }
 
 
@@ -405,32 +364,6 @@ void PreconditionBlock<MATRIX,inverse_type>::backward_step (
 }
 
 
-template <class MATRIX, typename inverse_type>
-const FullMatrix<inverse_type>&
-PreconditionBlock<MATRIX,inverse_type>::inverse(unsigned int i) const
-{
-  if (same_diagonal())
-    return var_inverse[0];
-  
-  Assert (i < var_inverse.size(), ExcIndexRange(i,0,var_inverse.size()));
-  return var_inverse[i];
-}
-
-
-template <class MATRIX, typename inverse_type>
-const FullMatrix<inverse_type>&
-PreconditionBlock<MATRIX,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];
-}
-
-
 template <class MATRIX, typename inverse_type>
 unsigned int PreconditionBlock<MATRIX,inverse_type>::block_size() const
 {
@@ -438,31 +371,6 @@ unsigned int PreconditionBlock<MATRIX,inverse_type>::block_size() const
 }
 
 
-template <class MATRIX, typename inverse_type>
-void
-PreconditionBlock<MATRIX,inverse_type>::set_same_diagonal()
-{
-  Assert(var_inverse.size()==0, ExcInverseMatricesAlreadyExist());
-  var_same_diagonal = true;
-}
-
-
-template <class MATRIX, typename inverse_type>
-bool
-PreconditionBlock<MATRIX,inverse_type>::same_diagonal() const
-{
-  return var_same_diagonal;
-}
-
-
-template <class MATRIX, typename inverse_type>
-bool
-PreconditionBlock<MATRIX,inverse_type>::inverses_ready() const
-{
-  return (var_inverse.size() != 0);
-}
-
-
 template <class MATRIX, typename inverse_type>
 void PreconditionBlock<MATRIX,inverse_type>::invert_diagblocks()
 {
@@ -470,21 +378,13 @@ void PreconditionBlock<MATRIX,inverse_type>::invert_diagblocks()
   Assert (blocksize!=0, ExcNotInitialized());
 
   const MATRIX &M=*A;
-  Assert (var_inverse.size()==0, ExcInverseMatricesAlreadyExist());
-
+  Assert (this->inverses_ready()==0, ExcInverseMatricesAlreadyExist());
+  
   FullMatrix<inverse_type> M_cell(blocksize);
 
-  if (same_diagonal())
+  if (this->same_diagonal())
     {
       deallog << "PreconditionBlock uses only one diagonal block" << std::endl;
-                                      // Invert only the first block
-                                      // This is a copy of the code in the
-                                      // 'else' part, stripped of the outer loop
-      if (store_diagonals)
-       var_diagonal.resize(1);
-      var_inverse.resize(1);
-      var_inverse[0].reinit(blocksize, blocksize);
-
       for (unsigned int row_cell=0; row_cell<blocksize; ++row_cell)
        {
          typename MATRIX::const_iterator entry = M.begin(row_cell);
@@ -496,45 +396,12 @@ void PreconditionBlock<MATRIX,inverse_type>::invert_diagblocks()
              ++entry;
            }
        }
-      if (store_diagonals)
-       var_diagonal[0] = M_cell;
-      var_inverse[0].invert(M_cell);
+      if (this->store_diagonals())
+       this->diagonal(0) = M_cell;
+      this->inverse(0).invert(M_cell);
     }
   else
     {
-                                      // 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.
-
-                                      // set the @p var_inverse array to the right
-                                      // size. we could do it like this:
-                                      // var_inverse = vector<>(nblocks,FullMatrix<>())
-                                      // but this would involve copying many
-                                      // FullMatrix objects.
-                                      //
-                                      // the following is a neat trick which
-                                      // avoids copying
-      if (store_diagonals)
-       {
-         std::vector<FullMatrix<inverse_type> >
-           tmp(nblocks, FullMatrix<inverse_type>(blocksize));
-         var_diagonal.swap (tmp);
-       }
-
-      if (true)
-       {
-         std::vector<FullMatrix<inverse_type> >
-           tmp(nblocks, FullMatrix<inverse_type>(blocksize));
-         var_inverse.swap (tmp);
-                                          // make sure the tmp object
-                                          // goes out of scope as
-                                          // soon as possible
-       };
-
       M_cell = 0;
       
       for (unsigned int cell=0; cell<nblocks; ++cell)
@@ -558,11 +425,12 @@ void PreconditionBlock<MATRIX,inverse_type>::invert_diagblocks()
                }
            }
 
-         if (store_diagonals)
-           var_diagonal[cell] = M_cell;
-         var_inverse[cell].invert(M_cell);
+         if (this->store_diagonals())
+           this->diagonal(cell) = M_cell;
+         this->inverse(cell).invert(M_cell);
        }
     }
+  this->inverses_computed(true);
 }
 
 
@@ -587,11 +455,8 @@ template <class MATRIX, typename inverse_type>
 unsigned int
 PreconditionBlock<MATRIX,inverse_type>::memory_consumption () const
 {
-  unsigned int mem = sizeof(*this);
-  for (unsigned int i=0; i<var_inverse.size(); ++i)
-    mem += MemoryConsumption::memory_consumption(var_inverse[i]);
-  for (unsigned int i=0; i<var_diagonal.size(); ++i)
-    mem += MemoryConsumption::memory_consumption(var_diagonal[i]);
+  unsigned int mem = sizeof(*this) - sizeof(PreconditionBlockBase<inverse_type>);
+  mem += PreconditionBlockBase<inverse_type>::memory_consumption();
   return mem;
 }
 
@@ -768,6 +633,18 @@ void PreconditionBlockJacobi<MATRIX,inverse_type>
 /*--------------------- PreconditionBlockSOR -----------------------*/
 
 
+template <class MATRIX, typename inverse_type>
+PreconditionBlockSOR<MATRIX,inverse_type>::PreconditionBlockSOR ()
+               : PreconditionBlock<MATRIX,inverse_type> (false)
+               
+{}
+
+template <class MATRIX, typename inverse_type>
+PreconditionBlockSOR<MATRIX,inverse_type>::PreconditionBlockSOR (bool store)
+               : PreconditionBlock<MATRIX,inverse_type> (store)
+               
+{}
+
 template <class MATRIX, typename inverse_type>
 template <typename number2>
 void PreconditionBlockSOR<MATRIX,inverse_type>::forward (
@@ -1059,9 +936,9 @@ void PreconditionBlockSOR<MATRIX,inverse_type>
 
 template <class MATRIX, typename inverse_type>
 PreconditionBlockSSOR<MATRIX,inverse_type>::PreconditionBlockSSOR ()
-{
-  this->store_diagonals = 1;
-}
+               : PreconditionBlockSOR<MATRIX,inverse_type> (true)
+               
+{}
 
 
 template <class MATRIX, typename inverse_type>
diff --git a/deal.II/lac/include/lac/precondition_block_base.h b/deal.II/lac/include/lac/precondition_block_base.h
new file mode 100644 (file)
index 0000000..e064730
--- /dev/null
@@ -0,0 +1,412 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//---------------------------------------------------------------------------
+#ifndef __deal2__precondition_block_base_h
+#define __deal2__precondition_block_base_h
+
+
+#include <base/config.h>
+#include <base/exceptions.h>
+#include <base/subscriptor.h>
+#include <base/smartpointer.h>
+#include <base/memory_consumption.h>
+
+#include <vector>
+
+DEAL_II_NAMESPACE_OPEN
+
+template <typename number> class FullMatrix;
+template <typename number> class Vector;
+
+/**
+ * A class storing the inverse diagonal blocks for block
+ * preconditioners and block relaxation methods.
+ *
+ * This class does the book keeping for preconditioners and relaxation
+ * methods based on inverting blocks on the diagonal of a matrix.
+ * It allows us to either store all diagonal blocks and their
+ * inverses or the same block for each entry, and it keeps track of
+ * the choice. Thus, after initializing it and filling the inverse
+ * diagonal blocks correctly, a derived class can use inverse() with
+ * an integer argument referring to the block number.
+ *
+ * Additionally, it allows the storage of the original diagonal
+ * blocks, not only the inverses. These are for instance used in the
+ * intermediate step of the SSOR preconditioner.
+ *
+ * @author Guido Kanschat
+ * @date 2010
+ */
+template <typename number>
+class PreconditionBlockBase
+{
+  public:
+                                    /**
+                                     * Constructor initializing
+                                     * default values.
+                                     */
+    PreconditionBlockBase(bool store_diagonals);
+
+                                    /**
+                                     * Deletes the inverse diagonal
+                                     * block matrices if existent hence
+                                     * leaves the class in the state
+                                     * that it had directly after
+                                     * calling the constructor.
+                                     */
+    void clear();
+
+                                    /**
+                                     * Resize to this number of
+                                     * diagonal blocks with the given
+                                     * block size. If
+                                     * <tt>compress</tt> is true,
+                                     * then only one block will be
+                                     * stored.
+                                     */
+    void reinit(unsigned int nblocks, unsigned int blocksize, bool compress);
+    
+                                    /**
+                                     * Tell the class that inverses
+                                     * are computed.
+                                     */
+    void inverses_computed(bool are_they);
+    
+                                    /**
+                                     * Use only the inverse of the
+                                     * first diagonal block to save
+                                     * memory and computation time.
+                                     *
+                                     * Possible applications:
+                                     * computing on a cartesian grid,
+                                     * all diagonal blocks are the
+                                     * same or all diagonal blocks
+                                     * are at least similar and
+                                     * inversion of one of them still
+                                     * yields a preconditioner.
+                                     */
+    void set_same_diagonal ();
+
+                                    /**
+                                     * Does the matrix use only one
+                                     * diagonal block?
+                                     */
+    bool same_diagonal () const;
+
+                                    /**
+                                     * Check, whether diagonal blocks
+                                     * (not their inverses)
+                                     * should be stored.
+                                     */
+    bool store_diagonals() const;
+
+                                    /**
+                                     * Return true, if inverses are
+                                     * ready for use.
+                                     */
+    bool inverses_ready () const;
+    
+                                    /**
+                                     * Checks whether the object is empty.
+                                     */
+    bool empty () const;
+
+                                    /**
+                                     * The number of blocks.
+                                     */
+    unsigned int size() const;
+    
+                                    /**
+                                     * Read-only access to entries.
+                                     * This function is only possible
+                                     * if the inverse diagonal blocks
+                                     * are stored.
+                                     */
+    number el(unsigned int i, unsigned int j) const;
+    
+                                    /**
+                                     * Access to the inverse diagonal
+                                     * blocks.
+                                     */
+    FullMatrix<number>& inverse (unsigned int i);
+    
+                                    /**
+                                     * Access to the inverse diagonal
+                                     * blocks.
+                                     */
+    const FullMatrix<number>& inverse (unsigned int i) const;
+    
+                                    /**
+                                     * Access to the diagonal
+                                     * blocks.
+                                     */
+    FullMatrix<number>& diagonal (unsigned int i);
+    
+                                    /**
+                                     * Access to the diagonal
+                                     * blocks.
+                                     */
+    const FullMatrix<number>& diagonal (unsigned int i) const;
+
+                                    /**
+                                     * Determine an estimate for the
+                                     * memory consumption (in bytes)
+                                     * of this object.
+                                     */
+    unsigned int memory_consumption () const;
+    
+                                    /**
+                                     * You are trying to access a
+                                     * diagonal block (not its
+                                     * inverse), but you decided not
+                                     * to store the diagonal blocks.
+                                     */
+    DeclException0 (ExcDiagonalsNotStored);
+
+  private:
+                                    /**
+                                     * The number of (inverse)
+                                     * diagonal blocks, if only one
+                                     * is stored.
+                                     */
+    unsigned int n_diagonal_blocks;
+    
+                                    /**
+                                     * Storage of the inverse
+                                     * matrices of the diagonal
+                                     * blocks matrices as
+                                     * <tt>FullMatrix<number></tt>
+                                     * matrices. Using
+                                     * <tt>number=float</tt> saves
+                                     * memory in comparison with
+                                     * <tt>number=double</tt>.
+                                     */
+    std::vector<FullMatrix<number> > var_inverse;
+
+                                    /**
+                                     * Storage of the original diagonal blocks.
+                                     *
+                                     * Used by the blocked SSOR method.
+                                     */
+    std::vector<FullMatrix<number> > var_diagonal;
+
+
+                                    /**
+                                     * This is true, if the field
+                                     * #var_diagonal is to be used.
+                                     */
+    bool var_store_diagonals;
+
+                                    /**
+                                     * This is true, if only one
+                                     * inverse is stored.
+                                     */
+    bool var_same_diagonal;
+
+                                    /**
+                                     * The inverse matrices are
+                                     * usable. Set by the parent
+                                     * class via inverses_computed().
+                                     */
+    bool var_inverses_ready;
+};
+
+//----------------------------------------------------------------------//
+
+template <typename number>
+inline
+PreconditionBlockBase<number>::PreconditionBlockBase(bool store)
+               :
+               n_diagonal_blocks(0),
+               var_store_diagonals(store),
+               var_same_diagonal(false),
+               var_inverses_ready(false)
+{}
+
+
+template <typename number>
+inline
+void
+PreconditionBlockBase<number>::clear()
+{
+    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());
+  var_same_diagonal = false;
+  var_inverses_ready = false;
+  n_diagonal_blocks = 0;
+}
+
+template <typename number>
+inline
+void
+PreconditionBlockBase<number>::reinit(unsigned int n, unsigned int b, bool compress)
+{
+  var_same_diagonal = compress;
+  var_inverses_ready = false;
+  n_diagonal_blocks = n;
+  
+  if (compress)
+    {
+      var_inverse.resize(1);
+      var_inverse[0].reinit(b,b);
+      if (store_diagonals())
+       {
+         var_diagonal.resize(1);
+         var_diagonal[0].reinit(b,b);
+       }
+    }
+  else
+    {
+                                      // set the arrays to the right
+                                      // size. we could do it like this:
+                                      // var_inverse = vector<>(nblocks,FullMatrix<>())
+                                      // but this would involve copying many
+                                      // FullMatrix objects.
+                                      //
+                                      // the following is a neat trick which
+                                      // avoids copying
+      if (store_diagonals())
+       {
+         std::vector<FullMatrix<number> >
+           tmp(n, FullMatrix<number>(b));
+         var_diagonal.swap (tmp);
+       }
+      
+      if (true)
+       {
+         std::vector<FullMatrix<number> >
+           tmp(n, FullMatrix<number>(b));
+         var_inverse.swap (tmp);
+                                          // make sure the tmp object
+                                          // goes out of scope as
+                                          // soon as possible
+       }
+    }
+}
+
+
+template <typename number>
+inline
+unsigned int
+PreconditionBlockBase<number>::size() const
+{
+  return n_diagonal_blocks;
+}
+
+
+template <typename number>
+inline
+const FullMatrix<number>&
+PreconditionBlockBase<number>::inverse(unsigned int i) const
+{
+  if (same_diagonal())
+    return var_inverse[0];
+  
+  Assert (i < var_inverse.size(), ExcIndexRange(i,0,var_inverse.size()));
+  return var_inverse[i];
+}
+
+
+template <typename number>
+inline
+const FullMatrix<number>&
+PreconditionBlockBase<number>::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];
+}
+
+
+template <typename number>
+inline
+FullMatrix<number>&
+PreconditionBlockBase<number>::inverse(unsigned int i)
+{
+  if (same_diagonal())
+    return var_inverse[0];
+  
+  Assert (i < var_inverse.size(), ExcIndexRange(i,0,var_inverse.size()));
+  return var_inverse[i];
+}
+
+
+template <typename number>
+inline
+FullMatrix<number>&
+PreconditionBlockBase<number>::diagonal(unsigned int i)
+{
+  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];
+}
+
+
+template <typename number>
+inline bool
+PreconditionBlockBase<number>::same_diagonal() const
+{
+  return var_same_diagonal;
+}
+
+
+template <typename number>
+inline bool
+PreconditionBlockBase<number>::store_diagonals() const
+{
+  return var_store_diagonals;
+}
+
+
+template <typename number>
+inline void
+PreconditionBlockBase<number>::inverses_computed(bool x)
+{
+  var_inverses_ready = x;
+}
+
+
+template <typename number>
+inline bool
+PreconditionBlockBase<number>::inverses_ready() const
+{
+  return var_inverses_ready;
+}
+
+
+template <typename number>
+inline
+unsigned int
+PreconditionBlockBase<number>::memory_consumption () const
+{
+  unsigned int mem = sizeof(*this);
+  for (unsigned int i=0; i<var_inverse.size(); ++i)
+    mem += MemoryConsumption::memory_consumption(var_inverse[i]);
+  for (unsigned int i=0; i<var_diagonal.size(); ++i)
+    mem += MemoryConsumption::memory_consumption(var_diagonal[i]);
+  return mem;
+}
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif

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.