tensor_product_matrix.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2017 - 2018 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#ifndef dealii_tensor_product_matrix_h
#define dealii_tensor_product_matrix_h


#include <deal.II/base/array_view.h>
#include <deal.II/base/config.h>
#include <deal.II/base/thread_management.h>
#include <deal.II/lac/lapack_full_matrix.h>
#include <deal.II/matrix_free/tensor_product_kernels.h>

DEAL_II_NAMESPACE_OPEN

template <typename> class Vector;
template <typename> class FullMatrix;
template <typename> class VectorizedArray;

template <int dim, typename Number, int size = -1>
class TensorProductMatrixSymmetricSumBase
{
public:
  unsigned int m () const;

  unsigned int n () const;

  void vmult (const ArrayView<Number> &dst,
              const ArrayView<const Number> &src) const;

  void apply_inverse (const ArrayView<Number> &dst,
                      const ArrayView<const Number> &src) const;

protected:
  TensorProductMatrixSymmetricSumBase () = default;

  std::array<Table<2,Number>,dim> mass_matrix;

  std::array<Table<2,Number>,dim> derivative_matrix;

  std::array<AlignedVector<Number>,dim> eigenvalues;

  std::array<Table<2,Number>,dim> eigenvectors;

private:
  mutable AlignedVector<Number> tmp_array;

  mutable Threads::Mutex mutex;
};



template <int dim, typename Number, int size = -1>
class TensorProductMatrixSymmetricSum
  : public TensorProductMatrixSymmetricSumBase<dim,Number,size>
{
public:
  TensorProductMatrixSymmetricSum () = default;

  TensorProductMatrixSymmetricSum (const std::array<Table<2,Number>,dim> &mass_matrix,
                                   const std::array<Table<2,Number>,dim> &derivative_matrix);

  TensorProductMatrixSymmetricSum (const std::array<FullMatrix<Number>,dim> &mass_matrix,
                                   const std::array<FullMatrix<Number>,dim> &derivative_matrix);

  TensorProductMatrixSymmetricSum (const Table<2,Number> &mass_matrix,
                                   const Table<2,Number> &derivative_matrix);

  void reinit (const std::array<Table<2,Number>,dim> &mass_matrix,
               const std::array<Table<2,Number>,dim> &derivative_matrix);

  void reinit (const std::array<FullMatrix<Number>,dim> &mass_matrix,
               const std::array<FullMatrix<Number>,dim> &derivative_matrix);

  void reinit (const Table<2,Number> &mass_matrix,
               const Table<2,Number> &derivative_matrix);

private:
  template <typename MatrixArray>
  void reinit_impl (MatrixArray &&mass_matrix,
                    MatrixArray &&derivative_matrix);
};



template <int dim, typename Number, int size>
class TensorProductMatrixSymmetricSum<dim,VectorizedArray<Number>,size>
  : public TensorProductMatrixSymmetricSumBase<dim,VectorizedArray<Number>,size>
{
public:
  TensorProductMatrixSymmetricSum () = default;

  TensorProductMatrixSymmetricSum (const std::array<Table<2,VectorizedArray<Number> >,dim> &mass_matrix,
                                   const std::array<Table<2,VectorizedArray<Number> >,dim> &derivative_matrix);

  TensorProductMatrixSymmetricSum (const Table<2,VectorizedArray<Number> > &mass_matrix,
                                   const Table<2,VectorizedArray<Number> > &derivative_matrix);

  void reinit (const std::array<Table<2,VectorizedArray<Number> >,dim> &mass_matrix,
               const std::array<Table<2,VectorizedArray<Number> >,dim> &derivative_matrix);

  void reinit (const Table<2,VectorizedArray<Number> > &mass_matrix,
               const Table<2,VectorizedArray<Number> > &derivative_matrix);

private:
  template <typename MatrixArray>
  void reinit_impl (MatrixArray &&mass_matrix,
                    MatrixArray &&derivative_matrix);
};


