mtl block-matrix data-structure and some block-solvers added

dune/amdis/linear_algebra/mtl/BlockMTLMatrix.hpp

+/** \file BlockMTLMatrix.h */
+
+#pragma once
+
+#include <boost/numeric/mtl/matrices.hpp>
+
+#include <dune/amdis/Basic.hpp>
+#include <dune/amdis/Loops.hpp>
+
+namespace AMDiS
+{
+  /// A wrapper for AMDiS::SolverMatrix to be used in MTL/ITL solvers
+  template <class MTLMatrix, size_t _N, size_t _M>
+  class BlockMTLMatrix
+    : public std::array<std::array<MTLMatrix, _M>, _N>
+  {
+    using Self = BlockMTLMatrix;
+    
+  public:
+    using size_type  = typename MTLMatrix::size_type;
+    using value_type = typename MTLMatrix::value_type;
+
+    template <size_t R, size_t C>
+    auto& operator()(const index_<R>, const index_<C>)
+    {
+      static_assert(R < N() && C < M(), "Indices out of range [0,N)x[0,M)");
+      return std::get<C>(std::get<R>(*this));
+    }
+    
+    /// Return a shared pointer to the i'th underlaying base vector
+    template <size_t R, size_t C>
+    auto const& operator()(const index_<R>, const index_<C>) const
+    {
+      static_assert(R < N() && C < M(), "Indices out of range [0,N)x[0,M)");
+      return std::get<C>(std::get<R>(*this));
+    }
+    
+    static constexpr size_t N() { return _N; }
+    static constexpr size_t M() { return _M; }
+
+    /// perform blockwise multiplication A*b -> x
+    template <class VectorIn, class VectorOut, class Assign>
+    void mult(VectorIn const& b, VectorOut& x, Assign) const
+    {
+      // create iranges to access array blocks
+      std::array<mtl::irange, _N> r_rows;
+      std::array<mtl::irange, _M> r_cols;
+      getRanges(r_rows, r_cols);
+      
+      For<0, _N>::loop([this, &r_rows, &r_cols, &b, &x](const auto _r) {
+        bool first = true;
+        VectorOut x_i(x[r_rows[_r]]);
+        For<0, _M>::loop([this, &b, &x_i, &r_cols, &first, _r](const auto _c) {
+          auto const& block = this->operator()(_r, _c);
+          if (num_rows(block) > 0) {
+            const VectorIn b_j(b[r_cols[_c]]);
+            if (first) {
+              Assign::first_update(x_i, block * b_j);
+              first = false;
+            }
+            else {
+              Assign::update(x_i, block * b_j);
+            }
+          }
+        });
+      });
+    }
+
+    template <class VectorIn>
+    mtl::vector::mat_cvec_multiplier<Self, VectorIn> 
+    operator*(VectorIn const& v) const
+    {
+      return {*this, v};
+    }
+    
+    void getRowRanges(std::array<mtl::irange, _N>& r_rows) const
+    {      
+      size_t start = 0;
+      For<0, _N>::loop([&](const auto _r) {
+        size_t finish = start + num_rows((*this)(_r,index_<0>()));
+        r_rows[_r].set(start, finish);
+        start = finish;
+      });
+    }
+    
+    void getColRanges(std::array<mtl::irange, _M>& r_cols) const
+    {      
+      size_t start = 0;
+      For<0, _M>::loop([&](const auto _c) {
+        size_t finish = start + num_cols((*this)(index_<0>(),_c));
+        r_cols[_c].set(start, finish);
+        start = finish;
+      });
+    }
+    
+    void getRanges(std::array<mtl::irange, _N>& r_rows, std::array<mtl::irange, _M>& r_cols) const
+    {      
+      getRowRanges(r_rows);
+      getColRanges(r_cols);
+    }
+  };
+
+  template <class MTLMatrix, size_t _N, size_t _M>
+  inline size_t num_rows(BlockMTLMatrix<MTLMatrix, _N, _M> const& A)
+  {
+    size_t nRows = 0;
+    For<0, _N>::loop([&](const auto _r) {
+      nRows += num_rows(A(_r,index_<0>()));
+    });
+    return nRows;
+  }
+
+  template <class MTLMatrix, size_t _N, size_t _M>
+  inline size_t num_cols(BlockMTLMatrix<MTLMatrix, _N, _M> const& A)
+  {
+    size_t nCols = 0;
+    For<0, _M>::loop([&](const auto _c) {
+      nCols += num_cols(A(index_<0>(),_c));
+    });
+    return nCols;
+  }
+
+  template <class MTLMatrix, size_t _N, size_t _M>
+  inline size_t size(BlockMTLMatrix<MTLMatrix, _N, _M> const& A)
+  {
+    return num_rows(A) * num_cols(A);
+  }
+
+  template <class MTLMatrix, size_t _N, size_t _M>
+  inline void set_to_zero(BlockMTLMatrix<MTLMatrix, _N, _M>& A) 
+  {
+    For<0, _N>::loop([&](const auto _r) {
+      For<0, _M>::loop([&](const auto _c) {
+        set_to_zero(A(_r,_c));
+      });
+    });
+  }
+
+} // end namespace AMDiS
+
+namespace mtl
+{
+  template <class MTLMatrix, size_t _N, size_t _M>
+  struct Collection<AMDiS::BlockMTLMatrix<MTLMatrix, _N, _M>>
+  {
+    using value_type = typename MTLMatrix::value_type;
+    using size_type  = typename MTLMatrix::size_type;
+  };
+
+  namespace ashape
+  {
+    template <class MTLMatrix, size_t _N, size_t _M>
+    struct ashape_aux<AMDiS::BlockMTLMatrix<MTLMatrix, _N, _M>>
+    {
+      using type = nonscal;
+    };
+
+  } // end namespace ashape
+
+} // end namespace mtl

