cleanup

git-svn-id: https://svn.dealii.org/trunk@22871 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/include/deal.II/lac/precondition_block.h b/deal.II/include/deal.II/lac/precondition_block.h
index 5d75ce1ae3441d15424eae194e5dd1bb78017140..dfb125d5d2040097dd2d857420ecec3f586ec0aa 100644 (file)
--- a/deal.II/include/deal.II/lac/precondition_block.h
+++ b/deal.II/include/deal.II/lac/precondition_block.h
@@ -126,6 +126,22 @@ class PreconditionBlock
                                          * are equal to save memory.
                                          */
        bool same_diagonal;
+                                        /**
+                                         * Choose the inversion
+                                         * method for the blocks.
+                                         */
+       typename PreconditionBlockBase<inverse_type>::Inversion inversion;
+
+                                        /**
+                                         * The if #inversion is SVD,
+                                         * the threshold below which
+                                         * a singular value will be
+                                         * considered zero and thus
+                                         * not inverted. This
+                                         * parameter is used in the
+                                         * call to LAPACKFullMatrix::compute_inverse_svd().
+                                         */
+       double threshold;
     };
 
 
@@ -309,6 +325,9 @@ class PreconditionBlock
     unsigned int block_size () const;
 
                                     /**
+                                     * @deprecated Use size()
+                                     * instead.
+                                     *
                                      * The number of blocks of the
                                      * matrix.
                                      */
@@ -371,11 +390,6 @@ class PreconditionBlock
                                      */
     double relaxation;
     
-                                    /**
-                                     * Number of blocks.
-                                     */
-    unsigned int nblocks;
-
                                     /**
                                      * The permutation vector
                                      */
@@ -547,53 +561,21 @@ class PreconditionBlockJacobi : public virtual Subscriptor,
         Accessor accessor;
     };
 
-                                    /**
-                                     * Make initialization function
-                                     * publicly available.
-                                     */
     PreconditionBlock<MATRIX, inverse_type>::initialize;
-
-                                    /**
-                                     * Make function of base class public again.
-                                     */
     PreconditionBlock<MATRIX, inverse_type>::clear;
-
-                                    /**
-                                     * Make function of base class public again.
-                                     */
     PreconditionBlock<MATRIX, inverse_type>::empty;
-
-                                    /**
-                                     * Make function of base class public again.
-                                     */
     PreconditionBlock<MATRIX, inverse_type>::el;
-
-                                    /**
-                                     * Make function of base class public again.
-                                     */
     PreconditionBlock<MATRIX, inverse_type>::set_same_diagonal;
-
-                                    /**
-                                     * Make function public.
-                                     */
     PreconditionBlock<MATRIX, inverse_type>::invert_diagblocks;
-
-                                    /**
-                                     * Make function public.
-                                     */
     PreconditionBlock<MATRIX, inverse_type>::block_size;
-
+    PreconditionBlockBase<inverse_type>::size;
+    PreconditionBlockBase<inverse_type>::inverse;
+    PreconditionBlockBase<inverse_type>::inverse_householder;
+    PreconditionBlockBase<inverse_type>::inverse_svd;
                                     /**
-                                     * Make function public.
+                                     * @deprecated Use size() instead
                                      */
     PreconditionBlock<MATRIX, inverse_type>::n_blocks;
-                                    /**
-                                     * Make function accessible to iterator.
-                                     */
-    PreconditionBlock<MATRIX, inverse_type>::inverse;
-                                    /**
-                                     * Make function accessible.
-                                     */
     PreconditionBlock<MATRIX, inverse_type>::set_permutation;
 
                                     /**
@@ -749,42 +731,18 @@ class PreconditionBlockSOR : public virtual Subscriptor,
                                     /**
                                      * Make type publicly available.
                                      */
