Compute block-residuals, instead of scalar ones

This transverses the matrix less, and uses the cache better.  On an example
with 6x6 blocks I can measure a speedup in the range of 30% for a full
multi-grid iteration (including transfer operators, too).

dune/solvers/iterationsteps/trustregiongsstep.cc

@@ -7,70 +7,49 @@ VectorType TrustRegionGSStep<MatrixType, VectorType>::getSol()
 
 template<class MatrixType, class VectorType>
 inline
-typename VectorType::field_type TrustRegionGSStep<MatrixType, VectorType>::
-residual(unsigned int index, int blockIndex) const
+void TrustRegionGSStep<MatrixType, VectorType>::iterate()
 {
     const MatrixType& mat = *this->mat_;
+    VectorType& x = *this->x_;
+    const std::vector<BoxConstraint<field_type,blocksize> >& obstacles = *this->obstacles_;
 
-    typedef typename MatrixType::row_type RowType;
-    typedef typename RowType::ConstIterator ColumnIterator;
-
-    /* The following loop computes
-     * \f[ sum_i = \sum_{j \ne i} A_{ij}x_j \f]
-     */
-    field_type r = (*this->rhs_)[index][blockIndex];
-    ColumnIterator cIt    = mat[index].template begin();
-    ColumnIterator cEndIt = mat[index].template end();
-
-    for (; cIt!=cEndIt; ++cIt) {
-        // r_i -= A_ij x_j
-        // r -= (*cIt)[blockIndex] * (*this->x_)[cIt.index()];
-        for (int i=0; i<blocksize; i++) {
-
-            // Omit diagonal element
-            if (cIt.index()==index && i==blockIndex)
-                continue;
+    for (size_t i=0; i<x.size(); i++) {
 
-            r -= (*cIt)[blockIndex][i] * (*this->x_)[cIt.index()][i];
+        // Dirichlet nodes in this block
+        std::bitset<blocksize> ignoreNodes(0);
+        if (this->ignoreNodes_)
+            ignoreNodes = (*this->ignoreNodes_)[i];
 
-        }
+        if (ignoreNodes.count() == blocksize)
+            continue;
 
-    }
+        VectorBlock blockResidual;
+        BlockGSStep<MatrixType,VectorType>::residual(i, blockResidual);
 
-    return r;
-}
-
-template<class MatrixType, class VectorType>
-inline
-void TrustRegionGSStep<MatrixType, VectorType>::iterate()
-{
-    const MatrixType& mat = *this->mat_;
-    const std::vector<BoxConstraint<field_type,blocksize> >& obstacles = *this->obstacles_;
-
-    for (size_t i=0; i<this->x_->size(); i++) {
+        mat[i][i].umv(x[i], blockResidual);
 
         // scalar gauss-seidel
         for (int j=0; j<blocksize; j++) {
 
-            if (this->ignoreNodes_ and (*this->ignoreNodes_)[i][j])
+            if (ignoreNodes[j])
                 continue;
 
             // Compute residual
-            field_type r = residual(i, j);
+            field_type r = blockResidual[j] - mat[i][i][j] * x[i];
+            r += mat[i][i][j][j] * x[i][j];
+
             const field_type& diag = mat[i][i][j][j];
 
             // Find local minimum
-            field_type& x = (*this->x_)[i][j];
-
             if (diag > 0) {
                 // 1d problem is convex
-                x = r / diag;
+                x[i][j] = r / diag;
 
-                if ( x < obstacles[i].lower(j)) {
-                    x = obstacles[i].lower(j);
+                if ( x[i][j] < obstacles[i].lower(j)) {
+                    x[i][j] = obstacles[i].lower(j);
                 }
-                else if ( x > obstacles[i].upper(j)) {
-                    x = obstacles[i].upper(j);
+                else if ( x[i][j] > obstacles[i].upper(j)) {
+                    x[i][j] = obstacles[i].upper(j);
                 }
 
             } else  {
@@ -87,7 +66,7 @@ void TrustRegionGSStep<MatrixType, VectorType>::iterate()
                 field_type lValue = 0.5*diag*lBound*lBound - r*lBound;
                 field_type uValue = 0.5*diag*uBound*uBound - r*uBound;
 
-                x = (lValue < uValue) ? lBound : uBound;
+                x[i][j] = (lValue < uValue) ? lBound : uBound;
 
             }
 

dune/solvers/iterationsteps/trustregiongsstep.hh

@@ -37,8 +37,6 @@
 
         virtual VectorType getSol();
 
-        field_type residual(unsigned int index, int blockIndex) const;
-
     };
 
 #include "trustregiongsstep.cc"