/*----------------------- Inline functions ----------------------------------*/

#ifndef DOXYGEN

namespace internal
{
  namespace TensorProductMatrix
  {
    template <typename Number>
    void
    spectral_assembly (const Number *mass_matrix,
                       const Number *derivative_matrix,
                       const unsigned int n_rows,
                       const unsigned int n_cols,
                       Number *eigenvalues,
                       Number *eigenvectors)
    {
      Assert (n_rows == n_cols, ExcNotImplemented());

      auto &&transpose_fill_nm
        = [](Number *out, const Number *in, const unsigned int n, const unsigned int m)
      {
        for (unsigned int mm = 0; mm < m; ++mm)
          for (unsigned int nn = 0; nn < n; ++nn)
            out[mm+nn*m] = *(in++);
      };

      std::vector<::Vector<Number> > eigenvecs(n_rows);
      LAPACKFullMatrix<Number> mass_copy(n_rows, n_cols);
      LAPACKFullMatrix<Number> deriv_copy(n_rows, n_cols);

      transpose_fill_nm (&(mass_copy(0,0)), mass_matrix, n_rows, n_cols);
      transpose_fill_nm (&(deriv_copy(0,0)), derivative_matrix, n_rows, n_cols);

      deriv_copy.compute_generalized_eigenvalues_symmetric (mass_copy, eigenvecs);
      AssertDimension (eigenvecs.size(), n_rows);
      for (unsigned int i=0; i<n_rows; ++i)
        for (unsigned int j=0; j<n_cols; ++j, ++eigenvectors)
          *eigenvectors = eigenvecs[j][i];

      for (unsigned int i=0; i<n_rows; ++i, ++eigenvalues)
        *eigenvalues = deriv_copy.eigenvalue(i).real();
    }
  }
}


template <int dim, typename Number, int size>
inline
unsigned int
TensorProductMatrixSymmetricSumBase<dim,Number,size>::m() const
{
  unsigned int m = mass_matrix[0].n_rows();
  for (unsigned int d = 1; d < dim; ++d)
    m *= mass_matrix[d].n_rows();
  return m;
}



template <int dim, typename Number, int size>
inline
unsigned int
TensorProductMatrixSymmetricSumBase<dim,Number,size>::n() const
{
  unsigned int n = mass_matrix[0].n_cols();
  for (unsigned int d = 1; d < dim; ++d)
    n *= mass_matrix[d].n_cols();
  return n;
}



template <int dim, typename Number, int size>
inline
void
TensorProductMatrixSymmetricSumBase<dim,Number,size>
::vmult (const ArrayView<Number> &dst_view,
         const ArrayView<const Number> &src_view) const
{
  AssertDimension (dst_view.size(), this->m());
  AssertDimension (src_view.size(), this->n());
  Threads::Mutex::ScopedLock lock(this->mutex);
  const unsigned int n = Utilities::fixed_power<dim>(size > 0 ? size : eigenvalues[0].size());
  tmp_array.resize_fast(n*2);
  constexpr int kernel_size = size > 0 ? size : 0;
  internal::EvaluatorTensorProduct<internal::evaluate_general,dim,kernel_size,kernel_size,Number>
  eval(AlignedVector<Number> {}, AlignedVector<Number> {},
       AlignedVector<Number> {}, mass_matrix[0].n_rows(), mass_matrix[0].n_rows());
  Number *t = tmp_array.begin();
  const Number *src = src_view.begin();
  Number *dst = &(dst_view[0]);

  if (dim == 1)
    {
      const Number *A = &derivative_matrix[0](0,0);
      eval.template apply<0, false, false> (A, src, dst);
    }

  else if (dim == 2)
    {
      const Number *A0 = &derivative_matrix[0](0,0);
      const Number *M0 = &mass_matrix[0](0,0);
      const Number *A1 = &derivative_matrix[1](0,0);
      const Number *M1 = &mass_matrix[1](0,0);
      eval.template apply<0, false, false> (M0, src, t);
      eval.template apply<1, false, false> (A1, t, dst);
      eval.template apply<0, false, false> (A0, src, t);
      eval.template apply<1, false, true>  (M1, t, dst);
    }

  else if (dim == 3)
    {
      const Number *A0 = &derivative_matrix[0](0,0);
      const Number *M0 = &mass_matrix[0](0,0);
      const Number *A1 = &derivative_matrix[1](0,0);
      const Number *M1 = &mass_matrix[1](0,0);
      const Number *A2 = &derivative_matrix[2](0,0);
      const Number *M2 = &mass_matrix[2](0,0);
      eval.template apply<0, false, false> (M0, src, t+n);
      eval.template apply<1, false, false> (M1, t+n, t);
      eval.template apply<2, false, false> (A2, t, dst);
      eval.template apply<1, false, false> (A1, t+n, t);
      eval.template apply<0, false, false> (A0, src, t+n);
      eval.template apply<1, false, true> (M1, t+n, t);
      eval.template apply<2, false, true> (M2, t, dst);
    }

  else
    AssertThrow(false, ExcNotImplemented());
}



