My previously committed version did not work for SparseMatrixEZ. This should now...

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

diff --git a/deal.II/lac/include/lac/precondition.h b/deal.II/lac/include/lac/precondition.h
index d0194995c4d8e19e1654887e237295bcb3a388be..21c48d756e96c1c660daf13e07ccfb6b1d0c74a1 100644 (file)
--- a/deal.II/lac/include/lac/precondition.h
+++ b/deal.II/lac/include/lac/precondition.h
@@ -18,6 +18,7 @@
 #include <base/config.h>
 #include <lac/vector_memory.h>
 #include <base/smartpointer.h>
+#include <base/template_constraints.h>
 
 DEAL_II_NAMESPACE_OPEN
 
@@ -928,16 +929,24 @@ PreconditionSSOR<MATRIX>::initialize (const MATRIX &rA,
 {
   this->PreconditionRelaxation<MATRIX>::initialize (rA, parameters);
 
+                                  // in case we have a SparseMatrix class,
+                                  // we can extract information about the
+                                  // diagonal.
+  const SparseMatrix<typename MATRIX::value_type> * mat = 
+    dynamic_cast<const SparseMatrix<typename MATRIX::value_type> *>(&*this->A);   
+
                                   // 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)
+  if (mat != 0)
     {
+      const std::size_t  * rowstart_ptr = 
+       mat->get_sparsity_pattern().get_rowstart_indices();
+      const unsigned int * const colnums = 
+       mat->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
@@ -945,11 +954,12 @@ PreconditionSSOR<MATRIX>::initialize (const MATRIX &rA,
                                       // 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;
+         pos_right_of_diagonal[row] = 
+           std::lower_bound (&colnums[*rowstart_ptr+1],
+                             &colnums[*(rowstart_ptr+1)],
+                             row)
+           - colnums;
+       }
     }
 }
 
@@ -971,7 +981,7 @@ inline void
 PreconditionSSOR<MATRIX>::Tvmult (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);
 }
 
 

diff --git a/deal.II/lac/include/lac/sparse_matrix.templates.h b/deal.II/lac/include/lac/sparse_matrix.templates.h
index a9e77117399b54792e04192b994e99ca96424221..dc3b73c97c7da4dfdf1a4371e24a24e073cab11d 100644 (file)
--- a/deal.II/lac/include/lac/sparse_matrix.templates.h
+++ b/deal.II/lac/include/lac/sparse_matrix.templates.h
@@ -1192,7 +1192,7 @@ SparseMatrix<number>::precondition_SSOR (Vector<somenumber>              &dst,
            pos_right_of_diagonal[row];
          Assert (first_right_of_diagonal_index <= *(rowstart_ptr+1),
                  ExcInternalError());
-         double s = *dst_ptr;
+         number s = *dst_ptr;
          for (unsigned int j=(*rowstart_ptr)+1; j<first_right_of_diagonal_index; ++j)
            s -= val[j] * dst(cols->colnums[j]);
 
@@ -1213,7 +1213,7 @@ SparseMatrix<number>::precondition_SSOR (Vector<somenumber>              &dst,
          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;
+         number s = *dst_ptr;
          for (unsigned int j=first_right_of_diagonal_index; j<end_row; ++j)
            s -= val[j] * dst(cols->colnums[j]);
       
@@ -1243,7 +1243,7 @@ SparseMatrix<number>::precondition_SSOR (Vector<somenumber>              &dst,
           -
           &cols->colnums[0]);
 
-      double s = *dst_ptr;
+      number s = *dst_ptr;
       for (unsigned int j=(*rowstart_ptr)+1; j<first_right_of_diagonal_index; ++j)
        s -= val[j] * dst(cols->colnums[j]);
 
@@ -1267,10 +1267,9 @@ SparseMatrix<number>::precondition_SSOR (Vector<somenumber>              &dst,
                                                                   &cols->colnums[end_row],
                                                                   static_cast<unsigned int>(row)) -
           &cols->colnums[0]);
-      double s = *dst_ptr;
+      number 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];
     };
 }

