vanka.hpp

// FEAT3: Finite Element Analysis Toolbox, Version 3
// Copyright (C) 2010 by Stefan Turek & the FEAT group
// FEAT3 is released under the GNU General Public License version 3,
// see the file 'copyright.txt' in the top level directory for details.
 
#pragma once
 
// includes, FEAT
#include <kernel/solver/base.hpp>
#include <kernel/adjacency/graph.hpp>
#include <kernel/lafem/dense_vector.hpp>
#include <kernel/lafem/dense_vector_blocked.hpp>
#include <kernel/lafem/power_vector.hpp>
#include <kernel/lafem/tuple_vector.hpp>
#include <kernel/lafem/sparse_matrix_csr.hpp>
#include <kernel/lafem/sparse_matrix_bcsr.hpp>
#include <kernel/lafem/power_diag_matrix.hpp>
#include <kernel/lafem/power_full_matrix.hpp>
#include <kernel/lafem/power_row_matrix.hpp>
#include <kernel/lafem/power_col_matrix.hpp>
#include <kernel/lafem/power_full_matrix.hpp>
#include <kernel/lafem/saddle_point_matrix.hpp>
#include <kernel/util/stop_watch.hpp>
 
// includes, system
#include <map>
#include <set>
#include <vector>
 
namespace FEAT
{
  namespace Solver
  {
    namespace Intern
    {
      template<typename Vector_>
      class VankaVector;
 
      template<typename Matrix_>
      class VankaMatrix;
 
      template<typename DT_, typename IT_>
      class VankaVector<LAFEM::DenseVector<DT_, IT_>>
      {
      public:
        typedef LAFEM::DenseVector<DT_, IT_> VectorType;
        static constexpr int dim = 1;
 
      protected:
        DT_* _vec_cor;
 
      public:
        explicit VankaVector(VectorType& vec_cor) :
          _vec_cor(vec_cor.elements())
        {
        }
 
        IT_ gather_def(DT_* x, const VectorType& vec, const IT_* idx, const IT_ n, const IT_ off) const
        {
          const DT_* vdef = vec.elements();
 
          for(IT_ i(0); i < n; ++i)
          {
            x[off+i] = vdef[idx[i]];
          }
          return off+n;
        }
 
        IT_ scatter_cor(const DT_ omega, const DT_* x, const IT_* idx, const IT_ n, const IT_ off)
        {
          for(IT_ i(0); i < n; ++i)
          {
            _vec_cor[idx[i]] += omega * x[off+i];
          }
          return off+n;
        }
      };
 
      template<typename DT_, typename IT_, int dim_>
      class VankaVector<LAFEM::DenseVectorBlocked<DT_, IT_, dim_>>
      {
      public:
        typedef LAFEM::DenseVectorBlocked<DT_, IT_, dim_> VectorType;
        typedef typename VectorType::ValueType ValueType;
        static constexpr int dim = dim_;
 
      protected:
        ValueType* _vec_cor;
 
      public:
        explicit VankaVector(VectorType& vec_cor) :
          _vec_cor(vec_cor.elements())
        {
        }
 
        IT_ gather_def(DT_* x, const VectorType& vec, const IT_* idx, const IT_ n, const IT_ off) const
        {
          const ValueType* vdef = vec.elements();
 
          for(IT_ i(0); i < n; ++i)
          {
            const ValueType& vi = vdef[idx[i]];
            for(int j(0); j < dim; ++j)
            {
              x[off + i*IT_(dim) + IT_(j)] = vi[j];
            }
          }
          return off + IT_(dim)*n;
        }
 
        IT_ scatter_cor(const DT_ omega, const DT_* x, const IT_* idx, const IT_ n, const IT_ off)
        {
          for(IT_ i(0); i < n; ++i)
          {
            ValueType& vi = _vec_cor[idx[i]];
            for(int j(0); j < dim; ++j)
            {
              vi[j] += omega * x[off + i*IT_(dim) + IT_(j)];
            }
          }
          return off + IT_(dim)*n;
        }
      };
 
      template<typename SubVector_, int dim_>
      class VankaVector<LAFEM::PowerVector<SubVector_, dim_>>
      {
      public:
        typedef LAFEM::PowerVector<SubVector_, dim_> VectorType;
        typedef VankaVector<SubVector_> FirstClass;
        typedef VankaVector<LAFEM::PowerVector<SubVector_, dim_-1>> RestClass;
 
        static constexpr int dim = FirstClass::dim + RestClass::dim;
 
      protected:
        FirstClass _first;
        RestClass _rest;
 
      public:
        explicit VankaVector(VectorType& vec_cor) :
          _first(vec_cor.first()),
          _rest(vec_cor.rest())
        {
        }
 
        template<typename DT_, typename IT_>
        IT_ gather_def(DT_* x, const VectorType& vec, const IT_* idx, const IT_ n, const IT_ off) const
        {
          IT_ noff = _first.gather_def(x, vec.first(), idx, n, off);
          return _rest.gather_def(x, vec.rest(), idx, n, noff);
        }
 
        template<typename DT_, typename IT_>
        IT_ scatter_cor(const DT_ omega, const DT_* x, const IT_* idx, const IT_ n, const IT_ off)
        {
          IT_ noff = _first.scatter_cor(omega, x, idx, n, off);
          return _rest.scatter_cor(omega, x, idx, n, noff);
        }
      };
 
      template<typename SubVector_>
      class VankaVector<LAFEM::PowerVector<SubVector_, 1>>
      {
      public:
        typedef LAFEM::PowerVector<SubVector_, 1> VectorType;
        typedef VankaVector<SubVector_> FirstClass;
 
        static constexpr int dim = FirstClass::dim;
 
      protected:
        FirstClass _first;
 
      public:
        explicit VankaVector(VectorType& vec_cor) :
          _first(vec_cor.first())
        {
        }
 
        template<typename DT_, typename IT_>
        IT_ gather_def(DT_* x, const VectorType& vec, const IT_* idx, const IT_ n, const IT_ off) const
        {
          return _first.gather_def(x, vec.first(), idx, n, off);
        }
 
        template<typename DT_, typename IT_>
        IT_ scatter_cor(const DT_ omega, const DT_* x, const IT_* idx, const IT_ n, const IT_ off)
        {
          return _first.scatter_cor(omega, x, idx, n, off);
        }
      };
 
      template<typename DT_, typename IT_>
      class VankaMatrix<LAFEM::SparseMatrixCSR<DT_, IT_>>
      {
      public:
        typedef LAFEM::SparseMatrixCSR<DT_, IT_> MatrixType;
        typedef LAFEM::DenseVector<DT_, IT_> VectorTypeR;
 
        static constexpr int row_dim = 1;
        static constexpr int col_dim = 1;
 
      protected:
        const IT_* _row_ptr;
        const IT_* _col_idx;
        const DT_* _mat_val;
 
      public:
        explicit VankaMatrix(const MatrixType& matrix) :
          _row_ptr(matrix.row_ptr()),
          _col_idx(matrix.col_ind()),
          _mat_val(matrix.val())
        {
        }
 
        std::pair<IT_, IT_> gather_full(
          DT_* data, const IT_* ridx, const IT_* cidx,
          const IT_ m, const IT_ n, const IT_ stride,
          const IT_ mo = IT_(0), const IT_ no = IT_(0)) const
        {
          // empty matrix?
          if(_mat_val == nullptr)
            return std::make_pair(mo+m, no+n);
 
          // loop over all local rows
          for(IT_ i(0); i < m; ++i)
          {
            // get row index
            const IT_ ri = ridx[i];
 
            // initialize loop variable for local columns
            IT_ j = IT_(0);
 
            // loop over the ri row of our matrix
            for(IT_ ai = _row_ptr[ri]; ai < _row_ptr[ri+1]; ++ai)
            {
              // get column index
              const IT_ rj = _col_idx[ai];
 
              // now search its position in our local matrix
              while((j < n) && (cidx[j] < rj))
              {
                ++j;
              }
              // did we find our local entry?
              if((j < n) && (cidx[j] == rj))
              {
                // found, so store it in our local matrix
                data[(mo + i) * stride + no + j] = _mat_val[ai];
              }
            }
          }
          return std::make_pair(mo+m, no+n);
        }
 
