Reference documentation for deal.II version 9.2.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
tensor_product_matrix.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2017 - 2020 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 
75 template <int dim, typename Number, int n_rows_1d = -1>
76 class TensorProductMatrixSymmetricSumBase
77 {
78 public:
83  using value_type = Number;
84 
89  static constexpr int n_rows_1d_static = n_rows_1d;
90 
96  unsigned int
97  m() const;
98 
104  unsigned int
105  n() const;
106 
113  void
114  vmult(const ArrayView<Number> &dst, const ArrayView<const Number> &src) const;
115 
122  void
123  apply_inverse(const ArrayView<Number> & dst,
124  const ArrayView<const Number> &src) const;
125 
126 protected:
130  TensorProductMatrixSymmetricSumBase() = default;
131 
135  std::array<Table<2, Number>, dim> mass_matrix;
136 
140  std::array<Table<2, Number>, dim> derivative_matrix;
141 
146  std::array<AlignedVector<Number>, dim> eigenvalues;
147 
152  std::array<Table<2, Number>, dim> eigenvectors;
153 
154 private:
158  mutable AlignedVector<Number> tmp_array;
159 
163  mutable Threads::Mutex mutex;
164 };
165 
166 
167 
241 template <int dim, typename Number, int n_rows_1d = -1>
242 class TensorProductMatrixSymmetricSum
243  : public TensorProductMatrixSymmetricSumBase<dim, Number, n_rows_1d>
244 {
245 public:
249  TensorProductMatrixSymmetricSum() = default;
250 
257  TensorProductMatrixSymmetricSum(
258  const std::array<Table<2, Number>, dim> &mass_matrix,
259  const std::array<Table<2, Number>, dim> &derivative_matrix);
260 
267  TensorProductMatrixSymmetricSum(
268  const std::array<FullMatrix<Number>, dim> &mass_matrix,
269  const std::array<FullMatrix<Number>, dim> &derivative_matrix);
270 
275  TensorProductMatrixSymmetricSum(const Table<2, Number> &mass_matrix,
276  const Table<2, Number> &derivative_matrix);
277 
289  void
290  reinit(const std::array<Table<2, Number>, dim> &mass_matrix,
291  const std::array<Table<2, Number>, dim> &derivative_matrix);
292 
298  void
299  reinit(const std::array<FullMatrix<Number>, dim> &mass_matrix,
300  const std::array<FullMatrix<Number>, dim> &derivative_matrix);
301 
307  void
308  reinit(const Table<2, Number> &mass_matrix,
309  const Table<2, Number> &derivative_matrix);
310 
311 private:
320  template <typename MatrixArray>
321  void
322  reinit_impl(MatrixArray &&mass_matrix, MatrixArray &&derivative_matrix);
323 };
324 
325 
326 
335 template <int dim, typename Number, int n_rows_1d>
336 class TensorProductMatrixSymmetricSum<dim, VectorizedArray<Number>, n_rows_1d>
337  : public TensorProductMatrixSymmetricSumBase<dim,
338  VectorizedArray<Number>,
339  n_rows_1d>
340 {
341 public:
345  TensorProductMatrixSymmetricSum() = default;
346 
353  TensorProductMatrixSymmetricSum(
354  const std::array<Table<2, VectorizedArray<Number>>, dim> &mass_matrix,
355  const std::array<Table<2, VectorizedArray<Number>>, dim>
356  &derivative_matrix);
357 
364  TensorProductMatrixSymmetricSum(
365  const Table<2, VectorizedArray<Number>> &mass_matrix,
366  const Table<2, VectorizedArray<Number>> &derivative_matrix);
367 
379  void
380  reinit(const std::array<Table<2, VectorizedArray<Number>>, dim> &mass_matrix,
381  const std::array<Table<2, VectorizedArray<Number>>, dim>
382  &derivative_matrix);
383 
389  void
390  reinit(const Table<2, VectorizedArray<Number>> &mass_matrix,
391  const Table<2, VectorizedArray<Number>> &derivative_matrix);
392 
393 private:
402  template <typename MatrixArray>
403  void
404  reinit_impl(MatrixArray &&mass_matrix, MatrixArray &&derivative_matrix);
405 };
406 
407 
408 /*----------------------- Inline functions ----------------------------------*/
409 
410 #ifndef DOXYGEN
411 
412 namespace internal
413 {
414  namespace TensorProductMatrix
415  {
424  template <typename Number>
425  void
426  spectral_assembly(const Number * mass_matrix,
427  const Number * derivative_matrix,
428  const unsigned int n_rows,
429  const unsigned int n_cols,
430  Number * eigenvalues,
431  Number * eigenvectors)
432  {
433  Assert(n_rows == n_cols, ExcNotImplemented());
434 
435  auto &&transpose_fill_nm = [](Number * out,
436  const Number * in,
437  const unsigned int n,
438  const unsigned int m) {
439  for (unsigned int mm = 0; mm < m; ++mm)
440  for (unsigned int nn = 0; nn < n; ++nn)
441  out[mm + nn * m] = *(in++);
442  };
443 
444  std::vector<::Vector<Number>> eigenvecs(n_rows);
445  LAPACKFullMatrix<Number> mass_copy(n_rows, n_cols);
446  LAPACKFullMatrix<Number> deriv_copy(n_rows, n_cols);
447 
448  transpose_fill_nm(&(mass_copy(0, 0)), mass_matrix, n_rows, n_cols);
449  transpose_fill_nm(&(deriv_copy(0, 0)), derivative_matrix, n_rows, n_cols);
450 
451  deriv_copy.compute_generalized_eigenvalues_symmetric(mass_copy,
452  eigenvecs);
453  AssertDimension(eigenvecs.size(), n_rows);
454  for (unsigned int i = 0; i < n_rows; ++i)
455  for (unsigned int j = 0; j < n_cols; ++j, ++eigenvectors)
456  *eigenvectors = eigenvecs[j][i];
457 
458  for (unsigned int i = 0; i < n_rows; ++i, ++eigenvalues)
459  *eigenvalues = deriv_copy.eigenvalue(i).real();
460  }
461  } // namespace TensorProductMatrix
462 } // namespace internal
463 
464 
465 template <int dim, typename Number, int n_rows_1d>
466 inline unsigned int
467 TensorProductMatrixSymmetricSumBase<dim, Number, n_rows_1d>::m() const
468 {
469  unsigned int m = mass_matrix[0].n_rows();
470  for (unsigned int d = 1; d < dim; ++d)
471  m *= mass_matrix[d].n_rows();
472  return m;
473 }
474 
475 
476 
477 template <int dim, typename Number, int n_rows_1d>
478 inline unsigned int
479 TensorProductMatrixSymmetricSumBase<dim, Number, n_rows_1d>::n() const
480 {
481  unsigned int n = mass_matrix[0].n_cols();
482  for (unsigned int d = 1; d < dim; ++d)
483  n *= mass_matrix[d].n_cols();
484  return n;
485 }
486 
487 
488 
489 template <int dim, typename Number, int n_rows_1d>
490 inline void
491 TensorProductMatrixSymmetricSumBase<dim, Number, n_rows_1d>::vmult(
492  const ArrayView<Number> & dst_view,
493  const ArrayView<const Number> &src_view) const
494 {
495  AssertDimension(dst_view.size(), this->m());
496  AssertDimension(src_view.size(), this->n());
497  std::lock_guard<std::mutex> lock(this->mutex);
498  const unsigned int n = Utilities::fixed_power<dim>(
499  n_rows_1d > 0 ? n_rows_1d : eigenvalues[0].size());
500  tmp_array.resize_fast(n * 2);
501  constexpr int kernel_size = n_rows_1d > 0 ? n_rows_1d : 0;
502  internal::EvaluatorTensorProduct<internal::evaluate_general,
503  dim,
504  kernel_size,
505  kernel_size,
506  Number>
507  eval(AlignedVector<Number>{},
508  AlignedVector<Number>{},
509  AlignedVector<Number>{},
510  mass_matrix[0].n_rows(),
511  mass_matrix[0].n_rows());
512  Number * t = tmp_array.begin();
513  const Number *src = src_view.begin();
514  Number * dst = dst_view.data();
515 
516  if (dim == 1)
517  {
518  const Number *A = &derivative_matrix[0](0, 0);
519  eval.template apply<0, false, false>(A, src, dst);
520  }
521 
522  else if (dim == 2)
523  {
524  const Number *A0 = &derivative_matrix[0](0, 0);
525  const Number *M0 = &mass_matrix[0](0, 0);
526  const Number *A1 = &derivative_matrix[1](0, 0);
527  const Number *M1 = &mass_matrix[1](0, 0);
528  eval.template apply<0, false, false>(M0, src, t);
529  eval.template apply<1, false, false>(A1, t, dst);
530  eval.template apply<0, false, false>(A0, src, t);
531  eval.template apply<1, false, true>(M1, t, dst);
532  }
533 
534  else if (dim == 3)
535  {
536  const Number *A0 = &derivative_matrix[0](0, 0);
537  const Number *M0 = &mass_matrix[0](0, 0);
538  const Number *A1 = &derivative_matrix[1](0, 0);
539  const Number *M1 = &mass_matrix[1](0, 0);
540  const Number *A2 = &derivative_matrix[2](0, 0);
541  const Number *M2 = &mass_matrix[2](0, 0);
542  eval.template apply<0, false, false>(M0, src, t + n);
543  eval.template apply<1, false, false>(M1, t + n, t);
544  eval.template apply<2, false, false>(A2, t, dst);
545  eval.template apply<1, false, false>(A1, t + n, t);
546  eval.template apply<0, false, false>(A0, src, t + n);
547  eval.template apply<1, false, true>(M1, t + n, t);
548  eval.template apply<2, false, true>(M2, t, dst);
549  }
550 
551  else
552  AssertThrow(false, ExcNotImplemented());
553 }
554 
555 
556 
557 template <int dim, typename Number, int n_rows_1d>
558 inline void
559 TensorProductMatrixSymmetricSumBase<dim, Number, n_rows_1d>::apply_inverse(
560  const ArrayView<Number> & dst_view,
561  const ArrayView<const Number> &src_view) const
562 {
563  AssertDimension(dst_view.size(), this->n());
564  AssertDimension(src_view.size(), this->m());
565  std::lock_guard<std::mutex> lock(this->mutex);
566  const unsigned int n = n_rows_1d > 0 ? n_rows_1d : eigenvalues[0].size();
567  tmp_array.resize_fast(Utilities::fixed_power<dim>(n));
568  constexpr int kernel_size = n_rows_1d > 0 ? n_rows_1d : 0;
569  internal::EvaluatorTensorProduct<internal::evaluate_general,
570  dim,
571  kernel_size,
572  kernel_size,
573  Number>
574  eval(AlignedVector<Number>(),
575  AlignedVector<Number>(),
576  AlignedVector<Number>(),
577  mass_matrix[0].n_rows(),
578  mass_matrix[0].n_rows());
579  Number * t = tmp_array.begin();
580  const Number *src = src_view.data();
581  Number * dst = dst_view.data();
582 
583  // NOTE: dof_to_quad has to be interpreted as 'dof to eigenvalue index'
584  // --> apply<.,true,.> (S,src,dst) calculates dst = S^T * src,
585  // --> apply<.,false,.> (S,src,dst) calculates dst = S * src,
586  // while the eigenvectors are stored column-wise in S, i.e.
587  // rows correspond to dofs whereas columns to eigenvalue indices!
588  if (dim == 1)
589  {
590  const Number *S = &eigenvectors[0](0, 0);
591  eval.template apply<0, true, false>(S, src, t);
592  for (unsigned int i = 0; i < n; ++i)
593  t[i] /= eigenvalues[0][i];
594  eval.template apply<0, false, false>(S, t, dst);
595  }
596 
597  else if (dim == 2)
598  {
599  const Number *S0 = &(eigenvectors[0](0, 0));
600  const Number *S1 = &(eigenvectors[1](0, 0));
601  eval.template apply<0, true, false>(S0, src, t);
602  eval.template apply<1, true, false>(S1, t, dst);
603  for (unsigned int i1 = 0, c = 0; i1 < n; ++i1)
604  for (unsigned int i0 = 0; i0 < n; ++i0, ++c)
605  dst[c] /= (eigenvalues[1][i1] + eigenvalues[0][i0]);
606  eval.template apply<0, false, false>(S0, dst, t);
607  eval.template apply<1, false, false>(S1, t, dst);
608  }
609 
610  else if (dim == 3)
611  {
612  const Number *S0 = &eigenvectors[0](0, 0);
613  const Number *S1 = &eigenvectors[1](0, 0);
614  const Number *S2 = &eigenvectors[2](0, 0);
615  eval.template apply<0, true, false>(S0, src, t);
616  eval.template apply<1, true, false>(S1, t, dst);
617  eval.template apply<2, true, false>(S2, dst, t);
618  for (unsigned int i2 = 0, c = 0; i2 < n; ++i2)
619  for (unsigned int i1 = 0; i1 < n; ++i1)
620  for (unsigned int i0 = 0; i0 < n; ++i0, ++c)
621  t[c] /=
622  (eigenvalues[2][i2] + eigenvalues[1][i1] + eigenvalues[0][i0]);
623  eval.template apply<0, false, false>(S0, t, dst);
624  eval.template apply<1, false, false>(S1, dst, t);
625  eval.template apply<2, false, false>(S2, t, dst);
626  }
627 
628  else
629  Assert(false, ExcNotImplemented());
630 }
631 
632 
633 //---------------------- TensorProductMatrixSymmetricSum ----------------------
634 
635 template <int dim, typename Number, int n_rows_1d>
636 inline TensorProductMatrixSymmetricSum<dim, Number, n_rows_1d>::
637  TensorProductMatrixSymmetricSum(
638  const std::array<Table<2, Number>, dim> &mass_matrix,
639  const std::array<Table<2, Number>, dim> &derivative_matrix)
640 {
641  reinit(mass_matrix, derivative_matrix);
642 }
643 
644 
645 
646 template <int dim, typename Number, int n_rows_1d>
647 inline TensorProductMatrixSymmetricSum<dim, Number, n_rows_1d>::
648  TensorProductMatrixSymmetricSum(
649  const std::array<FullMatrix<Number>, dim> &mass_matrix,
650  const std::array<FullMatrix<Number>, dim> &derivative_matrix)
651 {
652  reinit(mass_matrix, derivative_matrix);
653 }
654 
655 
656 
657 template <int dim, typename Number, int n_rows_1d>
658 inline TensorProductMatrixSymmetricSum<dim, Number, n_rows_1d>::
659  TensorProductMatrixSymmetricSum(const Table<2, Number> &mass_matrix,
660  const Table<2, Number> &derivative_matrix)
661 {
662  reinit(mass_matrix, derivative_matrix);
663 }
664 
665 
666 
667 template <int dim, typename Number, int n_rows_1d>
668 template <typename MatrixArray>
669 inline void
670 TensorProductMatrixSymmetricSum<dim, Number, n_rows_1d>::reinit_impl(
671  MatrixArray &&mass_matrices_,
672  MatrixArray &&derivative_matrices_)
673 {
674  auto &&mass_matrices = std::forward<MatrixArray>(mass_matrices_);
675  auto &&derivative_matrices = std::forward<MatrixArray>(derivative_matrices_);
676  this->mass_matrix = mass_matrices;
677  this->derivative_matrix = derivative_matrices;
678 
679  for (int dir = 0; dir < dim; ++dir)
680  {
681  Assert(n_rows_1d == -1 ||
682  (n_rows_1d > 0 && static_cast<unsigned int>(n_rows_1d) ==
683  mass_matrices[dir].n_rows()),
684  ExcDimensionMismatch(n_rows_1d, mass_matrices[dir].n_rows()));
685  AssertDimension(mass_matrices[dir].n_rows(), mass_matrices[dir].n_cols());
686  AssertDimension(mass_matrices[dir].n_rows(),
687  derivative_matrices[dir].n_rows());
688  AssertDimension(mass_matrices[dir].n_rows(),
689  derivative_matrices[dir].n_cols());
690 
691  this->eigenvectors[dir].reinit(mass_matrices[dir].n_cols(),
692  mass_matrices[dir].n_rows());
693  this->eigenvalues[dir].resize(mass_matrices[dir].n_cols());
694  internal::TensorProductMatrix::spectral_assembly<Number>(
695  &(mass_matrices[dir](0, 0)),
696  &(derivative_matrices[dir](0, 0)),
697  mass_matrices[dir].n_rows(),
698  mass_matrices[dir].n_cols(),
699  this->eigenvalues[dir].begin(),
700  &(this->eigenvectors[dir](0, 0)));
701  }
702 }
703 
704 
705 
706 template <int dim, typename Number, int n_rows_1d>
707 inline void
708 TensorProductMatrixSymmetricSum<dim, Number, n_rows_1d>::reinit(
709  const std::array<Table<2, Number>, dim> &mass_matrix,
710  const std::array<Table<2, Number>, dim> &derivative_matrix)
711 {
712  reinit_impl(mass_matrix, derivative_matrix);
713 }
714 
715 
716 
717 template <int dim, typename Number, int n_rows_1d>
718 inline void
719 TensorProductMatrixSymmetricSum<dim, Number, n_rows_1d>::reinit(
720  const std::array<FullMatrix<Number>, dim> &mass_matrix,
721  const std::array<FullMatrix<Number>, dim> &derivative_matrix)
722 {
723  std::array<Table<2, Number>, dim> mass_copy;
724  std::array<Table<2, Number>, dim> deriv_copy;
725 
726  std::transform(mass_matrix.cbegin(),
727  mass_matrix.cend(),
728  mass_copy.begin(),
729  [](const FullMatrix<Number> &m) -> Table<2, Number> {
730  return m;
731  });
732  std::transform(derivative_matrix.cbegin(),
733  derivative_matrix.cend(),
734  deriv_copy.begin(),
735  [](const FullMatrix<Number> &m) -> Table<2, Number> {
736  return m;
737  });
738 
739  reinit_impl(std::move(mass_copy), std::move(deriv_copy));
740 }
741 
742 
743 
744 template <int dim, typename Number, int n_rows_1d>
745 inline void
746 TensorProductMatrixSymmetricSum<dim, Number, n_rows_1d>::reinit(
747  const Table<2, Number> &mass_matrix,
748  const Table<2, Number> &derivative_matrix)
749 {
750  std::array<Table<2, Number>, dim> mass_matrices;
751  std::array<Table<2, Number>, dim> derivative_matrices;
752 
753  std::fill(mass_matrices.begin(), mass_matrices.end(), mass_matrix);
754  std::fill(derivative_matrices.begin(),
755  derivative_matrices.end(),
756  derivative_matrix);
757 
758  reinit_impl(std::move(mass_matrices), std::move(derivative_matrices));
759 }
760 
761 
762 
763 //------------- vectorized spec.: TensorProductMatrixSymmetricSum -------------
764 
765 template <int dim, typename Number, int n_rows_1d>
766 inline TensorProductMatrixSymmetricSum<dim,
767  VectorizedArray<Number>,
768  n_rows_1d>::
769  TensorProductMatrixSymmetricSum(
770  const std::array<Table<2, VectorizedArray<Number>>, dim> &mass_matrix,
771  const std::array<Table<2, VectorizedArray<Number>>, dim> &derivative_matrix)
772 {
773  reinit(mass_matrix, derivative_matrix);
774 }
775 
776 
777 
778 template <int dim, typename Number, int n_rows_1d>
779 inline TensorProductMatrixSymmetricSum<dim,
780  VectorizedArray<Number>,
781  n_rows_1d>::
782  TensorProductMatrixSymmetricSum(
783  const Table<2, VectorizedArray<Number>> &mass_matrix,
784  const Table<2, VectorizedArray<Number>> &derivative_matrix)
785 {
786  reinit(mass_matrix, derivative_matrix);
787 }
788 
789 
790 
791 template <int dim, typename Number, int n_rows_1d>
792 template <typename MatrixArray>
793 inline void
794 TensorProductMatrixSymmetricSum<dim, VectorizedArray<Number>, n_rows_1d>::
795  reinit_impl(MatrixArray &&mass_matrices_, MatrixArray &&derivative_matrices_)
796 {
797  auto &&mass_matrix = std::forward<MatrixArray>(mass_matrices_);
798  auto &&derivative_matrix = std::forward<MatrixArray>(derivative_matrices_);
799  this->mass_matrix = mass_matrix;
800  this->derivative_matrix = derivative_matrix;
801 
802  constexpr unsigned int macro_size = VectorizedArray<Number>::size();
803  std::size_t n_rows_max = (n_rows_1d > 0) ? n_rows_1d : 0;
804  if (n_rows_1d == -1)
805  for (unsigned int d = 0; d < dim; ++d)
806  n_rows_max = std::max(n_rows_max, mass_matrix[d].n_rows());
807  const std::size_t nm_flat_size_max = n_rows_max * n_rows_max * macro_size;
808  const std::size_t n_flat_size_max = n_rows_max * macro_size;
809 
810  std::vector<Number> mass_matrix_flat;
811  std::vector<Number> deriv_matrix_flat;
812  std::vector<Number> eigenvalues_flat;
813  std::vector<Number> eigenvectors_flat;
814  mass_matrix_flat.resize(nm_flat_size_max);
815  deriv_matrix_flat.resize(nm_flat_size_max);
816  eigenvalues_flat.resize(n_flat_size_max);
817  eigenvectors_flat.resize(nm_flat_size_max);
818  std::array<unsigned int, macro_size> offsets_nm;
819  std::array<unsigned int, macro_size> offsets_n;
820  for (int dir = 0; dir < dim; ++dir)
821  {
822  Assert(n_rows_1d == -1 ||
823  (n_rows_1d > 0 && static_cast<unsigned int>(n_rows_1d) ==
824  mass_matrix[dir].n_rows()),
825  ExcDimensionMismatch(n_rows_1d, mass_matrix[dir].n_rows()));
826  AssertDimension(mass_matrix[dir].n_rows(), mass_matrix[dir].n_cols());
827  AssertDimension(mass_matrix[dir].n_rows(),
828  derivative_matrix[dir].n_rows());
829  AssertDimension(mass_matrix[dir].n_rows(),
830  derivative_matrix[dir].n_cols());
831 
832  const unsigned int n_rows = mass_matrix[dir].n_rows();
833  const unsigned int n_cols = mass_matrix[dir].n_cols();
834  const unsigned int nm = n_rows * n_cols;
835  for (unsigned int vv = 0; vv < macro_size; ++vv)
836  offsets_nm[vv] = nm * vv;
837 
838  vectorized_transpose_and_store(false,
839  nm,
840  &(mass_matrix[dir](0, 0)),
841  offsets_nm.cbegin(),
842  mass_matrix_flat.data());
843  vectorized_transpose_and_store(false,
844  nm,
845  &(derivative_matrix[dir](0, 0)),
846  offsets_nm.cbegin(),
847  deriv_matrix_flat.data());
848 
849  const Number *mass_cbegin = mass_matrix_flat.data();
850  const Number *deriv_cbegin = deriv_matrix_flat.data();
851  Number * eigenvec_begin = eigenvectors_flat.data();
852  Number * eigenval_begin = eigenvalues_flat.data();
853  for (unsigned int lane = 0; lane < macro_size; ++lane)
854  internal::TensorProductMatrix::spectral_assembly<Number>(
855  mass_cbegin + nm * lane,
856  deriv_cbegin + nm * lane,
857  n_rows,
858  n_cols,
859  eigenval_begin + n_rows * lane,
860  eigenvec_begin + nm * lane);
861 
862  this->eigenvalues[dir].resize(n_rows);
863  this->eigenvectors[dir].reinit(n_rows, n_cols);
864  for (unsigned int vv = 0; vv < macro_size; ++vv)
865  offsets_n[vv] = n_rows * vv;
866  vectorized_load_and_transpose(n_rows,
867  eigenvalues_flat.data(),
868  offsets_n.cbegin(),
869  this->eigenvalues[dir].begin());
870  vectorized_load_and_transpose(nm,
871  eigenvectors_flat.data(),
872  offsets_nm.cbegin(),
873  &(this->eigenvectors[dir](0, 0)));
874  }
875 }
876 
877 
878 
879 template <int dim, typename Number, int n_rows_1d>
880 inline void
881 TensorProductMatrixSymmetricSum<dim, VectorizedArray<Number>, n_rows_1d>::
882  reinit(
883  const std::array<Table<2, VectorizedArray<Number>>, dim> &mass_matrix,
884  const std::array<Table<2, VectorizedArray<Number>>, dim> &derivative_matrix)
885 {
886  reinit_impl(mass_matrix, derivative_matrix);
887 }
888 
889 
890 
891 template <int dim, typename Number, int n_rows_1d>
892 inline void
893 TensorProductMatrixSymmetricSum<dim, VectorizedArray<Number>, n_rows_1d>::
894  reinit(const Table<2, VectorizedArray<Number>> &mass_matrix,
895  const Table<2, VectorizedArray<Number>> &derivative_matrix)
896 {
897  std::array<Table<2, VectorizedArray<Number>>, dim> mass_matrices;
898  std::array<Table<2, VectorizedArray<Number>>, dim> derivative_matrices;
899 
900  std::fill(mass_matrices.begin(), mass_matrices.end(), mass_matrix);
901  std::fill(derivative_matrices.begin(),
902  derivative_matrices.end(),
903  derivative_matrix);
904 
905  reinit_impl(std::move(mass_matrices), std::move(derivative_matrices));
906 }
907 
908 
909 
910 #endif
911 
912 DEAL_II_NAMESPACE_CLOSE
913 
914 #endif
array_view.h
ArrayView::begin
iterator begin() const
Definition: array_view.h:493
TensorProductMatrixSymmetricSumBase::vmult
void vmult(const ArrayView< Number > &dst, const ArrayView< const Number > &src) const
LinearAlgebraDealII::Vector
Vector< double > Vector
Definition: generic_linear_algebra.h:43
LocalIntegrators::L2::mass_matrix
void mass_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const double factor=1.)
Definition: l2.h:63
tensor_product_kernels.h
TensorProductMatrixSymmetricSumBase::eigenvalues
std::array< AlignedVector< Number >, dim > eigenvalues
Definition: tensor_product_matrix.h:146
Threads::Mutex
Definition: thread_management.h:91
TensorProductMatrixSymmetricSumBase::tmp_array
AlignedVector< Number > tmp_array
Definition: tensor_product_matrix.h:158
internal::EvaluatorTensorProduct
Definition: tensor_product_kernels.h:96
VectorizedArray::vectorized_transpose_and_store
void vectorized_transpose_and_store(const bool add_into, const unsigned int n_entries, const VectorizedArray< Number, width > *in, const unsigned int *offsets, Number *out)
Definition: vectorization.h:882
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
TensorProductMatrixSymmetricSumBase::TensorProductMatrixSymmetricSumBase
TensorProductMatrixSymmetricSumBase()=default
ArrayView
Definition: array_view.h:77
parallel::transform
void transform(const InputIterator &begin_in, const InputIterator &end_in, OutputIterator out, const Predicate &predicate, const unsigned int grainsize)
Definition: parallel.h:211
TensorProductMatrixSymmetricSumBase::n_rows_1d_static
static constexpr int n_rows_1d_static
Definition: tensor_product_matrix.h:89
SymmetricTensor::eigenvectors
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)
TensorProductMatrixSymmetricSumBase::apply_inverse
void apply_inverse(const ArrayView< Number > &dst, const ArrayView< const Number > &src) const
VectorizedArray< Number >
TensorProductMatrixSymmetricSum::reinit
void reinit(const std::array< Table< 2, Number >, dim > &mass_matrix, const std::array< Table< 2, Number >, dim > &derivative_matrix)
ArrayView::size
std::size_t size() const
Definition: array_view.h:484
Physics::Elasticity::Kinematics::d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
ArrayView::data
value_type * data() const noexcept
Definition: array_view.h:461
LAPACKFullMatrix< Number >
TensorProductMatrixSymmetricSumBase::n
unsigned int n() const
Table< 2, Number >
VectorizedArray::vectorized_load_and_transpose
void vectorized_load_and_transpose(const unsigned int n_entries, const Number *in, const unsigned int *offsets, VectorizedArray< Number, width > *out)
Definition: vectorization.h:807
thread_management.h
TrilinosWrappers::internal::begin
VectorType::value_type * begin(VectorType &V)
Definition: trilinos_sparse_matrix.cc:51
LAPACKSupport::eigenvalues
@ eigenvalues
Eigenvalue vector is filled.
Definition: lapack_support.h:68
internal::reinit
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:621
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
TensorProductMatrixSymmetricSumBase::mutex
Threads::Mutex mutex
Definition: tensor_product_matrix.h:163
VectorizedArrayBase< VectorizedArray< Number, width >, 1 >::size
static constexpr std::size_t size()
Definition: vectorization.h:261
AlignedVector< Number >
internal::evaluate_general
@ evaluate_general
Definition: tensor_product_kernels.h:42
lapack_full_matrix.h
AssertDimension
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1579
LAPACKSupport::A
static const char A
Definition: lapack_support.h:155
TensorProductMatrixSymmetricSum::TensorProductMatrixSymmetricSum
TensorProductMatrixSymmetricSum()=default
TensorProductMatrixSymmetricSumBase::derivative_matrix
std::array< Table< 2, Number >, dim > derivative_matrix
Definition: tensor_product_matrix.h:140
TensorProductMatrixSymmetricSumBase::m
unsigned int m() const
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1419
TensorProductMatrixSymmetricSumBase::eigenvectors
std::array< Table< 2, Number >, dim > eigenvectors
Definition: tensor_product_matrix.h:152
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
SymmetricTensor::eigenvalues
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
config.h
TensorProductMatrixSymmetricSum
Definition: tensor_product_matrix.h:242
FullMatrix
Definition: full_matrix.h:71
internal
Definition: aligned_vector.h:369
TensorProductMatrixSymmetricSumBase::mass_matrix
std::array< Table< 2, Number >, dim > mass_matrix
Definition: tensor_product_matrix.h:135
TensorProductMatrixSymmetricSumBase
Definition: tensor_product_matrix.h:76
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
TensorProductMatrixSymmetricSum::reinit_impl
void reinit_impl(MatrixArray &&mass_matrix, MatrixArray &&derivative_matrix)
AssertThrow
#define AssertThrow(cond, exc)
Definition: exceptions.h:1531
Utilities::MPI::max
T max(const T &t, const MPI_Comm &mpi_communicator)
eigenvectors
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)
Generated by   doxygen 1.8.17