template <int dim, typename Number, int size>
inline
void
TensorProductMatrixSymmetricSumBase<dim,Number,size>
::apply_inverse (const ArrayView<Number> &dst_view,
                 const ArrayView<const Number> &src_view) const
{
  AssertDimension (dst_view.size(), this->n());
  AssertDimension (src_view.size(), this->m());
  Threads::Mutex::ScopedLock lock(this->mutex);
  const unsigned int n = size > 0 ? size : eigenvalues[0].size();
  tmp_array.resize_fast (Utilities::fixed_power<dim>(n));
  constexpr int kernel_size = size > 0 ? size : 0;
  internal::EvaluatorTensorProduct<internal::evaluate_general,dim,kernel_size,kernel_size,Number>
  eval(AlignedVector<Number>(), AlignedVector<Number>(),
       AlignedVector<Number>(), mass_matrix[0].n_rows(), mass_matrix[0].n_rows());
  Number *t = tmp_array.begin();
  const Number *src = src_view.data();
  Number *dst = &(dst_view[0]);

  // NOTE: dof_to_quad has to be interpreted as 'dof to eigenvalue index'
  //       --> apply<.,true,.> (S,src,dst) calculates dst = S^T * src,
  //       --> apply<.,false,.> (S,src,dst) calculates dst = S * src,
  //       while the eigenvectors are stored column-wise in S, i.e.
  //       rows correspond to dofs whereas columns to eigenvalue indices!
  if (dim == 1)
    {
      const Number *S = &eigenvectors[0](0,0);
      eval.template apply<0, true, false> (S, src, t);
      for (unsigned int i=0; i<n; ++i)
        t[i] /= eigenvalues[0][i];
      eval.template apply<0, false, false> (S, t, dst);
    }

  else if (dim == 2)
    {
      const Number *S0 = &(eigenvectors[0](0,0));
      const Number *S1 = &(eigenvectors[1](0,0));
      eval.template apply<0, true, false> (S0, src, t);
      eval.template apply<1, true, false> (S1, t, dst);
      for (unsigned int i1=0, c=0; i1<n; ++i1)
        for (unsigned int i0=0; i0<n; ++i0, ++c)
          dst[c] /= (eigenvalues[1][i1] + eigenvalues[0][i0]);
      eval.template apply<0, false, false> (S0, dst, t);
      eval.template apply<1, false, false> (S1, t, dst);
    }

  else if (dim == 3)
    {
      const Number *S0 = &eigenvectors[0](0,0);
      const Number *S1 = &eigenvectors[1](0,0);
      const Number *S2 = &eigenvectors[2](0,0);
      eval.template apply<0, true, false> (S0, src, t);
      eval.template apply<1, true, false> (S1, t, dst);
      eval.template apply<2, true, false> (S2, dst, t);
      for (unsigned int i2=0, c=0; i2<n; ++i2)
        for (unsigned int i1=0; i1<n; ++i1)
          for (unsigned int i0=0; i0<n; ++i0, ++c)
            t[c] /= (eigenvalues[2][i2] + eigenvalues[1][i1] + eigenvalues[0][i0]);
      eval.template apply<0, false, false> (S0, t, dst);
      eval.template apply<1, false, false> (S1, dst, t);
      eval.template apply<2, false, false> (S2, t, dst);
    }

  else
    Assert(false, ExcNotImplemented());
}


