Merge branch 'feature/fastamg-nonblocked-support' into 'master'

Make FastAMG work with non-blocked matrices

See merge request !600

dune/istl/paamg/aggregates.hh

@@ -17,6 +17,7 @@
 #include <dune/common/sllist.hh>
 #include <dune/common/ftraits.hh>
 #include <dune/common/scalarmatrixview.hh>
+#include <dune/common/typetraits.hh>
 
 #include <utility>
 #include <set>
@@ -25,6 +26,7 @@
 #include <limits>
 #include <ostream>
 #include <tuple>
+#include <cmath>
 
 namespace Dune
 {
@@ -470,9 +472,13 @@ namespace Dune
        * @param m The matrix row to compute the norm of.
        */
       template<class M>
-      typename FieldTraits<typename M::field_type>::real_type operator()(const M& m) const
+      auto operator()(const M& m) const
       {
-        return m.infinity_norm();
+        using std::abs;
+        if constexpr(Dune::IsNumber<M>::value)
+          return abs(m);
+        else
+          return m.infinity_norm();
       }
     };
 

dune/istl/paamg/fastamg.hh

@@ -475,7 +475,12 @@ namespace Dune
           }
         }
         if(isDirichlet && hasDiagonal)
-          diag->solve(x[row.index()], b[row.index()]);
+        {
+          if constexpr (Dune::IsNumber<Block>::value)
+            x[row.index()] = b[row.index()]/(*diag);
+          else
+            diag->solve(x[row.index()], b[row.index()]);
+        }
       }
       if (verbosity_>0)
         std::cout<<" Preprocessing Dirichlet took "<<watch1.elapsed()<<std::endl;

dune/istl/paamg/fastamgsmoother.hh

@@ -6,6 +6,7 @@
 #define DUNE_ISTL_FASTAMGSMOOTHER_HH
 
 #include <cstddef>
+#include <dune/common/typetraits.hh>
 
 namespace Dune
 {
@@ -33,17 +34,31 @@ namespace Dune
           *dIter = *bIter;
 
           for (; col.index()<row.index(); ++col)
-            (*col).mmv(x[col.index()],*dIter);     // rhs -= sum_{j<i} a_ij * xnew_j
+          {
+            if constexpr (Dune::IsNumber<std::decay_t<decltype(*col)>>::value)
+              *dIter -= (*col)*x[col.index()];
+            else
+              (*col).mmv(x[col.index()],*dIter);     // rhs -= sum_{j<i} a_ij * xnew_j
+          }
           assert(row.index()==col.index());
           ColIterator diag=col;              // upper diagonal matrix not needed as x was 0 before.
 
           // Not recursive yet. Just solve with the diagonal
-          diag->solve(*xIter,*dIter);
+          if constexpr (Dune::IsNumber<std::decay_t<decltype(*diag)>>::value)
+            *xIter = (*dIter)/(*diag);
+          else
+            diag->solve(*xIter,*dIter);
+
           *dIter=0;   //as r=v
 
           // Update residual for the symmetric case
           for(col=(*row).begin(); col.index()<row.index(); ++col)
-            col->mmv(*xIter, d[col.index()]);     //d_j-=A_ij x_i
+          {
+            if constexpr (Dune::IsNumber<std::decay_t<decltype(*col)>>::value)
+              d[col.index()] -= (*col)*(*xIter);
+            else
+              col->mmv(*xIter, d[col.index()]);     //d_j-=A_ij x_i
+          }
         }
       }
     };
@@ -71,16 +86,29 @@ namespace Dune
           *dIter = *bIter;
 
           for (; col.index()>row.index(); --col)
-            (*col).mmv(x[col.index()],*dIter);     // rhs -= sum_{i>j} a_ij * xnew_j
+          {
+            if constexpr (Dune::IsNumber<std::decay_t<decltype(*col)>>::value)
+              *dIter -= (*col)*x[col.index()];
+            else
+              (*col).mmv(x[col.index()],*dIter);     // rhs -= sum_{j<i} a_ij * xnew_j
+          }
           assert(row.index()==col.index());
           ColIterator diag=col;
           YBlock v=*dIter;
           // upper diagonal matrix
           for (--col; col!=endCol; --col)
-            (*col).mmv(x[col.index()],v);     // v -= sum_{j<i} a_ij * xold_j
+          {
+            if constexpr (Dune::IsNumber<std::decay_t<decltype(*col)>>::value)
+              v -= (*col)*x[col.index()];
+            else
+              (*col).mmv(x[col.index()],v);     // v -= sum_{j<i} a_ij * xold_j
+          }
 
           // Not recursive yet. Just solve with the diagonal
-          diag->solve(*xIter,v);
+          if constexpr (Dune::IsNumber<std::decay_t<decltype(*diag)>>::value)
+            *xIter = v/(*diag);
+          else
+            diag->solve(*xIter,v);
 
           *dIter-=v;
 

dune/istl/paamg/test/fastamg.cc

@@ -29,7 +29,7 @@ void randomize(const M& mat, V& b)
 {
   V x=b;
 
-  srand((unsigned)std::clock());
+  srand(20030317);
 
   typedef typename V::iterator iterator;
   for(iterator i=x.begin(); i != x.end(); ++i)
@@ -38,7 +38,7 @@ void randomize(const M& mat, V& b)
   mat.mv(static_cast<const V&>(x), b);
 }
 
-template <int BS>
+template <class MatrixBlock, class VectorBlock>
 void testAMG(int N, int coarsenTarget, int ml)
 {
   std::cout<<"N="<<N<<" coarsenTarget="<<coarsenTarget<<" maxlevel="<<ml<<std::endl;
@@ -46,9 +46,7 @@ void testAMG(int N, int coarsenTarget, int ml)
   typedef Dune::ParallelIndexSet<int,LocalIndex,512> ParallelIndexSet;
 
   ParallelIndexSet indices;
-  typedef Dune::FieldMatrix<double,BS,BS> MatrixBlock;
   typedef Dune::BCRSMatrix<MatrixBlock> BCRSMat;
-  typedef Dune::FieldVector<double,BS> VectorBlock;
   typedef Dune::BlockVector<VectorBlock> Vector;
   typedef Dune::MatrixAdapter<BCRSMat,Vector,Vector> Operator;
   typedef Dune::Communication<void*> Comm;
@@ -137,8 +135,23 @@ try
   if(argc>3)
     ml = atoi(argv[3]);
 
-  testAMG<1>(N, coarsenTarget, ml);
-  testAMG<2>(N, coarsenTarget, ml);
+  {
+    using MB = double;
+    using VB = double;
+    testAMG<MB, VB>(N, coarsenTarget, ml);
+  }
+
+  {
+    using MB = Dune::FieldMatrix<double,1,1>;
+    using VB = Dune::FieldVector<double,1>;
+    testAMG<MB, VB>(N, coarsenTarget, ml);
+  }
+
+  {
+    using MB = Dune::FieldMatrix<double,2,2>;
+    using VB = Dune::FieldVector<double,2>;
+    testAMG<MB, VB>(N, coarsenTarget, ml);
+  }
 
   return 0;
 }