        IT_ gather_diag(DT_* data, const IT_* idx, const IT_ m, const IT_ mo = IT_(0)) const
        {
          // loop over all local rows
          for(IT_ i(0); i < m; ++i)
          {
            // get row index
            const IT_ ri = idx[i];
 
            // find diagonal entry
            for(IT_ ai = _row_ptr[ri]; ai < _row_ptr[ri+1]; ++ai)
            {
              // is it our diagonal?
              if(_col_idx[ai] == ri)
              {
                data[mo+i] = _mat_val[ai];
                break;
              }
            }
          }
          return mo + m;
        }
 
        IT_ mult_cor(DT_* x, const DT_ alpha, const VectorTypeR& vec_cor, const IT_* idx, const IT_ m, const IT_ off) const
        {
          if(_mat_val == nullptr)
            return off + m;
 
          const DT_* v = vec_cor.elements();
 
          // loop over all local rows
          for(IT_ i(0); i < m; ++i)
          {
            const IT_ vi = idx[i];
 
            // loop over all columns
            DT_ r = DT_(0);
            for(IT_ k = _row_ptr[vi]; k < _row_ptr[vi+1]; ++k)
            {
              r += _mat_val[k] * v[_col_idx[k]];
            }
            x[off+i] += alpha * r;
          }
          return off + m;
        }
      };
 
      // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
      // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
      // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
      template<typename DT_, typename IT_, int row_dim_, int col_dim_>
      class VankaMatrix<LAFEM::SparseMatrixBCSR<DT_, IT_, row_dim_, col_dim_>>
      {
      public:
        typedef LAFEM::SparseMatrixBCSR<DT_, IT_, row_dim_, col_dim_> MatrixType;
 
        typedef typename MatrixType::ValueType MatVal;
 
        static constexpr int row_dim = row_dim_;
        static constexpr int col_dim = col_dim_;
 
      protected:
        const IT_* _row_ptr;
        const IT_* _col_idx;
        const MatVal* _mat_val;
 
      public:
        explicit VankaMatrix(const MatrixType& matrix) :
          _row_ptr(matrix.row_ptr()),
          _col_idx(matrix.col_ind()),
          _mat_val(matrix.val())
        {
        }
 
        std::pair<IT_, IT_> gather_full(
          DT_* data, const IT_* ridx, const IT_* cidx,
          const IT_ m, const IT_ n, const IT_ stride,
          const IT_ mo = IT_(0), const IT_ no = IT_(0)) const
        {
          // empty matrix?
          XASSERTM(_mat_val != nullptr, "Vanka: invalid empty BCSR matrix");
 
          // loop over all local rows
          for(IT_ i(0); i < m; ++i)
          {
            // get row index
            const IT_ ri = ridx[i];
 
            // initialize loop variable for local columns
            IT_ j = IT_(0);
 
            // loop over the ri row of our matrix
            for(IT_ ai = _row_ptr[ri]; ai < _row_ptr[ri+1]; ++ai)
            {
              // get column index
              const IT_ rj = _col_idx[ai];
 
              // now search its position in our local matrix
              while((j < n) && (cidx[j] < rj))
              {
                ++j;
              }
              // did we find our local entry?
              if((j < n) && (cidx[j] == rj))
              {
                // found, so store it in our local matrix
                const MatVal& mv = _mat_val[ai];
 
                // copy matrix block
                for(int ii(0); ii < row_dim; ++ii)
                {
                  // compute offset
                  const IT_ lro = (mo + i * IT_(row_dim) + IT_(ii)) * stride + no + j * IT_(col_dim);
 
                  // copy block row
                  for(int jj(0); jj < col_dim; ++jj)
                  {
                    data[lro + IT_(jj)] = mv[ii][jj];
                  }
                }
              }
            }
          }
          return std::make_pair(mo + IT_(row_dim)*m, no + IT_(col_dim)*n);
        }
 
        IT_ gather_diag(DT_* data, const IT_* idx, const IT_ m, const IT_ mo = IT_(0)) const
        {
          // loop over all local rows
          for(IT_ i(0); i < m; ++i)
          {
            // get row index
            const IT_ ri = idx[i];
 
            // find diagonal entry
            for(IT_ ai = _row_ptr[ri]; ai < _row_ptr[ri+1]; ++ai)
            {
              // is it our diagonal?
              if(_col_idx[ai] == ri)
              {
                const MatVal& mv = _mat_val[ai];
 
                // copy block diagonal
                for(int k(0); k < row_dim; ++k)
                {
                  data[mo + i*IT_(row_dim) + IT_(k)] = mv[k][k];
                }
                break;
              }
            }
          }
          return mo + m;
        }
 
        IT_ mult_cor(DT_* x, const DT_ alpha, const LAFEM::DenseVectorBlocked<DT_, IT_, col_dim_>& vec_cor,
          const IT_* idx, const IT_ m, const IT_ off) const
        {
          if(_mat_val == nullptr)
            return off + m;
 
          typedef LAFEM::DenseVectorBlocked<DT_, IT_, col_dim_> VectorType;
          typedef typename VectorType::ValueType VecVal;
          const VecVal* v = vec_cor.elements();
 
          // loop over all local rows
          for(IT_ i(0); i < m; ++i)
          {
            const IT_ vi = idx[i];
 
            // loop over all columns
            for(IT_ k = _row_ptr[vi]; k < _row_ptr[vi+1]; ++k)
            {
              // get matrix and vector values
              const MatVal& mv = _mat_val[k];
              const VecVal& vv = v[_col_idx[k]];
 
              // loop over all block rows
              for(int ii(0); ii < row_dim; ++ii)
              {
                DT_ r = DT_(0);
                // loop over all block columns
                for(int jj(0); jj < col_dim; ++jj)
                {
                  r += mv[ii][jj] * vv[jj];
                }
                x[off + i*IT_(row_dim) + IT_(ii)] += alpha * r;
              }
            }
          }
          return off + m* IT_(row_dim);
        }
 
        IT_ mult_cor(DT_* x, const DT_ alpha, const LAFEM::DenseVector<DT_, IT_>& vec_cor,
          const IT_* idx, const IT_ m, const IT_ off) const
        {
          if(_mat_val == nullptr)
            return off + m;
 
          const DT_* v = vec_cor.elements();
 
          // loop over all local rows
          for(IT_ i(0); i < m; ++i)
          {
            const IT_ vi = idx[i];
 
            // loop over all columns
            for(IT_ k = _row_ptr[vi]; k < _row_ptr[vi+1]; ++k)
            {
              // get matrix and vector values
              const MatVal& mv = _mat_val[k];
              //const VecVal& vv = v[_col_idx[k]];
 
              // compute vector offset
              const IT_ vo = _col_idx[k] * IT_(col_dim);
 
              // loop over all block rows
              for(int ii(0); ii < row_dim; ++ii)
              {
                DT_ r = DT_(0);
                // loop over all block columns
                for(int jj(0); jj < col_dim; ++jj)
                {
                  r += mv[ii][jj] * v[vo + IT_(jj)];
                }
                x[off + i*IT_(row_dim) + IT_(ii)] += alpha * r;
              }
            }
          }
          return off + m* IT_(row_dim);
        }
      };
 
      // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
      // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
      // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 
      template<typename SubMatrix_, int dim_>
      class VankaMatrix<LAFEM::PowerDiagMatrix<SubMatrix_, dim_>>
      {
      public:
        typedef LAFEM::PowerDiagMatrix<SubMatrix_, dim_> MatrixType;
        typedef typename MatrixType::VectorTypeR VectorTypeR;
        typedef VankaMatrix<SubMatrix_> FirstClass;
        typedef VankaMatrix<LAFEM::PowerDiagMatrix<SubMatrix_, dim_-1>> RestClass;
 
        static constexpr int row_dim = FirstClass::row_dim + RestClass::row_dim;
        static constexpr int col_dim = FirstClass::col_dim + RestClass::col_dim;
 
      protected:
        FirstClass _first;
        RestClass _rest;
 
      public:
        explicit VankaMatrix(const MatrixType& matrix) :
          _first(matrix.first()),
          _rest(matrix.rest())
        {
        }
 
        template<typename DT_, typename IT_>
        std::pair<IT_, IT_> gather_full(
          DT_* data, const IT_* ridx, const IT_* cidx,
          const IT_ m, const IT_ n, const IT_ stride,
          const IT_ mo = IT_(0), const IT_ no = IT_(0)) const
        {
          std::pair<IT_, IT_> mn = _first.gather_full(data, ridx, cidx, m, n, stride, mo, no);
          return _rest.gather_full(data, ridx, cidx, m, n, stride, mn.first, mn.second);
        }
 
        template<typename DT_, typename IT_>
        IT_ gather_diag(DT_* data, const IT_* idx, const IT_ m, const IT_ mo = IT_(0)) const
        {
          IT_ mno = _first.gather_diag(data, idx, m, mo);
          return _rest.gather_diag(data, idx, m, mno);
        }
 
        template<typename DT_, typename IT_>
        IT_ mult_cor(DT_* x, const DT_ alpha, const VectorTypeR& vec_cor, const IT_* idx, const IT_ m, const IT_ off) const
        {
          IT_ noff = _first.mult_cor(x, alpha, vec_cor.first(), idx, m, off);
          return _rest.mult_cor(x, alpha, vec_cor.rest(), idx, m, noff);
        }
      };
 
      template<typename SubMatrix_>
      class VankaMatrix<LAFEM::PowerDiagMatrix<SubMatrix_, 1>>
      {
      public:
        typedef LAFEM::PowerDiagMatrix<SubMatrix_, 1> MatrixType;
        typedef typename MatrixType::VectorTypeR VectorTypeR;
        typedef VankaMatrix<SubMatrix_> FirstClass;
 
        static constexpr int row_dim = FirstClass::row_dim;
        static constexpr int col_dim = FirstClass::col_dim;
 
      protected:
        FirstClass _first;
 
      public:
        explicit VankaMatrix(const MatrixType& matrix) :
          _first(matrix.first())
        {
        }
 
        template<typename DT_, typename IT_>
        std::pair<IT_, IT_> gather_full(
          DT_* data, const IT_* ridx, const IT_* cidx,
          const IT_ m, const IT_ n, const IT_ stride,
          const IT_ mo = IT_(0), const IT_ no = IT_(0)) const
        {
          return _first.gather_full(data, ridx, cidx, m, n, stride, mo, no);
        }
 
        template<typename DT_, typename IT_>
        IT_ gather_diag(DT_* data, const IT_* idx, const IT_ m, const IT_ mo = IT_(0)) const
        {
          return _first.gather_diag(data, idx, m, mo);
        }
 
        template<typename DT_, typename IT_>
        IT_ mult_cor(DT_* x, const DT_ alpha, const VectorTypeR& vec_cor, const IT_* idx, const IT_ m, const IT_ off) const
        {
          return _first.mult_cor(x, alpha, vec_cor.first(), idx, m, off);
        }
      };
 
      template<typename SubMatrix_, int dim_>
      class VankaMatrix<LAFEM::PowerColMatrix<SubMatrix_, dim_>>
      {
      public:
        typedef LAFEM::PowerColMatrix<SubMatrix_, dim_> MatrixType;
        typedef typename MatrixType::VectorTypeR VectorTypeR;
        typedef VankaMatrix<SubMatrix_> FirstClass;
        typedef VankaMatrix<LAFEM::PowerColMatrix<SubMatrix_, dim_-1>> RestClass;
 
        static constexpr int row_dim = FirstClass::row_dim + RestClass::row_dim;
        static constexpr int col_dim = FirstClass::col_dim;
 
      protected:
        FirstClass _first;
        RestClass _rest;
 
      public:
        explicit VankaMatrix(const MatrixType& matrix) :
          _first(matrix.first()),
          _rest(matrix.rest())
        {
        }
 
        template<typename DT_, typename IT_>
        std::pair<IT_, IT_> gather_full(
          DT_* data, const IT_* ridx, const IT_* cidx,
          const IT_ m, const IT_ n, const IT_ stride,
          const IT_ mo = IT_(0), const IT_ no = IT_(0)) const
        {
          std::pair<IT_, IT_> mno = _first.gather_full(data, ridx, cidx, m, n, stride, mo, no);
          return _rest.gather_full(data, ridx, cidx, m, n, stride, mno.first, no);
        }
 
        template<int ro_, int co_, typename DT_, typename IT_>
        IT_ gather_diag_pfm(DT_* data, const IT_* idx, const IT_ m, const IT_ mo = IT_(0)) const
        {
          IT_ mno = _first.template gather_diag_pfm<ro_, co_>(data, idx, m, mo);
          return _rest.template gather_diag_pfm<ro_+1, co_>(data, idx, m, mno);
        }
 
        template<typename DT_, typename IT_>
        IT_ mult_cor(DT_* x, const DT_ alpha, const VectorTypeR& vec_cor, const IT_* idx, const IT_ m, const IT_ off) const
        {
          IT_ noff = _first.mult_cor(x, alpha, vec_cor, idx, m, off);
          return _rest.mult_cor(x, alpha, vec_cor, idx, m, noff);
        }
      };
 
      template<typename SubMatrix_>
      class VankaMatrix<LAFEM::PowerColMatrix<SubMatrix_, 1>>
      {
      public:
        typedef LAFEM::PowerColMatrix<SubMatrix_, 1> MatrixType;
        typedef typename MatrixType::VectorTypeR VectorTypeR;
        typedef VankaMatrix<SubMatrix_> FirstClass;
 
        static constexpr int row_dim = FirstClass::row_dim;
        static constexpr int col_dim = FirstClass::col_dim;
 
      protected:
        VankaMatrix<SubMatrix_> _first;
 
      public:
        explicit VankaMatrix(const MatrixType& matrix) :
          _first(matrix.first())
        {
        }
 
        template<typename DT_, typename IT_>
        std::pair<IT_, IT_> gather_full(
          DT_* data, const IT_* ridx, const IT_* cidx,
          const IT_ m, const IT_ n, const IT_ stride,
          const IT_ mo = IT_(0), const IT_ no = IT_(0)) const
        {
          return _first.gather_full(data, ridx, cidx, m, n, stride, mo, no);
        }
 
        template<int ro_, int co_, typename DT_, typename IT_>
        IT_ gather_diag_pfm(DT_* data, const IT_* idx, const IT_ m, const IT_ mo = IT_(0)) const
        {
          return _first.template gather_diag_pfm<ro_, co_>(data, idx, m, mo);
        }
 
        template<typename DT_, typename IT_>
        IT_ mult_cor(DT_* x, const DT_ alpha, const VectorTypeR& vec_cor, const IT_* idx, const IT_ m, const IT_ off) const
        {
          return _first.mult_cor(x, alpha, vec_cor, idx, m, off);
        }
      };
 
      template<typename SubMatrix_, int dim_>
      class VankaMatrix<LAFEM::PowerRowMatrix<SubMatrix_, dim_>>
      {
      public:
        typedef LAFEM::PowerRowMatrix<SubMatrix_, dim_> MatrixType;
        typedef typename MatrixType::VectorTypeR VectorTypeR;
        typedef VankaMatrix<SubMatrix_> FirstClass;
        typedef VankaMatrix<LAFEM::PowerRowMatrix<SubMatrix_, dim_-1>> RestClass;
 
