I had introduced some subtle errors in the SSOR preconditioner implementation. This...

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

diff --git a/deal.II/lac/include/lac/sparse_matrix.templates.h b/deal.II/lac/include/lac/sparse_matrix.templates.h
index 3fe59f41ffe5c94d807b5c20cfc51c1dfb05c297..e5874f41e04d7e715e50c7dc9088cfafcd0b8eca 100644 (file)
--- a/deal.II/lac/include/lac/sparse_matrix.templates.h
+++ b/deal.II/lac/include/lac/sparse_matrix.templates.h
@@ -1198,12 +1198,13 @@ SparseMatrix<number>::precondition_SSOR (Vector<somenumber>              &dst,
            pos_right_of_diagonal[row];
          Assert (first_right_of_diagonal_index <= *(rowstart_ptr+1),
                  ExcInternalError());
-         number s = *dst_ptr;
+         number s = 0;
          for (unsigned int j=(*rowstart_ptr)+1; j<first_right_of_diagonal_index; ++j)
-           s -= val[j] * dst(cols->colnums[j]);
+           s += val[j] * dst(cols->colnums[j]);
 
                                   // divide by diagonal element
-         *dst_ptr = s * om / val[*rowstart_ptr];
+         *dst_ptr -= s * om;
+         *dst_ptr /= val[*rowstart_ptr];
        };
   
       rowstart_ptr = &cols->rowstart[0];
@@ -1219,11 +1220,12 @@ 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];
-         number s = *dst_ptr;
+         number s = 0;
          for (unsigned int j=first_right_of_diagonal_index; j<end_row; ++j)
-           s -= val[j] * dst(cols->colnums[j]);
+           s += val[j] * dst(cols->colnums[j]);
       
-         *dst_ptr = s * om / val[*rowstart_ptr];
+         *dst_ptr -= s * om;
+         *dst_ptr /= val[*rowstart_ptr];
        };
       return;
     }
@@ -1249,12 +1251,13 @@ SparseMatrix<number>::precondition_SSOR (Vector<somenumber>              &dst,
           -
           &cols->colnums[0]);
 
-      number s = *dst_ptr;
+      number s = 0;
       for (unsigned int j=(*rowstart_ptr)+1; j<first_right_of_diagonal_index; ++j)
-       s -= val[j] * dst(cols->colnums[j]);
+       s += val[j] * dst(cols->colnums[j]);
 
                                       // divide by diagonal element
-      *dst_ptr = s * om / val[*rowstart_ptr];
+      *dst_ptr -= s * om;
+      *dst_ptr /= val[*rowstart_ptr];
     };
   
   rowstart_ptr = &cols->rowstart[0];
@@ -1273,10 +1276,11 @@ SparseMatrix<number>::precondition_SSOR (Vector<somenumber>              &dst,
                                                                   &cols->colnums[end_row],
                                                                   static_cast<unsigned int>(row)) -
           &cols->colnums[0]);
-      number s = *dst_ptr;
+      number s = 0;
       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];
+       s += val[j] * dst(cols->colnums[j]);
+      *dst_ptr -= s * om;
+      *dst_ptr /= val[*rowstart_ptr];
     };
 }
 

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 9ff6264c6a723a3a9a70037613c5ccbdaf615e6c..64f553d6722f0431a65421270cda3e56c0014dbd 100644 (file)
--- a/deal.II/lac/include/lac/sparse_matrix_ez.templates.h
+++ b/deal.II/lac/include/lac/sparse_matrix_ez.templates.h
@@ -342,12 +342,13 @@ SparseMatrixEZ<number>::precondition_SSOR (Vector<somenumber>       &dst,
   for (ri = row_info.begin(); ri != end; ++dst_ptr, ++src_ptr, ++ri)
     {
       Assert (ri->diagonal != RowInfo::invalid_diagonal, ExcNoDiagonal());
-      number s = *src_ptr;
+      number s = 0;
       const unsigned int end_row = ri->start + ri->diagonal;
       for (unsigned int i=ri->start;i<end_row;++i)
-       s -= data[i].value * dst(data[i].column);
+       s += data[i].value * dst(data[i].column);
       
-      *dst_ptr = s * om / data[ri->start + ri->diagonal].value;
+      *dst_ptr = *src_ptr - s * om;
+      *dst_ptr /= data[ri->start + ri->diagonal].value;
     }
                                   // Diagonal
   dst_ptr = dst.begin();
@@ -362,11 +363,12 @@ 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;
+      number s = 0;
       for (unsigned int i=rri->start+rri->diagonal+1;i<end_row;++i)
-       s -= om * data[i].value * dst(data[i].column);
+       s += data[i].value * dst(data[i].column);
 
-      *dst_ptr = s * om / data[rri->start + rri->diagonal].value;
+      *dst_ptr -= s * om;
+      *dst_ptr /= data[rri->start + rri->diagonal].value;
      }
 }