Reference documentation for deal.II version Git cb0bd54b52 2019-09-21 16:31:22 -0400
\(\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 // Forward declarations
32 #ifndef DOXYGEN
33 template <typename>
34 class Vector;
35 template <typename>
36 class FullMatrix;
37 #endif
38 
73 template <int dim, typename Number, int size = -1>
75 {
76 public:
82  unsigned int
83  m() const;
84 
90  unsigned int
91  n() const;
92 
99  void
100  vmult(const ArrayView<Number> &dst, const ArrayView<const Number> &src) const;
101 
108  void
109  apply_inverse(const ArrayView<Number> & dst,
110  const ArrayView<const Number> &src) const;
111 
112 protected:
117 
121  std::array<Table<2, Number>, dim> mass_matrix;
122 
126  std::array<Table<2, Number>, dim> derivative_matrix;
127 
132  std::array<AlignedVector<Number>, dim> eigenvalues;
133 
138  std::array<Table<2, Number>, dim> eigenvectors;
139 
140 private:
145 
150 };
151 
152 
153 
225 template <int dim, typename Number, int size = -1>
227  : public TensorProductMatrixSymmetricSumBase<dim, Number, size>
228 {
229 public:
234 
242  const std::array<Table<2, Number>, dim> &mass_matrix,
243  const std::array<Table<2, Number>, dim> &derivative_matrix);
244 
252  const std::array<FullMatrix<Number>, dim> &mass_matrix,
253  const std::array<FullMatrix<Number>, dim> &derivative_matrix);
254 
261 
273  void
274  reinit(const std::array<Table<2, Number>, dim> &mass_matrix,
275  const std::array<Table<2, Number>, dim> &derivative_matrix);
276 
282  void
283  reinit(const std::array<FullMatrix<Number>, dim> &mass_matrix,
284  const std::array<FullMatrix<Number>, dim> &derivative_matrix);
285 
291  void
292  reinit(const Table<2, Number> &mass_matrix,
293  const Table<2, Number> &derivative_matrix);
294 
295 private:
304  template <typename MatrixArray>
305  void
306  reinit_impl(MatrixArray &&mass_matrix, MatrixArray &&derivative_matrix);
307 };
308 
309 
310 
319 template <int dim, typename Number, int size>
322  VectorizedArray<Number>,
323  size>
324 {
325 public:
330 
338  const std::array<Table<2, VectorizedArray<Number>>, dim> &mass_matrix,
339  const std::array<Table<2, VectorizedArray<Number>>, dim>
341 
351 
363  void
364  reinit(const std::array<Table<2, VectorizedArray<Number>>, dim> &mass_matrix,
365  const std::array<Table<2, VectorizedArray<Number>>, dim>
367 
373  void
374  reinit(const Table<2, VectorizedArray<Number>> &mass_matrix,
376 
377 private:
386  template <typename MatrixArray>
387  void
388  reinit_impl(MatrixArray &&mass_matrix, MatrixArray &&derivative_matrix);
389 };
390 
391 
392 /*----------------------- Inline functions ----------------------------------*/
393 
394 #ifndef DOXYGEN
395 
396 namespace internal
397 {
398  namespace TensorProductMatrix
399  {
408  template <typename Number>
409  void
410  spectral_assembly(const Number * mass_matrix,
411  const Number * derivative_matrix,
412  const unsigned int n_rows,
413  const unsigned int n_cols,
414  Number * eigenvalues,
415  Number * eigenvectors)
416  {
417  Assert(n_rows == n_cols, ExcNotImplemented());
418 
419  auto &&transpose_fill_nm = [](Number * out,
420  const Number * in,
421  const unsigned int n,
422  const unsigned int m) {
423  for (unsigned int mm = 0; mm < m; ++mm)
424  for (unsigned int nn = 0; nn < n; ++nn)
425  out[mm + nn * m] = *(in++);
426  };
427 
428  std::vector<::Vector<Number>> eigenvecs(n_rows);
429  LAPACKFullMatrix<Number> mass_copy(n_rows, n_cols);
430  LAPACKFullMatrix<Number> deriv_copy(n_rows, n_cols);
431 
432  transpose_fill_nm(&(mass_copy(0, 0)), mass_matrix, n_rows, n_cols);
433  transpose_fill_nm(&(deriv_copy(0, 0)), derivative_matrix, n_rows, n_cols);
434 
435  deriv_copy.compute_generalized_eigenvalues_symmetric(mass_copy,
436  eigenvecs);
437  AssertDimension(eigenvecs.size(), n_rows);
438  for (unsigned int i = 0; i < n_rows; ++i)
439  for (unsigned int j = 0; j < n_cols; ++j, ++eigenvectors)
440  *eigenvectors = eigenvecs[j][i];
441 
442  for (unsigned int i = 0; i < n_rows; ++i, ++eigenvalues)
443  *eigenvalues = deriv_copy.eigenvalue(i).real();
444  }
445  } // namespace TensorProductMatrix
446 } // namespace internal
447 
448 
449 template <int dim, typename Number, int size>
450 inline unsigned int
452 {
453  unsigned int m = mass_matrix[0].n_rows();
454  for (unsigned int d = 1; d < dim; ++d)
455  m *= mass_matrix[d].n_rows();
456  return m;
457 }
458 
459 
460 
461 template <int dim, typename Number, int size>
462 inline unsigned int
464 {
465  unsigned int n = mass_matrix[0].n_cols();
466  for (unsigned int d = 1; d < dim; ++d)
467  n *= mass_matrix[d].n_cols();
468  return n;
469 }
470 
471 
472 
473 template <int dim, typename Number, int size>
474 inline void
476  const ArrayView<Number> & dst_view,
477  const ArrayView<const Number> &src_view) const
478 {
479  AssertDimension(dst_view.size(), this->m());
480  AssertDimension(src_view.size(), this->n());
481  std::lock_guard<std::mutex> lock(this->mutex);
482  const unsigned int n =
483  Utilities::fixed_power<dim>(size > 0 ? size : eigenvalues[0].size());
484  tmp_array.resize_fast(n * 2);
485  constexpr int kernel_size = size > 0 ? size : 0;
487  dim,
488  kernel_size,
489  kernel_size,
490  Number>
491  eval(AlignedVector<Number>{},
494  mass_matrix[0].n_rows(),
495  mass_matrix[0].n_rows());
496  Number * t = tmp_array.begin();
497  const Number *src = src_view.begin();
498  Number * dst = dst_view.data();
499 
500  if (dim == 1)
501  {
502  const Number *A = &derivative_matrix[0](0, 0);
503  eval.template apply<0, false, false>(A, src, dst);
504  }
505 
506  else if (dim == 2)
507  {
508  const Number *A0 = &derivative_matrix[0](0, 0);
509  const Number *M0 = &mass_matrix[0](0, 0);
510  const Number *A1 = &derivative_matrix[1](0, 0);
511  const Number *M1 = &mass_matrix[1](0, 0);
512  eval.template apply<0, false, false>(M0, src, t);
513  eval.template apply<1, false, false>(A1, t, dst);
514  eval.template apply<0, false, false>(A0, src, t);
515  eval.template apply<1, false, true>(M1, t, dst);
516  }
517 
518  else if (dim == 3)
519  {
520  const Number *A0 = &derivative_matrix[0](0, 0);
521  const Number *M0 = &mass_matrix[0](0, 0);
522  const Number *A1 = &derivative_matrix[1](0, 0);
523  const Number *M1 = &mass_matrix[1](0, 0);
524  const Number *A2 = &derivative_matrix[2](0, 0);
525  const Number *M2 = &mass_matrix[2](0, 0);
526  eval.template apply<0, false, false>(M0, src, t + n);
527  eval.template apply<1, false, false>(M1, t + n, t);
528  eval.template apply<2, false, false>(A2, t, dst);
529  eval.template apply<1, false, false>(A1, t + n, t);
530  eval.template apply<0, false, false>(A0, src, t + n);
531  eval.template apply<1, false, true>(M1, t + n, t);
532  eval.template apply<2, false, true>(M2, t, dst);
533  }
534 
535  else
536  AssertThrow(false, ExcNotImplemented());
537 }
538 
539 
540 
541 template <int dim, typename Number, int size>
542 inline void
544  const ArrayView<Number> & dst_view,
545  const ArrayView<const Number> &src_view) const
546 {
547  AssertDimension(dst_view.size(), this->n());
548  AssertDimension(src_view.size(), this->m());
549  std::lock_guard<std::mutex> lock(this->mutex);
550  const unsigned int n = size > 0 ? size : eigenvalues[0].size();
551  tmp_array.resize_fast(Utilities::fixed_power<dim>(n));
552  constexpr int kernel_size = size > 0 ? size : 0;
554  dim,
555  kernel_size,
556  kernel_size,
557  Number>
558  eval(AlignedVector<Number>(),
561  mass_matrix[0].n_rows(),
562  mass_matrix[0].n_rows());
563  Number * t = tmp_array.begin();
564  const Number *src = src_view.data();
565  Number * dst = dst_view.data();
566 
567  // NOTE: dof_to_quad has to be interpreted as 'dof to eigenvalue index'
568  // --> apply<.,true,.> (S,src,dst) calculates dst = S^T * src,
569  // --> apply<.,false,.> (S,src,dst) calculates dst = S * src,
570  // while the eigenvectors are stored column-wise in S, i.e.
571  // rows correspond to dofs whereas columns to eigenvalue indices!
572  if (dim == 1)
573  {
574  const Number *S = &eigenvectors[0](0, 0);
575  eval.template apply<0, true, false>(S, src, t);
576  for (unsigned int i = 0; i < n; ++i)
577  t[i] /= eigenvalues[0][i];
578  eval.template apply<0, false, false>(S, t, dst);
579  }
580 
581  else if (dim == 2)
582  {
583  const Number *S0 = &(eigenvectors[0](0, 0));
584  const Number *S1 = &(eigenvectors[1](0, 0));
585  eval.template apply<0, true, false>(S0, src, t);
586  eval.template apply<1, true, false>(S1, t, dst);
587  for (unsigned int i1 = 0, c = 0; i1 < n; ++i1)
588  for (unsigned int i0 = 0; i0 < n; ++i0, ++c)
589  dst[c] /= (eigenvalues[1][i1] + eigenvalues[0][i0]);
590  eval.template apply<0, false, false>(S0, dst, t);
591  eval.template apply<1, false, false>(S1, t, dst);
592  }
593 
594  else if (dim == 3)
595  {
596  const Number *S0 = &eigenvectors[0](0, 0);
597  const Number *S1 = &eigenvectors[1](0, 0);
598  const Number *S2 = &eigenvectors[2](0, 0);
599  eval.template apply<0, true, false>(S0, src, t);
600  eval.template apply<1, true, false>(S1, t, dst);
601  eval.template apply<2, true, false>(S2, dst, t);
602  for (unsigned int i2 = 0, c = 0; i2 < n; ++i2)
603  for (unsigned int i1 = 0; i1 < n; ++i1)
604  for (unsigned int i0 = 0; i0 < n; ++i0, ++c)
605  t[c] /=
606  (eigenvalues[2][i2] + eigenvalues[1][i1] + eigenvalues[0][i0]);
607  eval.template apply<0, false, false>(S0, t, dst);
608  eval.template apply<1, false, false>(S1, dst, t);
609  eval.template apply<2, false, false>(S2, t, dst);
610  }
611 
612  else
613  Assert(false, ExcNotImplemented());
614 }
615 
616 
617 //---------------------- TensorProductMatrixSymmetricSum ----------------------
618 
619 template <int dim, typename Number, int size>
622  const std::array<Table<2, Number>, dim> &mass_matrix,
623  const std::array<Table<2, Number>, dim> &derivative_matrix)
624 {
626 }
627 
628 
629 
630 template <int dim, typename Number, int size>
633  const std::array<FullMatrix<Number>, dim> &mass_matrix,
634  const std::array<FullMatrix<Number>, dim> &derivative_matrix)
635 {
637 }
638 
639 
640 
641 template <int dim, typename Number, int size>
645 {
646  reinit(mass_matrix, derivative_matrix);
647 }
648 
649 
650 
651 template <int dim, typename Number, int size>
652 template <typename MatrixArray>
653 inline void
655  MatrixArray &&mass_matrices_,
656  MatrixArray &&derivative_matrices_)
657 {
658  auto &&mass_matrices = std::forward<MatrixArray>(mass_matrices_);
659  auto &&derivative_matrices = std::forward<MatrixArray>(derivative_matrices_);
660  this->mass_matrix = mass_matrices;
661  this->derivative_matrix = derivative_matrices;
662 
663  for (int dir = 0; dir < dim; ++dir)
664  {
665  Assert(size == -1 || (size > 0 && static_cast<unsigned int>(size) ==
666  mass_matrices[dir].n_rows()),
667  ExcDimensionMismatch(size, mass_matrices[dir].n_rows()));
668  AssertDimension(mass_matrices[dir].n_rows(), mass_matrices[dir].n_cols());
669  AssertDimension(mass_matrices[dir].n_rows(),
670  derivative_matrices[dir].n_rows());
671  AssertDimension(mass_matrices[dir].n_rows(),
672  derivative_matrices[dir].n_cols());
673 
674  this->eigenvectors[dir].reinit(mass_matrices[dir].n_cols(),
675  mass_matrices[dir].n_rows());
676  this->eigenvalues[dir].resize(mass_matrices[dir].n_cols());
677  internal::TensorProductMatrix ::spectral_assembly<Number>(
678  &(mass_matrices[dir](0, 0)),
679  &(derivative_matrices[dir](0, 0)),
680  mass_matrices[dir].n_rows(),
681  mass_matrices[dir].n_cols(),
682  this->eigenvalues[dir].begin(),
683  &(this->eigenvectors[dir](0, 0)));
684  }
685 }
686 
687 
688 
689 template <int dim, typename Number, int size>
690 inline void
692  const std::array<Table<2, Number>, dim> &mass_matrix,
693  const std::array<Table<2, Number>, dim> &derivative_matrix)
694 {
695  reinit_impl(mass_matrix, derivative_matrix);
696 }
697 
698 
699 
700 template <int dim, typename Number, int size>
701 inline void
703  const std::array<FullMatrix<Number>, dim> &mass_matrix,
704  const std::array<FullMatrix<Number>, dim> &derivative_matrix)
705 {
706  std::array<Table<2, Number>, dim> mass_copy;
707  std::array<Table<2, Number>, dim> deriv_copy;
708 
709  std::transform(mass_matrix.cbegin(),
710  mass_matrix.cend(),
711  mass_copy.begin(),
712  [](const FullMatrix<Number> &m) -> Table<2, Number> {
713  return m;
714  });
715  std::transform(derivative_matrix.cbegin(),
716  derivative_matrix.cend(),
717  deriv_copy.begin(),
718  [](const FullMatrix<Number> &m) -> Table<2, Number> {
719  return m;
720  });
721 
722  reinit_impl(std::move(mass_copy), std::move(deriv_copy));
723 }
724 
725 
726 
727 template <int dim, typename Number, int size>
728 inline void
730  const Table<2, Number> &mass_matrix,
731  const Table<2, Number> &derivative_matrix)
732 {
733  std::array<Table<2, Number>, dim> mass_matrices;
734  std::array<Table<2, Number>, dim> derivative_matrices;
735 
736  std::fill(mass_matrices.begin(), mass_matrices.end(), mass_matrix);
737  std::fill(derivative_matrices.begin(),
738  derivative_matrices.end(),
740 
741  reinit_impl(std::move(mass_matrices), std::move(derivative_matrices));
742 }
743 
744 
745 
746 //------------- vectorized spec.: TensorProductMatrixSymmetricSum -------------
747 
748 template <int dim, typename Number, int size>
750  TensorProductMatrixSymmetricSum(
751  const std::array<Table<2, VectorizedArray<Number>>, dim> &mass_matrix,
752  const std::array<Table<2, VectorizedArray<Number>>, dim> &derivative_matrix)
753 {
754  reinit(mass_matrix, derivative_matrix);
755 }
756 
757 
758 
759 template <int dim, typename Number, int size>
761  TensorProductMatrixSymmetricSum(
762  const Table<2, VectorizedArray<Number>> &mass_matrix,
763  const Table<2, VectorizedArray<Number>> &derivative_matrix)
764 {
765  reinit(mass_matrix, derivative_matrix);
766 }
767 
768 
769 
770 template <int dim, typename Number, int size>
771 template <typename MatrixArray>
772 inline void
774  reinit_impl(MatrixArray &&mass_matrices_, MatrixArray &&derivative_matrices_)
775 {
776  auto &&mass_matrix = std::forward<MatrixArray>(mass_matrices_);
777  auto &&derivative_matrix = std::forward<MatrixArray>(derivative_matrices_);
778  this->mass_matrix = mass_matrix;
779  this->derivative_matrix = derivative_matrix;
780 
781  constexpr unsigned int macro_size = VectorizedArray<Number>::n_array_elements;
782  std::size_t n_rows_max = (size > 0) ? size : 0;
783  if (size == -1)
784  for (unsigned int d = 0; d < dim; ++d)
785  n_rows_max = std::max(n_rows_max, mass_matrix[d].n_rows());
786  const std::size_t nm_flat_size_max = n_rows_max * n_rows_max * macro_size;
787  const std::size_t n_flat_size_max = n_rows_max * macro_size;
788 
789  std::vector<Number> mass_matrix_flat;
790  std::vector<Number> deriv_matrix_flat;
791  std::vector<Number> eigenvalues_flat;
792  std::vector<Number> eigenvectors_flat;
793  mass_matrix_flat.resize(nm_flat_size_max);
794  deriv_matrix_flat.resize(nm_flat_size_max);
795  eigenvalues_flat.resize(n_flat_size_max);
796  eigenvectors_flat.resize(nm_flat_size_max);
797  std::array<unsigned int, macro_size> offsets_nm;
798  std::array<unsigned int, macro_size> offsets_n;
799  for (int dir = 0; dir < dim; ++dir)
800  {
801  Assert(size == -1 || (size > 0 && static_cast<unsigned int>(size) ==
802  mass_matrix[dir].n_rows()),
803  ExcDimensionMismatch(size, mass_matrix[dir].n_rows()));
804  AssertDimension(mass_matrix[dir].n_rows(), mass_matrix[dir].n_cols());
805  AssertDimension(mass_matrix[dir].n_rows(),
806  derivative_matrix[dir].n_rows());
807  AssertDimension(mass_matrix[dir].n_rows(),
808  derivative_matrix[dir].n_cols());
809 
810  const unsigned int n_rows = mass_matrix[dir].n_rows();
811  const unsigned int n_cols = mass_matrix[dir].n_cols();
812  const unsigned int nm = n_rows * n_cols;
813  for (unsigned int vv = 0; vv < macro_size; ++vv)
814  offsets_nm[vv] = nm * vv;
815 
816  vectorized_transpose_and_store(false,
817  nm,
818  &(mass_matrix[dir](0, 0)),
819  offsets_nm.cbegin(),
820  mass_matrix_flat.data());
821  vectorized_transpose_and_store(false,
822  nm,
823  &(derivative_matrix[dir](0, 0)),
824  offsets_nm.cbegin(),
825  deriv_matrix_flat.data());
826 
827  const Number *mass_cbegin = mass_matrix_flat.data();
828  const Number *deriv_cbegin = deriv_matrix_flat.data();
829  Number * eigenvec_begin = eigenvectors_flat.data();
830  Number * eigenval_begin = eigenvalues_flat.data();
831  for (unsigned int lane = 0; lane < macro_size; ++lane)
832  internal::TensorProductMatrix ::spectral_assembly<Number>(
833  mass_cbegin + nm * lane,
834  deriv_cbegin + nm * lane,
835  n_rows,
836  n_cols,
837  eigenval_begin + n_rows * lane,
838  eigenvec_begin + nm * lane);
839 
840  this->eigenvalues[dir].resize(n_rows);
841  this->eigenvectors[dir].reinit(n_rows, n_cols);
842  for (unsigned int vv = 0; vv < macro_size; ++vv)
843  offsets_n[vv] = n_rows * vv;
844  vectorized_load_and_transpose(n_rows,
845  eigenvalues_flat.data(),
846  offsets_n.cbegin(),
847  this->eigenvalues[dir].begin());
848  vectorized_load_and_transpose(nm,
849  eigenvectors_flat.data(),
850  offsets_nm.cbegin(),
851  &(this->eigenvectors[dir](0, 0)));
852  }
853 }
854 
855 
856 
857 template <int dim, typename Number, int size>
858 inline void
860  const std::array<Table<2, VectorizedArray<Number>>, dim> &mass_matrix,
861  const std::array<Table<2, VectorizedArray<Number>>, dim> &derivative_matrix)
862 {
863  reinit_impl(mass_matrix, derivative_matrix);
864 }
865 
866 
867 
868 template <int dim, typename Number, int size>
869 inline void
871  const Table<2, VectorizedArray<Number>> &mass_matrix,
872  const Table<2, VectorizedArray<Number>> &derivative_matrix)
873 {
874  std::array<Table<2, VectorizedArray<Number>>, dim> mass_matrices;
875  std::array<Table<2, VectorizedArray<Number>>, dim> derivative_matrices;
876 
877  std::fill(mass_matrices.begin(), mass_matrices.end(), mass_matrix);
878  std::fill(derivative_matrices.begin(),
879  derivative_matrices.end(),
881 
882  reinit_impl(std::move(mass_matrices), std::move(derivative_matrices));
883 }
884 
885 
886 
887 #endif
888 
889 DEAL_II_NAMESPACE_CLOSE
890 
891 #endif
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1567
std::array< Table< 2, Number >, dim > derivative_matrix
std::array< AlignedVector< Number >, dim > eigenvalues
void reinit(const std::array< Table< 2, Number >, dim > &mass_matrix, const std::array< Table< 2, Number >, dim > &derivative_matrix)
const TableIndices< N > & size() const
std::complex< number > eigenvalue(const size_type i) const
#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)
#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
void resize_fast(const size_type size)
std::array< Table< 2, Number >, dim > eigenvectors
iterator begin()
void vmult(const ArrayView< Number > &dst, const ArrayView< const Number > &src) const
void compute_generalized_eigenvalues_symmetric(LAPACKFullMatrix< number > &B, const number lower_bound, const number upper_bound, const number abs_accuracy, Vector< number > &eigenvalues, std::vector< Vector< number >> &eigenvectors, const types::blas_int itype=1)
static ::ExceptionBase & ExcNotImplemented()
iterator begin() const
Definition: array_view.h:480
void apply_inverse(const ArrayView< Number > &dst, const ArrayView< const Number > &src) const