// ------------------------------   TensorProductMatrixSymmetricSum   ------------------------------

template <int dim, typename Number, int size>
inline
TensorProductMatrixSymmetricSum<dim,Number,size>
::TensorProductMatrixSymmetricSum (const std::array<Table<2,Number>, dim> &mass_matrix,
                                   const std::array<Table<2,Number>, dim> &derivative_matrix)
{
  reinit (mass_matrix, derivative_matrix);
}



template <int dim, typename Number, int size>
inline
TensorProductMatrixSymmetricSum<dim,Number,size>
::TensorProductMatrixSymmetricSum(const std::array<FullMatrix<Number>, dim> &mass_matrix,
                                  const std::array<FullMatrix<Number>, dim> &derivative_matrix)
{
  reinit (mass_matrix, derivative_matrix);
}



template <int dim, typename Number, int size>
inline
TensorProductMatrixSymmetricSum<dim,Number,size>
::TensorProductMatrixSymmetricSum (const Table<2,Number> &mass_matrix,
                                   const Table<2,Number> &derivative_matrix)
{
  reinit (mass_matrix, derivative_matrix);
}



template <int dim, typename Number, int size>
template <typename MatrixArray>
inline
void
TensorProductMatrixSymmetricSum<dim,Number,size>
::reinit_impl (MatrixArray &&mass_matrices_,
               MatrixArray &&derivative_matrices_)
{
  auto &&mass_matrices = std::forward<MatrixArray>(mass_matrices_);
  auto &&derivative_matrices = std::forward<MatrixArray>(derivative_matrices_);
  this->mass_matrix = mass_matrices;
  this->derivative_matrix = derivative_matrices;

  for (int dir = 0; dir < dim; ++dir)
    {
      Assert (size == -1 || (size > 0 && static_cast<unsigned int>(size) == mass_matrices[dir].n_rows()),
              ExcDimensionMismatch(size, mass_matrices[dir].n_rows()));
      AssertDimension (mass_matrices[dir].n_rows(), mass_matrices[dir].n_cols());
      AssertDimension (mass_matrices[dir].n_rows(), derivative_matrices[dir].n_rows());
      AssertDimension (mass_matrices[dir].n_rows(), derivative_matrices[dir].n_cols());

      this->eigenvectors[dir].reinit (mass_matrices[dir].n_cols(), mass_matrices[dir].n_rows());
      this->eigenvalues[dir].resize (mass_matrices[dir].n_cols());
      internal::TensorProductMatrix
      ::spectral_assembly<Number> (&(mass_matrices[dir](0,0)),
                                   &(derivative_matrices[dir](0,0)),
                                   mass_matrices[dir].n_rows(),
                                   mass_matrices[dir].n_cols(),
                                   this->eigenvalues[dir].begin(),
                                   &(this->eigenvectors[dir](0,0)));
    }
}



template <int dim, typename Number, int size>
inline
void
TensorProductMatrixSymmetricSum<dim,Number,size>
::reinit (const std::array<Table<2,Number>, dim> &mass_matrix,
          const std::array<Table<2,Number>, dim> &derivative_matrix)
{
  reinit_impl (mass_matrix, derivative_matrix);
}