diff --git a/deal.II/lac/include/lac/sparse_matrix_ez.h b/deal.II/lac/include/lac/sparse_matrix_ez.h
index d49976abe35b5bab1be2ff165562ea9cc1f210e3..f89d46f87e6695ade2b52573e36f312391687758 100644 (file)
--- a/deal.II/lac/include/lac/sparse_matrix_ez.h
+++ b/deal.II/lac/include/lac/sparse_matrix_ez.h
@@ -713,7 +713,8 @@ class SparseMatrixEZ : public Subscriptor
     template <typename somenumber>
     void precondition_SSOR (Vector<somenumber>       &dst,
                            const Vector<somenumber> &src,
-                           const number              om = 1.) const;
+                           const number              om = 1.,
+                           const std::vector<unsigned int> &pos_right_of_diagonal = std::vector<unsigned int>()) const;
 
                                     /**
                                      * Apply SOR preconditioning matrix to @p src.

diff --git a/deal.II/lac/include/lac/sparse_matrix_ez.templates.h b/deal.II/lac/include/lac/sparse_matrix_ez.templates.h
index 0a069ec10a702187da883734ca3d4a148dc40123..9ff6264c6a723a3a9a70037613c5ccbdaf615e6c 100644 (file)
--- a/deal.II/lac/include/lac/sparse_matrix_ez.templates.h
+++ b/deal.II/lac/include/lac/sparse_matrix_ez.templates.h
@@ -326,7 +326,8 @@ template <typename somenumber>
 void
 SparseMatrixEZ<number>::precondition_SSOR (Vector<somenumber>       &dst,
                                           const Vector<somenumber> &src,
-                                          const number              om) const
+                                          const number              om,
+                                          const std::vector<unsigned int> &) const
 {
   Assert (m() == n(), ExcNotQuadratic());
   Assert (dst.size() == n(), ExcDimensionMismatch (dst.size(), n()));
@@ -344,9 +345,9 @@ SparseMatrixEZ<number>::precondition_SSOR (Vector<somenumber>       &dst,
       number s = *src_ptr;
       const unsigned int end_row = ri->start + ri->diagonal;
       for (unsigned int i=ri->start;i<end_row;++i)
-       s -= om * data[i].value * dst(data[i].column);
+       s -= data[i].value * dst(data[i].column);
       
-      *dst_ptr = s / data[ri->start + ri->diagonal].value;
+      *dst_ptr = s * om / data[ri->start + ri->diagonal].value;
     }
                                   // Diagonal
   dst_ptr = dst.begin();
@@ -361,9 +362,11 @@ SparseMatrixEZ<number>::precondition_SSOR (Vector<somenumber>       &dst,
    for (rri = row_info.rbegin(); rri != rend; --dst_ptr, ++rri)
      {
       const unsigned int end_row = rri->start + rri->length;
+      number s = *dst_ptr;
       for (unsigned int i=rri->start+rri->diagonal+1;i<end_row;++i)
-       *dst_ptr -= om * data[i].value * dst(data[i].column);
-      *dst_ptr /= data[rri->start + rri->diagonal].value;
+       s -= om * data[i].value * dst(data[i].column);
+
+      *dst_ptr = s * om / data[rri->start + rri->diagonal].value;
      }
 }
 

diff --git a/deal.II/lac/source/sparse_matrix_ez.inst.in b/deal.II/lac/source/sparse_matrix_ez.inst.in
index 29cd69e257df8196659b419d1e7c2eda716389eb..db7689a5bdcbdfad0b2d095e627450751dea7145 100644 (file)
--- a/deal.II/lac/source/sparse_matrix_ez.inst.in
+++ b/deal.II/lac/source/sparse_matrix_ez.inst.in
@@ -36,7 +36,8 @@ for (S1, S2 : REAL_SCALARS)
     template
       void SparseMatrixEZ<S1>::precondition_SSOR<S2> (Vector<S2> &,
                                                      const Vector<S2> &,
-                                                     const S1) const;
+                                                     const S1,
+                                                     const std::vector<unsigned int>&) const;
     template
       void SparseMatrixEZ<S1>::precondition_SOR<S2> (Vector<S2> &,
                                                     const Vector<S2> &,