-    using PreconditionBlock<MATRIX,inverse_type>::AdditionalData;
-
-                                    /**
-                                     * Make initialization function
-                                     * publicly available.
-                                     */
-    using PreconditionBlock<MATRIX, inverse_type>::initialize;
-
-                                    /**
-                                     * Make function of base class public again.
-                                     */
-    using PreconditionBlock<MATRIX, inverse_type>::clear;
-
-                                    /**
-                                     * Make function of base class public again.
-                                     */
-    using PreconditionBlock<MATRIX, inverse_type>::empty;
-
-                                    /**
-                                     * Make function of base class public again.
-                                     */
-    using PreconditionBlock<MATRIX, inverse_type>::el;
-
-                                    /**
-                                     * Make function of base class public again.
-                                     */
-    using PreconditionBlock<MATRIX, inverse_type>::set_same_diagonal;
-
-                                    /**
-                                     * Make function of base class public again.
-                                     */
-    using PreconditionBlock<MATRIX, inverse_type>::invert_diagblocks;
-                                    /**
-                                     * Make function accessible.
-                                     */
-    using PreconditionBlock<MATRIX, inverse_type>::set_permutation;
+    PreconditionBlock<MATRIX,inverse_type>::AdditionalData;
+    PreconditionBlock<MATRIX, inverse_type>::initialize;
+    PreconditionBlock<MATRIX, inverse_type>::clear;
+    PreconditionBlock<MATRIX, inverse_type>::empty;
+    PreconditionBlockBase<inverse_type>::size;
+    PreconditionBlockBase<inverse_type>::inverse;
+    PreconditionBlockBase<inverse_type>::inverse_householder;
+    PreconditionBlockBase<inverse_type>::inverse_svd;
+    PreconditionBlock<MATRIX, inverse_type>::el;
+    PreconditionBlock<MATRIX, inverse_type>::set_same_diagonal;
+    PreconditionBlock<MATRIX, inverse_type>::invert_diagblocks;
+    PreconditionBlock<MATRIX, inverse_type>::set_permutation;
 
                                     /**
                                      * Execute block SOR
@@ -957,46 +915,28 @@ class PreconditionBlockSSOR : public virtual Subscriptor,
                                      */
     PreconditionBlockSSOR ();
 
-                                    /**
-                                     * Make type publicly available.
-                                     */
-    using PreconditionBlockSOR<MATRIX,inverse_type>::AdditionalData;
+                                    // Keep AdditionalData accessible
+    PreconditionBlockSOR<MATRIX,inverse_type>::AdditionalData;
 
+                                    // The following are the
+                                    // functions of the base classes
+                                    // which we want to keep
+                                    // accessible.
                                     /**
                                      * Make initialization function
                                      * publicly available.
                                      */
