]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Implement a version of SSOR which is 40 per cent faster than the old version.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Tue, 7 Jul 1998 14:13:59 +0000 (14:13 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Tue, 7 Jul 1998 14:13:59 +0000 (14:13 +0000)
git-svn-id: https://svn.dealii.org/trunk@438 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/source/dsmatrix.cc

index 1b9df3fadc30e08e34e5895da92d44febaef905a..100f9f3a9e030c514590100c5743303db4c399ab 100644 (file)
@@ -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

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.