]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
allow for QR method in PreconditionBlockBase
authorkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 26 Oct 2010 19:56:35 +0000 (19:56 +0000)
committerkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 26 Oct 2010 19:56:35 +0000 (19:56 +0000)
git-svn-id: https://svn.dealii.org/trunk@22508 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/include/deal.II/lac/householder.h
deal.II/include/deal.II/lac/precondition_block.h
deal.II/include/deal.II/lac/precondition_block.templates.h
deal.II/include/deal.II/lac/precondition_block_base.h
deal.II/include/deal.II/lac/relaxation_block.h
deal.II/include/deal.II/lac/relaxation_block.templates.h

index c9a90496d609cf585c0614c46b86c2d516634266..e41be1b9dbc2e1e41c9444f0e20cad322efcea5d 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2005, 2006, 2007, 2008 by the deal.II authors
+//    Copyright (C) 2005, 2006, 2007, 2008, 2010 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -101,6 +101,18 @@ class Householder : private FullMatrix<number>
     double least_squares (BlockVector<number2> &dst,
                          const BlockVector<number2> &src) const;
 
+                                    /**
+                                     * A wrapper to least_squares(),
+                                     * implementing the standard
+                                     * MATRIX interface.
+                                     */
+    template<class VECTOR>
+    void vmult (VECTOR &dst, const VECTOR &src) const;
+
+    template<class VECTOR>
+    void Tvmult (VECTOR &dst, const VECTOR &src) const;
+
+    
   private:
                                     /**
                                      * Storage for the diagonal
@@ -280,6 +292,25 @@ Householder<number>::least_squares (BlockVector<number2>& dst,
 }
 
 
+template <typename number>
+template <class VECTOR>
+void
+Householder<number>::vmult (VECTOR &dst, const VECTOR &src) const
+{
+  least_squares (dst, src);
+}
+
+
+template <typename number>
+template <class VECTOR>
+void
+Householder<number>::Tvmult (VECTOR &, const VECTOR &) const
+{
+  Assert(false, ExcNotImplemented());
+}
+
+
+
 #endif // DOXYGEN
 
 DEAL_II_NAMESPACE_CLOSE
index 631e718bc117b73f1348f2b1b5a48abc9d145bfa..5d75ce1ae3441d15424eae194e5dd1bb78017140 100644 (file)
@@ -325,7 +325,11 @@ class PreconditionBlock
                                      * @{ */
 
                                     /**
-                                     * Exception
+                                     * For non-overlapping block
+                                     * preconditioners, the block
+                                     * size must divide the matrix
+                                     * size. If not, this exception
+                                     * gets thrown.
                                      */
     DeclException2 (ExcWrongBlockSize,
                    int, int,
index 0fb19a7c70801e20d1c5adfdf1ee2c4dbb984733..76eb4cacd3630cf19c34114b049e7809ba2a3549 100644 (file)
@@ -130,7 +130,18 @@ void PreconditionBlock<MATRIX,inverse_type>::invert_permuted_diagblocks(
        }
       if (this->store_diagonals())
        this->diagonal(0) = M_cell;
-      this->inverse(0).invert(M_cell);
+      switch (this->inversion)
+       {
+         case PreconditionBlockBase<inverse_type>::gauss_jordan:
+               this->inverse(0).invert(M_cell);
+               break;
+         case PreconditionBlockBase<inverse_type>::householder:
+               this->inverse_householder(0).initialize(M_cell);
+               break;
+         default:
+               Assert(false, ExcNotImplemented());
+
+       }
     }
   else
     {
@@ -170,7 +181,18 @@ void PreconditionBlock<MATRIX,inverse_type>::invert_permuted_diagblocks(
          
          if (this->store_diagonals())
            this->diagonal(cell) = M_cell;
-         this->inverse(cell).invert(M_cell);
+         switch (this->inversion)
+           {
+             case PreconditionBlockBase<inverse_type>::gauss_jordan:
+                   this->inverse(cell).invert(M_cell);
+                   break;
+             case PreconditionBlockBase<inverse_type>::householder:
+                   this->inverse_householder(cell).initialize(M_cell);
+                   break;
+             default:
+                   Assert(false, ExcNotImplemented());
+                   
+           }
        }
     }
   this->inverses_computed(true);