-    using PreconditionBlockSOR<MATRIX,inverse_type>::initialize;
-
-                                    /**
-                                     * Make function of base class public again.
-                                     */
-    using PreconditionBlockSOR<MATRIX,inverse_type>::clear;
-
-                                    /**
-                                     * Make function of base class public again.
-                                     */
-    using PreconditionBlockSOR<MATRIX,inverse_type>::set_permutation;
-
-                                    /**
-                                     * Make function of base class public again.
-                                     */
-    using PreconditionBlockSOR<MATRIX, inverse_type>::empty;
-
-                                    /**
-                                     * Make function of base class public again.
-                                     */
-    using PreconditionBlockSOR<MATRIX, inverse_type>::el;
-
-                                    /**
-                                     * Make function of base class public again.
-                                     */
-    using PreconditionBlockSOR<MATRIX,inverse_type>::set_same_diagonal;
-
-                                    /**
-                                     * Make function of base class public again.
-                                     */
-    using PreconditionBlockSOR<MATRIX,inverse_type>::invert_diagblocks;
+    PreconditionBlockSOR<MATRIX,inverse_type>::initialize;
+    PreconditionBlockSOR<MATRIX,inverse_type>::clear;
+    PreconditionBlockBase<inverse_type>::size;
+    PreconditionBlockBase<inverse_type>::inverse;
+    PreconditionBlockBase<inverse_type>::inverse_householder;
+    PreconditionBlockBase<inverse_type>::inverse_svd;
+    PreconditionBlockSOR<MATRIX,inverse_type>::set_permutation;
+    PreconditionBlockSOR<MATRIX, inverse_type>::empty;
+    PreconditionBlockSOR<MATRIX, inverse_type>::el;
+    PreconditionBlockSOR<MATRIX,inverse_type>::set_same_diagonal;
+    PreconditionBlockSOR<MATRIX,inverse_type>::invert_diagblocks;
 
                                     /**
                                      * Execute block SSOR
@@ -1039,21 +979,6 @@ class PreconditionBlockSSOR : public virtual Subscriptor,
 
 #ifndef DOXYGEN
 
-template<class MATRIX, typename inverse_type>
-inline
-PreconditionBlock<MATRIX, inverse_type>::AdditionalData::
-AdditionalData (const unsigned int block_size,
-               const double relaxation,
-               const bool invert_diagonal,
-               const bool same_diagonal)
-               :
-               relaxation (relaxation),
-               block_size(block_size),
-               invert_diagonal(invert_diagonal),
-               same_diagonal(same_diagonal)
-{}
-
-
 template<class MATRIX, typename inverse_type>
 inline bool
 PreconditionBlock<MATRIX, inverse_type>::empty () const
@@ -1068,7 +993,7 @@ template<class MATRIX, typename inverse_type>
 inline unsigned int
 PreconditionBlock<MATRIX, inverse_type>::n_blocks () const
 {
-  return nblocks;
+  return this->size();
 }
 
 

diff --git a/deal.II/include/deal.II/lac/precondition_block.templates.h b/deal.II/include/deal.II/lac/precondition_block.templates.h
index 4e28fa9d23106cdfda5abb66b9426d60792a5df2..67b0ab6087bdab4ab5531bb40a7a07731c84dc10 100644 (file)
--- a/deal.II/include/deal.II/lac/precondition_block.templates.h
+++ b/deal.II/include/deal.II/lac/precondition_block.templates.h
@@ -26,6 +26,22 @@
 DEAL_II_NAMESPACE_OPEN
 
 
+template<class MATRIX, typename inverse_type>
+PreconditionBlock<MATRIX, inverse_type>::AdditionalData::
+AdditionalData (const unsigned int block_size,
+               const double relaxation,
+               const bool invert_diagonal,
+               const bool same_diagonal)
+               :
+               relaxation (relaxation),
+               block_size(block_size),
+               invert_diagonal(invert_diagonal),
+               same_diagonal(same_diagonal),
+               inversion(PreconditionBlockBase<inverse_type>::gauss_jordan),
+               threshold(0.)
+{}
+
+
 template <typename number>
 PreconditionBlockBase<number>::~PreconditionBlockBase ()
 {}
@@ -67,8 +83,9 @@ void PreconditionBlock<MATRIX,inverse_type>::initialize (
   Assert (A->m()%bsize==0, ExcWrongBlockSize(bsize, A->m()));
   blocksize=bsize;
   relaxation = parameters.relaxation;
-  nblocks = A->m()/bsize;
-  this->reinit(nblocks, blocksize, parameters.same_diagonal);
+  const unsigned int nblocks = A->m()/bsize;
+  this->reinit(nblocks, blocksize, parameters.same_diagonal,
+              parameters.inversion);
   
   if (parameters.invert_diagonal)
     invert_diagblocks();
@@ -92,7 +109,7 @@ void PreconditionBlock<MATRIX,inverse_type>::initialize (
   Assert (A->m()%bsize==0, ExcWrongBlockSize(bsize, A->m()));
   blocksize=bsize;
   relaxation = parameters.relaxation;
-  nblocks = A->m()/bsize;
+  const unsigned int nblocks = A->m()/bsize;
   this->reinit(nblocks, blocksize, parameters.same_diagonal);
   
   if (parameters.invert_diagonal)
@@ -158,7 +175,7 @@ void PreconditionBlock<MATRIX,inverse_type>::invert_permuted_diagblocks(
                                       // of the unkowns.
       M_cell = 0;
       
-      for (unsigned int cell=0; cell<nblocks; ++cell)
+      for (unsigned int cell=0; cell<this->size(); ++cell)
        {
          const unsigned int cell_start = cell*blocksize;
          for (unsigned int row_cell=0; row_cell<blocksize; ++row_cell)
@@ -193,10 +210,10 @@ void PreconditionBlock<MATRIX,inverse_type>::invert_permuted_diagblocks(
              case PreconditionBlockBase<inverse_type>::householder:
                    this->inverse_householder(cell).initialize(M_cell);
                    break;
-         case PreconditionBlockBase<inverse_type>::svd:
-               this->inverse_svd(cell) = M_cell;
-               this->inverse_svd(cell).compute_inverse_svd(0.);
-               break;
+             case PreconditionBlockBase<inverse_type>::svd:
+                   this->inverse_svd(cell) = M_cell;
+                   this->inverse_svd(cell).compute_inverse_svd(0.);
+                   break;
              default:
                    Assert(false, ExcNotImplemented());
                    
@@ -241,7 +258,7 @@ void PreconditionBlock<MATRIX,inverse_type>::forward_step (
   FullMatrix<number> M_cell(this->blocksize);
 
                                   // Loop over all blocks
-  for (unsigned int cell=0; cell < this->nblocks; ++cell)
+  for (unsigned int cell=0; cell < this->size(); ++cell)
     {
       const unsigned int permuted_block_start = permuted
                                                ? permutation[block_start]
@@ -328,11 +345,11 @@ void PreconditionBlock<MATRIX,inverse_type>::backward_step (
                                       // row, column are the global numbering
                                       // of the unkowns.
   unsigned int row, row_cell;
-  unsigned int block_end=this->blocksize * this->nblocks;
+  unsigned int block_end=this->blocksize * this->size();
   number2 b_cell_row;
 
   FullMatrix<number> M_cell(this->blocksize);
-  for (unsigned int cell=this->nblocks; cell!=0 ;)
+  for (unsigned int cell=this->size(); cell!=0 ;)
     {
       --cell;
       const unsigned int block_start = block_end - this->blocksize;
@@ -433,13 +450,27 @@ void PreconditionBlock<MATRIX,inverse_type>::invert_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;
+         case PreconditionBlockBase<inverse_type>::svd:
+               this->inverse_svd(0) = M_cell;
+               this->inverse_svd(0).compute_inverse_svd(1.e-12);
+               break;
+         default:
+               Assert(false, ExcNotImplemented());
+       }
     }
   else
     {
       M_cell = 0;
       
-      for (unsigned int cell=0; cell<nblocks; ++cell)
+      for (unsigned int cell=0; cell<this->size(); ++cell)
        {
          const unsigned int cell_start = cell*blocksize;
          for (unsigned int row_cell=0; row_cell<blocksize; ++row_cell)
@@ -462,7 +493,21 @@ void PreconditionBlock<MATRIX,inverse_type>::invert_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;
+             case PreconditionBlockBase<inverse_type>::svd:
+                   this->inverse_svd(cell) = M_cell;
+                   this->inverse_svd(cell).compute_inverse_svd(1.e-12);
+                   break;
+             default:
+                   Assert(false, ExcNotImplemented());
+           }
        }
     }
   this->inverses_computed(true);
@@ -543,7 +588,7 @@ void PreconditionBlockJacobi<MATRIX,inverse_type>
   if (!this->inverses_ready())
     {
       FullMatrix<number> M_cell(this->blocksize);
-      for (unsigned int cell=0; cell < this->nblocks; ++cell)
+      for (unsigned int cell=0; cell < this->size(); ++cell)
        {
          for (row=cell*this->blocksize, row_cell=0;
               row_cell<this->blocksize;
@@ -569,7 +614,7 @@ void PreconditionBlockJacobi<MATRIX,inverse_type>
        }
     }
   else
-    for (unsigned int cell=0; cell < this->nblocks; ++cell)
+    for (unsigned int cell=0; cell < this->size(); ++cell)
       {
        for (row=cell*this->blocksize, row_cell=0;
             row_cell<this->blocksize;
@@ -730,7 +775,7 @@ void PreconditionBlockSOR<MATRIX,inverse_type>::forward (
                                   // inverses were not precomputed
   FullMatrix<number> M_cell(this->blocksize);
   
-  for (unsigned int cell=0; cell < this->nblocks; ++cell)
+  for (unsigned int cell=0; cell < this->size(); ++cell)
     {
       const unsigned int permuted_block_start = permuted
                                                ? this->permutation[block_start]
@@ -833,11 +878,11 @@ void PreconditionBlockSOR<MATRIX,inverse_type>::backward (
                                       // row, column are the global numbering
                                       // of the unkowns.
   unsigned int row, row_cell;
-  unsigned int block_end=this->blocksize * this->nblocks;
+  unsigned int block_end=this->blocksize * this->size();
   number2 b_cell_row;
 
   FullMatrix<number> M_cell(this->blocksize);
-  for (unsigned int cell=this->nblocks; cell!=0 ;)
+  for (unsigned int cell=this->size(); cell!=0 ;)
     {
       --cell;
       const unsigned int block_start = block_end - this->blocksize;
@@ -991,7 +1036,7 @@ void PreconditionBlockSSOR<MATRIX,inverse_type>::vmult (Vector<number2>       &d
   const double scaling = (2.-this->relaxation)/this->relaxation;
   
                                   // Multiply with diagonal blocks
-  for (unsigned int cell=0; cell < this->nblocks; ++cell)
+  for (unsigned int cell=0; cell < this->size(); ++cell)
     {
       unsigned int row = cell*this->blocksize;
       
@@ -1022,7 +1067,7 @@ void PreconditionBlockSSOR<MATRIX,inverse_type>::Tvmult (Vector<number2>       &
   const double scaling = (2.-this->relaxation)/this->relaxation;
   
                                   // Multiply with diagonal blocks
-  for (unsigned int cell=0; cell < this->nblocks; ++cell)
+  for (unsigned int cell=0; cell < this->size(); ++cell)
     {
       unsigned int row = cell*this->blocksize;