]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Provide a faster implementation of precondition_SSOR.
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 21 Apr 2009 08:02:28 +0000 (08:02 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 21 Apr 2009 08:02:28 +0000 (08:02 +0000)
git-svn-id: https://svn.dealii.org/trunk@18670 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/precondition.h
deal.II/lac/include/lac/sparse_decomposition.templates.h
deal.II/lac/include/lac/sparse_ilu.h
deal.II/lac/include/lac/sparse_ilu.templates.h
deal.II/lac/include/lac/sparse_matrix.h
deal.II/lac/include/lac/sparse_matrix.templates.h
deal.II/lac/source/sparse_matrix.inst.in

index 57227b9a593f862f3cdf4311134d8c800207e592..d0194995c4d8e19e1654887e237295bcb3a388be 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006 by the deal.II authors
+//    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2009 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -459,6 +459,26 @@ template <class MATRIX = SparseMatrix<double> >
 class PreconditionSSOR : public PreconditionRelaxation<MATRIX>
 {
   public:
+
+                                    /**
+                                     * A typedef to the base class.
+                                     */
+    typedef PreconditionRelaxation<MATRIX> BaseClass;
+
+       
+                                    /**
+                                     * Initialize matrix and
+                                     * relaxation parameter. The
+                                     * matrix is just stored in the
+                                     * preconditioner object. The
+                                     * relaxation parameter should be
+                                     * larger than zero and smaller
+                                     * than 2 for numerical
+                                     * reasons. It defaults to 1.
+                                     */
+    void initialize (const MATRIX &A,
+                    typename BaseClass::AdditionalData parameters = BaseClass::AdditionalData());
+
                                     /**
                                      * Apply preconditioner.
                                      */
@@ -474,6 +494,14 @@ class PreconditionSSOR : public PreconditionRelaxation<MATRIX>
                                      */
     template<class VECTOR>
     void Tvmult (VECTOR&, const VECTOR&) const;
+
+  private:
+                                    /**
+                                     * An array that stores for each matrix
+                                     * row where the first position after
+                                     * the diagonal is located.
+                                     */
+    std::vector<unsigned int> pos_right_of_diagonal;
 };
 
 
@@ -893,13 +921,46 @@ PreconditionSOR<MATRIX>::Tvmult (VECTOR &dst, const VECTOR &src) const
 
 //---------------------------------------------------------------------------
 
+template <class MATRIX>
+inline void
+PreconditionSSOR<MATRIX>::initialize (const MATRIX &rA,
+                                     typename BaseClass::AdditionalData parameters)
+{
+  this->PreconditionRelaxation<MATRIX>::initialize (rA, parameters);
+
+                                  // calculate the positions first after
+                                  // the diagonal.
+  const std::size_t  * rowstart_ptr = 
+    this->A->get_sparsity_pattern().get_rowstart_indices();
+  const unsigned int * const colnums = 
+    this->A->get_sparsity_pattern().get_column_numbers();
+  const unsigned int n = this->A->n();
+  pos_right_of_diagonal.resize(n);
+  for (unsigned int row=0; row<n; ++row, ++rowstart_ptr)
+    {
+                                      // find the first element in this line
+                                      // which is on the right of the diagonal.
+                                      // we need to precondition with the
+                                      // elements on the left only.
+                                      // note: the first entry in each
+                                      // line denotes the diagonal element,
+                                      // which we need not check.
+      pos_right_of_diagonal[row] = 
+       std::lower_bound (&colnums[*rowstart_ptr+1],
+                         &colnums[*(rowstart_ptr+1)],
+                         row)
+       - colnums;
+    }
+}
+
+
 template <class MATRIX>
 template<class VECTOR>
 inline void
 PreconditionSSOR<MATRIX>::vmult (VECTOR &dst, const VECTOR &src) const
 {
   Assert (this->A!=0, ExcNotInitialized());
-  this->A->precondition_SSOR (dst, src, this->relaxation);
+  this->A->precondition_SSOR (dst, src, this->relaxation, pos_right_of_diagonal);
 }
 
 
index 42b0b536fa585980f78143c4ee1f22b08e758b75..3d58ed0ae97e65a0f8f4ed23bb2d937f580175f8 100644 (file)
@@ -1,5 +1,5 @@
 //---------------------------------------------------------------------------
-//    Copyright (C) 2002, 2003, 2004, 2005, 2006 by the deal.II authors
+//    Copyright (C) 2002, 2003, 2004, 2005, 2006, 2009 by the deal.II authors
 //    by the deal.II authors and Stephen "Cheffo" Kolaroff
 //
 //    This file is subject to QPL and may not be  distributed
index b5938cbcec7009b829f24daaad65a2fc57178a33..8c4f2ceaf39d3dfc6d9f35e47501212bfcb6161c 100644 (file)
@@ -1,5 +1,5 @@
 //---------------------------------------------------------------------------
-//    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2008 by the deal.II authors
+//    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2008, 2009 by the deal.II authors
 //    by the deal.II authors and Stephen "Cheffo" Kolaroff
 //
 //    This file is subject to QPL and may not be  distributed
@@ -46,7 +46,7 @@ DEAL_II_NAMESPACE_OPEN
  * @<double@></tt>; others can be generated in application programs (see the
  * section on @ref Instantiations in the manual).
  * 
- * @author Wolfgang Bangerth, 2008; unified interface: Ralf Hartmann
+ * @author Wolfgang Bangerth, 2008, 2009; unified interface: Ralf Hartmann
  */
 template <typename number>
 class SparseILU : public SparseLUDecomposition<number>
index c19807c119b36ebbe13ba1125b81993bad516ee3..01b4a7164d8949631358eeed81ffdd6448ff166a 100644 (file)
@@ -1,5 +1,5 @@
 //---------------------------------------------------------------------------
-//    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2008 by the deal.II authors
+//    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2008, 2009 by the deal.II authors
 //    by the deal.II authors and Stephen "Cheffo" Kolaroff
 //
 //    This file is subject to QPL and may not be  distributed
index b6d348061365f9af652ab60b3e006ba9f7382034..f3948e2ee0042c7fe3d7976e2e0f8b5d4b856c5e 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name:  $
 //
-//    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors
+//    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -714,17 +714,20 @@ class SparseMatrix : public virtual Subscriptor
 
                                     /**
                                      * Return the number of actually
-                                     * nonzero elements of this
-                                     * matrix.
+                                     * nonzero elements of this matrix. It
+                                     * is possible to specify the parameter
+                                     * <tt>threshold</tt> in order to count
+                                     * only the elements that have absolute
+                                     * value greater than the threshold.
                                      *
-                                     * Note, that this function does
-                                     * (in contrary to
-                                     * n_nonzero_elements()) not
-                                     * count all entries of the
-                                     * sparsity pattern but only the
-                                     * ones that are nonzero.
-                                     */
-    unsigned int n_actually_nonzero_elements () const;
+                                     * Note, that this function does (in
+                                     * contrary to n_nonzero_elements())
+                                     * not count all entries of the
+                                     * sparsity pattern but only the ones
+                                     * that are nonzero (or whose absolute
+                                     * value is greater than threshold).
+                                     */
+    unsigned int n_actually_nonzero_elements (const double threshold = 0.) const;
     
                                     /**
                                      * Return a (constant) reference
@@ -1472,9 +1475,10 @@ class SparseMatrix : public virtual Subscriptor
                                      * <tt>src</tt>.
                                      */
     template <typename somenumber>
-    void precondition_SSOR (Vector<somenumber>       &dst,
-                           const Vector<somenumber> &src,
-                           const number              om = 1.) const;
+    void precondition_SSOR (Vector<somenumber>             &dst,
+                           const Vector<somenumber>       &src,
+                           const number                    om = 1.,
+                           const std::vector<unsigned int>&pos_right_of_diagonal=std::vector<unsigned int>()) const;
 
                                     /**
                                      * Apply SOR preconditioning
index d81775ef711a0e818cb1c656d61bd77e826b9c77..a9e77117399b54792e04192b994e99ca96424221 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors
+//    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -1155,9 +1155,10 @@ SparseMatrix<number>::precondition_Jacobi (Vector<somenumber>       &dst,
 template <typename number>
 template <typename somenumber>
 void
-SparseMatrix<number>::precondition_SSOR (Vector<somenumber>       &dst,
-                                        const Vector<somenumber> &src,
-                                        const number              om) const
+SparseMatrix<number>::precondition_SSOR (Vector<somenumber>              &dst,
+                                        const Vector<somenumber>        &src,
+                                        const number                     om,
+                                        const std::vector<unsigned int> &pos_right_of_diagonal) const
 {
                                   // to understand how this function works
                                   // you may want to take a look at the CVS
@@ -1175,6 +1176,55 @@ SparseMatrix<number>::precondition_SSOR (Vector<somenumber>       &dst,
   const std::size_t  *rowstart_ptr = &cols->rowstart[0];
   somenumber         *dst_ptr      = &dst(0);
 
+                                  // case when we have stored the position
+                                  // just right of the diagonal (then we
+                                  // don't have to search for it).
+  if (pos_right_of_diagonal.size() != 0)
+    {
+      Assert (pos_right_of_diagonal.size() == dst.size(),
+             ExcDimensionMismatch (pos_right_of_diagonal.size(), dst.size()));
+
+                                  // forward sweep
+      for (unsigned int row=0; row<n; ++row, ++dst_ptr, ++rowstart_ptr)
+       {
+         *dst_ptr = src(row);
+         const unsigned int first_right_of_diagonal_index = 
+           pos_right_of_diagonal[row];
+         Assert (first_right_of_diagonal_index <= *(rowstart_ptr+1),
+                 ExcInternalError());
+         double s = *dst_ptr;
+         for (unsigned int j=(*rowstart_ptr)+1; j<first_right_of_diagonal_index; ++j)
+           s -= val[j] * dst(cols->colnums[j]);
+
+                                  // divide by diagonal element
+         *dst_ptr = s * om / val[*rowstart_ptr];
+       };
+  
+      rowstart_ptr = &cols->rowstart[0];
+      dst_ptr      = &dst(0);
+      for (unsigned int row=0; row<n; ++row, ++rowstart_ptr, ++dst_ptr)
+       *dst_ptr *= (2.-om)*val[*rowstart_ptr];
+
+                                  // backward sweep
+      rowstart_ptr = &cols->rowstart[n-1];
+      dst_ptr      = &dst(n-1);
+      for (int row=n-1; row>=0; --row, --rowstart_ptr, --dst_ptr)
+       {
+         const unsigned int end_row = *(rowstart_ptr+1);
+         const unsigned int first_right_of_diagonal_index
+           = pos_right_of_diagonal[row];
+         double s = *dst_ptr;
+         for (unsigned int j=first_right_of_diagonal_index; j<end_row; ++j)
+           s -= val[j] * dst(cols->colnums[j]);
+      
+         *dst_ptr = s * om / val[*rowstart_ptr];
+       };
+      return;
+    }
+
+                                  // case when we need to get the position
+                                  // of the first element right of the
+                                  // diagonal manually for each sweep.
                                   // forward sweep
   for (unsigned int row=0; row<n; ++row, ++dst_ptr, ++rowstart_ptr)
     {
index 7766b8c87c6590c3b5853751bdfd8bd8ead49dfd..0edd36c17dd626355e84f01bef96b07369da5c5a 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id: sparse_matrix_matrix.in.h 15011 2007-08-22 16:59:41Z kanschat $
 //    Version: $Name$
 //
-//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2006, 2007 by the deal.II authors
+//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2006, 2007, 2009 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -64,7 +64,8 @@ for (S1, S2 : REAL_SCALARS)
     template void SparseMatrix<S1>::
       precondition_SSOR<S2> (Vector<S2> &,
                             const Vector<S2> &,
-                            const S1) const;
+                            const S1,
+                            const std::vector<unsigned int>&) const;
 
     template void SparseMatrix<S1>::
       precondition_SOR<S2> (Vector<S2> &,

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.