dune/amdis/linear_algebra/mtl/LinearSolver.hpp

@@ -6,6 +6,10 @@
 
 #include <dune/amdis/AdaptInfo.hpp>
 
+#include <dune/amdis/linear_algebra/mtl/itl/block_diagonal.hpp>
+#include <dune/amdis/linear_algebra/mtl/itl/fgmres.hpp>
+#include <dune/amdis/linear_algebra/mtl/itl/gmres2.hpp>
+
 namespace AMDiS
 {
   // creators for preconditioners, depending on matrix-type
@@ -131,7 +135,7 @@ namespace AMDiS
     }
     
     
-    void solve(AdaptInfo& adaptInfo, Matrix const& A, VectorX& x, VectorB const& b)
+    void solve(AdaptInfo& adaptInfo, Matrix const& A, VectorX& _x, VectorB const& _b)
     {      
       // create preconditioner
 //       PreconCreator<Matrix, VectorX, VectorB> preconCreator(name + "->precon");
@@ -157,28 +161,57 @@ namespace AMDiS
       size_t printCycle = 100;
       Parameters::get(name + "->print cycle", printCycle);
       
+      using BaseVector = mtl::dense_vector<typename Matrix::value_type>; 
+      BaseVector x(num_cols(A));
+      assign(_x, x);
+      BaseVector b(num_rows(A));
+      assign(_b, b);
       
-      auto const& A0 = A[0][0];
-      auto const& b0 = b[0];
-      auto& x0 = x[0];
-      
-      using Matrix0 = std::decay_t< decltype(A0) >;
-      using Vector0 = std::decay_t< decltype(b0) >;
-      
-      Vector0 r(b0); r -= A0*x0;
-      itl::cyclic_iteration<double> iter(r, maxIter, rTol, 0, printCycle);
+      BaseVector r(b); r -= A*x;
+      itl::cyclic_iteration<double> iter(two_norm(r), maxIter, rTol, 0, printCycle);
       iter.set_quite(info == 0);
 
       // no preconditioner
-      itl::pc::diagonal<Matrix0, typename Matrix0::value_type> P(A0);
+//       itl::pc::diagonal<Matrix> P(A);
+      itl::pc::identity<Matrix> P(A);
+      itl::pc::identity<Matrix> Pid(A);
+      
+      int restart = 30;
+      Parameters::get(name + "->restart", restart);
         
       // solve the linear system
-      itl::cg(A0, x0, b0, P, iter);
+      itl::fgmres(A, x, b, P, Pid, iter, restart);
+      
+      assign_back(_x, x);
       
       adaptInfo.setSolverIterations(iter.iterations());
       adaptInfo.setSolverResidual(iter.resid());
     }
     
