Reference documentation for deal.II version GIT relicensing-399-g79d89019c5 2024-04-16 15:00:02+00:00
\(\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\}}\)
Loading...
Searching...
No Matches
tensor_product_matrix.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2017 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_tensor_product_matrix_h
16#define dealii_tensor_product_matrix_h
17
18
19#include <deal.II/base/config.h>
20
23#include <deal.II/base/mutex.h>
25
27
29
31
32// Forward declarations
33#ifndef DOXYGEN
34template <typename>
35class Vector;
36template <typename>
37class FullMatrix;
38#endif
39
111template <int dim, typename Number, int n_rows_1d = -1>
113{
114public:
119 using value_type = Number;
120
125 static constexpr int n_rows_1d_static = n_rows_1d;
126
131
136 template <typename T>
138 const T &derivative_matrix);
139
157 template <typename T>
158 void
160
166 unsigned int
167 m() const;
168
174 unsigned int
175 n() const;
176
190 void
191 vmult(const ArrayView<Number> &dst, const ArrayView<const Number> &src) const;
192
198 void
200 const ArrayView<const Number> &src,
201 AlignedVector<Number> &tmp) const;
202
209 void
211 const ArrayView<const Number> &src) const;
212
216 std::size_t
218
219protected:
223 std::array<Table<2, Number>, dim> mass_matrix;
224
228 std::array<Table<2, Number>, dim> derivative_matrix;
229
234 std::array<AlignedVector<Number>, dim> eigenvalues;
235
240 std::array<Table<2, Number>, dim> eigenvectors;
241
242private:
247
252};
253
254
255
256namespace internal
257{
259 {
260 template <typename Number>
262 {
266 static constexpr std::size_t width = VectorizedArrayTrait::width();
267
269 std::pair<std::bitset<width>,
270 std::pair<Table<2, Number>, Table<2, Number>>>;
271
273 : eps(std::sqrt(std::numeric_limits<ScalarNumber>::epsilon()))
274 {}
275
276 bool
277 operator()(const MatrixPairType &left, const MatrixPairType &right) const
278 {
279 const auto &M_0 = left.second.first;
280 const auto &K_0 = left.second.second;
281 const auto &M_1 = right.second.first;
282 const auto &K_1 = right.second.second;
283
284 std::bitset<width> mask;
285
286 for (unsigned int v = 0; v < width; ++v)
287 mask[v] = left.first[v] && right.first[v];
288
289 const FloatingPointComparator<Number> comparator(
290 eps, false /*use relative tolerance*/, mask);
291
292 if (comparator(M_0, M_1))
293 return true;
294 else if (comparator(M_1, M_0))
295 return false;
296 else if (comparator(K_0, K_1))
297 return true;
298 else
299 return false;
300 }
301
302 private:
304 };
305 } // namespace TensorProductMatrixSymmetricSum
306} // namespace internal
307
308
309
333template <int dim, typename Number, int n_rows_1d = -1>
335{
336 using MatrixPairType = std::pair<Table<2, Number>, Table<2, Number>>;
337
338 using MatrixPairTypeWithMask = std::pair<
339 std::bitset<::internal::VectorizedArrayTrait<Number>::width()>,
341
342public:
347 {
352 const bool precompute_inverse_diagonal = true);
353
358
363 };
364
369 const AdditionalData &additional_data = AdditionalData());
370
375 void
376 reserve(const unsigned int size);
377
383 template <typename T>
384 void
385 insert(const unsigned int index, const T &Ms, const T &Ks);
386
391 void
393
397 void
398 apply_inverse(const unsigned int index,
399 const ArrayView<Number> &dst_in,
400 const ArrayView<const Number> &src_in) const;
401
405 std::size_t
407
416 std::size_t
418
419private:
424
429
434 std::vector<MatrixPairType> mass_and_derivative_matrices;
435
440 std::map<
442 unsigned int,
445
453 std::vector<unsigned int> indices;
454
459
464
469
474
479
483 std::vector<unsigned int> vector_ptr;
484
488 std::vector<unsigned int> matrix_ptr;
489
493 std::vector<unsigned int> vector_n_rows_1d;
494};
495
496
497/*----------------------- Inline functions ----------------------------------*/
498
499#ifndef DOXYGEN
500
501namespace internal
502{
504 {
513 template <typename Number>
514 void
515 spectral_assembly(const Number *mass_matrix,
516 const Number *derivative_matrix,
517 const unsigned int n_rows,
518 const unsigned int n_cols,
519 Number *eigenvalues,
520 Number *eigenvectors)
521 {
522 Assert(n_rows == n_cols, ExcNotImplemented());
523
524 std::vector<bool> constrained_dofs(n_rows, false);
525
526 for (unsigned int i = 0; i < n_rows; ++i)
527 {
528 if (mass_matrix[i + i * n_rows] == 0.0)
529 {
530 Assert(derivative_matrix[i + i * n_rows] == 0.0,
532
533 for (unsigned int j = 0; j < n_rows; ++j)
534 {
535 Assert(derivative_matrix[i + j * n_rows] == 0,
537 Assert(derivative_matrix[j + i * n_rows] == 0,
539 }
540
541 constrained_dofs[i] = true;
542 }
543 }
544
545 const auto transpose_fill_nm = [&constrained_dofs](Number *out,
546 const Number *in,
547 const unsigned int n,
548 const unsigned int m) {
549 for (unsigned int mm = 0, c = 0; mm < m; ++mm)
550 for (unsigned int nn = 0; nn < n; ++nn, ++c)
551 out[mm + nn * m] =
552 (mm == nn && constrained_dofs[mm]) ? Number(1.0) : in[c];
553 };
554
555 std::vector<::Vector<Number>> eigenvecs(n_rows);
556 LAPACKFullMatrix<Number> mass_copy(n_rows, n_cols);
557 LAPACKFullMatrix<Number> deriv_copy(n_rows, n_cols);
558
559 transpose_fill_nm(&(mass_copy(0, 0)), mass_matrix, n_rows, n_cols);
560 transpose_fill_nm(&(deriv_copy(0, 0)), derivative_matrix, n_rows, n_cols);
561
563 eigenvecs);
564 AssertDimension(eigenvecs.size(), n_rows);
565 for (unsigned int i = 0, c = 0; i < n_rows; ++i)
566 for (unsigned int j = 0; j < n_cols; ++j, ++c)
567 if (constrained_dofs[i] == false)
568 eigenvectors[c] = eigenvecs[j][i];
569
570 for (unsigned int i = 0; i < n_rows; ++i, ++eigenvalues)
571 *eigenvalues = deriv_copy.eigenvalue(i).real();
572 }
573
574
575
576 template <std::size_t dim, typename Number>
577 inline void
578 setup(const std::array<Table<2, Number>, dim> &mass_matrix,
579 const std::array<Table<2, Number>, dim> &derivative_matrix,
580 std::array<Table<2, Number>, dim> &eigenvectors,
581 std::array<AlignedVector<Number>, dim> &eigenvalues)
582 {
583 const unsigned int n_rows_1d = mass_matrix[0].n_cols();
584 (void)n_rows_1d;
585
586 for (unsigned int dir = 0; dir < dim; ++dir)
587 {
588 AssertDimension(n_rows_1d, mass_matrix[dir].n_cols());
589 AssertDimension(mass_matrix[dir].n_rows(), mass_matrix[dir].n_cols());
590 AssertDimension(mass_matrix[dir].n_rows(),
591 derivative_matrix[dir].n_rows());
592 AssertDimension(mass_matrix[dir].n_rows(),
593 derivative_matrix[dir].n_cols());
594
595 eigenvectors[dir].reinit(mass_matrix[dir].n_cols(),
596 mass_matrix[dir].n_rows());
597 eigenvalues[dir].resize(mass_matrix[dir].n_cols());
598 internal::TensorProductMatrixSymmetricSum::spectral_assembly<Number>(
599 &(mass_matrix[dir](0, 0)),
600 &(derivative_matrix[dir](0, 0)),
601 mass_matrix[dir].n_rows(),
602 mass_matrix[dir].n_cols(),
603 eigenvalues[dir].begin(),
604 &(eigenvectors[dir](0, 0)));
605 }
606 }
607
608
609
610 template <std::size_t dim, typename Number, std::size_t n_lanes>
611 inline void
612 setup(
613 const std::array<Table<2, VectorizedArray<Number, n_lanes>>, dim>
614 &mass_matrix,
615 const std::array<Table<2, VectorizedArray<Number, n_lanes>>, dim>
616 &derivative_matrix,
620 {
621 const unsigned int n_rows_1d = mass_matrix[0].n_cols();
622 constexpr unsigned int macro_size =
624 const std::size_t nm_flat_size_max = n_rows_1d * n_rows_1d * macro_size;
625 const std::size_t n_flat_size_max = n_rows_1d * macro_size;
626
627 std::vector<Number> mass_matrix_flat;
628 std::vector<Number> deriv_matrix_flat;
629 std::vector<Number> eigenvalues_flat;
630 std::vector<Number> eigenvectors_flat;
631 mass_matrix_flat.resize(nm_flat_size_max);
632 deriv_matrix_flat.resize(nm_flat_size_max);
633 eigenvalues_flat.resize(n_flat_size_max);
634 eigenvectors_flat.resize(nm_flat_size_max);
635 std::array<unsigned int, macro_size> offsets_nm;
636 std::array<unsigned int, macro_size> offsets_n;
637 for (unsigned int dir = 0; dir < dim; ++dir)
638 {
639 AssertDimension(n_rows_1d, mass_matrix[dir].n_cols());
640 AssertDimension(mass_matrix[dir].n_rows(), mass_matrix[dir].n_cols());
641 AssertDimension(mass_matrix[dir].n_rows(),
642 derivative_matrix[dir].n_rows());
643 AssertDimension(mass_matrix[dir].n_rows(),
644 derivative_matrix[dir].n_cols());
645
646 const unsigned int n_rows = mass_matrix[dir].n_rows();
647 const unsigned int n_cols = mass_matrix[dir].n_cols();
648 const unsigned int nm = n_rows * n_cols;
649 for (unsigned int vv = 0; vv < macro_size; ++vv)
650 offsets_nm[vv] = nm * vv;
651
653 nm,
654 &(mass_matrix[dir](0, 0)),
655 offsets_nm.cbegin(),
656 mass_matrix_flat.data());
658 nm,
659 &(derivative_matrix[dir](0, 0)),
660 offsets_nm.cbegin(),
661 deriv_matrix_flat.data());
662
663 const Number *mass_cbegin = mass_matrix_flat.data();
664 const Number *deriv_cbegin = deriv_matrix_flat.data();
665 Number *eigenvec_begin = eigenvectors_flat.data();
666 Number *eigenval_begin = eigenvalues_flat.data();
667 for (unsigned int lane = 0; lane < macro_size; ++lane)
668 internal::TensorProductMatrixSymmetricSum::spectral_assembly<
669 Number>(mass_cbegin + nm * lane,
670 deriv_cbegin + nm * lane,
671 n_rows,
672 n_cols,
673 eigenval_begin + n_rows * lane,
674 eigenvec_begin + nm * lane);
675
676 eigenvalues[dir].resize(n_rows);
677 eigenvectors[dir].reinit(n_rows, n_cols);
678 for (unsigned int vv = 0; vv < macro_size; ++vv)
679 offsets_n[vv] = n_rows * vv;
681 eigenvalues_flat.data(),
682 offsets_n.cbegin(),
683 eigenvalues[dir].begin());
685 eigenvectors_flat.data(),
686 offsets_nm.cbegin(),
687 &(eigenvectors[dir](0, 0)));
688 }
689 }
690
691
692
693 template <std::size_t dim, typename Number>
694 inline std::array<Table<2, Number>, dim>
695 convert(const std::array<Table<2, Number>, dim> &mass_matrix)
696 {
697 return mass_matrix;
698 }
699
700
701
702 template <std::size_t dim, typename Number>
703 inline std::array<Table<2, Number>, dim>
704 convert(const std::array<FullMatrix<Number>, dim> &mass_matrix)
705 {
706 std::array<Table<2, Number>, dim> mass_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
715 return mass_copy;
716 }
717
718
719
720 template <std::size_t dim, typename Number>
721 inline std::array<Table<2, Number>, dim>
722 convert(const Table<2, Number> &matrix)
723 {
724 std::array<Table<2, Number>, dim> matrices;
725
726 std::fill(matrices.begin(), matrices.end(), matrix);
727
728 return matrices;
729 }
730
731
732
733 template <int n_rows_1d_templated, std::size_t dim, typename Number>
734 void
735 vmult(Number *dst,
736 const Number *src,
738 const unsigned int n_rows_1d_non_templated,
739 const std::array<const Number *, dim> &mass_matrix,
740 const std::array<const Number *, dim> &derivative_matrix)
741 {
742 const unsigned int n_rows_1d = n_rows_1d_templated == 0 ?
743 n_rows_1d_non_templated :
744 n_rows_1d_templated;
745 const unsigned int n = Utilities::fixed_power<dim>(n_rows_1d);
746
747 tmp.resize_fast(n * 2);
748 Number *t = tmp.begin();
749
751 dim,
752 n_rows_1d_templated,
753 n_rows_1d_templated,
754 Number>
755 eval({}, {}, {}, n_rows_1d, n_rows_1d);
756
757 if (dim == 1)
758 {
759 const Number *A = derivative_matrix[0];
760 eval.template apply<0, false, false>(A, src, dst);
761 }
762
763 else if (dim == 2)
764 {
765 const Number *A0 = derivative_matrix[0];
766 const Number *M0 = mass_matrix[0];
767 const Number *A1 = derivative_matrix[1];
768 const Number *M1 = mass_matrix[1];
769 eval.template apply<0, false, false>(M0, src, t);
770 eval.template apply<1, false, false>(A1, t, dst);
771 eval.template apply<0, false, false>(A0, src, t);
772 eval.template apply<1, false, true>(M1, t, dst);
773 }
774
775 else if (dim == 3)
776 {
777 const Number *A0 = derivative_matrix[0];
778 const Number *M0 = mass_matrix[0];
779 const Number *A1 = derivative_matrix[1];
780 const Number *M1 = mass_matrix[1];
781 const Number *A2 = derivative_matrix[2];
782 const Number *M2 = mass_matrix[2];
783 eval.template apply<0, false, false>(M0, src, t + n);
784 eval.template apply<1, false, false>(M1, t + n, t);
785 eval.template apply<2, false, false>(A2, t, dst);
786 eval.template apply<1, false, false>(A1, t + n, t);
787 eval.template apply<0, false, false>(A0, src, t + n);
788 eval.template apply<1, false, true>(M1, t + n, t);
789 eval.template apply<2, false, true>(M2, t, dst);
790 }
791
792 else
794 }
795
796
797
798 template <int n_rows_1d_templated, std::size_t dim, typename Number>
799 void
800 apply_inverse(Number *dst,
801 const Number *src,
802 const unsigned int n_rows_1d_non_templated,
803 const std::array<const Number *, dim> &eigenvectors,
804 const std::array<const Number *, dim> &eigenvalues,
805 const Number *inverted_eigenvalues = nullptr)
806 {
807 const unsigned int n_rows_1d = n_rows_1d_templated == 0 ?
808 n_rows_1d_non_templated :
809 n_rows_1d_templated;
810
812 dim,
813 n_rows_1d_templated,
814 n_rows_1d_templated,
815 Number>
816 eval({}, {}, {}, n_rows_1d, n_rows_1d);
817
818 // NOTE: dof_to_quad has to be interpreted as 'dof to eigenvalue index'
819 // --> apply<.,true,.> (S,src,dst) calculates dst = S^T * src,
820 // --> apply<.,false,.> (S,src,dst) calculates dst = S * src,
821 // while the eigenvectors are stored column-wise in S, i.e.
822 // rows correspond to dofs whereas columns to eigenvalue indices!
823 if (dim == 1)
824 {
825 const Number *S = eigenvectors[0];
826 eval.template apply<0, true, false>(S, src, dst);
827
828 for (unsigned int i = 0; i < n_rows_1d; ++i)
829 if (inverted_eigenvalues)
830 dst[i] *= inverted_eigenvalues[i];
831 else
832 dst[i] /= eigenvalues[0][i];
833
834 eval.template apply<0, false, false>(S, dst, dst);
835 }
836
837 else if (dim == 2)
838 {
839 const Number *S0 = eigenvectors[0];
840 const Number *S1 = eigenvectors[1];
841 eval.template apply<0, true, false>(S0, src, dst);
842 eval.template apply<1, true, false>(S1, dst, dst);
843
844 for (unsigned int i1 = 0, c = 0; i1 < n_rows_1d; ++i1)
845 for (unsigned int i0 = 0; i0 < n_rows_1d; ++i0, ++c)
846 if (inverted_eigenvalues)
847 dst[c] *= inverted_eigenvalues[c];
848 else
849 dst[c] /= (eigenvalues[1][i1] + eigenvalues[0][i0]);
850
851 eval.template apply<1, false, false>(S1, dst, dst);
852 eval.template apply<0, false, false>(S0, dst, dst);
853 }
854
855 else if (dim == 3)
856 {
857 const Number *S0 = eigenvectors[0];
858 const Number *S1 = eigenvectors[1];
859 const Number *S2 = eigenvectors[2];
860 eval.template apply<0, true, false>(S0, src, dst);
861 eval.template apply<1, true, false>(S1, dst, dst);
862 eval.template apply<2, true, false>(S2, dst, dst);
863
864 for (unsigned int i2 = 0, c = 0; i2 < n_rows_1d; ++i2)
865 for (unsigned int i1 = 0; i1 < n_rows_1d; ++i1)
866 for (unsigned int i0 = 0; i0 < n_rows_1d; ++i0, ++c)
867 if (inverted_eigenvalues)
868 dst[c] *= inverted_eigenvalues[c];
869 else
870 dst[c] /= (eigenvalues[2][i2] + eigenvalues[1][i1] +
871 eigenvalues[0][i0]);
872
873 eval.template apply<2, false, false>(S2, dst, dst);
874 eval.template apply<1, false, false>(S1, dst, dst);
875 eval.template apply<0, false, false>(S0, dst, dst);
876 }
877
878 else
880 }
881
882
883
884 template <int n_rows_1d_templated, std::size_t dim, typename Number>
885 void
886 select_vmult(Number *dst,
887 const Number *src,
889 const unsigned int n_rows_1d,
890 const std::array<const Number *, dim> &mass_matrix,
891 const std::array<const Number *, dim> &derivative_matrix);
892
893
894
895 template <int n_rows_1d_templated, std::size_t dim, typename Number>
896 void
897 select_apply_inverse(Number *dst,
898 const Number *src,
899 const unsigned int n_rows_1d,
900 const std::array<const Number *, dim> &eigenvectors,
901 const std::array<const Number *, dim> &eigenvalues,
902 const Number *inverted_eigenvalues = nullptr);
903 } // namespace TensorProductMatrixSymmetricSum
904} // namespace internal
905
906
907template <int dim, typename Number, int n_rows_1d>
908inline unsigned int
910{
911 unsigned int m = mass_matrix[0].n_rows();
912 for (unsigned int d = 1; d < dim; ++d)
913 m *= mass_matrix[d].n_rows();
914 return m;
915}
916
917
918
919template <int dim, typename Number, int n_rows_1d>
920inline unsigned int
922{
923 unsigned int n = mass_matrix[0].n_cols();
924 for (unsigned int d = 1; d < dim; ++d)
925 n *= mass_matrix[d].n_cols();
926 return n;
927}
928
929
930
931template <int dim, typename Number, int n_rows_1d>
932inline void
934 const ArrayView<Number> &dst_view,
935 const ArrayView<const Number> &src_view) const
936{
937 std::lock_guard<std::mutex> lock(this->mutex);
938 this->vmult(dst_view, src_view, this->tmp_array);
939}
940
941
942
943template <int dim, typename Number, int n_rows_1d>
944inline void
946 const ArrayView<Number> &dst_view,
947 const ArrayView<const Number> &src_view,
948 AlignedVector<Number> &tmp_array) const
949{
950 AssertDimension(dst_view.size(), this->m());
951 AssertDimension(src_view.size(), this->n());
952
953 Number *dst = dst_view.begin();
954 const Number *src = src_view.begin();
955
956 std::array<const Number *, dim> mass_matrix, derivative_matrix;
957
958 for (unsigned int d = 0; d < dim; ++d)
959 {
960 mass_matrix[d] = &this->mass_matrix[d](0, 0);
961 derivative_matrix[d] = &this->derivative_matrix[d](0, 0);
962 }
963
964 const unsigned int n_rows_1d_non_templated = this->mass_matrix[0].n_rows();
965
966 if (n_rows_1d != -1)
967 internal::TensorProductMatrixSymmetricSum::vmult<
968 n_rows_1d == -1 ? 0 : n_rows_1d>(dst,
969 src,
970 tmp_array,
971 n_rows_1d_non_templated,
973 derivative_matrix);
974 else
975 internal::TensorProductMatrixSymmetricSum::select_vmult<1>(
976 dst,
977 src,
978 tmp_array,
979 n_rows_1d_non_templated,
980 mass_matrix,
981 derivative_matrix);
982}
983
984
985
986template <int dim, typename Number, int n_rows_1d>
987inline void
989 const ArrayView<Number> &dst_view,
990 const ArrayView<const Number> &src_view) const
991{
992 AssertDimension(dst_view.size(), this->n());
993 AssertDimension(src_view.size(), this->m());
994
995 Number *dst = dst_view.begin();
996 const Number *src = src_view.begin();
997
998 std::array<const Number *, dim> eigenvectors, eigenvalues;
999
1000 for (unsigned int d = 0; d < dim; ++d)
1001 {
1002 eigenvectors[d] = &this->eigenvectors[d](0, 0);
1003 eigenvalues[d] = this->eigenvalues[d].data();
1004 }
1005
1006 const unsigned int n_rows_1d_non_templated = this->mass_matrix[0].n_rows();
1007
1008 if (n_rows_1d != -1)
1009 internal::TensorProductMatrixSymmetricSum::apply_inverse<
1010 n_rows_1d == -1 ? 0 : n_rows_1d>(
1011 dst, src, n_rows_1d_non_templated, eigenvectors, eigenvalues);
1012 else
1013 internal::TensorProductMatrixSymmetricSum::select_apply_inverse<1>(
1014 dst, src, n_rows_1d_non_templated, eigenvectors, eigenvalues);
1015}
1016
1017
1018
1019template <int dim, typename Number, int n_rows_1d>
1020std::size_t
1022 const
1023{
1024 return MemoryConsumption::memory_consumption(mass_matrix) +
1025 MemoryConsumption::memory_consumption(derivative_matrix) +
1029}
1030
1031
1032
1033template <int dim, typename Number, int n_rows_1d>
1034template <typename T>
1036 TensorProductMatrixSymmetricSum(const T &mass_matrix,
1037 const T &derivative_matrix)
1038{
1039 reinit(mass_matrix, derivative_matrix);
1040}
1041
1042
1043
1044template <int dim, typename Number, int n_rows_1d>
1045template <typename T>
1046inline void
1048 const T &mass_matrix,
1049 const T &derivative_matrix)
1050{
1051 this->mass_matrix =
1052 internal::TensorProductMatrixSymmetricSum::convert<dim>(mass_matrix);
1053 this->derivative_matrix =
1054 internal::TensorProductMatrixSymmetricSum::convert<dim>(derivative_matrix);
1055
1056 internal::TensorProductMatrixSymmetricSum::setup(this->mass_matrix,
1057 this->derivative_matrix,
1058 this->eigenvectors,
1059 this->eigenvalues);
1060}
1061
1062
1063
1064template <int dim, typename Number, int n_rows_1d>
1066 AdditionalData::AdditionalData(const bool compress_matrices,
1067 const bool precompute_inverse_diagonal)
1068 : compress_matrices(compress_matrices)
1069 , precompute_inverse_diagonal(precompute_inverse_diagonal)
1070{}
1071
1072
1073
1074template <int dim, typename Number, int n_rows_1d>
1077 const AdditionalData &additional_data)
1078 : compress_matrices(additional_data.compress_matrices)
1079 , precompute_inverse_diagonal(additional_data.precompute_inverse_diagonal)
1080{}
1081
1082
1083
1084template <int dim, typename Number, int n_rows_1d>
1085void
1087 const unsigned int size)
1088{
1089 if (compress_matrices == false)
1090 mass_and_derivative_matrices.resize(size * dim);
1091 else
1092 indices.assign(size * dim, numbers::invalid_unsigned_int);
1093}
1094
1095
1096
1097template <int dim, typename Number, int n_rows_1d>
1098template <typename T>
1099void
1101 const unsigned int index,
1102 const T &Ms_in,
1103 const T &Ks_in)
1104{
1105 const auto Ms =
1106 internal::TensorProductMatrixSymmetricSum::convert<dim>(Ms_in);
1107 const auto Ks =
1108 internal::TensorProductMatrixSymmetricSum::convert<dim>(Ks_in);
1109
1110 for (unsigned int d = 0; d < dim; ++d)
1111 {
1112 if (compress_matrices == false)
1113 {
1114 const MatrixPairType matrix(Ms[d], Ks[d]);
1115 mass_and_derivative_matrices[index * dim + d] = matrix;
1116 }
1117 else
1118 {
1119 using VectorizedArrayTrait =
1121
1122 std::bitset<VectorizedArrayTrait::width()> mask;
1123
1124 for (unsigned int v = 0; v < VectorizedArrayTrait::width(); ++v)
1125 {
1126 typename VectorizedArrayTrait::value_type a = 0.0;
1127
1128 for (unsigned int i = 0; i < Ms[d].size(0); ++i)
1129 for (unsigned int j = 0; j < Ms[d].size(1); ++j)
1130 {
1131 a += std::abs(VectorizedArrayTrait::get(Ms[d][i][j], v));
1132 a += std::abs(VectorizedArrayTrait::get(Ks[d][i][j], v));
1133 }
1134
1135 mask[v] = (a != 0.0);
1136 }
1137
1138 const MatrixPairTypeWithMask matrix{mask, {Ms[d], Ks[d]}};
1139
1140 const auto ptr = cache.find(matrix);
1141
1142 if (ptr != cache.end())
1143 {
1144 const auto ptr_index = ptr->second;
1145 indices[index * dim + d] = ptr_index;
1146
1147 if ([&]() {
1148 for (unsigned int v = 0; v < VectorizedArrayTrait::width();
1149 ++v)
1150 if ((mask[v] == true) && (ptr->first.first[v] == false))
1151 return false;
1152
1153 return true;
1154 }())
1155 {
1156 // nothing to do
1157 }
1158 else
1159 {
1160 auto mask_new = ptr->first.first;
1161 auto Ms_new = ptr->first.second.first;
1162 auto Ks_new = ptr->first.second.second;
1163
1164 for (unsigned int v = 0; v < VectorizedArrayTrait::width();
1165 ++v)
1166 if (mask_new[v] == false && mask[v] == true)
1167 {
1168 mask_new[v] = true;
1169
1170 for (unsigned int i = 0; i < Ms_new.size(0); ++i)
1171 for (unsigned int j = 0; j < Ms_new.size(1); ++j)
1172 {
1173 VectorizedArrayTrait::get(Ms_new[i][j], v) =
1174 VectorizedArrayTrait::get(Ms[d][i][j], v);
1175 VectorizedArrayTrait::get(Ks_new[i][j], v) =
1176 VectorizedArrayTrait::get(Ks[d][i][j], v);
1177 }
1178 }
1179
1180 cache.erase(ptr);
1181
1182 const MatrixPairTypeWithMask entry_new{mask_new,
1183 {Ms_new, Ks_new}};
1184
1185 const auto ptr_ = cache.find(entry_new);
1186 AssertThrow(ptr_ == cache.end(), ExcNotImplemented());
1187
1188 cache[entry_new] = ptr_index;
1189 }
1190 }
1191 else
1192 {
1193 const auto size = cache.size();
1194 indices[index * dim + d] = size;
1195 cache[matrix] = size;
1196 }
1197 }
1198 }
1199}
1200
1201
1202
1203template <int dim, typename Number, int n_rows_1d>
1204void
1206{
1207 const auto store = [&](const unsigned int index,
1208 const MatrixPairType &M_and_K) {
1209 std::array<Table<2, Number>, 1> mass_matrix;
1210 mass_matrix[0] = M_and_K.first;
1211
1212 std::array<Table<2, Number>, 1> derivative_matrix;
1213 derivative_matrix[0] = M_and_K.second;
1214
1215 std::array<Table<2, Number>, 1> eigenvectors;
1216 std::array<AlignedVector<Number>, 1> eigenvalues;
1217
1218 internal::TensorProductMatrixSymmetricSum::setup(mass_matrix,
1219 derivative_matrix,
1221 eigenvalues);
1222
1223 for (unsigned int i = 0, m = matrix_ptr[index], v = vector_ptr[index];
1224 i < mass_matrix[0].n_rows();
1225 ++i, ++v)
1226 {
1227 for (unsigned int j = 0; j < mass_matrix[0].n_cols(); ++j, ++m)
1228 {
1229 this->mass_matrices[m] = mass_matrix[0][i][j];
1230 this->derivative_matrices[m] = derivative_matrix[0][i][j];
1231 this->eigenvectors[m] = eigenvectors[0][i][j];
1232 }
1233
1234 this->eigenvalues[v] = eigenvalues[0][i];
1235 }
1236 };
1237
1238 if (compress_matrices == false)
1239 {
1240 // case 1) no compression requested
1241
1242 AssertDimension(cache.size(), 0);
1243 AssertDimension(indices.size(), 0);
1244
1245 this->vector_ptr.resize(mass_and_derivative_matrices.size() + 1);
1246 this->matrix_ptr.resize(mass_and_derivative_matrices.size() + 1);
1247
1248 for (unsigned int i = 0; i < mass_and_derivative_matrices.size(); ++i)
1249 {
1250 const auto &M = mass_and_derivative_matrices[i].first;
1251
1252 this->vector_ptr[i + 1] = M.n_rows();
1253 this->matrix_ptr[i + 1] = M.n_rows() * M.n_cols();
1254 }
1255
1256 for (unsigned int i = 0; i < mass_and_derivative_matrices.size(); ++i)
1257 {
1258 this->vector_ptr[i + 1] += this->vector_ptr[i];
1259 this->matrix_ptr[i + 1] += this->matrix_ptr[i];
1260 }
1261
1262 this->mass_matrices.resize_fast(matrix_ptr.back());
1263 this->derivative_matrices.resize_fast(matrix_ptr.back());
1264 this->eigenvectors.resize_fast(matrix_ptr.back());
1265 this->eigenvalues.resize_fast(vector_ptr.back());
1266
1267 for (unsigned int i = 0; i < mass_and_derivative_matrices.size(); ++i)
1268 store(i, mass_and_derivative_matrices[i]);
1269
1270 mass_and_derivative_matrices.clear();
1271 }
1272 else if (cache.size() == indices.size())
1273 {
1274 // case 2) compression requested but none possible
1275
1276 this->vector_ptr.resize(cache.size() + 1);
1277 this->matrix_ptr.resize(cache.size() + 1);
1278
1279 std::map<unsigned int, MatrixPairType> inverted_cache;
1280
1281 for (const auto &i : cache)
1282 inverted_cache[i.second] = i.first.second;
1283
1284 for (unsigned int i = 0; i < indices.size(); ++i)
1285 {
1286 const auto &M = inverted_cache[indices[i]].first;
1287
1288 this->vector_ptr[i + 1] = M.n_rows();
1289 this->matrix_ptr[i + 1] = M.n_rows() * M.n_cols();
1290 }
1291
1292 for (unsigned int i = 0; i < cache.size(); ++i)
1293 {
1294 this->vector_ptr[i + 1] += this->vector_ptr[i];
1295 this->matrix_ptr[i + 1] += this->matrix_ptr[i];
1296 }
1297
1298 this->mass_matrices.resize_fast(matrix_ptr.back());
1299 this->derivative_matrices.resize_fast(matrix_ptr.back());
1300 this->eigenvectors.resize_fast(matrix_ptr.back());
1301 this->eigenvalues.resize_fast(vector_ptr.back());
1302
1303 for (unsigned int i = 0; i < indices.size(); ++i)
1304 store(i, inverted_cache[indices[i]]);
1305
1306 indices.clear();
1307 cache.clear();
1308 }
1309 else
1310 {
1311 // case 3) compress
1312
1313 this->vector_ptr.resize(cache.size() + 1);
1314 this->matrix_ptr.resize(cache.size() + 1);
1315
1316 for (const auto &i : cache)
1317 {
1318 const auto &M = i.first.second.first;
1319
1320 this->vector_ptr[i.second + 1] = M.n_rows();
1321 this->matrix_ptr[i.second + 1] = M.n_rows() * M.n_cols();
1322 }
1323
1324 for (unsigned int i = 0; i < cache.size(); ++i)
1325 {
1326 this->vector_ptr[i + 1] += this->vector_ptr[i];
1327 this->matrix_ptr[i + 1] += this->matrix_ptr[i];
1328 }
1329
1330 this->mass_matrices.resize_fast(matrix_ptr.back());
1331 this->derivative_matrices.resize_fast(matrix_ptr.back());
1332 this->eigenvectors.resize_fast(matrix_ptr.back());
1333 this->eigenvalues.resize_fast(vector_ptr.back());
1334
1335 for (const auto &i : cache)
1336 store(i.second, i.first.second);
1337
1338 cache.clear();
1339 }
1340
1341 if (precompute_inverse_diagonal)
1342 {
1343 if (dim == 1)
1344 {
1345 // 1D case: simply invert 1D eigenvalues
1346 for (unsigned int i = 0; i < this->eigenvalues.size(); ++i)
1347 this->eigenvalues[i] = Number(1.0) / this->eigenvalues[i];
1348 std::swap(this->inverted_eigenvalues, eigenvalues);
1349 }
1350 else
1351 {
1352 // 2D and 3D case: we have 2 or 3 1d eigenvalues so that we
1353 // need to combine these
1354
1355 // step 1) if eigenvalues/eigenvectors are compressed, we
1356 // need to compress the diagonal (the combination of ev
1357 // indices) as well. This is an optional step.
1358 std::vector<unsigned int> indices_ev;
1359
1360 if (indices.size() > 0)
1361 {
1362 // 1a) create cache (ev indics -> diag index)
1363 const unsigned int n_cells = indices.size() / dim;
1364 std::map<std::array<unsigned int, dim>, unsigned int> cache_ev;
1365 std::vector<unsigned int> cache_ev_idx(n_cells);
1366
1367 for (unsigned int i = 0, c = 0; i < n_cells; ++i)
1368 {
1369 std::array<unsigned int, dim> id;
1370
1371 for (unsigned int d = 0; d < dim; ++d, ++c)
1372 id[d] = indices[c];
1373
1374 const auto id_ptr = cache_ev.find(id);
1375
1376 if (id_ptr == cache_ev.end())
1377 {
1378 const auto size = cache_ev.size();
1379 cache_ev_idx[i] = size;
1380 cache_ev[id] = size;
1381 }
1382 else
1383 {
1384 cache_ev_idx[i] = id_ptr->second;
1385 }
1386 }
1387
1388 // 1b) store diagonal indices for each cell
1389 std::vector<unsigned int> new_indices;
1390 new_indices.reserve(indices.size() / dim * (dim + 1));
1391
1392 for (unsigned int i = 0, c = 0; i < n_cells; ++i)
1393 {
1394 for (unsigned int d = 0; d < dim; ++d, ++c)
1395 new_indices.push_back(indices[c]);
1396 new_indices.push_back(cache_ev_idx[i]);
1397 }
1398
1399 // 1c) transpose cache (diag index -> ev indices)
1400 indices_ev.resize(cache_ev.size() * dim);
1401 for (const auto &entry : cache_ev)
1402 for (unsigned int d = 0; d < dim; ++d)
1403 indices_ev[entry.second * dim + d] = entry.first[d];
1404
1405 std::swap(this->indices, new_indices);
1406 }
1407
1408 // step 2) allocate memory and set pointers
1409 const unsigned int n_diag =
1410 ((indices_ev.size() > 0) ? indices_ev.size() :
1411 (matrix_ptr.size() - 1)) /
1412 dim;
1413
1414 std::vector<unsigned int> new_vector_ptr(n_diag + 1, 0);
1415 std::vector<unsigned int> new_vector_n_rows_1d(n_diag, 0);
1416
1417 for (unsigned int i = 0; i < n_diag; ++i)
1418 {
1419 const unsigned int c = (indices_ev.size() > 0) ?
1420 indices_ev[dim * i + 0] :
1421 (dim * i + 0);
1422
1423 const unsigned int n_rows = vector_ptr[c + 1] - vector_ptr[c];
1424
1425 new_vector_n_rows_1d[i] = n_rows;
1426 new_vector_ptr[i + 1] = Utilities::pow(n_rows, dim);
1427 }
1428
1429 for (unsigned int i = 0; i < n_diag; ++i)
1430 new_vector_ptr[i + 1] += new_vector_ptr[i];
1431
1432 this->inverted_eigenvalues.resize(new_vector_ptr.back());
1433
1434 // step 3) loop over all unique diagonal entries and invert
1435 for (unsigned int i = 0; i < n_diag; ++i)
1436 {
1437 std::array<Number *, dim> evs;
1438
1439 for (unsigned int d = 0; d < dim; ++d)
1440 evs[d] =
1441 &this
1442 ->eigenvalues[this->vector_ptr[(indices_ev.size() > 0) ?
1443 indices_ev[dim * i + d] :
1444 (dim * i + d)]];
1445
1446 const unsigned int mm = new_vector_n_rows_1d[i];
1447 if (dim == 2)
1448 {
1449 for (unsigned int i1 = 0, c = 0; i1 < mm; ++i1)
1450 for (unsigned int i0 = 0; i0 < mm; ++i0, ++c)
1451 this->inverted_eigenvalues[new_vector_ptr[i] + c] =
1452 Number(1.0) / (evs[1][i1] + evs[0][i0]);
1453 }
1454 else
1455 {
1456 for (unsigned int i2 = 0, c = 0; i2 < mm; ++i2)
1457 for (unsigned int i1 = 0; i1 < mm; ++i1)
1458 for (unsigned int i0 = 0; i0 < mm; ++i0, ++c)
1459 this->inverted_eigenvalues[new_vector_ptr[i] + c] =
1460 Number(1.0) / (evs[2][i2] + evs[1][i1] + evs[0][i0]);
1461 }
1462 }
1463
1464 // step 4) clean up
1465 std::swap(this->vector_ptr, new_vector_ptr);
1466 std::swap(this->vector_n_rows_1d, new_vector_n_rows_1d);
1467 }
1468
1469 this->eigenvalues.clear();
1470 }
1471}
1472
1473
1474
1475template <int dim, typename Number, int n_rows_1d>
1476void
1478 apply_inverse(const unsigned int index,
1479 const ArrayView<Number> &dst_in,
1480 const ArrayView<const Number> &src_in) const
1481{
1482 Number *dst = dst_in.begin();
1483 const Number *src = src_in.begin();
1484
1485 if (this->eigenvalues.empty() == false)
1486 {
1487 std::array<const Number *, dim> eigenvectors;
1488 std::array<const Number *, dim> eigenvalues;
1489 unsigned int n_rows_1d_non_templated = 0;
1490
1491 for (unsigned int d = 0; d < dim; ++d)
1492 {
1493 const unsigned int translated_index =
1494 (indices.size() > 0) ? indices[dim * index + d] : (dim * index + d);
1495
1496 eigenvectors[d] =
1497 this->eigenvectors.data() + matrix_ptr[translated_index];
1498 eigenvalues[d] =
1499 this->eigenvalues.data() + vector_ptr[translated_index];
1500 n_rows_1d_non_templated =
1501 vector_ptr[translated_index + 1] - vector_ptr[translated_index];
1502 }
1503
1504 if (n_rows_1d != -1)
1505 internal::TensorProductMatrixSymmetricSum::apply_inverse<
1506 n_rows_1d == -1 ? 0 : n_rows_1d>(
1507 dst, src, n_rows_1d_non_templated, eigenvectors, eigenvalues);
1508 else
1509 internal::TensorProductMatrixSymmetricSum::select_apply_inverse<1>(
1510 dst, src, n_rows_1d_non_templated, eigenvectors, eigenvalues);
1511 }
1512 else
1513 {
1514 std::array<const Number *, dim> eigenvectors;
1515 const Number *inverted_eigenvalues = nullptr;
1516 unsigned int n_rows_1d_non_templated = 0;
1517
1518 for (unsigned int d = 0; d < dim; ++d)
1519 {
1520 const unsigned int translated_index =
1521 (indices.size() > 0) ?
1522 indices[((dim == 1) ? 1 : (dim + 1)) * index + d] :
1523 (dim * index + d);
1524
1525 eigenvectors[d] =
1526 this->eigenvectors.data() + matrix_ptr[translated_index];
1527 }
1528
1529 {
1530 const unsigned int translated_index =
1531 ((indices.size() > 0) && (dim != 1)) ?
1532 indices[(dim + 1) * index + dim] :
1533 index;
1534
1535 inverted_eigenvalues =
1536 this->inverted_eigenvalues.data() + vector_ptr[translated_index];
1537 n_rows_1d_non_templated =
1538 (dim == 1) ?
1539 (vector_ptr[translated_index + 1] - vector_ptr[translated_index]) :
1540 vector_n_rows_1d[translated_index];
1541 }
1542
1543 if (n_rows_1d != -1)
1544 internal::TensorProductMatrixSymmetricSum::apply_inverse<
1545 n_rows_1d == -1 ? 0 : n_rows_1d>(dst,
1546 src,
1547 n_rows_1d_non_templated,
1549 {},
1550 inverted_eigenvalues);
1551 else
1552 internal::TensorProductMatrixSymmetricSum::select_apply_inverse<1>(
1553 dst,
1554 src,
1555 n_rows_1d_non_templated,
1557 {},
1558 inverted_eigenvalues);
1559 }
1560}
1561
1562
1563
1564template <int dim, typename Number, int n_rows_1d>
1565std::size_t
1567 memory_consumption() const
1568{
1571 MemoryConsumption::memory_consumption(derivative_matrices) +
1576}
1577
1578
1579
1580template <int dim, typename Number, int n_rows_1d>
1581std::size_t
1583 storage_size() const
1584{
1585 if (matrix_ptr.empty())
1586 return 0; // if not initialized
1587
1588 return matrix_ptr.size() - 1;
1589}
1590
1591
1592
1593#endif
1594
1596
1597#endif
void resize_fast(const size_type new_size)
iterator begin()
iterator begin() const
Definition array_view.h:702
std::size_t size() const
Definition array_view.h:684
std::complex< typename numbers::NumberTraits< number >::real_type > eigenvalue(const size_type i) 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)
void apply_inverse(const unsigned int index, const ArrayView< Number > &dst_in, const ArrayView< const Number > &src_in) const
std::vector< MatrixPairType > mass_and_derivative_matrices
std::pair< std::bitset<::internal::VectorizedArrayTrait< Number >::width()>, MatrixPairType > MatrixPairTypeWithMask
std::pair< Table< 2, Number >, Table< 2, Number > > MatrixPairType
void reserve(const unsigned int size)
void insert(const unsigned int index, const T &Ms, const T &Ks)
std::map< MatrixPairTypeWithMask, unsigned int, internal::TensorProductMatrixSymmetricSum::MatrixPairComparator< Number > > cache
TensorProductMatrixSymmetricSumCollection(const AdditionalData &additional_data=AdditionalData())
void vmult(const ArrayView< Number > &dst, const ArrayView< const Number > &src, AlignedVector< Number > &tmp) const
std::array< Table< 2, Number >, dim > eigenvectors
std::array< Table< 2, Number >, dim > derivative_matrix
void reinit(const T &mass_matrix, const T &derivative_matrix)
void apply_inverse(const ArrayView< Number > &dst, const ArrayView< const Number > &src) const
void vmult(const ArrayView< Number > &dst, const ArrayView< const Number > &src) const
std::size_t memory_consumption() const
std::array< Table< 2, Number >, dim > mass_matrix
std::array< AlignedVector< Number >, dim > eigenvalues
TensorProductMatrixSymmetricSum(const T &mass_matrix, const T &derivative_matrix)
constexpr void clear()
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
Point< 2 > second
Definition grid_out.cc:4624
Point< 2 > first
Definition grid_out.cc:4623
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcInternalError()
#define AssertThrow(cond, exc)
#define DEAL_II_NOT_IMPLEMENTED()
for(unsigned int j=best_vertex+1;j< vertices.size();++j) if(vertices_to_use[j])
@ matrix
Contents is actually a matrix.
void mass_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const double factor=1.)
Definition l2.h:57
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
constexpr T pow(const T base, const int iexp)
Definition utilities.h:966
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition tria.cc:14951
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
static const unsigned int invalid_unsigned_int
Definition types.h:220
STL namespace.
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
AdditionalData(const bool compress_matrices=true, const bool precompute_inverse_diagonal=true)
bool operator()(const MatrixPairType &left, const MatrixPairType &right) const
std::pair< std::bitset< width >, std::pair< Table< 2, Number >, Table< 2, Number > > > MatrixPairType
static constexpr std::size_t width()
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)
void vectorized_load_and_transpose(const unsigned int n_entries, const Number *in, const unsigned int *offsets, VectorizedArray< Number, width > *out)
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)