template <int dim, typename Number, int size>
inline
void
TensorProductMatrixSymmetricSum<dim,Number,size>
::reinit (const std::array<FullMatrix<Number>, dim> &mass_matrix,
          const std::array<FullMatrix<Number>, dim> &derivative_matrix)
{
  std::array<Table<2,Number>,dim> mass_copy;
  std::array<Table<2,Number>,dim> deriv_copy;

  std::transform (mass_matrix.cbegin(), mass_matrix.cend(), mass_copy.begin(),
                  [] (const FullMatrix<Number> &m) ->Table<2,Number> {return m;});
  std::transform (derivative_matrix.cbegin(), derivative_matrix.cend(), deriv_copy.begin(),
                  [] (const FullMatrix<Number> &m) ->Table<2,Number> {return m;});

  reinit_impl (std::move(mass_copy), std::move(deriv_copy));
}



template <int dim, typename Number, int size>
inline
void
TensorProductMatrixSymmetricSum<dim,Number,size>
::reinit (const Table<2,Number> &mass_matrix,
          const Table<2,Number> &derivative_matrix)
{
  std::array<Table<2,Number>,dim> mass_matrices;
  std::array<Table<2,Number>,dim> derivative_matrices;

  std::fill (mass_matrices.begin(), mass_matrices.end(), mass_matrix);
  std::fill (derivative_matrices.begin(), derivative_matrices.end(), derivative_matrix);

  reinit_impl (std::move(mass_matrices), std::move(derivative_matrices));
}



// ------------------------------ vectorized spec.: TensorProductMatrixSymmetricSum   ------------------------------

template <int dim, typename Number, int size>
inline
TensorProductMatrixSymmetricSum<dim,VectorizedArray<Number>,size>
::TensorProductMatrixSymmetricSum (const std::array<Table<2,VectorizedArray<Number> >,dim> &mass_matrix,
                                   const std::array<Table<2,VectorizedArray<Number> >,dim> &derivative_matrix)
{
  reinit (mass_matrix, derivative_matrix);
}



template <int dim, typename Number, int size>
inline
TensorProductMatrixSymmetricSum<dim,VectorizedArray<Number>,size>
::TensorProductMatrixSymmetricSum (const Table<2,VectorizedArray<Number> > &mass_matrix,
                                   const Table<2,VectorizedArray<Number> > &derivative_matrix)
{
  reinit (mass_matrix, derivative_matrix);
}



