add option of inverting blocks with SVD

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

diff --git a/deal.II/include/deal.II/lac/matrix_block.h b/deal.II/include/deal.II/lac/matrix_block.h
index 1dba79fd2a72dc3a884f7a678432103ec90047e5..ca8d59c04a8092654be22c23410d8219d01d3d84 100644 (file)
--- 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 0b15a13d436af284411ffc7a6f576256bf2496b4..b3c3a6a004d78b515a1ccaf23bd49abeef8deeaa 100644 (file)
--- 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 8052b27c13bc1b1cfbf41c51a236ef35edcb10c4..d56119559329be3532986e943b7d9f0043693100 100644 (file)
--- 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 6298b4bfb3687c6720b47feee36b4ce2a7d972a8..dbc3a56980ba198e4d59d9dc34dec1fc2c4578d0 100644 (file)
--- 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());