        static constexpr int row_dim = FirstClass::row_dim;
        static constexpr int col_dim = FirstClass::col_dim + RestClass::col_dim;
 
      protected:
        FirstClass _first;
        RestClass _rest;
 
      public:
        explicit VankaMatrix(const MatrixType& matrix) :
          _first(matrix.first()),
          _rest(matrix.rest())
        {
        }
 
        template<typename DT_, typename IT_>
        std::pair<IT_, IT_> gather_full(
          DT_* data, const IT_* ridx, const IT_* cidx,
          const IT_ m, const IT_ n, const IT_ stride,
          const IT_ mo = IT_(0), const IT_ no = IT_(0)) const
        {
          std::pair<IT_, IT_> mno = _first.gather_full(data, ridx, cidx, m, n, stride, mo, no);
          return _rest.gather_full(data, ridx, cidx, m, n, stride, mo, mno.second);
        }
 
        template<int ro_, int co_, typename DT_, typename IT_>
        IT_ gather_diag_pfm(DT_* data, const IT_* idx, const IT_ m, const IT_ mo = IT_(0)) const
        {
          if(ro_ == co_)
            return _first.gather_diag(data, idx, m, mo);
          else
            return _rest.template gather_diag_pfm<ro_, co_+1>(data, idx, m, mo);
        }
 
        template<typename DT_, typename IT_>
        IT_ mult_cor(DT_* x, const DT_ alpha, const VectorTypeR& vec_cor, const IT_* idx, const IT_ m, const IT_ off) const
        {
          _first.mult_cor(x, alpha, vec_cor.first(), idx, m, off);
          return _rest.mult_cor(x, alpha, vec_cor.rest(), idx, m, off);
        }
      };
 
      template<typename SubMatrix_>
      class VankaMatrix<LAFEM::PowerRowMatrix<SubMatrix_, 1>>
      {
      public:
        typedef LAFEM::PowerRowMatrix<SubMatrix_, 1> MatrixType;
        typedef VankaMatrix<SubMatrix_> FirstClass;
        typedef typename MatrixType::VectorTypeR VectorTypeR;
 
        static constexpr int row_dim = FirstClass::row_dim;
        static constexpr int col_dim = FirstClass::col_dim;
 
      protected:
        FirstClass _first;
 
      public:
        explicit VankaMatrix(const MatrixType& matrix) :
          _first(matrix.first())
        {
        }
 
        template<typename DT_, typename IT_>
        std::pair<IT_, IT_> gather_full(
          DT_* data, const IT_* ridx, const IT_* cidx,
          const IT_ m, const IT_ n, const IT_ stride,
          const IT_ mo = IT_(0), const IT_ no = IT_(0)) const
        {
          return _first.gather_full(data, ridx, cidx, m, n, stride, mo, no);
        }
 
        template<int ro_, int co_, typename DT_, typename IT_>
        IT_ gather_diag_pfm(DT_* data, const IT_* idx, const IT_ m, const IT_ mo = IT_(0)) const
        {
          if(ro_ == co_)
            return _first.gather_diag(data, idx, m, mo);
          else
            return mo;
        }
 
        template<typename DT_, typename IT_>
        IT_ mult_cor(DT_* x, const DT_ alpha, const VectorTypeR& vec_cor, const IT_* idx, const IT_ m, const IT_ off) const
        {
          return _first.mult_cor(x, alpha, vec_cor.first(), idx, m, off);
        }
      };
 
      template<typename SubMatrix_, int dim_w_, int dim_h_>
      class VankaMatrix<LAFEM::PowerFullMatrix<SubMatrix_, dim_w_, dim_h_>>
      {
      public:
        typedef LAFEM::PowerFullMatrix<SubMatrix_, dim_w_, dim_h_> MatrixType;
        typedef typename MatrixType::VectorTypeR VectorTypeR;
        typedef VankaMatrix<typename MatrixType::ContClass> ContClass;
 
        static constexpr int row_dim = ContClass::row_dim;
        static constexpr int col_dim = ContClass::col_dim;
 
      protected:
        ContClass _cont;
 
      public:
        explicit VankaMatrix(const MatrixType& matrix) :
          _cont(matrix.get_container())
        {
        }
 
        template<typename DT_, typename IT_>
        std::pair<IT_, IT_> gather_full(
          DT_* data, const IT_* ridx, const IT_* cidx,
          const IT_ m, const IT_ n, const IT_ stride,
          const IT_ mo = IT_(0), const IT_ no = IT_(0)) const
        {
          return _cont.gather_full(data, ridx, cidx, m, n, stride, mo, no);
        }
 
        template<typename DT_, typename IT_>
        IT_ gather_diag(DT_* data, const IT_* idx, const IT_ m, const IT_ mo = IT_(0)) const
        {
          return _cont.template gather_diag_pfm<0,0>(data, idx, m, mo);
        }
 
        template<typename DT_, typename IT_>
        IT_ mult_cor(DT_* x, const DT_ alpha, const VectorTypeR& vec_cor, const IT_* idx, const IT_ m, const IT_ off) const
        {
          return _cont.mult_cor(x, alpha, vec_cor, idx, m, off);
        }
      };
 
      template<typename DT_, typename IT_>
      std::pair<const IT_*, const IT_*> vanka_graph(const LAFEM::SparseMatrixCSR<DT_, IT_>& matrix)
      {
        return std::make_pair(matrix.row_ptr(), matrix.col_ind());
      }
 
      template<typename DT_, typename IT_, int m_, int n_>
      std::pair<const IT_*, const IT_*> vanka_graph(const LAFEM::SparseMatrixBCSR<DT_, IT_, m_, n_>& matrix)
      {
        return std::make_pair(matrix.row_ptr(), matrix.col_ind());
      }
 
      template<typename DT_, typename IT_, int dim_>
      std::pair<const IT_*, const IT_*> vanka_graph(
        const LAFEM::PowerColMatrix<LAFEM::SparseMatrixCSR<DT_, IT_>, dim_>& matrix)
      {
        const auto& m = matrix.first();
        return std::make_pair(m.row_ptr(), m.col_ind());
      }
 
      template<typename DT_, typename IT_, int dim_>
      std::pair<const IT_*, const IT_*> vanka_graph(
        const LAFEM::PowerRowMatrix<LAFEM::SparseMatrixCSR<DT_, IT_>, dim_>& matrix)
      {
        const auto& m = matrix.first();
        return std::make_pair(m.row_ptr(), m.col_ind());
      }
 
      template<typename DT_, typename IT_, int dim_>
      std::pair<const IT_*, const IT_*> vanka_graph(
        const LAFEM::PowerDiagMatrix<LAFEM::SparseMatrixCSR<DT_, IT_>, dim_>& matrix)
      {
        const auto& m = matrix.first();
        return std::make_pair(m.row_ptr(), m.col_ind());
      }
 
      template<typename DT_, typename IT_, int row_dim_, int col_dim_>
      std::pair<const IT_*, const IT_*> vanka_graph(
        const LAFEM::PowerFullMatrix<LAFEM::SparseMatrixCSR<DT_, IT_>, row_dim_, col_dim_>& matrix)
      {
        const auto& m = matrix.template at<0,0>();
        return std::make_pair(m.row_ptr(), m.col_ind());
      }
 
    } // namespace Intern
 
    enum class VankaType
    {
      nodal_diag_mult = 0x000,
      nodal_full_mult = 0x001,
      block_diag_mult = 0x010,
      block_full_mult = 0x011,
      nodal_diag_add  = 0x100,
      nodal_full_add  = 0x101,
      block_diag_add  = 0x110,
      block_full_add  = 0x111
    };
 
