[!269] Implement preconditioners with scalar matrix entries

Merge branch 'implement-preconditioners-with-scalar-matrix-entries' into 'master'

See merge request [!269]

  [!269]: Nonecore/dune-istl/merge_requests/269

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/CMakeLists.txt

@@ -49,6 +49,8 @@ dune_add_test(SOURCES iotest.cc)
 
 dune_add_test(SOURCES inverseoperator2prectest.cc)
 
+dune_add_test(SOURCES preconditionerstest.cc)
+
 dune_add_test(SOURCES scaledidmatrixtest.cc)
 
 dune_add_test(SOURCES solvertest.cc)

dune/istl/test/preconditionerstest.cc

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+/** \file \brief Test the preconditioners in the file `preconditioners.hh`
+ */
+
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/operators.hh>
+#include <dune/istl/preconditioners.hh>
+#include <dune/istl/solvers.hh>
+#include <dune/istl/test/laplacian.hh>
+
+using namespace Dune;
+
+template <class Matrix, class Vector>
+void setupProblem(Matrix& matrix, Vector& b)
+{
+  int n=100;
+
+  setupLaplacian(matrix,n);
+
+  b.resize(n*n);
+  Vector x(n*n);
+
+  // Construct right-hand side such that '1' is the solution
+  b=0;
+  x=1;
+  matrix.mv(x, b);
+}
+
+template <class Matrix, class Vector, class Preconditioner>
+void testPreconditioner(const Matrix& matrix, const Vector& b, Vector& x, Preconditioner& prec)
+{
+  using Operator = MatrixAdapter<Matrix,Vector,Vector>;
+
+  Operator linearOperator(matrix);
+  InverseOperatorResult result;
+
+  Dune::LoopSolver<Vector> solver(linearOperator, prec, 1e-8,15,2);
+  auto residual = b;
+  solver.apply(x,residual,result);
+}
+
+template <class Matrix, class Vector>
+void testAllPreconditioners(const Matrix& matrix, const Vector& b)
+{
+  Vector x = b;  // Set the correct size
+  x = 0;
+  SeqSSOR<Matrix,Vector,Vector> seqSSOR(matrix, 1,1.0);
+  testPreconditioner(matrix, b, x, seqSSOR);
+
+  x = 0;
+  SeqSOR<Matrix,Vector,Vector> seqSOR(matrix, 1,1.0);
+  testPreconditioner(matrix, b, x, seqSOR);
+
+  x = 0;
+  SeqGS<Matrix,Vector,Vector> seqGS(matrix, 1,1.0);
+  testPreconditioner(matrix, b, x, seqGS);
+
+  x = 0;
+  SeqJac<Matrix,Vector,Vector> seqJac(matrix, 1,1.0);
+  testPreconditioner(matrix, b, x, seqJac);
+
+  x = 0;
+  SeqILU<Matrix,Vector,Vector> seqILU(matrix, 3, 1.2, true);
+  testPreconditioner(matrix, b, x, seqILU);
+
+  x = 0;
+  Richardson<Vector,Vector> richardson(1.5);
+  testPreconditioner(matrix, b, x, richardson);
+
+  x = 0;
+  SeqILDL<Matrix,Vector,Vector> seqILDL(matrix, 1.2);
+  testPreconditioner(matrix, b, x, seqILDL);
+}
+
+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> >;
+    Matrix matrix;
+    Vector b,x;
+    setupProblem(matrix, b);
+
+    testAllPreconditioners(matrix, b);
+  }
+
+  return 0;
+}
+catch (std::exception& e) {
+  std::cerr << e.what() << std::endl;
+  return 1;
+}