From: kanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Sun, 28 Nov 2010 23:21:29 +0000 (+0000)
Subject: cleanup
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=73c6e2399b2379e91c4eed666a8da616833a037a;p=dealii-svn.git

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 5d75ce1ae3..dfb125d5d2 100644
--- 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 4e28fa9d23..67b0ab6087 100644
--- 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;