+    template <class BlockVector, class BaseVector>
+    void assign(BlockVector const& bvec, BaseVector& vec)
+    {
+      size_t start = 0;
+      for (size_t r = 0; r < num_rows(bvec); ++r) {
+        size_t finish = start + num_rows(bvec[r]);
+        mtl::irange range(start, finish);
+        
+        vec[range] = bvec[r];
+      }
+    }
+    
+    template <class BlockVector, class BaseVector>
+    void assign_back(BlockVector& bvec, BaseVector const& vec)
+    {
+      size_t start = 0;
+      for (size_t r = 0; r < num_rows(bvec); ++r) {
+        size_t finish = start + num_rows(bvec[r]);
+        mtl::irange range(start, finish);
+        
+        bvec[r] = vec[range];
+      }
+    }
+    
   private:
     std::string name;
     std::string solverName;

dune/amdis/linear_algebra/mtl/SystemMatrix.hpp

@@ -11,6 +11,7 @@
 #include <dune/amdis/Basic.hpp>
 #include <dune/amdis/Log.hpp>
 
+#include <dune/amdis/linear_algebra/mtl/BlockMTLMatrix.hpp>
 #include <dune/amdis/linear_algebra/mtl/DOFMatrix.hpp>
 
 // for explicit instantiation
@@ -57,7 +58,7 @@ namespace AMDiS
           std::forward<RowFeSpace>(rowFeSpace), 
           std::get<C>(std::forward<ColFeSpaces>(colFeSpace)), 
           "matrix[" + std::to_string(R) + "][" + std::to_string(C) + "]",
-          multiMatrix[R][C] 
+          multiMatrix(index_<R>(), index_<C>()) 
           )...};
       }
     };
@@ -100,19 +101,13 @@ namespace AMDiS
     
     template <class RowFeSpace>
     using DOFMatrixGenerator = typename MultiMatrixRowGenerator<RowFeSpace>::DOFMatrices;
-      
-    template <size_t R, size_t C>
-    using Param = mtl::matrix::parameters<mtl::row_major, mtl::index::c_index, mtl::fixed::dimensions<R,C>>;
-      
-    template <class T, size_t S>
-    using FixedSizeMatrix = mtl::dense2D<T, Param<S,S>>;
     
   public:
     static constexpr size_t nComponent = std::tuple_size<FeSpaces>::value;
     AMDIS_STATIC_ASSERT( nComponent > 0 );
     
     using DOFMatrices  = MakeTuple_t<DOFMatrixGenerator, FeSpaces>;
-    using MultiMatrix  = FixedSizeMatrix<mtl::compressed2D<ValueType>, nComponent>;
+    using MultiMatrix  = BlockMTLMatrix<mtl::compressed2D<ValueType>, nComponent, nComponent>;
     
     
     /// Constructor that constructs new base-vectors
@@ -136,20 +131,20 @@ namespace AMDiS
     
     /// Return a shared pointer to the i'th underlaying base vector
     template <size_t R, size_t C>
-    auto& getMatrix(const index_<R> = index_<R>(), 
-                    const index_<C> = index_<C>())
+    auto& getMatrix(const index_<R> _r = index_<R>(), 
+                    const index_<C> _c = index_<C>())
     {
       static_assert(R < N() && C < M(), "Indices out of range [0,N)x[0,M)");
-      return matrix[R][C];
+      return matrix(_r, _c);
     }
     
     /// Return a shared pointer to the i'th underlaying base vector
     template <size_t R, size_t C>
-    auto const& getMatrix(const index_<R> = index_<R>(), 
-                          const index_<C> = index_<C>()) const
+    auto const& getMatrix(const index_<R> _r = index_<R>(), 
+                          const index_<C> _c = index_<C>()) const
     {
       static_assert(R < N() && C < M(), "Indices out of range [0,N)x[0,M)");
-      return matrix[R][C];
+      return matrix(_r, _c);
     }
     
     /// Return the underlying multi vector

dune/amdis/linear_algebra/mtl/itl/block_diagonal.hpp

