15#ifndef dealii_tensor_product_matrix_h
16#define dealii_tensor_product_matrix_h
111template <
int dim,
typename Number,
int n_rows_1d = -1>
136 template <
typename T>
157 template <
typename T>
260 template <
typename Number>
269 std::pair<std::bitset<width>,
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;
284 std::bitset<width>
mask;
286 for (
unsigned int v = 0; v <
width; ++v)
287 mask[v] = left.first[v] && right.first[v];
292 if (comparator(M_0, M_1))
294 else if (comparator(M_1, M_0))
296 else if (comparator(K_0, K_1))
333template <
int dim,
typename Number,
int n_rows_1d = -1>
339 std::bitset<::internal::VectorizedArrayTrait<Number>::width()>,
383 template <
typename T>
385 insert(
const unsigned int index,
const T &Ms,
const T &Ks);
513 template <
typename Number>
515 spectral_assembly(
const Number *mass_matrix,
516 const Number *derivative_matrix,
517 const unsigned int n_rows,
518 const unsigned int n_cols,
524 std::vector<bool> constrained_dofs(n_rows,
false);
526 for (
unsigned int i = 0; i < n_rows; ++i)
528 if (mass_matrix[i + i * n_rows] == 0.0)
530 Assert(derivative_matrix[i + i * n_rows] == 0.0,
533 for (
unsigned int j = 0; j < n_rows; ++j)
535 Assert(derivative_matrix[i + j * n_rows] == 0,
537 Assert(derivative_matrix[j + i * n_rows] == 0,
541 constrained_dofs[i] =
true;
545 const auto transpose_fill_nm = [&constrained_dofs](Number *out,
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)
552 (mm == nn && constrained_dofs[mm]) ? Number(1.0) : in[c];
555 std::vector<::Vector<Number>> eigenvecs(n_rows);
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);
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)
570 for (
unsigned int i = 0; i < n_rows; ++i, ++
eigenvalues)
576 template <std::
size_t dim,
typename Number>
583 const unsigned int n_rows_1d = mass_matrix[0].n_cols();
586 for (
unsigned int dir = 0; dir < dim; ++dir)
591 derivative_matrix[dir].n_rows());
593 derivative_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(),
610 template <std::
size_t dim,
typename Number, std::
size_t n_lanes>
621 const unsigned int n_rows_1d = mass_matrix[0].n_cols();
622 constexpr unsigned int macro_size =
624 const std::size_t nm_flat_size_max = n_rows_1d * n_rows_1d * macro_size;
625 const std::size_t n_flat_size_max = n_rows_1d * macro_size;
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)
642 derivative_matrix[dir].n_rows());
644 derivative_matrix[dir].n_cols());
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;
654 &(mass_matrix[dir](0, 0)),
656 mass_matrix_flat.data());
659 &(derivative_matrix[dir](0, 0)),
661 deriv_matrix_flat.data());
663 const Number *mass_cbegin = mass_matrix_flat.data();
664 const Number *deriv_cbegin = deriv_matrix_flat.data();
665 Number *eigenvec_begin = eigenvectors_flat.data();
666 Number *eigenval_begin = eigenvalues_flat.data();
667 for (
unsigned int lane = 0; lane < macro_size; ++lane)
668 internal::TensorProductMatrixSymmetricSum::spectral_assembly<
669 Number>(mass_cbegin + nm * lane,
670 deriv_cbegin + nm * lane,
673 eigenval_begin + n_rows * lane,
674 eigenvec_begin + nm * lane);
678 for (
unsigned int vv = 0; vv < macro_size; ++vv)
679 offsets_n[vv] = n_rows * vv;
681 eigenvalues_flat.data(),
685 eigenvectors_flat.data(),
693 template <std::
size_t dim,
typename Number>
694 inline std::array<Table<2, Number>, dim>
702 template <std::
size_t dim,
typename Number>
703 inline std::array<Table<2, Number>, dim>
706 std::array<Table<2, Number>, dim> mass_copy;
708 std::transform(mass_matrix.cbegin(),
720 template <std::
size_t dim,
typename Number>
721 inline std::array<Table<2, Number>, dim>
724 std::array<Table<2, Number>, dim> matrices;
726 std::fill(matrices.begin(), matrices.end(), matrix);
733 template <
int n_rows_1d_templated, std::
size_t dim,
typename Number>
738 const unsigned int n_rows_1d_non_templated,
739 const std::array<const Number *, dim> &mass_matrix,
740 const std::array<const Number *, dim> &derivative_matrix)
742 const unsigned int n_rows_1d = n_rows_1d_templated == 0 ?
743 n_rows_1d_non_templated :
745 const unsigned int n = Utilities::fixed_power<dim>(n_rows_1d);
748 Number *t = tmp.
begin();
755 eval({}, {}, {}, n_rows_1d, n_rows_1d);
759 const Number *A = derivative_matrix[0];
760 eval.template apply<0, false, false>(A, src, dst);
765 const Number *A0 = derivative_matrix[0];
766 const Number *M0 = mass_matrix[0];
767 const Number *A1 = derivative_matrix[1];
768 const Number *M1 = mass_matrix[1];
769 eval.template apply<0, false, false>(M0, src, t);
770 eval.template apply<1, false, false>(A1, t, dst);
771 eval.template apply<0, false, false>(A0, src, t);
772 eval.template apply<1, false, true>(M1, t, dst);
777 const Number *A0 = derivative_matrix[0];
778 const Number *M0 = mass_matrix[0];
779 const Number *A1 = derivative_matrix[1];
780 const Number *M1 = mass_matrix[1];
781 const Number *A2 = derivative_matrix[2];
782 const Number *M2 = mass_matrix[2];
783 eval.template apply<0, false, false>(M0, src, t + n);
784 eval.template apply<1, false, false>(M1, t + n, t);
785 eval.template apply<2, false, false>(A2, t, dst);
786 eval.template apply<1, false, false>(A1, t + n, t);
787 eval.template apply<0, false, false>(A0, src, t + n);
788 eval.template apply<1, false, true>(M1, t + n, t);
789 eval.template apply<2, false, true>(M2, t, dst);
798 template <
int n_rows_1d_templated, std::
size_t dim,
typename Number>
800 apply_inverse(Number *dst,
802 const unsigned int n_rows_1d_non_templated,
804 const std::array<const Number *, dim> &
eigenvalues,
805 const Number *inverted_eigenvalues =
nullptr)
807 const unsigned int n_rows_1d = n_rows_1d_templated == 0 ?
808 n_rows_1d_non_templated :
816 eval({}, {}, {}, n_rows_1d, n_rows_1d);
826 eval.template apply<0, true, false>(S, src, dst);
828 for (
unsigned int i = 0; i < n_rows_1d; ++i)
829 if (inverted_eigenvalues)
830 dst[i] *= inverted_eigenvalues[i];
834 eval.template apply<0, false, false>(S, dst, dst);
841 eval.template apply<0, true, false>(S0, src, dst);
842 eval.template apply<1, true, false>(S1, dst, dst);
844 for (
unsigned int i1 = 0, c = 0; i1 < n_rows_1d; ++i1)
845 for (
unsigned int i0 = 0; i0 < n_rows_1d; ++i0, ++c)
846 if (inverted_eigenvalues)
847 dst[c] *= inverted_eigenvalues[c];
851 eval.template apply<1, false, false>(S1, dst, dst);
852 eval.template apply<0, false, false>(S0, dst, dst);
860 eval.template apply<0, true, false>(S0, src, dst);
861 eval.template apply<1, true, false>(S1, dst, dst);
862 eval.template apply<2, true, false>(S2, dst, dst);
864 for (
unsigned int i2 = 0, c = 0; i2 < n_rows_1d; ++i2)
865 for (
unsigned int i1 = 0; i1 < n_rows_1d; ++i1)
866 for (
unsigned int i0 = 0; i0 < n_rows_1d; ++i0, ++c)
867 if (inverted_eigenvalues)
868 dst[c] *= inverted_eigenvalues[c];
873 eval.template apply<2, false, false>(S2, dst, dst);
874 eval.template apply<1, false, false>(S1, dst, dst);
875 eval.template apply<0, false, false>(S0, dst, dst);
884 template <
int n_rows_1d_templated, std::
size_t dim,
typename Number>
886 select_vmult(Number *dst,
889 const unsigned int n_rows_1d,
890 const std::array<const Number *, dim> &mass_matrix,
891 const std::array<const Number *, dim> &derivative_matrix);
895 template <
int n_rows_1d_templated, std::
size_t dim,
typename Number>
897 select_apply_inverse(Number *dst,
899 const unsigned int n_rows_1d,
901 const std::array<const Number *, dim> &
eigenvalues,
902 const Number *inverted_eigenvalues =
nullptr);
907template <
int dim,
typename Number,
int n_rows_1d>
912 for (
unsigned int d = 1;
d < dim; ++
d)
913 m *= mass_matrix[d].n_rows();
919template <
int dim,
typename Number,
int n_rows_1d>
924 for (
unsigned int d = 1;
d < dim; ++
d)
925 n *= mass_matrix[d].n_cols();
931template <
int dim,
typename Number,
int n_rows_1d>
937 std::lock_guard<std::mutex> lock(this->mutex);
938 this->vmult(dst_view, src_view, this->tmp_array);
943template <
int dim,
typename Number,
int n_rows_1d>
953 Number *dst = dst_view.
begin();
954 const Number *src = src_view.
begin();
956 std::array<const Number *, dim>
mass_matrix, derivative_matrix;
958 for (
unsigned int d = 0;
d < dim; ++
d)
961 derivative_matrix[
d] = &this->derivative_matrix[
d](0, 0);
964 const unsigned int n_rows_1d_non_templated = this->mass_matrix[0].n_rows();
967 internal::TensorProductMatrixSymmetricSum::vmult<
968 n_rows_1d == -1 ? 0 : n_rows_1d>(dst,
971 n_rows_1d_non_templated,
975 internal::TensorProductMatrixSymmetricSum::select_vmult<1>(
979 n_rows_1d_non_templated,
986template <
int dim,
typename Number,
int n_rows_1d>
995 Number *dst = dst_view.
begin();
996 const Number *src = src_view.
begin();
1000 for (
unsigned int d = 0;
d < dim; ++
d)
1006 const unsigned int n_rows_1d_non_templated = this->mass_matrix[0].n_rows();
1008 if (n_rows_1d != -1)
1009 internal::TensorProductMatrixSymmetricSum::apply_inverse<
1010 n_rows_1d == -1 ? 0 : n_rows_1d>(
1013 internal::TensorProductMatrixSymmetricSum::select_apply_inverse<1>(
1019template <
int dim,
typename Number,
int n_rows_1d>
1033template <
int dim,
typename Number,
int n_rows_1d>
1034template <
typename T>
1037 const T &derivative_matrix)
1039 reinit(mass_matrix, derivative_matrix);
1044template <
int dim,
typename Number,
int n_rows_1d>
1045template <
typename T>
1048 const T &mass_matrix,
1049 const T &derivative_matrix)
1052 internal::TensorProductMatrixSymmetricSum::convert<dim>(mass_matrix);
1053 this->derivative_matrix =
1054 internal::TensorProductMatrixSymmetricSum::convert<dim>(derivative_matrix);
1056 internal::TensorProductMatrixSymmetricSum::setup(this->mass_matrix,
1057 this->derivative_matrix,
1064template <
int dim,
typename Number,
int n_rows_1d>
1067 const bool precompute_inverse_diagonal)
1068 : compress_matrices(compress_matrices)
1069 , precompute_inverse_diagonal(precompute_inverse_diagonal)
1074template <
int dim,
typename Number,
int n_rows_1d>
1077 const AdditionalData &additional_data)
1078 : compress_matrices(additional_data.compress_matrices)
1079 , precompute_inverse_diagonal(additional_data.precompute_inverse_diagonal)
1084template <
int dim,
typename Number,
int n_rows_1d>
1087 const unsigned int size)
1089 if (compress_matrices ==
false)
1090 mass_and_derivative_matrices.resize(size * dim);
1097template <
int dim,
typename Number,
int n_rows_1d>
1098template <
typename T>
1101 const unsigned int index,
1106 internal::TensorProductMatrixSymmetricSum::convert<dim>(Ms_in);
1108 internal::TensorProductMatrixSymmetricSum::convert<dim>(Ks_in);
1110 for (
unsigned int d = 0;
d < dim; ++
d)
1112 if (compress_matrices ==
false)
1114 const MatrixPairType
matrix(Ms[d], Ks[d]);
1115 mass_and_derivative_matrices[
index * dim +
d] =
matrix;
1119 using VectorizedArrayTrait =
1122 std::bitset<VectorizedArrayTrait::width()>
mask;
1124 for (
unsigned int v = 0; v < VectorizedArrayTrait::width(); ++v)
1126 typename VectorizedArrayTrait::value_type a = 0.0;
1128 for (
unsigned int i = 0; i < Ms[
d].size(0); ++i)
1129 for (
unsigned int j = 0; j < Ms[
d].size(1); ++j)
1131 a +=
std::abs(VectorizedArrayTrait::get(Ms[d][i][j], v));
1132 a +=
std::abs(VectorizedArrayTrait::get(Ks[d][i][j], v));
1135 mask[v] = (a != 0.0);
1138 const MatrixPairTypeWithMask
matrix{
mask, {Ms[
d], Ks[
d]}};
1140 const auto ptr = cache.find(matrix);
1142 if (ptr != cache.end())
1144 const auto ptr_index = ptr->second;
1145 indices[
index * dim +
d] = ptr_index;
1148 for (
unsigned int v = 0; v < VectorizedArrayTrait::width();
1150 if ((mask[v] ==
true) && (ptr->first.first[v] ==
false))
1160 auto mask_new = ptr->first.first;
1161 auto Ms_new = ptr->first.second.first;
1162 auto Ks_new = ptr->first.second.second;
1164 for (
unsigned int v = 0; v < VectorizedArrayTrait::width();
1166 if (mask_new[v] ==
false && mask[v] ==
true)
1170 for (
unsigned int i = 0; i < Ms_new.size(0); ++i)
1171 for (
unsigned int j = 0; j < Ms_new.size(1); ++j)
1173 VectorizedArrayTrait::get(Ms_new[i][j], v) =
1174 VectorizedArrayTrait::get(Ms[d][i][j], v);
1175 VectorizedArrayTrait::get(Ks_new[i][j], v) =
1176 VectorizedArrayTrait::get(Ks[d][i][j], v);
1182 const MatrixPairTypeWithMask entry_new{mask_new,
1185 const auto ptr_ = cache.find(entry_new);
1188 cache[entry_new] = ptr_index;
1193 const auto size = cache.size();
1194 indices[
index * dim +
d] = size;
1203template <
int dim,
typename Number,
int n_rows_1d>
1207 const auto store = [&](
const unsigned int index,
1208 const MatrixPairType &M_and_K) {
1212 std::array<Table<2, Number>, 1> derivative_matrix;
1213 derivative_matrix[0] = M_and_K.second;
1218 internal::TensorProductMatrixSymmetricSum::setup(mass_matrix,
1223 for (
unsigned int i = 0, m = matrix_ptr[index], v = vector_ptr[index];
1227 for (
unsigned int j = 0; j <
mass_matrix[0].n_cols(); ++j, ++m)
1230 this->derivative_matrices[m] = derivative_matrix[0][i][j];
1238 if (compress_matrices ==
false)
1245 this->vector_ptr.resize(mass_and_derivative_matrices.size() + 1);
1246 this->matrix_ptr.resize(mass_and_derivative_matrices.size() + 1);
1248 for (
unsigned int i = 0; i < mass_and_derivative_matrices.size(); ++i)
1250 const auto &M = mass_and_derivative_matrices[i].first;
1252 this->vector_ptr[i + 1] = M.n_rows();
1253 this->matrix_ptr[i + 1] = M.n_rows() * M.n_cols();
1256 for (
unsigned int i = 0; i < mass_and_derivative_matrices.size(); ++i)
1258 this->vector_ptr[i + 1] += this->vector_ptr[i];
1259 this->matrix_ptr[i + 1] += this->matrix_ptr[i];
1262 this->mass_matrices.resize_fast(matrix_ptr.back());
1263 this->derivative_matrices.resize_fast(matrix_ptr.back());
1264 this->eigenvectors.resize_fast(matrix_ptr.back());
1265 this->eigenvalues.resize_fast(vector_ptr.back());
1267 for (
unsigned int i = 0; i < mass_and_derivative_matrices.size(); ++i)
1268 store(i, mass_and_derivative_matrices[i]);
1270 mass_and_derivative_matrices.clear();
1272 else if (cache.size() == indices.size())
1276 this->vector_ptr.resize(cache.size() + 1);
1277 this->matrix_ptr.resize(cache.size() + 1);
1279 std::map<unsigned int, MatrixPairType> inverted_cache;
1281 for (
const auto &i : cache)
1284 for (
unsigned int i = 0; i < indices.size(); ++i)
1286 const auto &M = inverted_cache[indices[i]].first;
1288 this->vector_ptr[i + 1] = M.n_rows();
1289 this->matrix_ptr[i + 1] = M.n_rows() * M.n_cols();
1292 for (
unsigned int i = 0; i < cache.size(); ++i)
1294 this->vector_ptr[i + 1] += this->vector_ptr[i];
1295 this->matrix_ptr[i + 1] += this->matrix_ptr[i];
1298 this->mass_matrices.resize_fast(matrix_ptr.back());
1299 this->derivative_matrices.resize_fast(matrix_ptr.back());
1300 this->eigenvectors.resize_fast(matrix_ptr.back());
1301 this->eigenvalues.resize_fast(vector_ptr.back());
1303 for (
unsigned int i = 0; i < indices.size(); ++i)
1304 store(i, inverted_cache[indices[i]]);
1313 this->vector_ptr.resize(cache.size() + 1);
1314 this->matrix_ptr.resize(cache.size() + 1);
1316 for (
const auto &i : cache)
1318 const auto &M = i.first.second.first;
1320 this->vector_ptr[i.second + 1] = M.n_rows();
1321 this->matrix_ptr[i.second + 1] = M.n_rows() * M.n_cols();
1324 for (
unsigned int i = 0; i < cache.size(); ++i)
1326 this->vector_ptr[i + 1] += this->vector_ptr[i];
1327 this->matrix_ptr[i + 1] += this->matrix_ptr[i];
1330 this->mass_matrices.resize_fast(matrix_ptr.back());
1331 this->derivative_matrices.resize_fast(matrix_ptr.back());
1332 this->eigenvectors.resize_fast(matrix_ptr.back());
1333 this->eigenvalues.resize_fast(vector_ptr.back());
1335 for (
const auto &i : cache)
1341 if (precompute_inverse_diagonal)
1346 for (
unsigned int i = 0; i < this->eigenvalues.size(); ++i)
1347 this->eigenvalues[i] = Number(1.0) / this->eigenvalues[i];
1348 std::swap(this->inverted_eigenvalues,
eigenvalues);
1358 std::vector<unsigned int> indices_ev;
1360 if (indices.size() > 0)
1363 const unsigned int n_cells = indices.size() / dim;
1364 std::map<std::array<unsigned int, dim>,
unsigned int> cache_ev;
1365 std::vector<unsigned int> cache_ev_idx(n_cells);
1367 for (
unsigned int i = 0, c = 0; i <
n_cells; ++i)
1369 std::array<unsigned int, dim> id;
1371 for (
unsigned int d = 0;
d < dim; ++
d, ++c)
1374 const auto id_ptr = cache_ev.find(
id);
1376 if (id_ptr == cache_ev.end())
1378 const auto size = cache_ev.size();
1379 cache_ev_idx[i] = size;
1380 cache_ev[id] = size;
1384 cache_ev_idx[i] = id_ptr->second;
1389 std::vector<unsigned int> new_indices;
1390 new_indices.reserve(indices.size() / dim * (dim + 1));
1392 for (
unsigned int i = 0, c = 0; i <
n_cells; ++i)
1394 for (
unsigned int d = 0;
d < dim; ++
d, ++c)
1395 new_indices.push_back(indices[c]);
1396 new_indices.push_back(cache_ev_idx[i]);
1400 indices_ev.resize(cache_ev.size() * dim);
1401 for (
const auto &entry : cache_ev)
1403 indices_ev[entry.second * dim + d] = entry.first[d];
1405 std::swap(this->indices, new_indices);
1409 const unsigned int n_diag =
1410 ((indices_ev.size() > 0) ? indices_ev.size() :
1411 (matrix_ptr.size() - 1)) /
1414 std::vector<unsigned int> new_vector_ptr(n_diag + 1, 0);
1415 std::vector<unsigned int> new_vector_n_rows_1d(n_diag, 0);
1417 for (
unsigned int i = 0; i < n_diag; ++i)
1419 const unsigned int c = (indices_ev.size() > 0) ?
1420 indices_ev[dim * i + 0] :
1423 const unsigned int n_rows = vector_ptr[c + 1] - vector_ptr[c];
1425 new_vector_n_rows_1d[i] = n_rows;
1429 for (
unsigned int i = 0; i < n_diag; ++i)
1430 new_vector_ptr[i + 1] += new_vector_ptr[i];
1432 this->inverted_eigenvalues.resize(new_vector_ptr.back());
1435 for (
unsigned int i = 0; i < n_diag; ++i)
1437 std::array<Number *, dim> evs;
1439 for (
unsigned int d = 0;
d < dim; ++
d)
1442 ->
eigenvalues[this->vector_ptr[(indices_ev.size() > 0) ?
1443 indices_ev[dim * i +
d] :
1446 const unsigned int mm = new_vector_n_rows_1d[i];
1449 for (
unsigned int i1 = 0, c = 0; i1 < mm; ++i1)
1450 for (
unsigned int i0 = 0; i0 < mm; ++i0, ++c)
1451 this->inverted_eigenvalues[new_vector_ptr[i] + c] =
1452 Number(1.0) / (evs[1][i1] + evs[0][i0]);
1456 for (
unsigned int i2 = 0, c = 0; i2 < mm; ++i2)
1457 for (
unsigned int i1 = 0; i1 < mm; ++i1)
1458 for (
unsigned int i0 = 0; i0 < mm; ++i0, ++c)
1459 this->inverted_eigenvalues[new_vector_ptr[i] + c] =
1460 Number(1.0) / (evs[2][i2] + evs[1][i1] + evs[0][i0]);
1465 std::swap(this->vector_ptr, new_vector_ptr);
1466 std::swap(this->vector_n_rows_1d, new_vector_n_rows_1d);
1469 this->eigenvalues.clear();
1475template <
int dim,
typename Number,
int n_rows_1d>
1482 Number *dst = dst_in.
begin();
1483 const Number *src = src_in.
begin();
1485 if (this->eigenvalues.empty() ==
false)
1489 unsigned int n_rows_1d_non_templated = 0;
1491 for (
unsigned int d = 0;
d < dim; ++
d)
1493 const unsigned int translated_index =
1494 (indices.size() > 0) ? indices[dim * index + d] : (dim *
index +
d);
1497 this->eigenvectors.data() + matrix_ptr[translated_index];
1499 this->eigenvalues.data() + vector_ptr[translated_index];
1500 n_rows_1d_non_templated =
1501 vector_ptr[translated_index + 1] - vector_ptr[translated_index];
1504 if (n_rows_1d != -1)
1505 internal::TensorProductMatrixSymmetricSum::apply_inverse<
1506 n_rows_1d == -1 ? 0 : n_rows_1d>(
1509 internal::TensorProductMatrixSymmetricSum::select_apply_inverse<1>(
1515 const Number *inverted_eigenvalues =
nullptr;
1516 unsigned int n_rows_1d_non_templated = 0;
1518 for (
unsigned int d = 0;
d < dim; ++
d)
1520 const unsigned int translated_index =
1521 (indices.size() > 0) ?
1522 indices[((dim == 1) ? 1 : (dim + 1)) *
index +
d] :
1526 this->eigenvectors.data() + matrix_ptr[translated_index];
1530 const unsigned int translated_index =
1531 ((indices.size() > 0) && (dim != 1)) ?
1532 indices[(dim + 1) *
index + dim] :
1535 inverted_eigenvalues =
1536 this->inverted_eigenvalues.data() + vector_ptr[translated_index];
1537 n_rows_1d_non_templated =
1539 (vector_ptr[translated_index + 1] - vector_ptr[translated_index]) :
1540 vector_n_rows_1d[translated_index];
1543 if (n_rows_1d != -1)
1544 internal::TensorProductMatrixSymmetricSum::apply_inverse<
1545 n_rows_1d == -1 ? 0 : n_rows_1d>(dst,
1547 n_rows_1d_non_templated,
1550 inverted_eigenvalues);
1552 internal::TensorProductMatrixSymmetricSum::select_apply_inverse<1>(
1555 n_rows_1d_non_templated,
1558 inverted_eigenvalues);
1564template <
int dim,
typename Number,
int n_rows_1d>
1580template <
int dim,
typename Number,
int n_rows_1d>
1585 if (matrix_ptr.empty())
1588 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)
static const 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 > >, std::integral_constant< int, dim >::value > eigenvectors(const SymmetricTensor< 2, dim, Number > &T, const SymmetricTensorEigenvectorMethod method=SymmetricTensorEigenvectorMethod::ql_implicit_shifts)
void vectorized_load_and_transpose(const unsigned int n_entries, const Number *in, const unsigned int *offsets, VectorizedArray< Number, width > *out)
void vectorized_transpose_and_store(const bool add_into, const unsigned int n_entries, const VectorizedArray< Number, width > *in, const unsigned int *offsets, Number *out)