template <int dim, typename Number, int size>
template <typename MatrixArray>
inline
void
TensorProductMatrixSymmetricSum<dim,VectorizedArray<Number>,size>
::reinit_impl (MatrixArray &&mass_matrices_,
               MatrixArray &&derivative_matrices_)
{
  auto &&mass_matrix = std::forward<MatrixArray>(mass_matrices_);
  auto &&derivative_matrix = std::forward<MatrixArray>(derivative_matrices_);
  this->mass_matrix = mass_matrix;
  this->derivative_matrix = derivative_matrix;

  constexpr unsigned int macro_size = VectorizedArray<Number>::n_array_elements;
  std::size_t n_rows_max = (size > 0) ? size : 0 ;
  if (size == -1)
    for (unsigned int d = 0; d < dim; ++d)
      n_rows_max = std::max(n_rows_max, mass_matrix[d].n_rows());
  const std::size_t nm_flat_size_max = n_rows_max * n_rows_max * macro_size;
  const std::size_t n_flat_size_max  = n_rows_max * macro_size;

  std::vector<Number> mass_matrix_flat;
  std::vector<Number> deriv_matrix_flat;
  std::vector<Number> eigenvalues_flat;
  std::vector<Number> eigenvectors_flat;
  mass_matrix_flat.resize (nm_flat_size_max);
  deriv_matrix_flat.resize (nm_flat_size_max);
  eigenvalues_flat.resize (n_flat_size_max);
  eigenvectors_flat.resize (nm_flat_size_max);
  std::array<unsigned int,macro_size> offsets_nm;
  std::array<unsigned int,macro_size> offsets_n;
  for (int dir = 0; dir < dim; ++dir)
    {
      Assert (size == -1 ||
              (size > 0 && static_cast<unsigned int>(size) == mass_matrix[dir].n_rows()),
              ExcDimensionMismatch(size, mass_matrix[dir].n_rows()));
      AssertDimension (mass_matrix[dir].n_rows(), mass_matrix[dir].n_cols());
      AssertDimension (mass_matrix[dir].n_rows(), derivative_matrix[dir].n_rows());
      AssertDimension (mass_matrix[dir].n_rows(), derivative_matrix[dir].n_cols());

      const unsigned int n_rows = mass_matrix[dir].n_rows();
      const unsigned int n_cols = mass_matrix[dir].n_cols();
      const unsigned int nm = n_rows * n_cols;
      for (unsigned int vv=0; vv<macro_size; ++vv)
        offsets_nm[vv] = nm * vv;

      vectorized_transpose_and_store (false, nm, &(mass_matrix[dir](0,0)),
                                      offsets_nm.cbegin(), mass_matrix_flat.data());
      vectorized_transpose_and_store (false, nm, &(derivative_matrix[dir](0,0)),
                                      offsets_nm.cbegin(), deriv_matrix_flat.data());

      const Number *mass_cbegin = mass_matrix_flat.data();
      const Number *deriv_cbegin = deriv_matrix_flat.data();
      Number *eigenvec_begin = eigenvectors_flat.data();
      Number *eigenval_begin = eigenvalues_flat.data();
      for (unsigned int lane = 0; lane < macro_size; ++lane)
        internal::TensorProductMatrix
        ::spectral_assembly<Number> (mass_cbegin+nm*lane, deriv_cbegin+nm*lane, n_rows, n_cols,
                                     eigenval_begin+n_rows*lane, eigenvec_begin+nm*lane);

      this->eigenvalues[dir].resize (n_rows);
      this->eigenvectors[dir].reinit (n_rows, n_cols);
      for (unsigned int vv=0; vv<macro_size; ++vv)
        offsets_n[vv] = n_rows * vv;
      vectorized_load_and_transpose (n_rows, eigenvalues_flat.data(),
                                     offsets_n.cbegin(), this->eigenvalues[dir].begin());
      vectorized_load_and_transpose (nm, eigenvectors_flat.data(),
                                     offsets_nm.cbegin(), &(this->eigenvectors[dir](0,0)));
    }
}



template <int dim, typename Number, int size>
inline
void
TensorProductMatrixSymmetricSum<dim,VectorizedArray<Number>,size>
::reinit (const std::array<Table<2,VectorizedArray<Number> >,dim> &mass_matrix,
          const std::array<Table<2,VectorizedArray<Number> >,dim> &derivative_matrix)
{
  reinit_impl (mass_matrix, derivative_matrix);
}



template <int dim, typename Number, int size>
inline
void
TensorProductMatrixSymmetricSum<dim,VectorizedArray<Number>,size>
::reinit (const Table<2,VectorizedArray<Number> > &mass_matrix,
          const Table<2,VectorizedArray<Number> > &derivative_matrix)
{
  std::array<Table<2,VectorizedArray<Number> >,dim> mass_matrices;
  std::array<Table<2,VectorizedArray<Number> >,dim> derivative_matrices;

  std::fill (mass_matrices.begin(), mass_matrices.end(), mass_matrix);
  std::fill (derivative_matrices.begin(), derivative_matrices.end(), derivative_matrix);

  reinit_impl (std::move(mass_matrices), std::move(derivative_matrices));
}



#endif

DEAL_II_NAMESPACE_CLOSE

#endif