15#ifndef dealii_tensor_product_matrix_h
16#define dealii_tensor_product_matrix_h
114template <
int dim,
typename Number,
int n_rows_1d = -1>
139 template <
typename T>
160 template <
typename T>
263 template <
typename Number>
272 std::pair<std::bitset<width>,
282 const auto &M_0 = left.second.first;
283 const auto &K_0 = left.second.second;
284 const auto &M_1 = right.second.first;
285 const auto &K_1 = right.second.second;
287 std::bitset<width>
mask;
289 for (
unsigned int v = 0; v <
width; ++v)
290 mask[v] = left.first[v] && right.first[v];
295 if (comparator(M_0, M_1))
297 else if (comparator(M_1, M_0))
299 else if (comparator(K_0, K_1))
336template <
int dim,
typename Number,
int n_rows_1d = -1>
342 std::bitset<::internal::VectorizedArrayTrait<Number>::width()>,
386 template <
typename T>
388 insert(
const unsigned int index,
const T &Ms,
const T &Ks);
516 template <
typename Number>
518 spectral_assembly(
const Number *mass_matrix,
519 const Number *derivative_matrix,
520 const unsigned int n_rows,
521 const unsigned int n_cols,
527 std::vector<bool> constrained_dofs(n_rows,
false);
529 for (
unsigned int i = 0; i < n_rows; ++i)
531 if (mass_matrix[i + i * n_rows] == 0.0)
533 Assert(derivative_matrix[i + i * n_rows] == 0.0,
536 for (
unsigned int j = 0; j < n_rows; ++j)
538 Assert(derivative_matrix[i + j * n_rows] == 0,
540 Assert(derivative_matrix[j + i * n_rows] == 0,
544 constrained_dofs[i] =
true;
548 const auto transpose_fill_nm = [&constrained_dofs](Number *out,
550 const unsigned int n,
551 const unsigned int m) {
552 for (
unsigned int mm = 0, c = 0; mm < m; ++mm)
553 for (
unsigned int nn = 0; nn < n; ++nn, ++c)
555 (mm == nn && constrained_dofs[mm]) ? Number(1.0) : in[c];
558 std::vector<::Vector<Number>> eigenvecs(n_rows);
562 transpose_fill_nm(&(mass_copy(0, 0)), mass_matrix, n_rows, n_cols);
563 transpose_fill_nm(&(deriv_copy(0, 0)), derivative_matrix, n_rows, n_cols);
568 for (
unsigned int i = 0, c = 0; i < n_rows; ++i)
569 for (
unsigned int j = 0; j < n_cols; ++j, ++c)
570 if (constrained_dofs[i] ==
false)
573 for (
unsigned int i = 0; i < n_rows; ++i, ++
eigenvalues)
579 template <std::
size_t dim,
typename Number>
586 const unsigned int n_rows_1d = mass_matrix[0].n_cols();
589 for (
unsigned int dir = 0; dir < dim; ++dir)
594 derivative_matrix[dir].n_rows());
596 derivative_matrix[dir].n_cols());
599 mass_matrix[dir].n_rows());
600 eigenvalues[dir].resize(mass_matrix[dir].n_cols());
601 internal::TensorProductMatrixSymmetricSum::spectral_assembly<Number>(
602 &(mass_matrix[dir](0, 0)),
603 &(derivative_matrix[dir](0, 0)),
604 mass_matrix[dir].n_rows(),
605 mass_matrix[dir].n_cols(),
613 template <std::
size_t dim,
typename Number, std::
size_t n_lanes>
624 const unsigned int n_rows_1d = mass_matrix[0].n_cols();
625 constexpr unsigned int macro_size =
627 const std::size_t nm_flat_size_max = n_rows_1d * n_rows_1d * macro_size;
628 const std::size_t n_flat_size_max = n_rows_1d * macro_size;
630 std::vector<Number> mass_matrix_flat;
631 std::vector<Number> deriv_matrix_flat;
632 std::vector<Number> eigenvalues_flat;
633 std::vector<Number> eigenvectors_flat;
634 mass_matrix_flat.resize(nm_flat_size_max);
635 deriv_matrix_flat.resize(nm_flat_size_max);
636 eigenvalues_flat.resize(n_flat_size_max);
637 eigenvectors_flat.resize(nm_flat_size_max);
638 std::array<unsigned int, macro_size> offsets_nm;
639 std::array<unsigned int, macro_size> offsets_n;
640 for (
unsigned int dir = 0; dir < dim; ++dir)
645 derivative_matrix[dir].n_rows());
647 derivative_matrix[dir].n_cols());
649 const unsigned int n_rows = mass_matrix[dir].n_rows();
650 const unsigned int n_cols = mass_matrix[dir].n_cols();
651 const unsigned int nm = n_rows * n_cols;
652 for (
unsigned int vv = 0; vv < macro_size; ++vv)
653 offsets_nm[vv] = nm * vv;
655 vectorized_transpose_and_store<Number, n_lanes>(
658 &(mass_matrix[dir](0, 0)),
660 mass_matrix_flat.data());
661 vectorized_transpose_and_store<Number, n_lanes>(
664 &(derivative_matrix[dir](0, 0)),
666 deriv_matrix_flat.data());
668 const Number *mass_cbegin = mass_matrix_flat.data();
669 const Number *deriv_cbegin = deriv_matrix_flat.data();
670 Number *eigenvec_begin = eigenvectors_flat.data();
671 Number *eigenval_begin = eigenvalues_flat.data();
672 for (
unsigned int lane = 0; lane < macro_size; ++lane)
673 internal::TensorProductMatrixSymmetricSum::spectral_assembly<
674 Number>(mass_cbegin + nm * lane,
675 deriv_cbegin + nm * lane,
678 eigenval_begin + n_rows * lane,
679 eigenvec_begin + nm * lane);
683 for (
unsigned int vv = 0; vv < macro_size; ++vv)
684 offsets_n[vv] = n_rows * vv;
685 vectorized_load_and_transpose<Number, n_lanes>(
687 eigenvalues_flat.data(),
690 vectorized_load_and_transpose<Number, n_lanes>(
692 eigenvectors_flat.data(),
700 template <std::
size_t dim,
typename Number>
701 inline std::array<Table<2, Number>, dim>
709 template <std::
size_t dim,
typename Number>
710 inline std::array<Table<2, Number>, dim>
713 std::array<Table<2, Number>, dim> mass_copy;
715 std::transform(mass_matrix.cbegin(),
727 template <std::
size_t dim,
typename Number>
728 inline std::array<Table<2, Number>, dim>
731 std::array<Table<2, Number>, dim> matrices;
733 std::fill(matrices.begin(), matrices.end(), matrix);
740 template <
int n_rows_1d_templated, std::
size_t dim,
typename Number>
745 const unsigned int n_rows_1d_non_templated,
746 const std::array<const Number *, dim> &mass_matrix,
747 const std::array<const Number *, dim> &derivative_matrix)
749 const unsigned int n_rows_1d = n_rows_1d_templated == 0 ?
750 n_rows_1d_non_templated :
752 const unsigned int n = Utilities::fixed_power<dim>(n_rows_1d);
755 Number *t = tmp.
begin();
762 eval({}, {}, {}, n_rows_1d, n_rows_1d);
766 const Number *A = derivative_matrix[0];
767 eval.template apply<0, false, false>(A, src, dst);
772 const Number *A0 = derivative_matrix[0];
773 const Number *M0 = mass_matrix[0];
774 const Number *A1 = derivative_matrix[1];
775 const Number *M1 = mass_matrix[1];
776 eval.template apply<0, false, false>(M0, src, t);
777 eval.template apply<1, false, false>(A1, t, dst);
778 eval.template apply<0, false, false>(A0, src, t);
779 eval.template apply<1, false, true>(M1, t, dst);
784 const Number *A0 = derivative_matrix[0];
785 const Number *M0 = mass_matrix[0];
786 const Number *A1 = derivative_matrix[1];
787 const Number *M1 = mass_matrix[1];
788 const Number *A2 = derivative_matrix[2];
789 const Number *M2 = mass_matrix[2];
790 eval.template apply<0, false, false>(M0, src, t + n);
791 eval.template apply<1, false, false>(M1, t + n, t);
792 eval.template apply<2, false, false>(A2, t, dst);
793 eval.template apply<1, false, false>(A1, t + n, t);
794 eval.template apply<0, false, false>(A0, src, t + n);
795 eval.template apply<1, false, true>(M1, t + n, t);
796 eval.template apply<2, false, true>(M2, t, dst);
805 template <
int n_rows_1d_templated, std::
size_t dim,
typename Number>
807 apply_inverse(Number *dst,
809 const unsigned int n_rows_1d_non_templated,
811 const std::array<const Number *, dim> &
eigenvalues,
812 const Number *inverted_eigenvalues =
nullptr)
814 const unsigned int n_rows_1d = n_rows_1d_templated == 0 ?
815 n_rows_1d_non_templated :
823 eval({}, {}, {}, n_rows_1d, n_rows_1d);
833 eval.template apply<0, true, false>(S, src, dst);
835 for (
unsigned int i = 0; i < n_rows_1d; ++i)
836 if (inverted_eigenvalues)
837 dst[i] *= inverted_eigenvalues[i];
841 eval.template apply<0, false, false>(S, dst, dst);
848 eval.template apply<0, true, false>(S0, src, dst);
849 eval.template apply<1, true, false>(S1, dst, dst);
851 for (
unsigned int i1 = 0, c = 0; i1 < n_rows_1d; ++i1)
852 for (
unsigned int i0 = 0; i0 < n_rows_1d; ++i0, ++c)
853 if (inverted_eigenvalues)
854 dst[c] *= inverted_eigenvalues[c];
858 eval.template apply<1, false, false>(S1, dst, dst);
859 eval.template apply<0, false, false>(S0, dst, dst);
867 eval.template apply<0, true, false>(S0, src, dst);
868 eval.template apply<1, true, false>(S1, dst, dst);
869 eval.template apply<2, true, false>(S2, dst, dst);
871 for (
unsigned int i2 = 0, c = 0; i2 < n_rows_1d; ++i2)
872 for (
unsigned int i1 = 0; i1 < n_rows_1d; ++i1)
873 for (
unsigned int i0 = 0; i0 < n_rows_1d; ++i0, ++c)
874 if (inverted_eigenvalues)
875 dst[c] *= inverted_eigenvalues[c];
880 eval.template apply<2, false, false>(S2, dst, dst);
881 eval.template apply<1, false, false>(S1, dst, dst);
882 eval.template apply<0, false, false>(S0, dst, dst);
891 template <
int n_rows_1d_templated, std::
size_t dim,
typename Number>
893 select_vmult(Number *dst,
896 const unsigned int n_rows_1d,
897 const std::array<const Number *, dim> &mass_matrix,
898 const std::array<const Number *, dim> &derivative_matrix);
902 template <
int n_rows_1d_templated, std::
size_t dim,
typename Number>
904 select_apply_inverse(Number *dst,
906 const unsigned int n_rows_1d,
908 const std::array<const Number *, dim> &
eigenvalues,
909 const Number *inverted_eigenvalues =
nullptr);
914template <
int dim,
typename Number,
int n_rows_1d>
919 for (
unsigned int d = 1;
d < dim; ++
d)
920 m *= mass_matrix[d].n_rows();
926template <
int dim,
typename Number,
int n_rows_1d>
931 for (
unsigned int d = 1;
d < dim; ++
d)
932 n *= mass_matrix[d].n_cols();
938template <
int dim,
typename Number,
int n_rows_1d>
944 std::lock_guard<std::mutex> lock(this->mutex);
945 this->vmult(dst_view, src_view, this->tmp_array);
950template <
int dim,
typename Number,
int n_rows_1d>
960 Number *dst = dst_view.
begin();
961 const Number *src = src_view.
begin();
963 std::array<const Number *, dim>
mass_matrix, derivative_matrix;
965 for (
unsigned int d = 0;
d < dim; ++
d)
968 derivative_matrix[
d] = &this->derivative_matrix[
d](0, 0);
971 const unsigned int n_rows_1d_non_templated = this->mass_matrix[0].n_rows();
974 internal::TensorProductMatrixSymmetricSum::vmult<
975 n_rows_1d == -1 ? 0 : n_rows_1d>(dst,
978 n_rows_1d_non_templated,
982 internal::TensorProductMatrixSymmetricSum::select_vmult<1>(
986 n_rows_1d_non_templated,
993template <
int dim,
typename Number,
int n_rows_1d>
1002 Number *dst = dst_view.
begin();
1003 const Number *src = src_view.
begin();
1007 for (
unsigned int d = 0;
d < dim; ++
d)
1013 const unsigned int n_rows_1d_non_templated = this->mass_matrix[0].n_rows();
1015 if (n_rows_1d != -1)
1016 internal::TensorProductMatrixSymmetricSum::apply_inverse<
1017 n_rows_1d == -1 ? 0 : n_rows_1d>(
1020 internal::TensorProductMatrixSymmetricSum::select_apply_inverse<1>(
1026template <
int dim,
typename Number,
int n_rows_1d>
1040template <
int dim,
typename Number,
int n_rows_1d>
1041template <
typename T>
1044 const T &derivative_matrix)
1046 reinit(mass_matrix, derivative_matrix);
1051template <
int dim,
typename Number,
int n_rows_1d>
1052template <
typename T>
1055 const T &mass_matrix,
1056 const T &derivative_matrix)
1059 internal::TensorProductMatrixSymmetricSum::convert<dim>(mass_matrix);
1060 this->derivative_matrix =
1061 internal::TensorProductMatrixSymmetricSum::convert<dim>(derivative_matrix);
1063 internal::TensorProductMatrixSymmetricSum::setup(this->mass_matrix,
1064 this->derivative_matrix,
1071template <
int dim,
typename Number,
int n_rows_1d>
1074 const bool precompute_inverse_diagonal)
1075 : compress_matrices(compress_matrices)
1076 , precompute_inverse_diagonal(precompute_inverse_diagonal)
1081template <
int dim,
typename Number,
int n_rows_1d>
1084 const AdditionalData &additional_data)
1085 : compress_matrices(additional_data.compress_matrices)
1086 , precompute_inverse_diagonal(additional_data.precompute_inverse_diagonal)
1091template <
int dim,
typename Number,
int n_rows_1d>
1094 const unsigned int size)
1096 if (compress_matrices ==
false)
1097 mass_and_derivative_matrices.resize(
size * dim);
1104template <
int dim,
typename Number,
int n_rows_1d>
1105template <
typename T>
1108 const unsigned int index,
1113 internal::TensorProductMatrixSymmetricSum::convert<dim>(Ms_in);
1115 internal::TensorProductMatrixSymmetricSum::convert<dim>(Ks_in);
1117 for (
unsigned int d = 0;
d < dim; ++
d)
1119 if (compress_matrices ==
false)
1121 const MatrixPairType
matrix(Ms[d], Ks[d]);
1122 mass_and_derivative_matrices[
index * dim +
d] =
matrix;
1126 using VectorizedArrayTrait =
1129 std::bitset<VectorizedArrayTrait::width()>
mask;
1131 for (
unsigned int v = 0; v < VectorizedArrayTrait::width(); ++v)
1133 typename VectorizedArrayTrait::value_type a = 0.0;
1135 for (
unsigned int i = 0; i < Ms[
d].size(0); ++i)
1136 for (
unsigned int j = 0; j < Ms[
d].size(1); ++j)
1138 a +=
std::abs(VectorizedArrayTrait::get(Ms[d][i][j], v));
1139 a +=
std::abs(VectorizedArrayTrait::get(Ks[d][i][j], v));
1142 mask[v] = (a != 0.0);
1145 const MatrixPairTypeWithMask
matrix{
mask, {Ms[
d], Ks[
d]}};
1147 const auto ptr = cache.find(matrix);
1149 if (ptr != cache.end())
1151 const auto ptr_index = ptr->second;
1152 indices[
index * dim +
d] = ptr_index;
1155 for (
unsigned int v = 0; v < VectorizedArrayTrait::width();
1157 if ((mask[v] ==
true) && (ptr->first.first[v] ==
false))
1167 auto mask_new = ptr->first.first;
1168 auto Ms_new = ptr->first.second.first;
1169 auto Ks_new = ptr->first.second.second;
1171 for (
unsigned int v = 0; v < VectorizedArrayTrait::width();
1173 if (mask_new[v] ==
false && mask[v] ==
true)
1177 for (
unsigned int i = 0; i < Ms_new.size(0); ++i)
1178 for (
unsigned int j = 0; j < Ms_new.size(1); ++j)
1180 VectorizedArrayTrait::get(Ms_new[i][j], v) =
1181 VectorizedArrayTrait::get(Ms[d][i][j], v);
1182 VectorizedArrayTrait::get(Ks_new[i][j], v) =
1183 VectorizedArrayTrait::get(Ks[d][i][j], v);
1189 const MatrixPairTypeWithMask entry_new{mask_new,
1192 const auto ptr_ = cache.find(entry_new);
1195 cache[entry_new] = ptr_index;
1200 const auto size = cache.size();
1210template <
int dim,
typename Number,
int n_rows_1d>
1214 const auto store = [&](
const unsigned int index,
1215 const MatrixPairType &M_and_K) {
1219 std::array<Table<2, Number>, 1> derivative_matrix;
1220 derivative_matrix[0] = M_and_K.second;
1225 internal::TensorProductMatrixSymmetricSum::setup(mass_matrix,
1230 for (
unsigned int i = 0, m = matrix_ptr[index], v = vector_ptr[index];
1234 for (
unsigned int j = 0; j <
mass_matrix[0].n_cols(); ++j, ++m)
1237 this->derivative_matrices[m] = derivative_matrix[0][i][j];
1245 if (compress_matrices ==
false)
1252 this->vector_ptr.resize(mass_and_derivative_matrices.size() + 1);
1253 this->matrix_ptr.resize(mass_and_derivative_matrices.size() + 1);
1255 for (
unsigned int i = 0; i < mass_and_derivative_matrices.size(); ++i)
1257 const auto &M = mass_and_derivative_matrices[i].first;
1259 this->vector_ptr[i + 1] = M.n_rows();
1260 this->matrix_ptr[i + 1] = M.n_rows() * M.n_cols();
1263 for (
unsigned int i = 0; i < mass_and_derivative_matrices.size(); ++i)
1265 this->vector_ptr[i + 1] += this->vector_ptr[i];
1266 this->matrix_ptr[i + 1] += this->matrix_ptr[i];
1269 this->mass_matrices.resize_fast(matrix_ptr.back());
1270 this->derivative_matrices.resize_fast(matrix_ptr.back());
1271 this->eigenvectors.resize_fast(matrix_ptr.back());
1272 this->eigenvalues.resize_fast(vector_ptr.back());
1274 for (
unsigned int i = 0; i < mass_and_derivative_matrices.size(); ++i)
1275 store(i, mass_and_derivative_matrices[i]);
1277 mass_and_derivative_matrices.clear();
1279 else if (cache.size() == indices.size())
1283 this->vector_ptr.resize(cache.size() + 1);
1284 this->matrix_ptr.resize(cache.size() + 1);
1286 std::map<unsigned int, MatrixPairType> inverted_cache;
1288 for (
const auto &i : cache)
1291 for (
unsigned int i = 0; i < indices.size(); ++i)
1293 const auto &M = inverted_cache[indices[i]].first;
1295 this->vector_ptr[i + 1] = M.n_rows();
1296 this->matrix_ptr[i + 1] = M.n_rows() * M.n_cols();
1299 for (
unsigned int i = 0; i < cache.size(); ++i)
1301 this->vector_ptr[i + 1] += this->vector_ptr[i];
1302 this->matrix_ptr[i + 1] += this->matrix_ptr[i];
1305 this->mass_matrices.resize_fast(matrix_ptr.back());
1306 this->derivative_matrices.resize_fast(matrix_ptr.back());
1307 this->eigenvectors.resize_fast(matrix_ptr.back());
1308 this->eigenvalues.resize_fast(vector_ptr.back());
1310 for (
unsigned int i = 0; i < indices.size(); ++i)
1311 store(i, inverted_cache[indices[i]]);
1320 this->vector_ptr.resize(cache.size() + 1);
1321 this->matrix_ptr.resize(cache.size() + 1);
1323 for (
const auto &i : cache)
1325 const auto &M = i.first.second.first;
1327 this->vector_ptr[i.second + 1] = M.n_rows();
1328 this->matrix_ptr[i.second + 1] = M.n_rows() * M.n_cols();
1331 for (
unsigned int i = 0; i < cache.size(); ++i)
1333 this->vector_ptr[i + 1] += this->vector_ptr[i];
1334 this->matrix_ptr[i + 1] += this->matrix_ptr[i];
1337 this->mass_matrices.resize_fast(matrix_ptr.back());
1338 this->derivative_matrices.resize_fast(matrix_ptr.back());
1339 this->eigenvectors.resize_fast(matrix_ptr.back());
1340 this->eigenvalues.resize_fast(vector_ptr.back());
1342 for (
const auto &i : cache)
1348 if (precompute_inverse_diagonal)
1353 for (
unsigned int i = 0; i < this->eigenvalues.size(); ++i)
1354 this->eigenvalues[i] = Number(1.0) / this->eigenvalues[i];
1355 std::swap(this->inverted_eigenvalues,
eigenvalues);
1365 std::vector<unsigned int> indices_ev;
1367 if (indices.size() > 0)
1370 const unsigned int n_cells = indices.size() / dim;
1371 std::map<std::array<unsigned int, dim>,
unsigned int> cache_ev;
1372 std::vector<unsigned int> cache_ev_idx(n_cells);
1374 for (
unsigned int i = 0, c = 0; i <
n_cells; ++i)
1376 std::array<unsigned int, dim> id;
1378 for (
unsigned int d = 0;
d < dim; ++
d, ++c)
1381 const auto id_ptr = cache_ev.find(
id);
1383 if (id_ptr == cache_ev.end())
1385 const auto size = cache_ev.size();
1386 cache_ev_idx[i] =
size;
1387 cache_ev[id] =
size;
1391 cache_ev_idx[i] = id_ptr->second;
1396 std::vector<unsigned int> new_indices;
1397 new_indices.reserve(indices.size() / dim * (dim + 1));
1399 for (
unsigned int i = 0, c = 0; i <
n_cells; ++i)
1401 for (
unsigned int d = 0;
d < dim; ++
d, ++c)
1402 new_indices.push_back(indices[c]);
1403 new_indices.push_back(cache_ev_idx[i]);
1407 indices_ev.resize(cache_ev.size() * dim);
1408 for (
const auto &entry : cache_ev)
1409 for (unsigned
int d = 0;
d < dim; ++
d)
1410 indices_ev[entry.second * dim + d] = entry.first[d];
1412 std::swap(this->indices, new_indices);
1416 const unsigned int n_diag =
1417 ((indices_ev.size() > 0) ? indices_ev.size() :
1418 (matrix_ptr.size() - 1)) /
1421 std::vector<unsigned int> new_vector_ptr(n_diag + 1, 0);
1422 std::vector<unsigned int> new_vector_n_rows_1d(n_diag, 0);
1424 for (
unsigned int i = 0; i < n_diag; ++i)
1426 const unsigned int c = (indices_ev.size() > 0) ?
1427 indices_ev[dim * i + 0] :
1430 const unsigned int n_rows = vector_ptr[c + 1] - vector_ptr[c];
1432 new_vector_n_rows_1d[i] = n_rows;
1436 for (
unsigned int i = 0; i < n_diag; ++i)
1437 new_vector_ptr[i + 1] += new_vector_ptr[i];
1439 this->inverted_eigenvalues.resize(new_vector_ptr.back());
1442 for (
unsigned int i = 0; i < n_diag; ++i)
1444 std::array<Number *, dim> evs;
1446 for (
unsigned int d = 0;
d < dim; ++
d)
1449 ->
eigenvalues[this->vector_ptr[(indices_ev.size() > 0) ?
1450 indices_ev[dim * i +
d] :
1453 const unsigned int mm = new_vector_n_rows_1d[i];
1456 for (
unsigned int i1 = 0, c = 0; i1 < mm; ++i1)
1457 for (
unsigned int i0 = 0; i0 < mm; ++i0, ++c)
1458 this->inverted_eigenvalues[new_vector_ptr[i] + c] =
1459 Number(1.0) / (evs[1][i1] + evs[0][i0]);
1463 for (
unsigned int i2 = 0, c = 0; i2 < mm; ++i2)
1464 for (
unsigned int i1 = 0; i1 < mm; ++i1)
1465 for (
unsigned int i0 = 0; i0 < mm; ++i0, ++c)
1466 this->inverted_eigenvalues[new_vector_ptr[i] + c] =
1467 Number(1.0) / (evs[2][i2] + evs[1][i1] + evs[0][i0]);
1472 std::swap(this->vector_ptr, new_vector_ptr);
1473 std::swap(this->vector_n_rows_1d, new_vector_n_rows_1d);
1476 this->eigenvalues.clear();
1482template <
int dim,
typename Number,
int n_rows_1d>
1489 Number *dst = dst_in.
begin();
1490 const Number *src = src_in.
begin();
1492 if (this->eigenvalues.empty() ==
false)
1496 unsigned int n_rows_1d_non_templated = 0;
1498 for (
unsigned int d = 0;
d < dim; ++
d)
1500 const unsigned int translated_index =
1501 (indices.size() > 0) ? indices[dim * index + d] : (dim *
index +
d);
1504 this->eigenvectors.data() + matrix_ptr[translated_index];
1506 this->eigenvalues.data() + vector_ptr[translated_index];
1507 n_rows_1d_non_templated =
1508 vector_ptr[translated_index + 1] - vector_ptr[translated_index];
1511 if (n_rows_1d != -1)
1512 internal::TensorProductMatrixSymmetricSum::apply_inverse<
1513 n_rows_1d == -1 ? 0 : n_rows_1d>(
1516 internal::TensorProductMatrixSymmetricSum::select_apply_inverse<1>(
1522 const Number *inverted_eigenvalues =
nullptr;
1523 unsigned int n_rows_1d_non_templated = 0;
1525 for (
unsigned int d = 0;
d < dim; ++
d)
1527 const unsigned int translated_index =
1528 (indices.size() > 0) ?
1529 indices[((dim == 1) ? 1 : (dim + 1)) *
index +
d] :
1533 this->eigenvectors.data() + matrix_ptr[translated_index];
1537 const unsigned int translated_index =
1538 ((indices.size() > 0) && (dim != 1)) ?
1539 indices[(dim + 1) *
index + dim] :
1542 inverted_eigenvalues =
1543 this->inverted_eigenvalues.data() + vector_ptr[translated_index];
1544 n_rows_1d_non_templated =
1546 (vector_ptr[translated_index + 1] - vector_ptr[translated_index]) :
1547 vector_n_rows_1d[translated_index];
1550 if (n_rows_1d != -1)
1551 internal::TensorProductMatrixSymmetricSum::apply_inverse<
1552 n_rows_1d == -1 ? 0 : n_rows_1d>(dst,
1554 n_rows_1d_non_templated,
1557 inverted_eigenvalues);
1559 internal::TensorProductMatrixSymmetricSum::select_apply_inverse<1>(
1562 n_rows_1d_non_templated,
1565 inverted_eigenvalues);
1571template <
int dim,
typename Number,
int n_rows_1d>
1587template <
int dim,
typename Number,
int n_rows_1d>
1592 if (matrix_ptr.empty())
1595 return matrix_ptr.size() - 1;
void resize_fast(const size_type new_size)
std::complex< typename numbers::NumberTraits< number >::real_type > eigenvalue(const size_type i) const
void compute_generalized_eigenvalues_symmetric(LAPACKFullMatrix< number > &B, const number lower_bound, const number upper_bound, const number abs_accuracy, Vector< number > &eigenvalues, std::vector< Vector< number > > &eigenvectors, const types::blas_int itype=1)
std::size_t storage_size() const
AlignedVector< Number > mass_matrices
const bool precompute_inverse_diagonal
std::size_t memory_consumption() const
std::vector< unsigned int > matrix_ptr
void apply_inverse(const unsigned int index, const ArrayView< Number > &dst_in, const ArrayView< const Number > &src_in) const
AlignedVector< Number > eigenvectors
std::vector< unsigned int > vector_ptr
const bool compress_matrices
std::vector< MatrixPairType > mass_and_derivative_matrices
std::pair< std::bitset<::internal::VectorizedArrayTrait< Number >::width()>, MatrixPairType > MatrixPairTypeWithMask
std::pair< Table< 2, Number >, Table< 2, Number > > MatrixPairType
void reserve(const unsigned int size)
void insert(const unsigned int index, const T &Ms, const T &Ks)
AlignedVector< Number > inverted_eigenvalues
std::vector< unsigned int > indices
std::vector< unsigned int > vector_n_rows_1d
std::map< MatrixPairTypeWithMask, unsigned int, internal::TensorProductMatrixSymmetricSum::MatrixPairComparator< Number > > cache
AlignedVector< Number > eigenvalues
AlignedVector< Number > derivative_matrices
TensorProductMatrixSymmetricSumCollection(const AdditionalData &additional_data=AdditionalData())
void vmult(const ArrayView< Number > &dst, const ArrayView< const Number > &src, AlignedVector< Number > &tmp) const
std::array< Table< 2, Number >, dim > eigenvectors
std::array< Table< 2, Number >, dim > derivative_matrix
void reinit(const T &mass_matrix, const T &derivative_matrix)
void apply_inverse(const ArrayView< Number > &dst, const ArrayView< const Number > &src) const
void vmult(const ArrayView< Number > &dst, const ArrayView< const Number > &src) const
AlignedVector< Number > tmp_array
static constexpr int n_rows_1d_static
std::size_t memory_consumption() const
TensorProductMatrixSymmetricSum()=default
std::array< Table< 2, Number >, dim > mass_matrix
std::array< AlignedVector< Number >, dim > eigenvalues
TensorProductMatrixSymmetricSum(const T &mass_matrix, const T &derivative_matrix)
static constexpr std::size_t size()
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcInternalError()
#define AssertThrow(cond, exc)
#define DEAL_II_NOT_IMPLEMENTED()
@ matrix
Contents is actually a matrix.
void mass_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const double factor=1.)
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
constexpr T pow(const T base, const int iexp)
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
constexpr unsigned int invalid_unsigned_int
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
bool precompute_inverse_diagonal
AdditionalData(const bool compress_matrices=true, const bool precompute_inverse_diagonal=true)
typename VectorizedArrayTrait::value_type ScalarNumber
bool operator()(const MatrixPairType &left, const MatrixPairType &right) const
std::pair< std::bitset< width >, std::pair< Table< 2, Number >, Table< 2, Number > > > MatrixPairType
static constexpr std::size_t width
static constexpr std::size_t width()
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
std::array< std::pair< Number, Tensor< 1, dim, Number > >, dim > eigenvectors(const SymmetricTensor< 2, dim, Number > &T, const SymmetricTensorEigenvectorMethod method=SymmetricTensorEigenvectorMethod::ql_implicit_shifts)