Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
tensor_product_matrix.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2017 - 2019 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_tensor_product_matrix_h
17 #define dealii_tensor_product_matrix_h
18 
19 
20 #include <deal.II/base/config.h>
21 
22 #include <deal.II/base/array_view.h>
23 #include <deal.II/base/thread_management.h>
24 
25 #include <deal.II/lac/lapack_full_matrix.h>
26 
27 #include <deal.II/matrix_free/tensor_product_kernels.h>
28 
29 DEAL_II_NAMESPACE_OPEN
30 
31 template <typename>
32 class Vector;
33 template <typename>
34 class FullMatrix;
35 template <typename>
36 class VectorizedArray;
37 
72 template <int dim, typename Number, int size = -1>
74 {
75 public:
81  unsigned int
82  m() const;
83 
89  unsigned int
90  n() const;
91 
98  void
99  vmult(const ArrayView<Number> &dst, const ArrayView<const Number> &src) const;
100 
107  void
108  apply_inverse(const ArrayView<Number> & dst,
109  const ArrayView<const Number> &src) const;
110 
111 protected:
116 
120  std::array<Table<2, Number>, dim> mass_matrix;
121 
125  std::array<Table<2, Number>, dim> derivative_matrix;
126 
131  std::array<AlignedVector<Number>, dim> eigenvalues;
132 
137  std::array<Table<2, Number>, dim> eigenvectors;
138 
139 private:
144 
149 };
150 
151 
152 
224 template <int dim, typename Number, int size = -1>
226  : public TensorProductMatrixSymmetricSumBase<dim, Number, size>
227 {
228 public:
233 
241  const std::array<Table<2, Number>, dim> &mass_matrix,
242  const std::array<Table<2, Number>, dim> &derivative_matrix);
243 
251  const std::array<FullMatrix<Number>, dim> &mass_matrix,
252  const std::array<FullMatrix<Number>, dim> &derivative_matrix);
253 
260 
272  void
273  reinit(const std::array<Table<2, Number>, dim> &mass_matrix,
274  const std::array<Table<2, Number>, dim> &derivative_matrix);
275 
281  void
282  reinit(const std::array<FullMatrix<Number>, dim> &mass_matrix,
283  const std::array<FullMatrix<Number>, dim> &derivative_matrix);
284 
290  void
293 
294 private:
303  template <typename MatrixArray>
304  void
305  reinit_impl(MatrixArray &&mass_matrix, MatrixArray &&derivative_matrix);
306 };
307 
308 
309 
318 template <int dim, typename Number, int size>
321  VectorizedArray<Number>,
322  size>
323 {
324 public:
329 
337  const std::array<Table<2, VectorizedArray<Number>>, dim> &mass_matrix,
338  const std::array<Table<2, VectorizedArray<Number>>, dim>
340 
350 
362  void
363  reinit(const std::array<Table<2, VectorizedArray<Number>>, dim> &mass_matrix,
364  const std::array<Table<2, VectorizedArray<Number>>, dim>
366 
372  void
375 
376 private:
385  template <typename MatrixArray>
386  void
387  reinit_impl(MatrixArray &&mass_matrix, MatrixArray &&derivative_matrix);
388 };
389 
390 
391 /*----------------------- Inline functions ----------------------------------*/
392 
393 #ifndef DOXYGEN
394 
395 namespace internal
396 {
397  namespace TensorProductMatrix
398  {
407  template <typename Number>
408  void
409  spectral_assembly(const Number * mass_matrix,
410  const Number * derivative_matrix,
411  const unsigned int n_rows,
412  const unsigned int n_cols,
413  Number * eigenvalues,
414  Number * eigenvectors)
415  {
416  Assert(n_rows == n_cols, ExcNotImplemented());
417 
418  auto &&transpose_fill_nm = [](Number * out,
419  const Number * in,
420  const unsigned int n,
421  const unsigned int m) {
422  for (unsigned int mm = 0; mm < m; ++mm)
423  for (unsigned int nn = 0; nn < n; ++nn)
424  out[mm + nn * m] = *(in++);
425  };
426 
427  std::vector<::Vector<Number>> eigenvecs(n_rows);
428  LAPACKFullMatrix<Number> mass_copy(n_rows, n_cols);
429  LAPACKFullMatrix<Number> deriv_copy(n_rows, n_cols);
430 
431  transpose_fill_nm(&(mass_copy(0, 0)), mass_matrix, n_rows, n_cols);
432  transpose_fill_nm(&(deriv_copy(0, 0)), derivative_matrix, n_rows, n_cols);
433 
434  deriv_copy.compute_generalized_eigenvalues_symmetric(mass_copy,
435  eigenvecs);
436  AssertDimension(eigenvecs.size(), n_rows);
437  for (unsigned int i = 0; i < n_rows; ++i)
438  for (unsigned int j = 0; j < n_cols; ++j, ++eigenvectors)
439  *eigenvectors = eigenvecs[j][i];
440 
441  for (unsigned int i = 0; i < n_rows; ++i, ++eigenvalues)
442  *eigenvalues = deriv_copy.eigenvalue(i).real();
443  }
444  } // namespace TensorProductMatrix
445 } // namespace internal
446 
447 
448 template <int dim, typename Number, int size>
449 inline unsigned int
451 {
452  unsigned int m = mass_matrix[0].n_rows();
453  for (unsigned int d = 1; d < dim; ++d)
454  m *= mass_matrix[d].n_rows();
455  return m;
456 }
457 
458 
459 
460 template <int dim, typename Number, int size>
461 inline unsigned int
463 {
464  unsigned int n = mass_matrix[0].n_cols();
465  for (unsigned int d = 1; d < dim; ++d)
466  n *= mass_matrix[d].n_cols();
467  return n;
468 }
469 
470 
471 
472 template <int dim, typename Number, int size>
473 inline void
475  const ArrayView<Number> & dst_view,
476  const ArrayView<const Number> &src_view) const
477 {
478  AssertDimension(dst_view.size(), this->m());
479  AssertDimension(src_view.size(), this->n());
480  std::lock_guard<std::mutex> lock(this->mutex);
481  const unsigned int n =
482  Utilities::fixed_power<dim>(size > 0 ? size : eigenvalues[0].size());
483  tmp_array.resize_fast(n * 2);
484  constexpr int kernel_size = size > 0 ? size : 0;
486  dim,
487  kernel_size,
488  kernel_size,
489  Number>
490  eval(AlignedVector<Number>{},
493  mass_matrix[0].n_rows(),
494  mass_matrix[0].n_rows());
495  Number * t = tmp_array.begin();
496  const Number *src = src_view.begin();
497  Number * dst = dst_view.data();
498 
499  if (dim == 1)
500  {
501  const Number *A = &derivative_matrix[0](0, 0);
502  eval.template apply<0, false, false>(A, src, dst);
503  }
504 
505  else if (dim == 2)
506  {
507  const Number *A0 = &derivative_matrix[0](0, 0);
508  const Number *M0 = &mass_matrix[0](0, 0);
509  const Number *A1 = &derivative_matrix[1](0, 0);
510  const Number *M1 = &mass_matrix[1](0, 0);
511  eval.template apply<0, false, false>(M0, src, t);
512  eval.template apply<1, false, false>(A1, t, dst);
513  eval.template apply<0, false, false>(A0, src, t);
514  eval.template apply<1, false, true>(M1, t, dst);
515  }
516 
517  else if (dim == 3)
518  {
519  const Number *A0 = &derivative_matrix[0](0, 0);
520  const Number *M0 = &mass_matrix[0](0, 0);
521  const Number *A1 = &derivative_matrix[1](0, 0);
522  const Number *M1 = &mass_matrix[1](0, 0);
523  const Number *A2 = &derivative_matrix[2](0, 0);
524  const Number *M2 = &mass_matrix[2](0, 0);
525  eval.template apply<0, false, false>(M0, src, t + n);
526  eval.template apply<1, false, false>(M1, t + n, t);
527  eval.template apply<2, false, false>(A2, t, dst);
528  eval.template apply<1, false, false>(A1, t + n, t);
529  eval.template apply<0, false, false>(A0, src, t + n);
530  eval.template apply<1, false, true>(M1, t + n, t);
531  eval.template apply<2, false, true>(M2, t, dst);
532  }
533 
534  else
535  AssertThrow(false, ExcNotImplemented());
536 }
537 
538 
539 
540 template <int dim, typename Number, int size>
541 inline void
543  const ArrayView<Number> & dst_view,
544  const ArrayView<const Number> &src_view) const
545 {
546  AssertDimension(dst_view.size(), this->n());
547  AssertDimension(src_view.size(), this->m());
548  std::lock_guard<std::mutex> lock(this->mutex);
549  const unsigned int n = size > 0 ? size : eigenvalues[0].size();
550  tmp_array.resize_fast(Utilities::fixed_power<dim>(n));
551  constexpr int kernel_size = size > 0 ? size : 0;
553  dim,
554  kernel_size,
555  kernel_size,
556  Number>
557  eval(AlignedVector<Number>(),
560  mass_matrix[0].n_rows(),
561  mass_matrix[0].n_rows());
562  Number * t = tmp_array.begin();
563  const Number *src = src_view.data();
564  Number * dst = dst_view.data();
565 
566  // NOTE: dof_to_quad has to be interpreted as 'dof to eigenvalue index'
567  // --> apply<.,true,.> (S,src,dst) calculates dst = S^T * src,
568  // --> apply<.,false,.> (S,src,dst) calculates dst = S * src,
569  // while the eigenvectors are stored column-wise in S, i.e.
570  // rows correspond to dofs whereas columns to eigenvalue indices!
571  if (dim == 1)
572  {
573  const Number *S = &eigenvectors[0](0, 0);
574  eval.template apply<0, true, false>(S, src, t);
575  for (unsigned int i = 0; i < n; ++i)
576  t[i] /= eigenvalues[0][i];
577  eval.template apply<0, false, false>(S, t, dst);
578  }
579 
580  else if (dim == 2)
581  {
582  const Number *S0 = &(eigenvectors[0](0, 0));
583  const Number *S1 = &(eigenvectors[1](0, 0));
584  eval.template apply<0, true, false>(S0, src, t);
585  eval.template apply<1, true, false>(S1, t, dst);
586  for (unsigned int i1 = 0, c = 0; i1 < n; ++i1)
587  for (unsigned int i0 = 0; i0 < n; ++i0, ++c)
588  dst[c] /= (eigenvalues[1][i1] + eigenvalues[0][i0]);
589  eval.template apply<0, false, false>(S0, dst, t);
590  eval.template apply<1, false, false>(S1, t, dst);
591  }
592 
593  else if (dim == 3)
594  {
595  const Number *S0 = &eigenvectors[0](0, 0);
596  const Number *S1 = &eigenvectors[1](0, 0);
597  const Number *S2 = &eigenvectors[2](0, 0);
598  eval.template apply<0, true, false>(S0, src, t);
599  eval.template apply<1, true, false>(S1, t, dst);
600  eval.template apply<2, true, false>(S2, dst, t);
601  for (unsigned int i2 = 0, c = 0; i2 < n; ++i2)
602  for (unsigned int i1 = 0; i1 < n; ++i1)
603  for (unsigned int i0 = 0; i0 < n; ++i0, ++c)
604  t[c] /=
605  (eigenvalues[2][i2] + eigenvalues[1][i1] + eigenvalues[0][i0]);
606  eval.template apply<0, false, false>(S0, t, dst);
607  eval.template apply<1, false, false>(S1, dst, t);
608  eval.template apply<2, false, false>(S2, t, dst);
609  }
610 
611  else
612  Assert(false, ExcNotImplemented());
613 }
614 
615 
616 //---------------------- TensorProductMatrixSymmetricSum ----------------------
617 
618 template <int dim, typename Number, int size>
621  const std::array<Table<2, Number>, dim> &mass_matrix,
622  const std::array<Table<2, Number>, dim> &derivative_matrix)
623 {
624  reinit(mass_matrix, derivative_matrix);
625 }
626 
627 
628 
629 template <int dim, typename Number, int size>
632  const std::array<FullMatrix<Number>, dim> &mass_matrix,
633  const std::array<FullMatrix<Number>, dim> &derivative_matrix)
634 {
635  reinit(mass_matrix, derivative_matrix);
636 }
637 
638 
639 
640 template <int dim, typename Number, int size>
643  const Table<2, Number> &derivative_matrix)
644 {
645  reinit(mass_matrix, derivative_matrix);
646 }
647 
648 
649 
650 template <int dim, typename Number, int size>
651 template <typename MatrixArray>
652 inline void
654  MatrixArray &&mass_matrices_,
655  MatrixArray &&derivative_matrices_)
656 {
657  auto &&mass_matrices = std::forward<MatrixArray>(mass_matrices_);
658  auto &&derivative_matrices = std::forward<MatrixArray>(derivative_matrices_);
659  this->mass_matrix = mass_matrices;
660  this->derivative_matrix = derivative_matrices;
661 
662  for (int dir = 0; dir < dim; ++dir)
663  {
664  Assert(size == -1 || (size > 0 && static_cast<unsigned int>(size) ==
665  mass_matrices[dir].n_rows()),
666  ExcDimensionMismatch(size, mass_matrices[dir].n_rows()));
667  AssertDimension(mass_matrices[dir].n_rows(), mass_matrices[dir].n_cols());
668  AssertDimension(mass_matrices[dir].n_rows(),
669  derivative_matrices[dir].n_rows());
670  AssertDimension(mass_matrices[dir].n_rows(),
671  derivative_matrices[dir].n_cols());
672 
673  this->eigenvectors[dir].reinit(mass_matrices[dir].n_cols(),
674  mass_matrices[dir].n_rows());
675  this->eigenvalues[dir].resize(mass_matrices[dir].n_cols());
676  internal::TensorProductMatrix ::spectral_assembly<Number>(
677  &(mass_matrices[dir](0, 0)),
678  &(derivative_matrices[dir](0, 0)),
679  mass_matrices[dir].n_rows(),
680  mass_matrices[dir].n_cols(),
681  this->eigenvalues[dir].begin(),
682  &(this->eigenvectors[dir](0, 0)));
683  }
684 }
685 
686 
687 
688 template <int dim, typename Number, int size>
689 inline void
691  const std::array<Table<2, Number>, dim> &mass_matrix,
692  const std::array<Table<2, Number>, dim> &derivative_matrix)
693 {
694  reinit_impl(mass_matrix, derivative_matrix);
695 }
696 
697 
698 
699 template <int dim, typename Number, int size>
700 inline void
702  const std::array<FullMatrix<Number>, dim> &mass_matrix,
703  const std::array<FullMatrix<Number>, dim> &derivative_matrix)
704 {
705  std::array<Table<2, Number>, dim> mass_copy;
706  std::array<Table<2, Number>, dim> deriv_copy;
707 
708  std::transform(mass_matrix.cbegin(),
709  mass_matrix.cend(),
710  mass_copy.begin(),
711  [](const FullMatrix<Number> &m) -> Table<2, Number> {
712  return m;
713  });
714  std::transform(derivative_matrix.cbegin(),
715  derivative_matrix.cend(),
716  deriv_copy.begin(),
717  [](const FullMatrix<Number> &m) -> Table<2, Number> {
718  return m;
719  });
720 
721  reinit_impl(std::move(mass_copy), std::move(deriv_copy));
722 }
723 
724 
725 
726 template <int dim, typename Number, int size>
727 inline void
729  const Table<2, Number> &mass_matrix,
730  const Table<2, Number> &derivative_matrix)
731 {
732  std::array<Table<2, Number>, dim> mass_matrices;
733  std::array<Table<2, Number>, dim> derivative_matrices;
734 
735  std::fill(mass_matrices.begin(), mass_matrices.end(), mass_matrix);
736  std::fill(derivative_matrices.begin(),
737  derivative_matrices.end(),
738  derivative_matrix);
739 
740  reinit_impl(std::move(mass_matrices), std::move(derivative_matrices));
741 }
742 
743 
744 
745 //------------- vectorized spec.: TensorProductMatrixSymmetricSum -------------
746 
747 template <int dim, typename Number, int size>
749  TensorProductMatrixSymmetricSum(
750  const std::array<Table<2, VectorizedArray<Number>>, dim> &mass_matrix,
751  const std::array<Table<2, VectorizedArray<Number>>, dim> &derivative_matrix)
752 {
753  reinit(mass_matrix, derivative_matrix);
754 }
755 
756 
757 
758 template <int dim, typename Number, int size>
760  TensorProductMatrixSymmetricSum(
761  const Table<2, VectorizedArray<Number>> &mass_matrix,
762  const Table<2, VectorizedArray<Number>> &derivative_matrix)
763 {
764  reinit(mass_matrix, derivative_matrix);
765 }
766 
767 
768 
769 template <int dim, typename Number, int size>
770 template <typename MatrixArray>
771 inline void
773  reinit_impl(MatrixArray &&mass_matrices_, MatrixArray &&derivative_matrices_)
774 {
775  auto &&mass_matrix = std::forward<MatrixArray>(mass_matrices_);
776  auto &&derivative_matrix = std::forward<MatrixArray>(derivative_matrices_);
777  this->mass_matrix = mass_matrix;
778  this->derivative_matrix = derivative_matrix;
779 
780  constexpr unsigned int macro_size = VectorizedArray<Number>::n_array_elements;
781  std::size_t n_rows_max = (size > 0) ? size : 0;
782  if (size == -1)
783  for (unsigned int d = 0; d < dim; ++d)
784  n_rows_max = std::max(n_rows_max, mass_matrix[d].n_rows());
785  const std::size_t nm_flat_size_max = n_rows_max * n_rows_max * macro_size;
786  const std::size_t n_flat_size_max = n_rows_max * macro_size;
787 
788  std::vector<Number> mass_matrix_flat;
789  std::vector<Number> deriv_matrix_flat;
790  std::vector<Number> eigenvalues_flat;
791  std::vector<Number> eigenvectors_flat;
792  mass_matrix_flat.resize(nm_flat_size_max);
793  deriv_matrix_flat.resize(nm_flat_size_max);
794  eigenvalues_flat.resize(n_flat_size_max);
795  eigenvectors_flat.resize(nm_flat_size_max);
796  std::array<unsigned int, macro_size> offsets_nm;
797  std::array<unsigned int, macro_size> offsets_n;
798  for (int dir = 0; dir < dim; ++dir)
799  {
800  Assert(size == -1 || (size > 0 && static_cast<unsigned int>(size) ==
801  mass_matrix[dir].n_rows()),
802  ExcDimensionMismatch(size, mass_matrix[dir].n_rows()));
803  AssertDimension(mass_matrix[dir].n_rows(), mass_matrix[dir].n_cols());
804  AssertDimension(mass_matrix[dir].n_rows(),
805  derivative_matrix[dir].n_rows());
806  AssertDimension(mass_matrix[dir].n_rows(),
807  derivative_matrix[dir].n_cols());
808 
809  const unsigned int n_rows = mass_matrix[dir].n_rows();
810  const unsigned int n_cols = mass_matrix[dir].n_cols();
811  const unsigned int nm = n_rows * n_cols;
812  for (unsigned int vv = 0; vv < macro_size; ++vv)
813  offsets_nm[vv] = nm * vv;
814 
816  nm,
817  &(mass_matrix[dir](0, 0)),
818  offsets_nm.cbegin(),
819  mass_matrix_flat.data());
821  nm,
822  &(derivative_matrix[dir](0, 0)),
823  offsets_nm.cbegin(),
824  deriv_matrix_flat.data());
825 
826  const Number *mass_cbegin = mass_matrix_flat.data();
827  const Number *deriv_cbegin = deriv_matrix_flat.data();
828  Number * eigenvec_begin = eigenvectors_flat.data();
829  Number * eigenval_begin = eigenvalues_flat.data();
830  for (unsigned int lane = 0; lane < macro_size; ++lane)
831  internal::TensorProductMatrix ::spectral_assembly<Number>(
832  mass_cbegin + nm * lane,
833  deriv_cbegin + nm * lane,
834  n_rows,
835  n_cols,
836  eigenval_begin + n_rows * lane,
837  eigenvec_begin + nm * lane);
838 
839  this->eigenvalues[dir].resize(n_rows);
840  this->eigenvectors[dir].reinit(n_rows, n_cols);
841  for (unsigned int vv = 0; vv < macro_size; ++vv)
842  offsets_n[vv] = n_rows * vv;
844  eigenvalues_flat.data(),
845  offsets_n.cbegin(),
846  this->eigenvalues[dir].begin());
848  eigenvectors_flat.data(),
849  offsets_nm.cbegin(),
850  &(this->eigenvectors[dir](0, 0)));
851  }
852 }
853 
854 
855 
856 template <int dim, typename Number, int size>
857 inline void
859  const std::array<Table<2, VectorizedArray<Number>>, dim> &mass_matrix,
860  const std::array<Table<2, VectorizedArray<Number>>, dim> &derivative_matrix)
861 {
862  reinit_impl(mass_matrix, derivative_matrix);
863 }
864 
865 
866 
867 template <int dim, typename Number, int size>
868 inline void
870  const Table<2, VectorizedArray<Number>> &mass_matrix,
871  const Table<2, VectorizedArray<Number>> &derivative_matrix)
872 {
873  std::array<Table<2, VectorizedArray<Number>>, dim> mass_matrices;
874  std::array<Table<2, VectorizedArray<Number>>, dim> derivative_matrices;
875 
876  std::fill(mass_matrices.begin(), mass_matrices.end(), mass_matrix);
877  std::fill(derivative_matrices.begin(),
878  derivative_matrices.end(),
879  derivative_matrix);
880 
881  reinit_impl(std::move(mass_matrices), std::move(derivative_matrices));
882 }
883 
884 
885 
886 #endif
887 
888 DEAL_II_NAMESPACE_CLOSE
889 
890 #endif
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1567
std::array< Table< 2, Number >, dim > derivative_matrix
std::array< AlignedVector< Number >, dim > eigenvalues
std::array< std::pair< Number, Tensor< 1, dim, Number > >, std::integral_constant< int, dim >::value > eigenvectors(const SymmetricTensor< 2, dim, Number > &T, const SymmetricTensorEigenvectorMethod method=SymmetricTensorEigenvectorMethod::ql_implicit_shifts)
void reinit(const std::array< Table< 2, Number >, dim > &mass_matrix, const std::array< Table< 2, Number >, dim > &derivative_matrix)
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
std::array< std::pair< Number, Tensor< 1, dim, Number > >, std::integral_constant< int, dim >::value > eigenvectors(const SymmetricTensor< 2, dim, Number > &T, const SymmetricTensorEigenvectorMethod method=SymmetricTensorEigenvectorMethod::ql_implicit_shifts)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
std::array< Table< 2, Number >, dim > mass_matrix
LinearAlgebra::distributed::Vector< Number > Vector
std::size_t size() const
Definition: array_view.h:471
void reinit_impl(MatrixArray &&mass_matrix, MatrixArray &&derivative_matrix)
void vectorized_transpose_and_store(const bool add_into, const unsigned int n_entries, const VectorizedArray< Number > *in, const unsigned int *offsets, Number *out)
void mass_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const double factor=1.)
Definition: l2.h:63
#define Assert(cond, exc)
Definition: exceptions.h:1407
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
value_type * data() const noexcept
Definition: array_view.h:448
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
void vectorized_load_and_transpose(const unsigned int n_entries, const Number *in, const unsigned int *offsets, VectorizedArray< Number > *out)
std::array< Table< 2, Number >, dim > eigenvectors
void vmult(const ArrayView< Number > &dst, const ArrayView< const Number > &src) const
static ::ExceptionBase & ExcNotImplemented()
iterator begin() const
Definition: array_view.h:480
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
void apply_inverse(const ArrayView< Number > &dst, const ArrayView< const Number > &src) const