    class VankaFactorError :
      public SolverException
    {
    public:
      VankaFactorError() : SolverException("Vanka Factorization Error") {}
    };
 
    template<typename Matrix_, typename Filter_>
    class Vanka;
 
    template<typename MatrixA_, typename MatrixB_, typename MatrixD_, typename Filter_>
    class Vanka<LAFEM::SaddlePointMatrix<MatrixA_, MatrixB_, MatrixD_>, Filter_> :
      public SolverBase<LAFEM::TupleVector<typename MatrixB_::VectorTypeL, typename MatrixD_::VectorTypeL>>
    {
    public:
      typedef LAFEM::SaddlePointMatrix<MatrixA_, MatrixB_, MatrixD_> MatrixType;
      typedef Filter_ FilterType;
 
      typedef typename MatrixType::VectorTypeR VectorType;
      typedef typename MatrixType::DataType DataType;
      typedef typename MatrixType::IndexType IndexType;
 
      typedef typename MatrixD_::VectorTypeR VectorV;
      typedef typename MatrixB_::VectorTypeR VectorP;
 
    protected:
      // our Vanka matrix types
      typedef Intern::VankaMatrix<MatrixA_> VankaMatrixA;
      typedef Intern::VankaMatrix<MatrixB_> VankaMatrixB;
      typedef Intern::VankaMatrix<MatrixD_> VankaMatrixD;
 
      // dimension sanity checks
      static_assert(VankaMatrixA::row_dim == VankaMatrixA::col_dim, "Matrix A has invalid dimensions");
      static_assert(VankaMatrixA::row_dim == VankaMatrixB::row_dim, "Matrices A and B have incompatible dimensions");
      static_assert(VankaMatrixA::col_dim == VankaMatrixD::col_dim, "Matrices A and D have incompatible dimensions");
      static_assert(VankaMatrixB::col_dim == VankaMatrixD::row_dim, "Matrices B and D have incompatible dimensions");
 
      // for now, B and D cannot have more than 1 pressure dimension...
      static_assert(VankaMatrixB::col_dim == 1, "Invalid pressure space dimension");
 
      const MatrixType& _matrix;
      const FilterType& _filter;
      VankaType _type;
      DataType _omega;
      Index _num_iter;
      IndexType _degree_v;
      IndexType _degree_p;
      std::vector<IndexType> _block_v_ptr, _block_v_idx;
      std::vector<IndexType> _block_p_ptr, _block_p_idx;
      std::vector<DataType> _data;
      std::vector<DataType> _vdef, _vcor;
      VectorType _vec_scale, _vec_tmp1, _vec_tmp2;
 
      // stop watch for symbolic factorization
      StopWatch watch_init_symbolic;
      // stop watch for numeric factorization
      StopWatch watch_init_numeric;
      // stop watch for apply time
      StopWatch watch_apply;
 
    public:
      explicit Vanka(const MatrixType& matrix, const FilterType& filter, VankaType type,
        DataType omega = DataType(1), Index num_iter = Index(1)) :
        _matrix(matrix),
        _filter(filter),
        _type(type),
        _omega(omega),
        _num_iter(num_iter)
      {
      }
 
      virtual String name() const override
      {
        return "Vanka";
      }
 
      virtual void init_symbolic() override
      {
        watch_init_symbolic.start();
 
        bool block = ((int(_type) & 0x010) != 0);
        bool multi = ((int(_type) & 0x100) == 0);
 
        // compute pressure block graph
        if(block)
        {
          this->_build_p_block();
        }
        else
        {
          this->_build_p_nodal();
        }
 
        // compute velocity block graph
        this->_build_v_block();
 
        // allocate memory for numerical factorization
        this->_alloc_data();
 
        // allocate temporary vector for additive
        if(!multi)
        {
          this->_vec_scale = this->_matrix.create_vector_r();
          this->_vec_tmp1 = this->_matrix.create_vector_r();
          this->_vec_tmp2 = this->_matrix.create_vector_r();
        }
 
        watch_init_symbolic.stop();
      }
 
      virtual void done_symbolic() override
      {
        _vec_tmp2.clear();
        _vec_tmp1.clear();
        _vec_scale.clear();
        _vdef.clear();
        _vcor.clear();
        _data.clear();
        _block_v_ptr.clear();
        _block_v_idx.clear();
        _block_p_ptr.clear();
        _block_p_idx.clear();
      }
 
      virtual void init_numeric() override
      {
        watch_init_numeric.start();
 
        bool full = ((int(_type) & 0x001) != 0);
        bool multi = ((int(_type) & 0x100) == 0);
 
        if(full)
        {
          this->_factor_full();
        }
        else
        {
          this->_factor_diag();
        }
 
        if(!multi)
        {
          this->_calc_scale();
        }
 
        watch_init_numeric.stop();
      }
 
      virtual Status apply(VectorType& vec_cor, const VectorType& vec_def) override
      {
        watch_apply.start();
 
        bool full = ((int(_type) & 0x001) != 0);
        if(full)
        {
          this->_apply_full(vec_cor, vec_def);
        }
        else
        {
          this->_apply_diag(vec_cor, vec_def);
        }
 
        watch_apply.stop();
 
        return Status::success;
      }
 
      std::size_t data_size() const
      {
        return _data.size();
      }
 
      std::size_t bytes() const
      {
        return sizeof(IndexType) * (_block_v_ptr.size() + _block_v_idx.size() + _block_p_ptr.size() + _block_p_idx.size())
          + sizeof(DataType) * (_data.size() + _vdef.size() + _vcor.size())
          + _vec_scale.bytes() + _vec_tmp1.bytes() + _vec_tmp2.bytes();
      }
 
      void reset_timings()
      {
        watch_init_symbolic.reset();
        watch_init_numeric.reset();
        watch_apply.reset();
      }
 
      double time_init_symbolic() const
      {
        return watch_init_symbolic.elapsed();
      }
 
      double time_init_numeric() const
      {
        return watch_init_numeric.elapsed();
      }
 
      double time_apply() const
      {
        return watch_apply.elapsed();
      }
 
    protected:
      void _build_p_block()
      {
        // fetch matrix dimensions
        const IndexType m = IndexType(_matrix.block_d().rows());
 
        // fetch the matrix arrays
        auto graph_b = Intern::vanka_graph(_matrix.block_b());
        const IndexType* row_ptr_b = graph_b.first;
        const IndexType* col_idx_b = graph_b.second;
        auto graph_d = Intern::vanka_graph(_matrix.block_d());
        const IndexType* row_ptr_d = graph_d.first;
        const IndexType* col_idx_d = graph_d.second;
 
        // clear block arrays
        _block_p_ptr.clear();
        _block_p_idx.clear();
        _block_p_ptr.push_back(IndexType(0));
 
        // local map
        std::map<IndexType, int> map_s;
 
        // allocate mask vector
        std::vector<int> mask(std::size_t(m), 0);
 
        // loop over all pressure nodes
        for(Index i(0); i < m; ++i)
        {
          // is this node already processed?
          if(mask[i] != 0)
            continue;
 
          // clear the local map
          map_s.clear();
 
          // okay, loop over all entries of D
          for(IndexType kd = row_ptr_d[i]; kd < row_ptr_d[i+1]; ++kd)
          {
            // fetch the column index of D
            const IndexType col_d = col_idx_d[kd];
 
            // loop over the row of B
            for(IndexType kb = row_ptr_b[col_d]; kb < row_ptr_b[col_d+1]; ++kb)
            {
              // fetch the column index of B
              const IndexType col_b = col_idx_b[kb];
 
              // insert into map
              auto ib = map_s.emplace(col_b, 1);
              if(!ib.second)
              {
                // node already exists, increment counter then
                ++(ib.first->second);
              }
            }
          }
 
          // compute the map degree
          int ideg = 0;
          for(auto it = map_s.begin(); it != map_s.end(); ++it)
            ideg = Math::max(ideg, it->second);
 
          // loop over all map entries with maximum degree
          for(auto it = map_s.begin(); it != map_s.end(); ++it)
          {
            if(ideg == it->second)
            {
              // insert into map
              _block_p_idx.push_back(it->first);
              mask[it->first] = 1;
            }
          }
 
          // push row
          _block_p_ptr.push_back(IndexType(_block_p_idx.size()));
        }
      }
 
      void _build_p_nodal()
      {
        const IndexType m = IndexType(_matrix.block_d().rows());
 
        // clear block arrays
        _block_p_ptr.clear();
        _block_p_idx.clear();
        _block_p_ptr.reserve(m+1);
        _block_p_idx.reserve(m);
        for(IndexType i(0); i < m; ++i)
        {
          _block_p_ptr.push_back(i);
          _block_p_idx.push_back(i);
        }
        _block_p_ptr.push_back(m);
      }
 
      void _build_v_block()
      {
        // fetch the matrix arrays
        auto graph_d = Intern::vanka_graph(_matrix.block_d());
        const IndexType* row_ptr_d = graph_d.first;
        const IndexType* col_idx_d = graph_d.second;
 
        // fetch number of pressure blocks
        const IndexType m = IndexType(_block_p_ptr.size()-1);
 
        // clear block arrays
        _block_v_ptr.clear();
        _block_v_idx.clear();
        _block_v_ptr.reserve(m+1);
        _block_v_ptr.push_back(IndexType(0));
 
        std::set<IndexType> set_v;
 
        // loop over all pressure blocks
        for(IndexType i(0); i < m; ++i)
        {
          // loop over all pressure dofs in the current block
          for(IndexType j = _block_p_ptr[i]; j < _block_p_ptr[i+1]; ++j)
          {
            // get the pressure dof index
            const IndexType pix = _block_p_idx[j];
 
            // now loop over the corresponding row of D
            for(IndexType k = row_ptr_d[pix]; k < row_ptr_d[pix+1]; ++k)
            {
              // insert velocity dof into block set
              set_v.insert(col_idx_d[k]);
            }
          }
 
          // push velocity block
          for(auto it = set_v.begin(); it != set_v.end(); ++it)
          {
            _block_v_idx.push_back(*it);
          }
 
          // update pointer
          _block_v_ptr.push_back(IndexType(_block_v_idx.size()));
 
          // clear local set
          set_v.clear();
        }
      }
 
      void _alloc_data()
      {
        // get the dimension of our system
        const IndexType dim = IndexType(Intern::VankaMatrix<MatrixA_>::row_dim);
 
        // use diagonal?
        bool diag = ((int(_type) & 1) == 0);
 
        // number of non-zero entries
        IndexType nze = IndexType(0);
 
        // reset degrees
        _degree_v = _degree_p = IndexType(0);
 
        // loop over all blocks
        const IndexType m = IndexType(_block_p_ptr.size()-1);
        for(IndexType i(0); i < m; ++i)
        {
          // get number of velocity and pressure dofs
          IndexType nv = _block_v_ptr[i+1] - _block_v_ptr[i];
          IndexType np = _block_p_ptr[i+1] - _block_p_ptr[i];
 
          // update degrees
          _degree_v = Math::max(nv, _degree_v);
          _degree_p = Math::max(np, _degree_p);
 
          // update count
          if(diag)
          {
            // main diagonal of A
            nze += dim * nv;
            // matrices B and D
            nze += IndexType(2) * dim * nv * np;
            // Schur-Complement matrix S
            nze += np * np;
          }
          else
          {
            // dense local systems
            nze += Math::sqr(dim*nv + np);
          }
        }
 
        // finally, allocate our memory
        _data.resize(nze);
        _vdef.resize(dim*_degree_v + _degree_p);
        _vcor.resize(dim*_degree_v + _degree_p);
      }
 
      void _calc_scale()
      {
        // create a Vanka Vector
        this->_vec_scale.format();
        Intern::VankaVector<VectorV> vanka_v(this->_vec_scale.template at<0>());
        Intern::VankaVector<VectorP> vanka_p(this->_vec_scale.template at<1>());
 
        // get our data arrays
        const IndexType* vptr = _block_v_ptr.data();
        const IndexType* vidx = _block_v_idx.data();
        const IndexType* pptr = _block_p_ptr.data();
        const IndexType* pidx = _block_p_idx.data();
 
        // create a vector of ones
        std::vector<DataType> vec_one(vanka_v.dim*_degree_v + _degree_p, DataType(1));
        const DataType* vone = vec_one.data();
 
        // loop over all blocks
        const IndexType nblocks = IndexType(_block_v_ptr.size() - 1);
        for(IndexType iblock(0); iblock < nblocks; ++iblock)
        {
          // get number of velocity and pressure dofs
          const IndexType nv = vptr[iblock+1] - vptr[iblock];
          const IndexType np = pptr[iblock+1] - pptr[iblock];
 
          // scatter ones
          vanka_v.scatter_cor(DataType(1), vone, &vidx[vptr[iblock]], nv, IndexType(0));
          vanka_p.scatter_cor(DataType(1), vone, &pidx[pptr[iblock]], np, IndexType(0));
        }
 
        // invert components
        this->_vec_scale.component_invert(this->_vec_scale);
      }
 
      void _factor_full()
      {
        Intern::VankaMatrix<MatrixA_> vanka_a(_matrix.block_a());
        Intern::VankaMatrix<MatrixB_> vanka_b(_matrix.block_b());
        Intern::VankaMatrix<MatrixD_> vanka_d(_matrix.block_d());
 
        // get the dimension of our system
        const IndexType dim = IndexType(vanka_a.row_dim);
 
        // get our data arrays
        DataType* data = _data.data();
        const IndexType* vptr = _block_v_ptr.data();
        const IndexType* vidx = _block_v_idx.data();
        const IndexType* pptr = _block_p_ptr.data();
        const IndexType* pidx = _block_p_idx.data();
 
        // format block data
        ::memset(data, 0, sizeof(DataType) * _data.size());
 
        // allocate pivot array
        std::vector<IndexType> pivot(3*(dim*_degree_v + _degree_p));
 
        // current block offset
        IndexType block_offset = IndexType(0);
 
        // loop over all blocks
        const IndexType nblocks = IndexType(_block_v_ptr.size() - 1);
        for(IndexType iblock(0); iblock < nblocks; ++iblock)
        {
          // get number of velocity and pressure dofs
          const IndexType nv = vptr[iblock+1] - vptr[iblock];
          const IndexType np = pptr[iblock+1] - pptr[iblock];
 
          // get our local indices
          const IndexType* loc_vidx = &vidx[vptr[iblock]];
          const IndexType* loc_pidx = &pidx[pptr[iblock]];
 
          // compute matrix stride
          const IndexType n = dim*nv + np;
          const IndexType block_size = n*n;
 
          // get our block data array pointer
          DataType* block_data = &data[block_offset];
 
          // gather matrix a
          std::pair<IndexType,IndexType> ao =
            vanka_a.gather_full(block_data, loc_vidx, loc_vidx, nv, nv, n);
          vanka_b.gather_full(block_data, loc_vidx, loc_pidx, nv, np, n, IndexType(0), ao.second);
          vanka_d.gather_full(block_data, loc_pidx, loc_vidx, np, nv, n, ao.first, IndexType(0));
 
          // invert local matrix block
          /*DataType ret =*/ Math::invert_matrix(n, n, block_data, pivot.data());
          //if(!Math::isnormal(ret))
          //{
            // Don't throw any exceptions, as this often leads to false alerts
            //throw VankaFactorError();
          //}
 
          // increment block data offset
          block_offset += block_size;
        }
      }
 
      void _apply_full(VectorType& vec_cor, const VectorType& vec_def)
      {
        // format correction vector
        vec_cor.format();
 
        // do we use the multiplicative variant?
        const bool multi = ((int(_type) & 0x100) == 0);
 
        // additive?
        if(!multi)
        {
          this->_vec_tmp1.copy(vec_def);
          this->_vec_tmp2.format();
        }
 
        // get our sub-vectors
        VectorV& vec_cv = (multi ? vec_cor.template at<0>() : this->_vec_tmp2.template at<0>());
        VectorP& vec_cp = (multi ? vec_cor.template at<1>() : this->_vec_tmp2.template at<1>());
        const VectorV& vec_dv = (multi ? vec_def.template at<0>() : this->_vec_tmp1.template at<0>());
        const VectorP& vec_dp = (multi ? vec_def.template at<1>() : this->_vec_tmp1.template at<1>());
 
        // create our vanka matrix and vector objects
        Intern::VankaMatrix<MatrixA_> vanka_a(_matrix.block_a());
        Intern::VankaMatrix<MatrixB_> vanka_b(_matrix.block_b());
        Intern::VankaMatrix<MatrixD_> vanka_d(_matrix.block_d());
        Intern::VankaVector<VectorV> vanka_v(vec_cv);
        Intern::VankaVector<VectorP> vanka_p(vec_cp);
 
        // get velocity vector dimension
        const IndexType velo_dim = IndexType(vanka_v.dim);
 
        // get block count
        const IndexType num_blocks = IndexType(_block_v_ptr.size()-1);
 
        // get block data arrays
        const IndexType* vptr = _block_v_ptr.data();
        const IndexType* vidx = _block_v_idx.data();
        const IndexType* pptr = _block_p_ptr.data();
        const IndexType* pidx = _block_p_idx.data();
 
        // get local data and vectors
        const DataType* lmat_data = _data.data();
        DataType* lcor = _vcor.data();
        DataType* ldef = _vdef.data();
 
        Index flops(0);
 
        // iterate
        for(IndexType iter(0); iter < _num_iter; ++iter)
        {
          // additive variant?
          if((!multi) && (iter > IndexType(0)))
          {
            // compute current defect
            this->_matrix.apply(_vec_tmp1, vec_cor, vec_def, -DataType(1));
            _vec_tmp2.format();
          }
 
          TimeStamp stamp_kernel;
 
          // reset block offset
          IndexType block_offset = IndexType(0);
 
          // loop over all blocks
          for(IndexType iblock(0); iblock < num_blocks; ++iblock)
          {
            // get local sizes
            const IndexType nv = vptr[iblock+1] - vptr[iblock];
            const IndexType np = pptr[iblock+1] - pptr[iblock];
 
            // get our local indices
            const IndexType* loc_vidx = &vidx[vptr[iblock]];
            const IndexType* loc_pidx = &pidx[pptr[iblock]];
 
            // compute local degree
            const IndexType n = velo_dim * nv + np;
 
            // get our local matrix data
            const DataType* lmat = &lmat_data[block_offset];
 
            // get local vectors
            DataType* lcor_v = lcor;
            DataType* lcor_p = &lcor[velo_dim * nv];
            DataType* ldef_v = ldef;
            DataType* ldef_p = &ldef[velo_dim * nv];
 
            // gather local defect
            vanka_v.gather_def(ldef_v, vec_dv, loc_vidx, nv, IndexType(0));
            vanka_p.gather_def(ldef_p, vec_dp, loc_pidx, np, IndexType(0));
 
            // multiplicative variants only:
            if(multi)
            {
              // subtract A*x
              vanka_a.mult_cor(ldef_v, -DataType(1), vec_cv, loc_vidx, nv, IndexType(0));
              vanka_b.mult_cor(ldef_v, -DataType(1), vec_cp, loc_vidx, nv, IndexType(0));
              vanka_d.mult_cor(ldef_p, -DataType(1), vec_cv, loc_pidx, np, IndexType(0));
            }
 
            // solve local system
            for(IndexType i(0); i < n; ++i)
            {
              lcor[i] = DataType(0);
              for(IndexType j(0); j < n; ++j)
              {
                lcor[i] += lmat[i*n + j] * ldef[j];
              }
            }
            flops += 2*n*n;
 
            // scatter result
            vanka_v.scatter_cor(_omega, lcor_v, loc_vidx, nv, IndexType(0));
            vanka_p.scatter_cor(_omega, lcor_p, loc_pidx, np, IndexType(0));
 
            // update block offset
            block_offset += n*n;
          }
 
          Statistics::add_time_precon(stamp_kernel.elapsed_now());
 
          // additive variant?
          if(!multi)
          {
            // multiply by scaling vector
            this->_vec_tmp2.component_product(this->_vec_tmp2, this->_vec_scale);
 
            // update correction vector
            vec_cor.axpy(this->_vec_tmp2);
          }
 
          // apply filter
          _filter.filter_cor(vec_cor);
        }
 
        Statistics::add_flops(flops);
      }
 
      void _factor_diag()
      {
        Intern::VankaMatrix<MatrixA_> vanka_a(_matrix.block_a());
        Intern::VankaMatrix<MatrixB_> vanka_b(_matrix.block_b());
        Intern::VankaMatrix<MatrixD_> vanka_d(_matrix.block_d());
 
        // get the dimension of our system
        const IndexType dim = IndexType(vanka_a.row_dim);
 
        // get our data arrays
        DataType* data = _data.data();
        const IndexType* vptr = _block_v_ptr.data();
        const IndexType* vidx = _block_v_idx.data();
        const IndexType* pptr = _block_p_ptr.data();
        const IndexType* pidx = _block_p_idx.data();
 
        // format block data
        ::memset(data, 0, sizeof(DataType) * _data.size());
 
        // allocate pivot array (pressure only)
        std::vector<IndexType> pivot(3*_degree_p);
 
        // current block offset
        IndexType block_offset = IndexType(0);
 
        // loop over all blocks
        const IndexType nblocks = IndexType(_block_v_ptr.size() - 1);
        for(IndexType iblock(0); iblock < nblocks; ++iblock)
        {
          // get number of velocity and pressure dofs
          const IndexType nv = vptr[iblock+1] - vptr[iblock];
          const IndexType np = pptr[iblock+1] - pptr[iblock];
          const IndexType dnv = dim*nv;
 
          // get our local indices
          const IndexType* loc_vidx = &vidx[vptr[iblock]];
          const IndexType* loc_pidx = &pidx[pptr[iblock]];
 
          // compute local block size
          const IndexType block_size = dnv + IndexType(2) * dnv * np + np*np;
 
          // get our block data array pointer
          DataType* block_data = &data[block_offset];
 
          // get our local matrices
          DataType* loc_a =  block_data;
          DataType* loc_b = &block_data[dnv];
          DataType* loc_d = &block_data[dnv + dnv*np];
          DataType* loc_s = &block_data[dnv + IndexType(2)*dnv*np];
 
          // gather diag(A)
          vanka_a.gather_diag(loc_a, loc_vidx, nv);
          // gather B and D
          vanka_b.gather_full(loc_b, loc_vidx, loc_pidx, nv, np, np);
          vanka_d.gather_full(loc_d, loc_pidx, loc_vidx, np, nv, dnv);
 
          // invert diag(A) and pre-multiply D by diag(A)^{-1}
          for(IndexType i(0); i < dnv; ++i)
          {
            // invert a_ii
            loc_a[i] = DataType(1) / loc_a[i];
 
            // make sure we have a normal value
            if(!Math::isnormal(loc_a[i]))
            {
              throw VankaFactorError();
            }
 
            // pre-multiply D by diag(A)^{-1}
            for(IndexType j(0); j < np; ++j)
            {
              loc_d[j*dnv + i] *= loc_a[i];
            }
          }
 
          // calculate Schur-complement of A:
          // S := -D * diag(A)^{-1} * B
          for(IndexType i(0); i < np; ++i)
          {
            for(IndexType j(0); j < np; ++j)
            {
              DataType s = DataType(0);
              for(IndexType k(0); k < dnv; ++k)
              {
                // Note: D is already pre-multiplied by diag(A)^{-1}
                s += loc_d[i*dnv + k] * loc_b[k*np + j];
              }
              loc_s[i*np + j] = -s;
            }
          }
 
          // invert local matrix S
          DataType det = Math::invert_matrix(np, np, loc_s, pivot.data());
          if(!Math::isnormal(det))
          {
            throw VankaFactorError();
          }
 
          // increment block data offset
          block_offset += block_size;
        }
      }
 
      void _apply_diag(VectorType& vec_cor, const VectorType& vec_def)
      {
        // format correction vector
        vec_cor.format();
 
        // do we use the multiplicative variant?
        const bool multi = ((int(_type) & 0x100) == 0);
 
        // additive?
        if(!multi)
        {
          this->_vec_tmp1.copy(vec_def);
          this->_vec_tmp2.format();
        }
 
        // get our sub-vectors
        VectorV& vec_cv = (multi ? vec_cor.template at<0>() : this->_vec_tmp2.template at<0>());
        VectorP& vec_cp = (multi ? vec_cor.template at<1>() : this->_vec_tmp2.template at<1>());
        const VectorV& vec_dv = (multi ? vec_def.template at<0>() : this->_vec_tmp1.template at<0>());
        const VectorP& vec_dp = (multi ? vec_def.template at<1>() : this->_vec_tmp1.template at<1>());
 
        // create our vanka matrix and vector objects
        Intern::VankaMatrix<MatrixA_> vanka_a(_matrix.block_a());
        Intern::VankaMatrix<MatrixB_> vanka_b(_matrix.block_b());
        Intern::VankaMatrix<MatrixD_> vanka_d(_matrix.block_d());
        Intern::VankaVector<VectorV> vanka_v(vec_cv);
        Intern::VankaVector<VectorP> vanka_p(vec_cp);
 
        // get velocity vector dimension
        const IndexType velo_dim = IndexType(vanka_v.dim);
 
        // get block count
        const IndexType num_blocks = IndexType(_block_v_ptr.size()-1);
 
        // get block data arrays
        const IndexType* vptr = _block_v_ptr.data();
        const IndexType* vidx = _block_v_idx.data();
        const IndexType* pptr = _block_p_ptr.data();
        const IndexType* pidx = _block_p_idx.data();
 
        // get local data and vectors
        const DataType* data = _data.data();
        DataType* lcor = _vcor.data();
        DataType* ldef = _vdef.data();
 
        Index flops(0);
 
        // iterate
        for(IndexType iter(0); iter < _num_iter; ++iter)
        {
          // additive variant?
          if((!multi) && (iter > IndexType(0)))
          {
            // compute current defect
            this->_matrix.apply(_vec_tmp1, vec_cor, vec_def, -DataType(1));
            _vec_tmp2.format();
          }
 
          TimeStamp stamp_kernel;
 
          // reset block offset
          IndexType block_offset = IndexType(0);
 
          // loop over all blocks
          for(IndexType iblock(0); iblock < num_blocks; ++iblock)
          {
            // get local sizes
            const IndexType nv = vptr[iblock+1] - vptr[iblock];
            const IndexType np = pptr[iblock+1] - pptr[iblock];
            const IndexType dnv = velo_dim*nv;
 
            // get our local indices
            const IndexType* loc_vidx = &vidx[vptr[iblock]];
            const IndexType* loc_pidx = &pidx[pptr[iblock]];
 
            // compute local block size
            const IndexType block_size = dnv + IndexType(2) * dnv * np + np*np;
 
            // get our block data array pointer
            const DataType* block_data = &data[block_offset];
 
            // get our local matrices
            const DataType* loc_a =  block_data;
            const DataType* loc_b = &block_data[dnv];
            const DataType* loc_d = &block_data[dnv + dnv*np];
            const DataType* loc_s = &block_data[dnv + IndexType(2)*dnv*np];
 
            // get local vectors
            DataType* lcor_v = lcor;
            DataType* lcor_p = &lcor[dnv];
            DataType* ldef_v = ldef;
            DataType* ldef_p = &ldef[dnv];
 
            // gather local defect
            vanka_v.gather_def(ldef_v, vec_dv, loc_vidx, nv, IndexType(0));
            vanka_p.gather_def(ldef_p, vec_dp, loc_pidx, np, IndexType(0));
 
            // multiplicative variants only:
            if(multi)
            {
              // subtract A*x
              vanka_a.mult_cor(ldef_v, -DataType(1), vec_cv, loc_vidx, nv, IndexType(0));
              vanka_b.mult_cor(ldef_v, -DataType(1), vec_cp, loc_vidx, nv, IndexType(0));
              vanka_d.mult_cor(ldef_p, -DataType(1), vec_cv, loc_pidx, np, IndexType(0));
            }
 
            // update pressure RHS:
            // g_p := f_p - D * diag(A)^{-1} * f_u
            for(IndexType i(0); i < np; ++i)
            {
              DataType r = DataType(0);
              for(IndexType j(0); j < dnv; ++j)
              {
                // Note: D is already pre-multiplied by diag(A)^{-1}
                r += loc_d[i*dnv + j] * ldef_v[j];
              }
              ldef_p[i] -= r;
            }
            flops += 2*np*dnv;
 
            // solve pressure:
            // p := S^{-1} * g_p
            for(IndexType i(0); i < np; ++i)
            {
              DataType r = DataType(0);
              for(IndexType j(0); j < np; ++j)
              {
                r += loc_s[i*np + j] * ldef_p[j];
              }
              lcor_p[i] = r;
            }
            flops += 2*np*np;
 
            // update velocity RHS and solve velocity
            for(IndexType i(0); i < dnv; ++i)
            {
              DataType xb = DataType(0);
              for(IndexType j(0); j < np; ++j)
              {
                xb += loc_b[i*np + j] * lcor_p[j];
              }
              // solve: u := diag(A)^{-1} * (f_u - B*p)
              lcor_v[i] = loc_a[i] * (ldef_v[i] - xb);
            }
            flops += dnv * (2*np + 2);
 
            // scatter result
            vanka_v.scatter_cor(_omega, lcor_v, loc_vidx, nv, IndexType(0));
            vanka_p.scatter_cor(_omega, lcor_p, loc_pidx, np, IndexType(0));
 
            // update block offset
            block_offset += block_size;
          }
 
          Statistics::add_time_precon(stamp_kernel.elapsed_now());
 
          // additive variant?
          if(!multi)
          {
            // multiply by scaling vector
            this->_vec_tmp2.component_product(this->_vec_tmp2, this->_vec_scale);
 
            // update correction vector
            vec_cor.axpy(this->_vec_tmp2);
          }
 
          // apply filter
          _filter.filter_cor(vec_cor);
        }
 
        Statistics::add_flops(flops);
      }
    }; // class Vanka<...>
 
    template<typename MatrixA_, typename MatrixB_, typename MatrixD_, typename Filter_>
    std::shared_ptr<Vanka<LAFEM::SaddlePointMatrix<MatrixA_, MatrixB_, MatrixD_>, Filter_>> new_vanka(
      const LAFEM::SaddlePointMatrix<MatrixA_, MatrixB_, MatrixD_>& matrix,
      const Filter_& filter,
      VankaType type,
      typename MatrixA_::DataType omega = typename MatrixA_::DataType(1),
      Index num_iter = Index(1))
    {
      return std::make_shared<Vanka<LAFEM::SaddlePointMatrix<MatrixA_, MatrixB_, MatrixD_>, Filter_>>
        (matrix, filter, type, omega, num_iter);
    }
  } // namespace Solver
} // namespace FEAT