deal.II version GIT relicensing-2291-g668cd86249 2024-12-24 11:30:00+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
21#include <deal.II/base/array_view.h>
22#include <deal.II/base/floating_point_comparator.h>
23#include <deal.II/base/mutex.h>
24#include <deal.II/base/vectorization.h>
25
26#include <deal.II/lac/lapack_full_matrix.h>
27
28#include <deal.II/matrix_free/tensor_product_kernels.h>
29
30DEAL_II_NAMESPACE_OPEN
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>
112class TensorProductMatrixSymmetricSum
113{
114public:
119 using value_type = Number;
120
125 static constexpr int n_rows_1d_static = n_rows_1d;
126
130 TensorProductMatrixSymmetricSum() = default;
131
136 template <typename T>
137 TensorProductMatrixSymmetricSum(const T &mass_matrix,
138 const T &derivative_matrix);
139
157 template <typename T>
158 void
159 reinit(const T &mass_matrix, const T &derivative_matrix);
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
199 vmult(const ArrayView<Number> &dst,
200 const ArrayView<const Number> &src,
201 AlignedVector<Number> &tmp) const;
202
209 void
210 apply_inverse(const ArrayView<Number> &dst,
211 const ArrayView<const Number> &src) const;
212
216 std::size_t
217 memory_consumption() const;
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:
246 mutable AlignedVector<Number> tmp_array;
247
251 mutable Threads::Mutex mutex;
252};
253
254
255
256namespace internal
257{
258 namespace TensorProductMatrixSymmetricSum
259 {
260 template <typename Number>
261 struct MatrixPairComparator
262 {
263 using VectorizedArrayTrait =
264 ::internal::VectorizedArrayTrait<Number>;
265 using ScalarNumber = typename VectorizedArrayTrait::value_type;
266 static constexpr std::size_t width = VectorizedArrayTrait::width();
267
268 using MatrixPairType =
269 std::pair<std::bitset<width>,
270 std::pair<Table<2, Number>, Table<2, Number>>>;
271
272 MatrixPairComparator()
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:
303 const ScalarNumber eps;
304 };
305 } // namespace TensorProductMatrixSymmetricSum
306} // namespace internal
307
308
309
333template <int dim, typename Number, int n_rows_1d = -1>
334class TensorProductMatrixSymmetricSumCollection
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()>,
340 MatrixPairType>;
341
342public:
346 struct AdditionalData
347 {
351 AdditionalData(const bool compress_matrices = true,
352 const bool precompute_inverse_diagonal = true);
353
357 bool compress_matrices;
358
362 bool precompute_inverse_diagonal;
363 };
364
368 TensorProductMatrixSymmetricSumCollection(
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
392 finalize();
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
406 memory_consumption() const;
407
416 std::size_t
417 storage_size() const;
418
419private:
423 const bool compress_matrices;
424
428 const bool precompute_inverse_diagonal;
429
434 std::vector<MatrixPairType> mass_and_derivative_matrices;
435
440 std::map<
441 MatrixPairTypeWithMask,
442 unsigned int,
443 internal::TensorProductMatrixSymmetricSum::MatrixPairComparator<Number>>
444 cache;
445
453 std::vector<unsigned int> indices;
454
458 AlignedVector<Number> mass_matrices;
459
463 AlignedVector<Number> derivative_matrices;
464
468 AlignedVector<Number> eigenvectors;
469
473 AlignedVector<Number> eigenvalues;
474
478 AlignedVector<Number> inverted_eigenvalues;
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{
503 namespace TensorProductMatrixSymmetricSum
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,
531 ExcInternalError());
532
533 for (unsigned int j = 0; j < n_rows; ++j)
534 {
535 Assert(derivative_matrix[i + j * n_rows] == 0,
536 ExcInternalError());
537 Assert(derivative_matrix[j + i * n_rows] == 0,
538 ExcInternalError());
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
562 deriv_copy.compute_generalized_eigenvalues_symmetric(mass_copy,
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,
617 std::array<Table<2, VectorizedArray<Number, n_lanes>>, dim> &eigenvectors,
618 std::array<AlignedVector<VectorizedArray<Number, n_lanes>>, dim>
619 &eigenvalues)
620 {
621 const unsigned int n_rows_1d = mass_matrix[0].n_cols();
622 constexpr unsigned int macro_size =
623 VectorizedArray<Number, n_lanes>::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
652 vectorized_transpose_and_store<Number, n_lanes>(
653 false,
654 nm,
655 &(mass_matrix[dir](0, 0)),
656 offsets_nm.data(),
657 mass_matrix_flat.data());
658 vectorized_transpose_and_store<Number, n_lanes>(
659 false,
660 nm,
661 &(derivative_matrix[dir](0, 0)),
662 offsets_nm.data(),
663 deriv_matrix_flat.data());
664
665 const Number *mass_cbegin = mass_matrix_flat.data();
666 const Number *deriv_cbegin = deriv_matrix_flat.data();
667 Number *eigenvec_begin = eigenvectors_flat.data();
668 Number *eigenval_begin = eigenvalues_flat.data();
669 for (unsigned int lane = 0; lane < macro_size; ++lane)
670 internal::TensorProductMatrixSymmetricSum::spectral_assembly<
671 Number>(mass_cbegin + nm * lane,
672 deriv_cbegin + nm * lane,
673 n_rows,
674 n_cols,
675 eigenval_begin + n_rows * lane,
676 eigenvec_begin + nm * lane);
677
678 eigenvalues[dir].resize(n_rows);
679 eigenvectors[dir].reinit(n_rows, n_cols);
680 for (unsigned int vv = 0; vv < macro_size; ++vv)
681 offsets_n[vv] = n_rows * vv;
682 vectorized_load_and_transpose<Number, n_lanes>(
683 n_rows,
684 eigenvalues_flat.data(),
685 offsets_n.data(),
686 eigenvalues[dir].begin());
687 vectorized_load_and_transpose<Number, n_lanes>(
688 nm,
689 eigenvectors_flat.data(),
690 offsets_nm.data(),
691 &(eigenvectors[dir](0, 0)));
692 }
693 }
694
695
696
697 template <std::size_t dim, typename Number>
698 inline std::array<Table<2, Number>, dim>
699 convert(const std::array<Table<2, Number>, dim> &mass_matrix)
700 {
701 return mass_matrix;
702 }
703
704
705
706 template <std::size_t dim, typename Number>
707 inline std::array<Table<2, Number>, dim>
708 convert(const std::array<FullMatrix<Number>, dim> &mass_matrix)
709 {
710 std::array<Table<2, Number>, dim> mass_copy;
711
712 std::transform(mass_matrix.cbegin(),
713 mass_matrix.cend(),
714 mass_copy.begin(),
715 [](const FullMatrix<Number> &m) -> Table<2, Number> {
716 return m;
717 });
718
719 return mass_copy;
720 }
721
722
723
724 template <std::size_t dim, typename Number>
725 inline std::array<Table<2, Number>, dim>
726 convert(const Table<2, Number> &matrix)
727 {
728 std::array<Table<2, Number>, dim> matrices;
729
730 std::fill(matrices.begin(), matrices.end(), matrix);
731
732 return matrices;
733 }
734
735
736
737 template <int n_rows_1d_templated, std::size_t dim, typename Number>
738 void
739 vmult(Number *dst,
740 const Number *src,
741 AlignedVector<Number> &tmp,
742 const unsigned int n_rows_1d_non_templated,
743 const std::array<const Number *, dim> &mass_matrix,
744 const std::array<const Number *, dim> &derivative_matrix)
745 {
746 const unsigned int n_rows_1d = n_rows_1d_templated == 0 ?
747 n_rows_1d_non_templated :
748 n_rows_1d_templated;
749 const unsigned int n = Utilities::fixed_power<dim>(n_rows_1d);
750
751 tmp.resize_fast(n * 2);
752 Number *t = tmp.begin();
753
754 internal::EvaluatorTensorProduct<internal::evaluate_general,
755 dim,
756 n_rows_1d_templated,
757 n_rows_1d_templated,
758 Number>
759 eval({}, {}, {}, n_rows_1d, n_rows_1d);
760
761 if (dim == 1)
762 {
763 const Number *A = derivative_matrix[0];
764 eval.template apply<0, false, false>(A, src, dst);
765 }
766
767 else if (dim == 2)
768 {
769 const Number *A0 = derivative_matrix[0];
770 const Number *M0 = mass_matrix[0];
771 const Number *A1 = derivative_matrix[1];
772 const Number *M1 = mass_matrix[1];
773 eval.template apply<0, false, false>(M0, src, t);
774 eval.template apply<1, false, false>(A1, t, dst);
775 eval.template apply<0, false, false>(A0, src, t);
776 eval.template apply<1, false, true>(M1, t, dst);
777 }
778
779 else if (dim == 3)
780 {
781 const Number *A0 = derivative_matrix[0];
782 const Number *M0 = mass_matrix[0];
783 const Number *A1 = derivative_matrix[1];
784 const Number *M1 = mass_matrix[1];
785 const Number *A2 = derivative_matrix[2];
786 const Number *M2 = mass_matrix[2];
787 eval.template apply<0, false, false>(M0, src, t + n);
788 eval.template apply<1, false, false>(M1, t + n, t);
789 eval.template apply<2, false, false>(A2, t, dst);
790 eval.template apply<1, false, false>(A1, t + n, t);
791 eval.template apply<0, false, false>(A0, src, t + n);
792 eval.template apply<1, false, true>(M1, t + n, t);
793 eval.template apply<2, false, true>(M2, t, dst);
794 }
795
796 else
797 AssertThrow(false, ExcNotImplemented());
798 }
799
800
801
802 template <int n_rows_1d_templated, std::size_t dim, typename Number>
803 void
804 apply_inverse(Number *dst,
805 const Number *src,
806 const unsigned int n_rows_1d_non_templated,
807 const std::array<const Number *, dim> &eigenvectors,
808 const std::array<const Number *, dim> &eigenvalues,
809 const Number *inverted_eigenvalues = nullptr)
810 {
811 const unsigned int n_rows_1d = n_rows_1d_templated == 0 ?
812 n_rows_1d_non_templated :
813 n_rows_1d_templated;
814
815 internal::EvaluatorTensorProduct<internal::evaluate_general,
816 dim,
817 n_rows_1d_templated,
818 n_rows_1d_templated,
819 Number>
820 eval({}, {}, {}, n_rows_1d, n_rows_1d);
821
822 // NOTE: dof_to_quad has to be interpreted as 'dof to eigenvalue index'
823 // --> apply<.,true,.> (S,src,dst) calculates dst = S^T * src,
824 // --> apply<.,false,.> (S,src,dst) calculates dst = S * src,
825 // while the eigenvectors are stored column-wise in S, i.e.
826 // rows correspond to dofs whereas columns to eigenvalue indices!
827 if (dim == 1)
828 {
829 const Number *S = eigenvectors[0];
830 eval.template apply<0, true, false>(S, src, dst);
831
832 for (unsigned int i = 0; i < n_rows_1d; ++i)
833 if (inverted_eigenvalues)
834 dst[i] *= inverted_eigenvalues[i];
835 else
836 dst[i] /= eigenvalues[0][i];
837
838 eval.template apply<0, false, false>(S, dst, dst);
839 }
840
841 else if (dim == 2)
842 {
843 const Number *S0 = eigenvectors[0];
844 const Number *S1 = eigenvectors[1];
845 eval.template apply<0, true, false>(S0, src, dst);
846 eval.template apply<1, true, false>(S1, dst, dst);
847
848 for (unsigned int i1 = 0, c = 0; i1 < n_rows_1d; ++i1)
849 for (unsigned int i0 = 0; i0 < n_rows_1d; ++i0, ++c)
850 if (inverted_eigenvalues)
851 dst[c] *= inverted_eigenvalues[c];
852 else
853 dst[c] /= (eigenvalues[1][i1] + eigenvalues[0][i0]);
854
855 eval.template apply<1, false, false>(S1, dst, dst);
856 eval.template apply<0, false, false>(S0, dst, dst);
857 }
858
859 else if (dim == 3)
860 {
861 const Number *S0 = eigenvectors[0];
862 const Number *S1 = eigenvectors[1];
863 const Number *S2 = eigenvectors[2];
864 eval.template apply<0, true, false>(S0, src, dst);
865 eval.template apply<1, true, false>(S1, dst, dst);
866 eval.template apply<2, true, false>(S2, dst, dst);
867
868 for (unsigned int i2 = 0, c = 0; i2 < n_rows_1d; ++i2)
869 for (unsigned int i1 = 0; i1 < n_rows_1d; ++i1)
870 for (unsigned int i0 = 0; i0 < n_rows_1d; ++i0, ++c)
871 if (inverted_eigenvalues)
872 dst[c] *= inverted_eigenvalues[c];
873 else
874 dst[c] /= (eigenvalues[2][i2] + eigenvalues[1][i1] +
875 eigenvalues[0][i0]);
876
877 eval.template apply<2, false, false>(S2, dst, dst);
878 eval.template apply<1, false, false>(S1, dst, dst);
879 eval.template apply<0, false, false>(S0, dst, dst);
880 }
881
882 else
883 DEAL_II_NOT_IMPLEMENTED();
884 }
885
886
887
888 template <int n_rows_1d_templated, std::size_t dim, typename Number>
889 void
890 select_vmult(Number *dst,
891 const Number *src,
892 AlignedVector<Number> &tmp,
893 const unsigned int n_rows_1d,
894 const std::array<const Number *, dim> &mass_matrix,
895 const std::array<const Number *, dim> &derivative_matrix);
896
897
898
899 template <int n_rows_1d_templated, std::size_t dim, typename Number>
900 void
901 select_apply_inverse(Number *dst,
902 const Number *src,
903 const unsigned int n_rows_1d,
904 const std::array<const Number *, dim> &eigenvectors,
905 const std::array<const Number *, dim> &eigenvalues,
906 const Number *inverted_eigenvalues = nullptr);
907 } // namespace TensorProductMatrixSymmetricSum
908} // namespace internal
909
910
911template <int dim, typename Number, int n_rows_1d>
912inline unsigned int
913TensorProductMatrixSymmetricSum<dim, Number, n_rows_1d>::m() const
914{
915 unsigned int m = mass_matrix[0].n_rows();
916 for (unsigned int d = 1; d < dim; ++d)
917 m *= mass_matrix[d].n_rows();
918 return m;
919}
920
921
922
923template <int dim, typename Number, int n_rows_1d>
924inline unsigned int
925TensorProductMatrixSymmetricSum<dim, Number, n_rows_1d>::n() const
926{
927 unsigned int n = mass_matrix[0].n_cols();
928 for (unsigned int d = 1; d < dim; ++d)
929 n *= mass_matrix[d].n_cols();
930 return n;
931}
932
933
934
935template <int dim, typename Number, int n_rows_1d>
936inline void
937TensorProductMatrixSymmetricSum<dim, Number, n_rows_1d>::vmult(
938 const ArrayView<Number> &dst_view,
939 const ArrayView<const Number> &src_view) const
940{
941 std::lock_guard<std::mutex> lock(this->mutex);
942 this->vmult(dst_view, src_view, this->tmp_array);
943}
944
945
946
947template <int dim, typename Number, int n_rows_1d>
948inline void
949TensorProductMatrixSymmetricSum<dim, Number, n_rows_1d>::vmult(
950 const ArrayView<Number> &dst_view,
951 const ArrayView<const Number> &src_view,
952 AlignedVector<Number> &tmp_array) const
953{
954 AssertDimension(dst_view.size(), this->m());
955 AssertDimension(src_view.size(), this->n());
956
957 Number *dst = dst_view.begin();
958 const Number *src = src_view.begin();
959
960 std::array<const Number *, dim> mass_matrix, derivative_matrix;
961
962 for (unsigned int d = 0; d < dim; ++d)
963 {
964 mass_matrix[d] = &this->mass_matrix[d](0, 0);
965 derivative_matrix[d] = &this->derivative_matrix[d](0, 0);
966 }
967
968 const unsigned int n_rows_1d_non_templated = this->mass_matrix[0].n_rows();
969
970 if (n_rows_1d != -1)
971 internal::TensorProductMatrixSymmetricSum::vmult<
972 n_rows_1d == -1 ? 0 : n_rows_1d>(dst,
973 src,
974 tmp_array,
975 n_rows_1d_non_templated,
976 mass_matrix,
977 derivative_matrix);
978 else
979 internal::TensorProductMatrixSymmetricSum::select_vmult<1>(
980 dst,
981 src,
982 tmp_array,
983 n_rows_1d_non_templated,
984 mass_matrix,
985 derivative_matrix);
986}
987
988
989
990template <int dim, typename Number, int n_rows_1d>
991inline void
992TensorProductMatrixSymmetricSum<dim, Number, n_rows_1d>::apply_inverse(
993 const ArrayView<Number> &dst_view,
994 const ArrayView<const Number> &src_view) const
995{
996 AssertDimension(dst_view.size(), this->n());
997 AssertDimension(src_view.size(), this->m());
998
999 Number *dst = dst_view.begin();
1000 const Number *src = src_view.begin();
1001
1002 std::array<const Number *, dim> eigenvectors, eigenvalues;
1003
1004 for (unsigned int d = 0; d < dim; ++d)
1005 {
1006 eigenvectors[d] = &this->eigenvectors[d](0, 0);
1007 eigenvalues[d] = this->eigenvalues[d].data();
1008 }
1009
1010 const unsigned int n_rows_1d_non_templated = this->mass_matrix[0].n_rows();
1011
1012 if (n_rows_1d != -1)
1013 internal::TensorProductMatrixSymmetricSum::apply_inverse<
1014 n_rows_1d == -1 ? 0 : n_rows_1d>(
1015 dst, src, n_rows_1d_non_templated, eigenvectors, eigenvalues);
1016 else
1017 internal::TensorProductMatrixSymmetricSum::select_apply_inverse<1>(
1018 dst, src, n_rows_1d_non_templated, eigenvectors, eigenvalues);
1019}
1020
1021
1022
1023template <int dim, typename Number, int n_rows_1d>
1024std::size_t
1025TensorProductMatrixSymmetricSum<dim, Number, n_rows_1d>::memory_consumption()
1026 const
1027{
1028 return MemoryConsumption::memory_consumption(mass_matrix) +
1029 MemoryConsumption::memory_consumption(derivative_matrix) +
1030 MemoryConsumption::memory_consumption(eigenvalues) +
1031 MemoryConsumption::memory_consumption(eigenvectors) +
1032 MemoryConsumption::memory_consumption(tmp_array);
1033}
1034
1035
1036
1037template <int dim, typename Number, int n_rows_1d>
1038template <typename T>
1039inline TensorProductMatrixSymmetricSum<dim, Number, n_rows_1d>::
1040 TensorProductMatrixSymmetricSum(const T &mass_matrix,
1041 const T &derivative_matrix)
1042{
1043 reinit(mass_matrix, derivative_matrix);
1044}
1045
1046
1047
1048template <int dim, typename Number, int n_rows_1d>
1049template <typename T>
1050inline void
1051TensorProductMatrixSymmetricSum<dim, Number, n_rows_1d>::reinit(
1052 const T &mass_matrix,
1053 const T &derivative_matrix)
1054{
1055 this->mass_matrix =
1056 internal::TensorProductMatrixSymmetricSum::convert<dim>(mass_matrix);
1057 this->derivative_matrix =
1058 internal::TensorProductMatrixSymmetricSum::convert<dim>(derivative_matrix);
1059
1060 internal::TensorProductMatrixSymmetricSum::setup(this->mass_matrix,
1061 this->derivative_matrix,
1062 this->eigenvectors,
1063 this->eigenvalues);
1064}
1065
1066
1067
1068template <int dim, typename Number, int n_rows_1d>
1069TensorProductMatrixSymmetricSumCollection<dim, Number, n_rows_1d>::
1070 AdditionalData::AdditionalData(const bool compress_matrices,
1071 const bool precompute_inverse_diagonal)
1072 : compress_matrices(compress_matrices)
1073 , precompute_inverse_diagonal(precompute_inverse_diagonal)
1074{}
1075
1076
1077
1078template <int dim, typename Number, int n_rows_1d>
1079TensorProductMatrixSymmetricSumCollection<dim, Number, n_rows_1d>::
1080 TensorProductMatrixSymmetricSumCollection(
1081 const AdditionalData &additional_data)
1082 : compress_matrices(additional_data.compress_matrices)
1083 , precompute_inverse_diagonal(additional_data.precompute_inverse_diagonal)
1084{}
1085
1086
1087
1088template <int dim, typename Number, int n_rows_1d>
1089void
1090TensorProductMatrixSymmetricSumCollection<dim, Number, n_rows_1d>::reserve(
1091 const unsigned int size)
1092{
1093 if (compress_matrices == false)
1094 mass_and_derivative_matrices.resize(size * dim);
1095 else
1096 indices.assign(size * dim, numbers::invalid_unsigned_int);
1097}
1098
1099
1100
1101template <int dim, typename Number, int n_rows_1d>
1102template <typename T>
1103void
1104TensorProductMatrixSymmetricSumCollection<dim, Number, n_rows_1d>::insert(
1105 const unsigned int index,
1106 const T &Ms_in,
1107 const T &Ks_in)
1108{
1109 const auto Ms =
1110 internal::TensorProductMatrixSymmetricSum::convert<dim>(Ms_in);
1111 const auto Ks =
1112 internal::TensorProductMatrixSymmetricSum::convert<dim>(Ks_in);
1113
1114 for (unsigned int d = 0; d < dim; ++d)
1115 {
1116 if (compress_matrices == false)
1117 {
1118 const MatrixPairType matrix(Ms[d], Ks[d]);
1119 mass_and_derivative_matrices[index * dim + d] = matrix;
1120 }
1121 else
1122 {
1123 using VectorizedArrayTrait =
1124 ::internal::VectorizedArrayTrait<Number>;
1125
1126 std::bitset<VectorizedArrayTrait::width()> mask;
1127
1128 for (unsigned int v = 0; v < VectorizedArrayTrait::width(); ++v)
1129 {
1130 typename VectorizedArrayTrait::value_type a = 0.0;
1131
1132 for (unsigned int i = 0; i < Ms[d].size(0); ++i)
1133 for (unsigned int j = 0; j < Ms[d].size(1); ++j)
1134 {
1135 a += std::abs(VectorizedArrayTrait::get(Ms[d][i][j], v));
1136 a += std::abs(VectorizedArrayTrait::get(Ks[d][i][j], v));
1137 }
1138
1139 mask[v] = (a != 0.0);
1140 }
1141
1142 const MatrixPairTypeWithMask matrix{mask, {Ms[d], Ks[d]}};
1143
1144 const auto ptr = cache.find(matrix);
1145
1146 if (ptr != cache.end())
1147 {
1148 const auto ptr_index = ptr->second;
1149 indices[index * dim + d] = ptr_index;
1150
1151 if ([&]() {
1152 for (unsigned int v = 0; v < VectorizedArrayTrait::width();
1153 ++v)
1154 if ((mask[v] == true) && (ptr->first.first[v] == false))
1155 return false;
1156
1157 return true;
1158 }())
1159 {
1160 // nothing to do
1161 }
1162 else
1163 {
1164 auto mask_new = ptr->first.first;
1165 auto Ms_new = ptr->first.second.first;
1166 auto Ks_new = ptr->first.second.second;
1167
1168 for (unsigned int v = 0; v < VectorizedArrayTrait::width();
1169 ++v)
1170 if (mask_new[v] == false && mask[v] == true)
1171 {
1172 mask_new[v] = true;
1173
1174 for (unsigned int i = 0; i < Ms_new.size(0); ++i)
1175 for (unsigned int j = 0; j < Ms_new.size(1); ++j)
1176 {
1177 VectorizedArrayTrait::get(Ms_new[i][j], v) =
1178 VectorizedArrayTrait::get(Ms[d][i][j], v);
1179 VectorizedArrayTrait::get(Ks_new[i][j], v) =
1180 VectorizedArrayTrait::get(Ks[d][i][j], v);
1181 }
1182 }
1183
1184 cache.erase(ptr);
1185
1186 const MatrixPairTypeWithMask entry_new{mask_new,
1187 {Ms_new, Ks_new}};
1188
1189 const auto ptr_ = cache.find(entry_new);
1190 AssertThrow(ptr_ == cache.end(), ExcNotImplemented());
1191
1192 cache[entry_new] = ptr_index;
1193 }
1194 }
1195 else
1196 {
1197 const auto size = cache.size();
1198 indices[index * dim + d] = size;
1199 cache[matrix] = size;
1200 }
1201 }
1202 }
1203}
1204
1205
1206
1207template <int dim, typename Number, int n_rows_1d>
1208void
1209TensorProductMatrixSymmetricSumCollection<dim, Number, n_rows_1d>::finalize()
1210{
1211 const auto store = [&](const unsigned int index,
1212 const MatrixPairType &M_and_K) {
1213 std::array<Table<2, Number>, 1> mass_matrix;
1214 mass_matrix[0] = M_and_K.first;
1215
1216 std::array<Table<2, Number>, 1> derivative_matrix;
1217 derivative_matrix[0] = M_and_K.second;
1218
1219 std::array<Table<2, Number>, 1> eigenvectors;
1220 std::array<AlignedVector<Number>, 1> eigenvalues;
1221
1222 internal::TensorProductMatrixSymmetricSum::setup(mass_matrix,
1223 derivative_matrix,
1224 eigenvectors,
1225 eigenvalues);
1226
1227 for (unsigned int i = 0, m = matrix_ptr[index], v = vector_ptr[index];
1228 i < mass_matrix[0].n_rows();
1229 ++i, ++v)
1230 {
1231 for (unsigned int j = 0; j < mass_matrix[0].n_cols(); ++j, ++m)
1232 {
1233 this->mass_matrices[m] = mass_matrix[0][i][j];
1234 this->derivative_matrices[m] = derivative_matrix[0][i][j];
1235 this->eigenvectors[m] = eigenvectors[0][i][j];
1236 }
1237
1238 this->eigenvalues[v] = eigenvalues[0][i];
1239 }
1240 };
1241
1242 if (compress_matrices == false)
1243 {
1244 // case 1) no compression requested
1245
1246 AssertDimension(cache.size(), 0);
1247 AssertDimension(indices.size(), 0);
1248
1249 this->vector_ptr.resize(mass_and_derivative_matrices.size() + 1);
1250 this->matrix_ptr.resize(mass_and_derivative_matrices.size() + 1);
1251
1252 for (unsigned int i = 0; i < mass_and_derivative_matrices.size(); ++i)
1253 {
1254 const auto &M = mass_and_derivative_matrices[i].first;
1255
1256 this->vector_ptr[i + 1] = M.n_rows();
1257 this->matrix_ptr[i + 1] = M.n_rows() * M.n_cols();
1258 }
1259
1260 for (unsigned int i = 0; i < mass_and_derivative_matrices.size(); ++i)
1261 {
1262 this->vector_ptr[i + 1] += this->vector_ptr[i];
1263 this->matrix_ptr[i + 1] += this->matrix_ptr[i];
1264 }
1265
1266 this->mass_matrices.resize_fast(matrix_ptr.back());
1267 this->derivative_matrices.resize_fast(matrix_ptr.back());
1268 this->eigenvectors.resize_fast(matrix_ptr.back());
1269 this->eigenvalues.resize_fast(vector_ptr.back());
1270
1271 for (unsigned int i = 0; i < mass_and_derivative_matrices.size(); ++i)
1272 store(i, mass_and_derivative_matrices[i]);
1273
1274 mass_and_derivative_matrices.clear();
1275 }
1276 else if (cache.size() == indices.size())
1277 {
1278 // case 2) compression requested but none possible
1279
1280 this->vector_ptr.resize(cache.size() + 1);
1281 this->matrix_ptr.resize(cache.size() + 1);
1282
1283 std::map<unsigned int, MatrixPairType> inverted_cache;
1284
1285 for (const auto &i : cache)
1286 inverted_cache[i.second] = i.first.second;
1287
1288 for (unsigned int i = 0; i < indices.size(); ++i)
1289 {
1290 const auto &M = inverted_cache[indices[i]].first;
1291
1292 this->vector_ptr[i + 1] = M.n_rows();
1293 this->matrix_ptr[i + 1] = M.n_rows() * M.n_cols();
1294 }
1295
1296 for (unsigned int i = 0; i < cache.size(); ++i)
1297 {
1298 this->vector_ptr[i + 1] += this->vector_ptr[i];
1299 this->matrix_ptr[i + 1] += this->matrix_ptr[i];
1300 }
1301
1302 this->mass_matrices.resize_fast(matrix_ptr.back());
1303 this->derivative_matrices.resize_fast(matrix_ptr.back());
1304 this->eigenvectors.resize_fast(matrix_ptr.back());
1305 this->eigenvalues.resize_fast(vector_ptr.back());
1306
1307 for (unsigned int i = 0; i < indices.size(); ++i)
1308 store(i, inverted_cache[indices[i]]);
1309
1310 indices.clear();
1311 cache.clear();
1312 }
1313 else
1314 {
1315 // case 3) compress
1316
1317 this->vector_ptr.resize(cache.size() + 1);
1318 this->matrix_ptr.resize(cache.size() + 1);
1319
1320 for (const auto &i : cache)
1321 {
1322 const auto &M = i.first.second.first;
1323
1324 this->vector_ptr[i.second + 1] = M.n_rows();
1325 this->matrix_ptr[i.second + 1] = M.n_rows() * M.n_cols();
1326 }
1327
1328 for (unsigned int i = 0; i < cache.size(); ++i)
1329 {
1330 this->vector_ptr[i + 1] += this->vector_ptr[i];
1331 this->matrix_ptr[i + 1] += this->matrix_ptr[i];
1332 }
1333
1334 this->mass_matrices.resize_fast(matrix_ptr.back());
1335 this->derivative_matrices.resize_fast(matrix_ptr.back());
1336 this->eigenvectors.resize_fast(matrix_ptr.back());
1337 this->eigenvalues.resize_fast(vector_ptr.back());
1338
1339 for (const auto &i : cache)
1340 store(i.second, i.first.second);
1341
1342 cache.clear();
1343 }
1344
1345 if (precompute_inverse_diagonal)
1346 {
1347 if (dim == 1)
1348 {
1349 // 1D case: simply invert 1D eigenvalues
1350 for (unsigned int i = 0; i < this->eigenvalues.size(); ++i)
1351 this->eigenvalues[i] = Number(1.0) / this->eigenvalues[i];
1352 std::swap(this->inverted_eigenvalues, eigenvalues);
1353 }
1354 else
1355 {
1356 // 2D and 3D case: we have 2 or 3 1d eigenvalues so that we
1357 // need to combine these
1358
1359 // step 1) if eigenvalues/eigenvectors are compressed, we
1360 // need to compress the diagonal (the combination of ev
1361 // indices) as well. This is an optional step.
1362 std::vector<unsigned int> indices_ev;
1363
1364 if (indices.size() > 0)
1365 {
1366 // 1a) create cache (ev indics -> diag index)
1367 const unsigned int n_cells = indices.size() / dim;
1368 std::map<std::array<unsigned int, dim>, unsigned int> cache_ev;
1369 std::vector<unsigned int> cache_ev_idx(n_cells);
1370
1371 for (unsigned int i = 0, c = 0; i < n_cells; ++i)
1372 {
1373 std::array<unsigned int, dim> id;
1374
1375 for (unsigned int d = 0; d < dim; ++d, ++c)
1376 id[d] = indices[c];
1377
1378 const auto id_ptr = cache_ev.find(id);
1379
1380 if (id_ptr == cache_ev.end())
1381 {
1382 const auto size = cache_ev.size();
1383 cache_ev_idx[i] = size;
1384 cache_ev[id] = size;
1385 }
1386 else
1387 {
1388 cache_ev_idx[i] = id_ptr->second;
1389 }
1390 }
1391
1392 // 1b) store diagonal indices for each cell
1393 std::vector<unsigned int> new_indices;
1394 new_indices.reserve(indices.size() / dim * (dim + 1));
1395
1396 for (unsigned int i = 0, c = 0; i < n_cells; ++i)
1397 {
1398 for (unsigned int d = 0; d < dim; ++d, ++c)
1399 new_indices.push_back(indices[c]);
1400 new_indices.push_back(cache_ev_idx[i]);
1401 }
1402
1403 // 1c) transpose cache (diag index -> ev indices)
1404 indices_ev.resize(cache_ev.size() * dim);
1405 for (const auto &entry : cache_ev)
1406 for (unsigned int d = 0; d < dim; ++d)
1407 indices_ev[entry.second * dim + d] = entry.first[d];
1408
1409 std::swap(this->indices, new_indices);
1410 }
1411
1412 // step 2) allocate memory and set pointers
1413 const unsigned int n_diag =
1414 ((indices_ev.size() > 0) ? indices_ev.size() :
1415 (matrix_ptr.size() - 1)) /
1416 dim;
1417
1418 std::vector<unsigned int> new_vector_ptr(n_diag + 1, 0);
1419 std::vector<unsigned int> new_vector_n_rows_1d(n_diag, 0);
1420
1421 for (unsigned int i = 0; i < n_diag; ++i)
1422 {
1423 const unsigned int c = (indices_ev.size() > 0) ?
1424 indices_ev[dim * i + 0] :
1425 (dim * i + 0);
1426
1427 const unsigned int n_rows = vector_ptr[c + 1] - vector_ptr[c];
1428
1429 new_vector_n_rows_1d[i] = n_rows;
1430 new_vector_ptr[i + 1] = Utilities::pow(n_rows, dim);
1431 }
1432
1433 for (unsigned int i = 0; i < n_diag; ++i)
1434 new_vector_ptr[i + 1] += new_vector_ptr[i];
1435
1436 this->inverted_eigenvalues.resize(new_vector_ptr.back());
1437
1438 // step 3) loop over all unique diagonal entries and invert
1439 for (unsigned int i = 0; i < n_diag; ++i)
1440 {
1441 std::array<Number *, dim> evs;
1442
1443 for (unsigned int d = 0; d < dim; ++d)
1444 evs[d] =
1445 &this
1446 ->eigenvalues[this->vector_ptr[(indices_ev.size() > 0) ?
1447 indices_ev[dim * i + d] :
1448 (dim * i + d)]];
1449
1450 const unsigned int mm = new_vector_n_rows_1d[i];
1451 if (dim == 2)
1452 {
1453 for (unsigned int i1 = 0, c = 0; i1 < mm; ++i1)
1454 for (unsigned int i0 = 0; i0 < mm; ++i0, ++c)
1455 this->inverted_eigenvalues[new_vector_ptr[i] + c] =
1456 Number(1.0) / (evs[1][i1] + evs[0][i0]);
1457 }
1458 else
1459 {
1460 for (unsigned int i2 = 0, c = 0; i2 < mm; ++i2)
1461 for (unsigned int i1 = 0; i1 < mm; ++i1)
1462 for (unsigned int i0 = 0; i0 < mm; ++i0, ++c)
1463 this->inverted_eigenvalues[new_vector_ptr[i] + c] =
1464 Number(1.0) / (evs[2][i2] + evs[1][i1] + evs[0][i0]);
1465 }
1466 }
1467
1468 // step 4) clean up
1469 std::swap(this->vector_ptr, new_vector_ptr);
1470 std::swap(this->vector_n_rows_1d, new_vector_n_rows_1d);
1471 }
1472
1473 this->eigenvalues.clear();
1474 }
1475}
1476
1477
1478
1479template <int dim, typename Number, int n_rows_1d>
1480void
1481TensorProductMatrixSymmetricSumCollection<dim, Number, n_rows_1d>::
1482 apply_inverse(const unsigned int index,
1483 const ArrayView<Number> &dst_in,
1484 const ArrayView<const Number> &src_in) const
1485{
1486 Number *dst = dst_in.begin();
1487 const Number *src = src_in.begin();
1488
1489 if (this->eigenvalues.empty() == false)
1490 {
1491 std::array<const Number *, dim> eigenvectors;
1492 std::array<const Number *, dim> eigenvalues;
1493 unsigned int n_rows_1d_non_templated = 0;
1494
1495 for (unsigned int d = 0; d < dim; ++d)
1496 {
1497 const unsigned int translated_index =
1498 (indices.size() > 0) ? indices[dim * index + d] : (dim * index + d);
1499
1500 eigenvectors[d] =
1501 this->eigenvectors.data() + matrix_ptr[translated_index];
1502 eigenvalues[d] =
1503 this->eigenvalues.data() + vector_ptr[translated_index];
1504 n_rows_1d_non_templated =
1505 vector_ptr[translated_index + 1] - vector_ptr[translated_index];
1506 }
1507
1508 if (n_rows_1d != -1)
1509 internal::TensorProductMatrixSymmetricSum::apply_inverse<
1510 n_rows_1d == -1 ? 0 : n_rows_1d>(
1511 dst, src, n_rows_1d_non_templated, eigenvectors, eigenvalues);
1512 else
1513 internal::TensorProductMatrixSymmetricSum::select_apply_inverse<1>(
1514 dst, src, n_rows_1d_non_templated, eigenvectors, eigenvalues);
1515 }
1516 else
1517 {
1518 std::array<const Number *, dim> eigenvectors;
1519 const Number *inverted_eigenvalues = nullptr;
1520 unsigned int n_rows_1d_non_templated = 0;
1521
1522 for (unsigned int d = 0; d < dim; ++d)
1523 {
1524 const unsigned int translated_index =
1525 (indices.size() > 0) ?
1526 indices[((dim == 1) ? 1 : (dim + 1)) * index + d] :
1527 (dim * index + d);
1528
1529 eigenvectors[d] =
1530 this->eigenvectors.data() + matrix_ptr[translated_index];
1531 }
1532
1533 {
1534 const unsigned int translated_index =
1535 ((indices.size() > 0) && (dim != 1)) ?
1536 indices[(dim + 1) * index + dim] :
1537 index;
1538
1539 inverted_eigenvalues =
1540 this->inverted_eigenvalues.data() + vector_ptr[translated_index];
1541 n_rows_1d_non_templated =
1542 (dim == 1) ?
1543 (vector_ptr[translated_index + 1] - vector_ptr[translated_index]) :
1544 vector_n_rows_1d[translated_index];
1545 }
1546
1547 if (n_rows_1d != -1)
1548 internal::TensorProductMatrixSymmetricSum::apply_inverse<
1549 n_rows_1d == -1 ? 0 : n_rows_1d>(dst,
1550 src,
1551 n_rows_1d_non_templated,
1552 eigenvectors,
1553 {},
1554 inverted_eigenvalues);
1555 else
1556 internal::TensorProductMatrixSymmetricSum::select_apply_inverse<1>(
1557 dst,
1558 src,
1559 n_rows_1d_non_templated,
1560 eigenvectors,
1561 {},
1562 inverted_eigenvalues);
1563 }
1564}
1565
1566
1567
1568template <int dim, typename Number, int n_rows_1d>
1569std::size_t
1570TensorProductMatrixSymmetricSumCollection<dim, Number, n_rows_1d>::
1571 memory_consumption() const
1572{
1573 return MemoryConsumption::memory_consumption(indices) +
1574 MemoryConsumption::memory_consumption(mass_matrices) +
1575 MemoryConsumption::memory_consumption(derivative_matrices) +
1576 MemoryConsumption::memory_consumption(eigenvectors) +
1577 MemoryConsumption::memory_consumption(eigenvalues) +
1578 MemoryConsumption::memory_consumption(matrix_ptr) +
1579 MemoryConsumption::memory_consumption(vector_ptr);
1580}
1581
1582
1583
1584template <int dim, typename Number, int n_rows_1d>
1585std::size_t
1586TensorProductMatrixSymmetricSumCollection<dim, Number, n_rows_1d>::
1587 storage_size() const
1588{
1589 if (matrix_ptr.empty())
1590 return 0; // if not initialized
1591
1592 return matrix_ptr.size() - 1;
1593}
1594
1595
1596
1597#endif
1598
1599DEAL_II_NAMESPACE_CLOSE
1600
1601#endif
array_view.h
AlignedVector::resize_fast
void resize_fast(const size_type new_size)
AlignedVector::begin
iterator begin()
ArrayView::begin
iterator begin() const
Definition array_view.h:702
ArrayView::size
std::size_t size() const
Definition array_view.h:684
LAPACKFullMatrix::eigenvalue
std::complex< typename numbers::NumberTraits< number >::real_type > eigenvalue(const size_type i) const
Definition lapack_full_matrix.h:1208
LAPACKFullMatrix::compute_generalized_eigenvalues_symmetric
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)
Definition lapack_full_matrix.cc:2235
TensorProductMatrixSymmetricSumCollection::apply_inverse
void apply_inverse(const unsigned int index, const ArrayView< Number > &dst_in, const ArrayView< const Number > &src_in) const
TensorProductMatrixSymmetricSumCollection::mass_and_derivative_matrices
std::vector< MatrixPairType > mass_and_derivative_matrices
Definition tensor_product_matrix.h:434
TensorProductMatrixSymmetricSumCollection::MatrixPairTypeWithMask
std::pair< std::bitset<::internal::VectorizedArrayTrait< Number >::width()>, MatrixPairType > MatrixPairTypeWithMask
Definition tensor_product_matrix.h:340
TensorProductMatrixSymmetricSumCollection::MatrixPairType
std::pair< Table< 2, Number >, Table< 2, Number > > MatrixPairType
Definition tensor_product_matrix.h:336
TensorProductMatrixSymmetricSumCollection::reserve
void reserve(const unsigned int size)
TensorProductMatrixSymmetricSumCollection::insert
void insert(const unsigned int index, const T &Ms, const T &Ks)
TensorProductMatrixSymmetricSumCollection::cache
std::map< MatrixPairTypeWithMask, unsigned int, internal::TensorProductMatrixSymmetricSum::MatrixPairComparator< Number > > cache
Definition tensor_product_matrix.h:444
TensorProductMatrixSymmetricSumCollection::TensorProductMatrixSymmetricSumCollection
TensorProductMatrixSymmetricSumCollection(const AdditionalData &additional_data=AdditionalData())
TensorProductMatrixSymmetricSum::vmult
void vmult(const ArrayView< Number > &dst, const ArrayView< const Number > &src, AlignedVector< Number > &tmp) const
TensorProductMatrixSymmetricSum::eigenvectors
std::array< Table< 2, Number >, dim > eigenvectors
Definition tensor_product_matrix.h:240
TensorProductMatrixSymmetricSum::derivative_matrix
std::array< Table< 2, Number >, dim > derivative_matrix
Definition tensor_product_matrix.h:228
TensorProductMatrixSymmetricSum::reinit
void reinit(const T &mass_matrix, const T &derivative_matrix)
TensorProductMatrixSymmetricSum::apply_inverse
void apply_inverse(const ArrayView< Number > &dst, const ArrayView< const Number > &src) const
TensorProductMatrixSymmetricSum::vmult
void vmult(const ArrayView< Number > &dst, const ArrayView< const Number > &src) const
TensorProductMatrixSymmetricSum::memory_consumption
std::size_t memory_consumption() const
TensorProductMatrixSymmetricSum::mass_matrix
std::array< Table< 2, Number >, dim > mass_matrix
Definition tensor_product_matrix.h:223
TensorProductMatrixSymmetricSum::eigenvalues
std::array< AlignedVector< Number >, dim > eigenvalues
Definition tensor_product_matrix.h:234
TensorProductMatrixSymmetricSum::TensorProductMatrixSymmetricSum
TensorProductMatrixSymmetricSum(const T &mass_matrix, const T &derivative_matrix)
Tensor::clear
constexpr void clear()
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:498
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
floating_point_comparator.h
second
Point< 2 > second
Definition grid_out.cc:4630
first
Point< 2 > first
Definition grid_out.cc:4629
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Assert
#define Assert(cond, exc)
Definition exceptions.h:1659
AssertDimension
#define AssertDimension(dim1, dim2)
Definition exceptions.h:2026
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
AssertThrow
#define AssertThrow(cond, exc)
Definition exceptions.h:1780
DEAL_II_NOT_IMPLEMENTED
#define DEAL_II_NOT_IMPLEMENTED()
Definition exceptions.h:1855
lapack_full_matrix.h
size
std::size_t size
Definition mpi.cc:734
LAPACKSupport::matrix
@ matrix
Contents is actually a matrix.
Definition lapack_support.h:57
LocalIntegrators::L2::mass_matrix
void mass_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const double factor=1.)
Definition l2.h:57
MemoryConsumption::memory_consumption
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
Definition memory_consumption.h:271
Physics::Elasticity::Kinematics::d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Utilities::pow
constexpr T pow(const T base, const int iexp)
Definition utilities.h:966
internal::TriangulationImplementation::n_cells
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition tria.cc:14874
internal::reinit
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition matrix_block.h:599
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition types.h:220
std
STL namespace.
std::abs
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6927
TensorProductMatrixSymmetricSumCollection::AdditionalData::AdditionalData
AdditionalData(const bool compress_matrices=true, const bool precompute_inverse_diagonal=true)
internal::TensorProductMatrixSymmetricSum::MatrixPairComparator::operator()
bool operator()(const MatrixPairType &left, const MatrixPairType &right) const
Definition tensor_product_matrix.h:277
internal::TensorProductMatrixSymmetricSum::MatrixPairComparator::MatrixPairType
std::pair< std::bitset< width >, std::pair< Table< 2, Number >, Table< 2, Number > > > MatrixPairType
Definition tensor_product_matrix.h:270
internal::VectorizedArrayTrait< Number >::width
static constexpr std::size_t width()
Definition vectorization.h:6415
eigenvalues
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
eigenvectors
std::array< std::pair< Number, Tensor< 1, dim, Number > >, dim > eigenvectors(const SymmetricTensor< 2, dim, Number > &T, const SymmetricTensorEigenvectorMethod method=SymmetricTensorEigenvectorMethod::ql_implicit_shifts)
tensor_product_kernels.h
Generated by doxygen 1.9.8