From: Wolfgang Bangerth <bangerth@math.tamu.edu>
Date: Tue, 7 Jul 1998 14:13:59 +0000 (+0000)
Subject: Implement a version of SSOR which is 40 per cent faster than the old version.
X-Git-Tag: v8.0.0~22816
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=7aca368b123b3de3eacae66f9b5b954ba31761d5;p=dealii.git

Implement a version of SSOR which is 40 per cent faster than the old version.


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

diff --git a/deal.II/lac/source/dsmatrix.cc b/deal.II/lac/source/dsmatrix.cc
index 1b9df3fadc..100f9f3a9e 100644
--- a/deal.II/lac/source/dsmatrix.cc
+++ b/deal.II/lac/source/dsmatrix.cc
@@ -586,31 +586,59 @@ void
 dSMatrix::precondition_SSOR (dVector& dst, const dVector& src,
 			     const double om) const
 {
+				   // to understand how this function works
+				   // you may want to take a look at the CVS
+				   // archives to see the original version
+				   // which is much clearer...
   Assert (cols != 0, ExcMatrixNotInitialized());
   Assert (val != 0, ExcMatrixNotInitialized());
   Assert (m() == n(), ExcMatrixNotSquare());
 
-  const unsigned int  n = src.size();
-  unsigned int  j;
+  const unsigned int  n            = src.size();
+  const unsigned int *rowstart_ptr = &cols->rowstart[0];
+  double             *dst_ptr      = &dst(0);
   
-  for (unsigned int row=0; row<n; ++row)
+  for (unsigned int row=0; row<n; ++row, ++dst_ptr, ++rowstart_ptr)
     {
-      dst(row) = src(row);
-      for (j=cols->rowstart[row]; j<cols->rowstart[row+1] ;j++)
-	if ((unsigned int)cols->colnums[j] < row)
-	  dst(row) -= om* val[j] * dst(cols->colnums[j]);
-      dst(row) /= val[cols->rowstart[row]];
-    }
-  for (unsigned int row=0; row<n; ++row)
-    dst(row) *= (2.-om)*val[cols->rowstart[row]];
+      *dst_ptr = src(row);
+				       // 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.
+      const unsigned int first_right_of_diagonal_index
+	= (lower_bound (&cols->colnums[*rowstart_ptr+1],
+			&cols->colnums[*(rowstart_ptr+1)],
+			static_cast<signed int>(row)) -
+	   &cols->colnums[0]);
+				       
+      for (unsigned int j=(*rowstart_ptr)+1; j<first_right_of_diagonal_index; ++j)
+	*dst_ptr -= om* val[j] * dst(cols->colnums[j]);
+      *dst_ptr /= val[*rowstart_ptr];
+    };
   
-  for (int row=n-1; row>=0; --row)
+  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];
+
+  rowstart_ptr = &cols->rowstart[n-1];
+  dst_ptr      = &dst(n-1);
+  for (int row=n-1; row>=0; --row, --rowstart_ptr, --dst_ptr)
     {
-      for (j=cols->rowstart[row];j<cols->rowstart[row+1];j++)
+      const unsigned int first_right_of_diagonal_index
+	= (lower_bound (&cols->colnums[*rowstart_ptr+1],
+			&cols->colnums[*(rowstart_ptr+1)],
+			static_cast<signed int>(row)) -
+	   &cols->colnums[0]);
+      for (unsigned int j=first_right_of_diagonal_index; j<*(rowstart_ptr+1); ++j)
 	if (cols->colnums[j] > row)
-	  dst(row) -= om* val[j] * dst(cols->colnums[j]);
-      dst(row) /= val[cols->rowstart[row]];
-    }
+	  *dst_ptr -= om* val[j] * dst(cols->colnums[j]);
+      
+      *dst_ptr /= val[*rowstart_ptr];
+    };
 }
 
 void