+#pragma once
+
+#include <boost/numeric/itl/pc/diagonal.hpp>
+
+#include <dune/amdis/linear_algebra/mtl/BlockMTLMatrix.hpp>
+
+namespace itl
+{
+  namespace pc
+  {
+    /// Diagonal Preconditioner
+    template <class MTLMatrix, size_t _N, size_t _M, class Value>
+    class diagonal<AMDiS::BlockMTLMatrix<MTLMatrix, _N, _M>, Value>
+    {
+      using Self       = diagonal;
+      using MatrixType = AMDiS::BlockMTLMatrix<MTLMatrix, _N, _M>;
+      
+    public:
+      using value_type = Value;
+      using size_type  = typename mtl::Collection<MatrixType>::size_type;
+
+      /// Constructor takes matrix reference
+      explicit diagonal(const MatrixType& A) : inv_diag(num_rows(A))
+      {
+        MTL_THROW_IF(num_rows(A) != num_cols(A), mtl::matrix_not_square());
+        using math::reciprocal;
+        
+        std::array<mtl::irange, _N> r_rows;
+        A.getRowRanges(r_rows);
+
+        for (size_t i = 0; i < _N; i++)
+          inv_diag[r_rows[i]] = mtl::matrix::diagonal(A[i][i]);
+
+        for (size_type i = 0; i < num_rows(A); ++i)
+          inv_diag[i]= reciprocal(inv_diag[i]);
+      }
+
+      /// Member function solve, better use free function solve
+      template <typename Vector>
+      Vector solve(const Vector& x) const
+      {
+        Vector y(resource(x));
+        solve(x, y);
+        return y;
+      }
+
+      template <typename VectorIn, typename VectorOut>
+      void solve(const VectorIn& x, VectorOut& y) const
+      {
+        y.checked_change_resource(x);
+        MTL_THROW_IF(size(x) != size(inv_diag), mtl::incompatible_size());
+        for (size_type i= 0; i < size(inv_diag); ++i)
+          y[i]= inv_diag[i] * x[i];
+      }
+
+      /// Member function for solving adjoint problem, better use free function adjoint_solve
+      template <typename Vector>
+      Vector adjoint_solve(const Vector& x) const
+      {
+        Vector y(resource(x));
+        adjoint_solve(x, y);
+        return y;
+      }
+
+      template <typename VectorIn, typename VectorOut>
+      void adjoint_solve(const VectorIn& x, VectorOut& y) const
+      {
+        using mtl::conj;
+        y.checked_change_resource(x);
+        MTL_THROW_IF(size(x) != size(inv_diag), mtl::incompatible_size());
+        for (size_type i= 0; i < size(inv_diag); ++i)
+          y[i]= conj(inv_diag[i]) * x[i];
+      }
+
+    protected:
+      mtl::vector::dense_vector<value_type>    inv_diag;
+    };
+
+
+  }
+} // namespace itl::pc

dune/amdis/linear_algebra/mtl/itl/details.hpp

+/******************************************************************************
+ *
+ * AMDiS - Adaptive multidimensional simulations
+ *
+ * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
+ * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
+ *
+ * Authors:
+ * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
+ *
+ * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+ * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+ *
+ *
+ * This file is part of AMDiS
+ *
+ * See also license.opensource.txt in the distribution.
+ *
+ ******************************************************************************/
+
+// Written by Simon Praetorius
+
+
+#ifndef ITL_DETAIL_INCLUDE
+#define ITL_DETAIL_INCLUDE
+
+#include <boost/math/special_functions/sign.hpp>
+
+namespace itl
+{
+
+  namespace details
+  {
+
+    /// Compute the Givens rotation matrix parameters for a and b.
+    //     template<typename T>
+    //     void rotmat(const T& a, const T& b , T& c, T& s)
+    //     {
+    //       using std::abs; using std::sqrt; using mtl::conj;
+    //
+    //       const T zero = math::zero(T());
+    //       if (a == zero) {
+    // 	c = 0.0;
+    // 	s = 1.0;
+    //       } else {
+    // 	double temp = abs(a) / sqrt( conj(a)*a + conj(b)*b );
+    // 	c = temp;
+    // 	s = temp * (b / a);
+    //       }
+    //     }
+
+    inline void rotmat(const double& a, const double& b , double& c, double& s)
+    {
+      using std::abs;
+      using std::sqrt;
+      if ( b == 0.0 )
+      {
+        c = 1.0;
+        s = 0.0;
+      }
+      else if ( abs(b) > abs(a) )
+      {
+        double temp = a / b;
+        s = 1.0 / sqrt( 1.0 + temp*temp );
+        c = temp * s;
+      }
+      else
+      {
+        double temp = b / a;
+        c = 1.0 / sqrt( 1.0 + temp*temp );
+        s = temp * c;
+      }
+    }
+  }
+}
+
+#endif // ITL_DETAIL_INCLUDE

dune/amdis/linear_algebra/mtl/itl/fgmres.hpp

