From 348e79a5b4af19f485cb55b65e030a2ba98ea73f Mon Sep 17 00:00:00 2001
From: Guido Kanschat <dr.guido.kanschat@gmail.com>
Date: Tue, 9 Nov 2010 15:01:03 +0000
Subject: [PATCH] add option of inverting blocks with SVD

git-svn-id: https://svn.dealii.org/trunk@22643 0785d39b-7218-0410-832d-ea1e28bc413d
---
 deal.II/include/deal.II/lac/matrix_block.h    |  1 +
 .../deal.II/lac/precondition_block_base.h     | 38 ++++++++++++++++
 .../include/deal.II/lac/relaxation_block.h    | 44 +++++++++++++++++++
 .../deal.II/lac/relaxation_block.templates.h  |  8 ++--
 4 files changed, 88 insertions(+), 3 deletions(-)

diff --git a/deal.II/include/deal.II/lac/matrix_block.h b/deal.II/include/deal.II/lac/matrix_block.h
index 1dba79fd2a..ca8d59c04a 100644
--- a/deal.II/include/deal.II/lac/matrix_block.h
+++ b/deal.II/include/deal.II/lac/matrix_block.h
@@ -465,6 +465,7 @@ class MatrixBlockVector
     NamedData<std_cxx1x::shared_ptr<value_type> >::subscribe;
     NamedData<std_cxx1x::shared_ptr<value_type> >::unsubscribe;
     NamedData<std_cxx1x::shared_ptr<value_type> >::size;
+    NamedData<std_cxx1x::shared_ptr<value_type> >::name;
 };
 
 
diff --git a/deal.II/include/deal.II/lac/precondition_block_base.h b/deal.II/include/deal.II/lac/precondition_block_base.h
index 0b15a13d43..b3c3a6a004 100644
--- a/deal.II/include/deal.II/lac/precondition_block_base.h
+++ b/deal.II/include/deal.II/lac/precondition_block_base.h
@@ -205,6 +205,18 @@ class PreconditionBlockBase
 				      */
     const FullMatrix<number>& inverse (unsigned int i) const;
     
+				     /**
+				      * Access to the inverse diagonal
+				      * blocks if Inversion is #householder.
+				      */
+    const Householder<number>& inverse_householder (unsigned int i) const;
+    
+				     /**
+				      * Access to the inverse diagonal
+				      * blocks if Inversion is #householder.
+				      */
+    const LAPACKFullMatrix<number>& inverse_svd (unsigned int i) const;
+    
 				     /**
 				      * Access to the diagonal
 				      * blocks.
@@ -507,6 +519,32 @@ PreconditionBlockBase<number>::inverse(unsigned int i) const
 }
 
 
+template <typename number>
+inline
+const Householder<number>&
+PreconditionBlockBase<number>::inverse_householder(unsigned int i) const
+{
+  if (same_diagonal())
+    return var_inverse_householder[0];
+  
+  AssertIndexRange (i, var_inverse_householder.size());
+  return var_inverse_householder[i];
+}
+
+
+template <typename number>
+inline
+const LAPACKFullMatrix<number>&
+PreconditionBlockBase<number>::inverse_svd(unsigned int i) const
+{
+  if (same_diagonal())
+    return var_inverse_svd[0];
+  
+  AssertIndexRange (i, var_inverse_svd.size());
+  return var_inverse_svd[i];
+}
+
+
 template <typename number>
 inline
 const FullMatrix<number>&
diff --git a/deal.II/include/deal.II/lac/relaxation_block.h b/deal.II/include/deal.II/lac/relaxation_block.h
index 8052b27c13..d561195593 100644
--- a/deal.II/include/deal.II/lac/relaxation_block.h
+++ b/deal.II/include/deal.II/lac/relaxation_block.h
@@ -106,6 +106,17 @@ class RelaxationBlock :
 					  * 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;
     };
 
 				     /**
@@ -267,6 +278,22 @@ class RelaxationBlockSOR : public virtual Subscriptor,
 				      * Make function of base class public again.
 				      */
     using RelaxationBlock<MATRIX, inverse_type>::empty;
+				     /**
+				      * Make function of base class public again.
+				      */
+    using RelaxationBlock<MATRIX, inverse_type>::size;
+				     /**
+				      * Make function of base class public again.
+				      */
+    using RelaxationBlock<MATRIX, inverse_type>::inverse;
+				     /**
+				      * Make function of base class public again.
+				      */
+    using RelaxationBlock<MATRIX, inverse_type>::inverse_householder;
+				     /**
+				      * Make function of base class public again.
+				      */
+    using RelaxationBlock<MATRIX, inverse_type>::inverse_svd;
 				     /**
 				      * Perform one step of the SOR
 				      * iteration.
@@ -335,6 +362,23 @@ class RelaxationBlockSSOR : public virtual Subscriptor,
 				      * Make function of base class public again.
 				      */
     using RelaxationBlock<MATRIX, inverse_type>::empty;
+
+				     /**
+				      * Make function of base class public again.
+				      */
+    using RelaxationBlock<MATRIX, inverse_type>::size;
+				     /**
+				      * Make function of base class public again.
+				      */
+    using RelaxationBlock<MATRIX, inverse_type>::inverse;
+				     /**
+				      * Make function of base class public again.
+				      */
+    using RelaxationBlock<MATRIX, inverse_type>::inverse_householder;
+				     /**
+				      * Make function of base class public again.
+				      */
+    using RelaxationBlock<MATRIX, inverse_type>::inverse_svd;
 				     /**
 				      * Perform one step of the SOR
 				      * iteration.
diff --git a/deal.II/include/deal.II/lac/relaxation_block.templates.h b/deal.II/include/deal.II/lac/relaxation_block.templates.h
index 6298b4bfb3..dbc3a56980 100644
--- a/deal.II/include/deal.II/lac/relaxation_block.templates.h
+++ b/deal.II/include/deal.II/lac/relaxation_block.templates.h
@@ -25,7 +25,8 @@ RelaxationBlock<MATRIX,inverse_type>::AdditionalData::AdditionalData ()
 		relaxation(1.),
 		invert_diagonal(true),
 		same_diagonal(false),
-		inversion(PreconditionBlockBase<inverse_type>::gauss_jordan)
+		inversion(PreconditionBlockBase<inverse_type>::gauss_jordan),
+		threshold(0.)
 {}
 
 
@@ -41,7 +42,8 @@ RelaxationBlock<MATRIX,inverse_type>::AdditionalData::AdditionalData (
 		relaxation(relaxation),
 		invert_diagonal(invert_diagonal),
 		same_diagonal(same_diagonal),
-		inversion(PreconditionBlockBase<inverse_type>::gauss_jordan)
+		inversion(PreconditionBlockBase<inverse_type>::gauss_jordan),
+		threshold(0.)
 {}
 
 
@@ -135,7 +137,7 @@ RelaxationBlock<MATRIX,inverse_type>::invert_diagblocks ()
 	      case PreconditionBlockBase<inverse_type>::svd:
 		    this->inverse_svd(block).reinit(bs, bs);
 		    this->inverse_svd(block) = M_cell;
-		    this->inverse_svd(block).compute_inverse_svd(0.);
+		    this->inverse_svd(block).compute_inverse_svd(additional_data.threshold);
 		    break;
 	      default:
 		    Assert(false, ExcNotImplemented());
-- 
2.39.5