@@ -250,9 +272,9 @@ void PreconditionBlock<MATRIX,inverse_type>::forward_step (
       if (this->inverses_ready())
        {
          if (transpose_diagonal)
-           this->inverse(cell).Tvmult(x_cell, b_cell);
+           this->inverse_Tvmult(cell, x_cell, b_cell);
          else
-           this->inverse(cell).vmult(x_cell, b_cell);
+           this->inverse_vmult(cell, x_cell, b_cell);
        }
       else
        {
@@ -347,9 +369,9 @@ void PreconditionBlock<MATRIX,inverse_type>::backward_step (
       if (this->inverses_ready())
        {
          if (transpose_diagonal)
-           this->inverse(cell).Tvmult(x_cell, b_cell);
+           this->inverse_Tvmult(cell, x_cell, b_cell);
          else
-           this->inverse(cell).vmult(x_cell, b_cell);
+           this->inverse_vmult(cell, x_cell, b_cell);
        }
       else
        {
@@ -547,7 +569,7 @@ void PreconditionBlockJacobi<MATRIX,inverse_type>
          {
            b_cell(row_cell)=src(row);
          }
-       this->inverse(cell).vmult(x_cell, b_cell);
+       this->inverse_vmult(cell, x_cell, b_cell);
                                         // distribute x_cell to dst
        for (row=cell*this->blocksize, row_cell=0;
             row_cell<this->blocksize;
@@ -737,9 +759,9 @@ void PreconditionBlockSOR<MATRIX,inverse_type>::forward (
       if (this->inverses_ready())
        {
          if (transpose_diagonal)
-           this->inverse(cell).Tvmult(x_cell, b_cell);
+           this->inverse_Tvmult(cell, x_cell, b_cell);
          else
-           this->inverse(cell).vmult(x_cell, b_cell);
+           this->inverse_vmult(cell, x_cell, b_cell);
        }
       else
        {
@@ -851,9 +873,9 @@ void PreconditionBlockSOR<MATRIX,inverse_type>::backward (
       if (this->inverses_ready())
        {
          if (transpose_diagonal)
-           this->inverse(cell).Tvmult(x_cell, b_cell);
+           this->inverse_Tvmult(cell, x_cell, b_cell);
          else
-           this->inverse(cell).vmult(x_cell, b_cell);
+           this->inverse_vmult(cell, x_cell, b_cell);
        }
       else
        {
index b7962bd7fb1ec5a9da62cc5d9e03dcc96086a68c..d8bc501f1bce0939009965a4ecac20125bc87d3b 100644 (file)
@@ -19,6 +19,7 @@
 #include <base/subscriptor.h>
 #include <base/smartpointer.h>
 #include <base/memory_consumption.h>
+#include <lac/householder.h>
 
 #include <vector>
 
@@ -50,11 +51,32 @@ template <typename number>
 class PreconditionBlockBase
 {
   public:
+                                    /**
+                                     * Choose a method for inverting
+                                     * the blocks, and thus the data
+                                     * type for the inverses.
+                                     */
+    enum Inversion
+    {
+                                          /**
+                                           * Use the standard
+                                           * Gauss-Jacobi method
+                                           * implemented in FullMatrix::inverse().
+                                           */
+         gauss_jordan,
+                                          /**
+                                           * Use QR decomposition of
+                                           * the Householder class.
+                                           */
+         householder
+    };
+    
                                     /**
                                      * Constructor initializing
                                      * default values.
                                      */
-    PreconditionBlockBase(bool store_diagonals = false);
+    PreconditionBlockBase(bool store_diagonals = false,
+                         Inversion method = gauss_jordan);
 
                                     /**
                                      * The virtual destructor
@@ -78,7 +100,8 @@ class PreconditionBlockBase
                                      * then only one block will be
                                      * stored.
                                      */
-    void reinit(unsigned int nblocks, unsigned int blocksize, bool compress);
+    void reinit(unsigned int nblocks, unsigned int blocksize, bool compress,
+               Inversion method = gauss_jordan);
     
                                     /**
                                      * Tell the class that inverses
@@ -137,13 +160,33 @@ class PreconditionBlockBase
                                      * are stored.
                                      */
     number el(unsigned int i, unsigned int j) const;
+
+                                    /**
+                                     * Multiply with the inverse
+                                     * block at position <tt>i</tt>.
+                                     */
+    template <typename number2>
+    void inverse_vmult(unsigned int i, Vector<number2>& dst, const Vector<number2>& src) const;
+    
+                                    /**
+                                     * Multiply with the transposed inverse
+                                     * block at position <tt>i</tt>.
+                                     */
+    template <typename number2>
+    void inverse_Tvmult(unsigned int i, Vector<number2>& dst, const Vector<number2>& src) const;
     
                                     /**
                                      * Access to the inverse diagonal
-                                     * blocks.
+                                     * blocks if Inversion is #gauss_jordan.
                                      */
     FullMatrix<number>& inverse (unsigned int i);
     
+                                    /**
+                                     * Access to the inverse diagonal
+                                     * blocks if Inversion is #householder.
+                                     */
+    Householder<number>& inverse_householder (unsigned int i);
+    
                                     /**
                                      * Access to the inverse diagonal
                                      * blocks.
@@ -177,6 +220,12 @@ class PreconditionBlockBase
                                      */
     DeclException0 (ExcDiagonalsNotStored);
 
+  protected:
+                                    /**
+                                     * The method used for inverting blocks.
+                                     */
+    Inversion inversion;
+    
   private:
                                     /**
                                      * The number of (inverse)
@@ -185,17 +234,31 @@ class PreconditionBlockBase
                                      */
     unsigned int n_diagonal_blocks;
     
-                                    /**
+                                    /**
                                      * Storage of the inverse
                                      * matrices of the diagonal
                                      * blocks matrices as
                                      * <tt>FullMatrix<number></tt>
-                                     * matrices. Using
+                                     * matrices, if Inversion
+                                     * #gauss_jordan is used. Using
                                      * <tt>number=float</tt> saves
                                      * memory in comparison with
                                      * <tt>number=double</tt>.
                                      */
-    std::vector<FullMatrix<number> > var_inverse;
+    std::vector<FullMatrix<number> > var_inverse_full;
+
+                                    /**
+                                     * Storage of the inverse
+                                     * matrices of the diagonal
+                                     * blocks matrices as
+                                     * <tt>Householder<number></tt>
+                                     * matrices if Inversion
+                                     * #householder is used. Using
+                                     * <tt>number=float</tt> saves
+                                     * memory in comparison with
+                                     * <tt>number=double</tt>.
+                                     */
+    std::vector<Householder<number> > var_inverse_householder;
 
                                     /**
                                      * Storage of the original diagonal blocks.
@@ -229,8 +292,10 @@ class PreconditionBlockBase
 
 template <typename number>
 inline
-PreconditionBlockBase<number>::PreconditionBlockBase(bool store)
+PreconditionBlockBase<number>::PreconditionBlockBase(
+  bool store, Inversion method)
                :
+               inversion(method),
                n_diagonal_blocks(0),
                var_store_diagonals(store),
                var_same_diagonal(false),
@@ -243,8 +308,10 @@ inline
 void
 PreconditionBlockBase<number>::clear()
 {
-    if (var_inverse.size()!=0)
-    var_inverse.erase(var_inverse.begin(), var_inverse.end());
+  if (var_inverse_full.size()!=0)
+    var_inverse_full.erase(var_inverse_full.begin(), var_inverse_full.end());
+  if (var_inverse_full.size()!=0)
+    var_inverse_householder.erase(var_inverse_householder.begin(), var_inverse_householder.end());
   if (var_diagonal.size()!=0)
     var_diagonal.erase(var_diagonal.begin(), var_diagonal.end());
   var_same_diagonal = false;
@@ -255,16 +322,29 @@ PreconditionBlockBase<number>::clear()
 template <typename number>
 inline
 void
-PreconditionBlockBase<number>::reinit(unsigned int n, unsigned int b, bool compress)
+PreconditionBlockBase<number>::reinit(unsigned int n, unsigned int b, bool compress,
+Inversion method)
 {
+  inversion = method;
   var_same_diagonal = compress;
   var_inverses_ready = false;
   n_diagonal_blocks = n;
   
   if (compress)
     {
-      var_inverse.resize(1);
-      var_inverse[0].reinit(b,b);
+      switch(inversion)
+       {
+         case gauss_jordan:
+               var_inverse_full.resize(1);
+               var_inverse_full[0].reinit(b,b);
+               break;
+         case householder:
+               var_inverse_householder.resize(1);
+               break;
+         default:
+               Assert(false, ExcNotImplemented());
+       }
+      
       if (store_diagonals())
        {
          var_diagonal.resize(1);
@@ -287,15 +367,21 @@ PreconditionBlockBase<number>::reinit(unsigned int n, unsigned int b, bool compr
            tmp(n, FullMatrix<number>(b));
          var_diagonal.swap (tmp);
        }
-      
-      if (true)
+
+      switch(inversion)
        {
-         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
+         case gauss_jordan:
+         {
+           std::vector<FullMatrix<number> >
+             tmp(n, FullMatrix<number>(b));
+           var_inverse_full.swap (tmp);
+           break;
+         }
+         case householder:
+               var_inverse_householder.resize(n);
+               break;
+         default:      
+               Assert(false, ExcNotImplemented());     
        }
     }
 }
@@ -309,6 +395,55 @@ PreconditionBlockBase<number>::size() const
   return n_diagonal_blocks;
 }
 
+template <typename number>
+template <typename number2>
+inline
+void
+PreconditionBlockBase<number>::inverse_vmult(
+  unsigned int i, Vector<number2>& dst, const Vector<number2>& src) const
+{
+  const unsigned int ii = same_diagonal() ? 0U : i;
+  
+  switch (inversion)
+    {
+      case gauss_jordan:
+           AssertIndexRange (ii, var_inverse_full.size());
+           var_inverse_full[ii].vmult(dst, src);
+           break;
+      case householder:
+           AssertIndexRange (ii, var_inverse_householder.size());
+           var_inverse_householder[ii].vmult(dst, src);
+           break;
+      default:
+           Assert(false, ExcNotImplemented());
+    }
+}
+
+
+template <typename number>
+template <typename number2>
+inline
+void
+PreconditionBlockBase<number>::inverse_Tvmult(
+  unsigned int i, Vector<number2>& dst, const Vector<number2>& src) const
+{
+  const unsigned int ii = same_diagonal() ? 0U : i;
+  
+  switch (inversion)
+    {
+      case gauss_jordan:
+           AssertIndexRange (ii, var_inverse_full.size());
+           var_inverse_full[ii].Tvmult(dst, src);
+           break;
+      case householder:
+           AssertIndexRange (ii, var_inverse_householder.size());
+           var_inverse_householder[ii].Tvmult(dst, src);
+           break;
+      default:
+           Assert(false, ExcNotImplemented());
+    }
+}
+
 
 template <typename number>
 inline
@@ -316,10 +451,10 @@ const FullMatrix<number>&
 PreconditionBlockBase<number>::inverse(unsigned int i) const
 {
   if (same_diagonal())
-    return var_inverse[0];
+    return var_inverse_full[0];
   
-  Assert (i < var_inverse.size(), ExcIndexRange(i,0,var_inverse.size()));
-  return var_inverse[i];
+  Assert (i < var_inverse_full.size(), ExcIndexRange(i,0,var_inverse_full.size()));
+  return var_inverse_full[i];
 }
 
 
@@ -344,10 +479,23 @@ FullMatrix<number>&
 PreconditionBlockBase<number>::inverse(unsigned int i)
 {
   if (same_diagonal())
-    return var_inverse[0];
+    return var_inverse_full[0];
+  
+  Assert (i < var_inverse_full.size(), ExcIndexRange(i,0,var_inverse_full.size()));
+  return var_inverse_full[i];
+}
+
+
+template <typename number>
+inline
+Householder<number>&
+PreconditionBlockBase<number>::inverse_householder(unsigned int i)
+{
+  if (same_diagonal())
+    return var_inverse_householder[0];
   
-  Assert (i < var_inverse.size(), ExcIndexRange(i,0,var_inverse.size()));
-  return var_inverse[i];
+  AssertIndexRange (i, var_inverse_householder.size());
+  return var_inverse_householder[i];
 }
 
 
@@ -404,8 +552,8 @@ 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_inverse_full.size(); ++i)
+    mem += MemoryConsumption::memory_consumption(var_inverse_full[i]);
   for (unsigned int i=0; i<var_diagonal.size(); ++i)
     mem += MemoryConsumption::memory_consumption(var_diagonal[i]);
   return mem;
index 61d26ba291e9c568f6a5103f2c55d22744e82432..8052b27c13bc1b1cfbf41c51a236ef35edcb10c4 100644 (file)
@@ -101,6 +101,11 @@ class RelaxationBlock :
                                          * are equal to save memory.
                                          */
        bool same_diagonal;
+                                        /**
+                                         * Choose the inversion
+                                         * method for the blocks.
+                                         */
+       typename PreconditionBlockBase<inverse_type>::Inversion inversion;
     };
 
                                     /**
index 4206ff170bef350046e9b107d7d62956a9bf7844..003ea60b6d9062439d9c3b03813c1e16b71667f6 100644 (file)
@@ -24,7 +24,8 @@ RelaxationBlock<MATRIX,inverse_type>::AdditionalData::AdditionalData ()
                :
                relaxation(1.),
                invert_diagonal(true),
-               same_diagonal(false)
+               same_diagonal(false),
+               inversion(PreconditionBlockBase<inverse_type>::gauss_jordan)
 {}
 
 
@@ -39,7 +40,8 @@ RelaxationBlock<MATRIX,inverse_type>::AdditionalData::AdditionalData (
                block_list(&bl),
                relaxation(relaxation),
                invert_diagonal(invert_diagonal),
-               same_diagonal(same_diagonal)
+               same_diagonal(same_diagonal),
+               inversion(PreconditionBlockBase<inverse_type>::gauss_jordan)
 {}
 
 
@@ -57,7 +59,8 @@ RelaxationBlock<MATRIX,inverse_type>::initialize (
   A = &M;
   additional_data = parameters;
   
-  this->reinit(additional_data.block_list->size(), 0, additional_data.same_diagonal);
+  this->reinit(additional_data.block_list->size(), 0, additional_data.same_diagonal,
+              additional_data.inversion);
   
   if (parameters.invert_diagonal)
     invert_diagblocks();
@@ -120,8 +123,18 @@ RelaxationBlock<MATRIX,inverse_type>::invert_diagblocks ()
              this->diagonal(block).reinit(bs, bs);
              this->diagonal(block) = M_cell;
            }
-         this->inverse(block).reinit(bs, bs);
-         this->inverse(block).invert(M_cell);
+         switch(this->inversion)
+           {
+             case PreconditionBlockBase<inverse_type>::gauss_jordan:
+                   this->inverse(block).reinit(bs, bs);
+                   this->inverse(block).invert(M_cell);
+                   break;
+             case PreconditionBlockBase<inverse_type>::householder:
+                   this->inverse_householder(block).initialize(M_cell);
+                   break;
+             default:
+                   Assert(false, ExcNotImplemented());
+           }
        }
     }
   this->inverses_computed(true);
@@ -162,7 +175,7 @@ RelaxationBlock<MATRIX,inverse_type>::do_step (
            b_cell(row_cell) -= entry->value() * prev(entry->column());
        }
                                       // Apply inverse diagonal
-      this->inverse(block).vmult(x_cell, b_cell);
+      this->inverse_vmult(block, x_cell, b_cell);
                                       // Store in result vector
       row=additional_data.block_list->begin(block);
       for (unsigned int row_cell=0; row_cell<bs; ++row_cell, ++row)

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.