+/******************************************************************************
+ *
+ * AMDiS - Adaptive multidimensional simulations
+ *
+ * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
+ * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
+ *
+ * Authors:
+ * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
+ *
+ * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+ * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+ *
+ *
+ * This file is part of AMDiS
+ *
+ * See also license.opensource.txt in the distribution.
+ *
+ ******************************************************************************/
+
+// Written by Simon Praetorius
+
+
+#ifndef ITL_FGMRES_INCLUDE
+#define ITL_FGMRES_INCLUDE
+
+#include <algorithm>
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/vector/dense_vector.hpp>
+#include <boost/numeric/mtl/matrix/dense2D.hpp>
+#include <boost/numeric/mtl/matrix/multi_vector.hpp>
+#include <boost/numeric/mtl/operation/givens.hpp>
+#include <boost/numeric/mtl/operation/two_norm.hpp>
+#include <boost/numeric/mtl/utility/exception.hpp>
+#include <boost/numeric/mtl/utility/irange.hpp>
+
+#include "details.hpp"
+
+namespace itl
+{
+
+  /// Flexible Generalized Minimal Residual method (without restart)
+  /// Cite: Youcef Saad, A Flexible Inner-Outer Preconditioned GMRES Algorithm, 1993
+  /** It computes at most kmax_in iterations (or size(x) depending on what is smaller)
+      regardless on whether the termination criterion is reached or not.   **/
+  template <typename Matrix, typename Vector, typename Preconditioner, typename Iteration>
+  int fgmres_full(const Matrix& A, Vector& x, const Vector& b,
+                  Preconditioner& P, Iteration& iter)
+  {
+    using mtl::irange;
+    using mtl::iall;
+    using std::abs;
+    using std::sqrt;
+    using math::reciprocal;
+    typedef typename mtl::Collection<Vector>::value_type Scalar;
+    typedef typename mtl::Collection<Vector>::size_type  Size;
+
+    if (size(b) == 0) throw mtl::logic_error("empty rhs vector");
+
+    const Scalar                zero= math::zero(Scalar()), breakdown_tol= 1.e-16, kappa = 10.0;
+    Scalar                      rho, bnrm2, temp, hr;
+    Size                        k, kmax(std::min(size(x), Size(iter.max_iterations() - iter.iterations())));
+    Vector                      w(resource(x)), r(b - A*x);
+    mtl::matrix::multi_vector<Vector>   V(Vector(resource(x), zero), kmax+1);
+    mtl::matrix::multi_vector<Vector>   Z(Vector(resource(x), zero), kmax+1);
+    mtl::matrix::dense2D<Scalar>        H(kmax+1, kmax);
+
+    mtl::vector::dense_vector<Scalar>   sn(kmax, zero), cs(kmax, zero), s(kmax+1, zero), y(kmax, zero);  // replicated in distributed solvers
+
+    bnrm2 = two_norm(b);
+    if (bnrm2 == zero)
+    {
+      set_to_zero(x);
+      // b == 0 => solution = 0
+      return iter.terminate(bnrm2);
+    }
+
+    rho = two_norm(r);				// norm of preconditioned residual
+    if (iter.finished(rho))			// initial guess is good enough solution
+      return iter;
+
+    V.vector(0) = r * reciprocal(rho);
+    H = zero;
+    s[0] = rho;
+
+    // FGMRES iteration
+    for (k= 0; k < kmax && !iter.finished(rho) ; ++k, ++iter)
+    {
+      Z.vector(k) = solve(P, V.vector(k));
+      w = A * Z.vector(k);
+      temp = two_norm(w);
+
+      for (Size j= 0; j < k+1; j++)
+      {
+        H[j][k] = dot(V.vector(j), w);
+        w -= H[j][k] * V.vector(j);
+      }
+      H[k+1][k]= two_norm(w);
+
+      // reorthogonalization, only if "heuristic" condition is fulfilled
+      if (H[k+1][k] < temp * reciprocal(kappa))
+      {
+        for (Size i= 0; i < k+1; i++)
+        {
+          hr = dot(w, V.vector(i));
+          H[i][k] += hr;
+          w -= hr * V.vector(i);
+        }
+        temp = two_norm(w);
+
+        if (temp < H[k+1][k] * reciprocal(kappa))
+        {
+          set_to_zero(w);
+          H[k+1][k] = 0.0;
+        }
+        else
+        {
+          H[k+1][k] = temp;