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