Implement the standard preconditioners with scalar vector entries

dune/istl/gsetc.hh

@@ -8,6 +8,9 @@
 #include <iostream>
 #include <iomanip>
 #include <string>
+
+#include <dune/common/hybridutilities.hh>
+
 #include "multitypeblockvector.hh"
 #include "multitypeblockmatrix.hh"
 
@@ -385,12 +388,23 @@ namespace Dune {
         rhs = b[i.index()];           // rhs = b_i
         coliterator endj=(*i).end();
         coliterator j=(*i).begin();
+        Hybrid::ifElse(IsNumber<typename M::block_type>(),
+          [&](auto id) {
+            for (; j.index()<i.index(); ++j)
+              rhs -= id(*j) * x[j.index()];
+            coliterator diag=j++;           // *diag = a_ii and increment coliterator j from a_ii to a_i+1,i to skip diagonal
+            for (; j != endj; ++j)
+              rhs -= id(*j) * x[j.index()];
+            x[i.index()] = rhs / id(*diag);
+          },
+          [&](auto id) {
         for (; j.index()<i.index(); ++j)           // iterate over a_ij with j < i
-          (*j).mmv(x[j.index()],rhs);               // rhs -= sum_{j<i} a_ij * xnew_j
+          id(*j).mmv(x[j.index()],rhs);               // rhs -= sum_{j<i} a_ij * xnew_j
         coliterator diag=j++;           // *diag = a_ii and increment coliterator j from a_ii to a_i+1,i to skip diagonal
         for (; j != endj; ++j)           // iterate over a_ij with j > i
-          (*j).mmv(x[j.index()],rhs);               // rhs -= sum_{j>i} a_ij * xold_j
+          id(*j).mmv(x[j.index()],rhs);               // rhs -= sum_{j>i} a_ij * xold_j
         algmeta_itsteps<I-1,typename M::block_type>::dbgs(*diag,x[i.index()],rhs,w);           // if I==1: xnew_i = rhs/a_ii
+        });
       }
       // next two lines: xnew_i = w / a_ii * (b_i - sum_{j<i} a_ij * xnew_j - sum_{j>=i} a_ij * xold_j) + (1-w)*xold;
       x *= w;
@@ -417,13 +431,25 @@ namespace Dune {
         rhs = b[i.index()];           // rhs = b_i
         coliterator endj=(*i).end();           // iterate over a_ij with j < i
         coliterator j=(*i).begin();
+        Hybrid::ifElse(IsNumber<typename M::block_type>(),
+          [&](auto id) {
+            for (; j.index()<i.index(); ++j)
+              rhs -= id(*j) * x[j.index()];               //  rhs -= sum_{j<i} a_ij * xnew_j
+            coliterator diag=j;           // *diag = a_ii
+            for (; j!=endj; ++j)
+              rhs -= id(*j) * x[j.index()];               // rhs -= sum_{j<i} a_ij * xnew_j
+            v = rhs / id(*diag);
+            id(x[i.index()]) += w*id(v);           // x_i = w / a_ii * (b_i - sum_{j<i} a_ij * xnew_j - sum_{j>=i} a_ij * xold_j)
+          },
+          [&](auto id) {
         for (; j.index()<i.index(); ++j)
-          (*j).mmv(x[j.index()],rhs);               //  rhs -= sum_{j<i} a_ij * xnew_j
+          id(*j).mmv(x[j.index()],rhs);               //  rhs -= sum_{j<i} a_ij * xnew_j
         coliterator diag=j;           // *diag = a_ii
         for (; j!=endj; ++j)
-          (*j).mmv(x[j.index()],rhs);               // rhs -= sum_{j<i} a_ij * xnew_j
+          id(*j).mmv(x[j.index()],rhs);               // rhs -= sum_{j<i} a_ij * xnew_j
         algmeta_itsteps<I-1,typename M::block_type>::bsorf(*diag,v,rhs,w);           // if blocksize I==1: v = rhs/a_ii
-        x[i.index()].axpy(w,v);           // x_i = w / a_ii * (b_i - sum_{j<i} a_ij * xnew_j - sum_{j>=i} a_ij * xold_j)
+        id(x[i.index()]).axpy(w,v);           // x_i = w / a_ii * (b_i - sum_{j<i} a_ij * xnew_j - sum_{j>=i} a_ij * xold_j)
+          });
       }
     }
 
@@ -447,13 +473,25 @@ namespace Dune {
         rhs = b[i.index()];
         coliterator endj=(*i).end();
         coliterator j=(*i).begin();
+        Hybrid::ifElse(IsNumber<typename M::block_type>(),
+          [&](auto id) {
+            for (; j.index()<i.index(); ++j)
+              rhs -= id(*j) * x[j.index()];
+            coliterator diag=j;
+            for (; j!=endj; ++j)
+              rhs -= id(*j) * x[j.index()];
+            v = rhs / id(*diag);
+            x[i.index()] += w*id(v);
+          },
+          [&](auto id) {
         for (; j.index()<i.index(); ++j)
-          (*j).mmv(x[j.index()],rhs);
+          id(*j).mmv(x[j.index()],rhs);
         coliterator diag=j;
         for (; j!=endj; ++j)
-          (*j).mmv(x[j.index()],rhs);
+          id(*j).mmv(x[j.index()],rhs);
         algmeta_itsteps<I-1,typename M::block_type>::bsorb(*diag,v,rhs,w);
-        x[i.index()].axpy(w,v);
+        id(x[i.index()]).axpy(w,v);
+          });
       }
     }
 
@@ -473,12 +511,23 @@ namespace Dune {
         rhs = b[i.index()];
         coliterator endj=(*i).end();
         coliterator j=(*i).begin();
+        Hybrid::ifElse(IsNumber<typename M::block_type>(),
+          [&](auto id) {
+            for (; j.index()<i.index(); ++j)
+              rhs -= id(*j) * x[j.index()];
+            coliterator diag=j;
+            for (; j!=endj; ++j)
+              rhs -= id(*j) * x[j.index()];
+            v[i.index()] = rhs / id(*diag);
+          },
+          [&](auto id) {
         for (; j.index()<i.index(); ++j)
-          (*j).mmv(x[j.index()],rhs);
+          id(*j).mmv(x[j.index()],rhs);
         coliterator diag=j;
         for (; j!=endj; ++j)
-          (*j).mmv(x[j.index()],rhs);
+          id(*j).mmv(x[j.index()],rhs);
         algmeta_itsteps<I-1,typename M::block_type>::dbjac(*diag,v[i.index()],rhs,w);
+          });
       }
       x.axpy(w,v);
     }

dune/istl/test/preconditionerstest.cc

@@ -82,6 +82,16 @@ void testAllPreconditioners(const Matrix& matrix, const Vector& b)
 
 int main() try
 {
+  {
+    using Matrix = BCRSMatrix<double>;
+    using Vector = BlockVector<double>;
+    Matrix matrix;
+    Vector b,x;
+    setupProblem(matrix, b);
+
+    testAllPreconditioners(matrix, b);
+  }
+
   {
     using Matrix = BCRSMatrix<FieldMatrix<double,1,1> >;
     using Vector = BlockVector<FieldVector<double,1> >;