16 #ifndef dealii_tensor_product_matrix_h
17 #define dealii_tensor_product_matrix_h
115 template <
int dim,
typename Number,
int n_rows_1d = -1>
140 template <
typename T>
161 template <
typename T>
281 template <
typename Number>
290 std::pair<std::bitset<width>,
300 const auto &M_0 = left.second.first;
301 const auto &K_0 = left.second.second;
302 const auto &M_1 = right.second.first;
303 const auto &K_1 = right.second.second;
305 std::bitset<width>
mask;
307 for (
unsigned int v = 0; v <
width; ++v)
308 mask[v] = left.first[v] && right.first[v];
313 if (comparator(M_0, M_1))
315 else if (comparator(M_1, M_0))
317 else if (comparator(K_0, K_1))
354 template <
int dim,
typename Number,
int n_rows_1d = -1>
360 std::bitset<::internal::VectorizedArrayTrait<Number>::width()>,
404 template <
typename T>
406 insert(
const unsigned int index,
const T &Ms,
const T &Ks);
535 template <
typename Number>
538 const Number * derivative_matrix,
539 const unsigned int n_rows,
540 const unsigned int n_cols,
546 std::vector<bool> constrained_dofs(n_rows,
false);
548 for (
unsigned int i = 0; i < n_rows; ++i)
552 Assert(derivative_matrix[i + i * n_rows] == 0.0,
555 for (
unsigned int j = 0; j < n_rows; ++j)
557 Assert(derivative_matrix[i + j * n_rows] == 0,
559 Assert(derivative_matrix[j + i * n_rows] == 0,
563 constrained_dofs[i] =
true;
567 const auto transpose_fill_nm = [&constrained_dofs](Number * out,
569 const unsigned int n,
570 const unsigned int m) {
571 for (
unsigned int mm = 0, c = 0; mm < m; ++mm)
572 for (
unsigned int nn = 0; nn < n; ++nn, ++c)
574 (mm == nn && constrained_dofs[mm]) ? Number(1.0) : in[c];
577 std::vector<::Vector<Number>> eigenvecs(n_rows);
581 transpose_fill_nm(&(mass_copy(0, 0)),
mass_matrix, n_rows, n_cols);
582 transpose_fill_nm(&(deriv_copy(0, 0)), derivative_matrix, n_rows, n_cols);
587 for (
unsigned int i = 0, c = 0; i < n_rows; ++i)
588 for (
unsigned int j = 0; j < n_cols; ++j, ++c)
589 if (constrained_dofs[i] ==
false)
592 for (
unsigned int i = 0; i < n_rows; ++i, ++
eigenvalues)
598 template <std::
size_t dim,
typename Number>
605 const unsigned int n_rows_1d =
mass_matrix[0].n_cols();
608 for (
unsigned int dir = 0; dir < dim; ++dir)
613 derivative_matrix[dir].n_rows());
615 derivative_matrix[dir].n_cols());
620 internal::TensorProductMatrixSymmetricSum::spectral_assembly<Number>(
622 &(derivative_matrix[dir](0, 0)),
632 template <std::
size_t dim,
typename Number, std::
size_t n_lanes>
643 const unsigned int n_rows_1d =
mass_matrix[0].n_cols();
644 constexpr
unsigned int macro_size =
646 const std::size_t nm_flat_size_max = n_rows_1d * n_rows_1d * macro_size;
647 const std::size_t n_flat_size_max = n_rows_1d * macro_size;
649 std::vector<Number> mass_matrix_flat;
650 std::vector<Number> deriv_matrix_flat;
651 std::vector<Number> eigenvalues_flat;
652 std::vector<Number> eigenvectors_flat;
653 mass_matrix_flat.resize(nm_flat_size_max);
654 deriv_matrix_flat.resize(nm_flat_size_max);
655 eigenvalues_flat.resize(n_flat_size_max);
656 eigenvectors_flat.resize(nm_flat_size_max);
657 std::array<unsigned int, macro_size> offsets_nm;
658 std::array<unsigned int, macro_size> offsets_n;
659 for (
unsigned int dir = 0; dir < dim; ++dir)
664 derivative_matrix[dir].n_rows());
666 derivative_matrix[dir].n_cols());
668 const unsigned int n_rows =
mass_matrix[dir].n_rows();
669 const unsigned int n_cols =
mass_matrix[dir].n_cols();
670 const unsigned int nm = n_rows * n_cols;
671 for (
unsigned int vv = 0; vv < macro_size; ++vv)
672 offsets_nm[vv] = nm * vv;
678 mass_matrix_flat.data());
681 &(derivative_matrix[dir](0, 0)),
683 deriv_matrix_flat.data());
685 const Number *mass_cbegin = mass_matrix_flat.data();
686 const Number *deriv_cbegin = deriv_matrix_flat.data();
687 Number * eigenvec_begin = eigenvectors_flat.data();
688 Number * eigenval_begin = eigenvalues_flat.data();
689 for (
unsigned int lane = 0; lane < macro_size; ++lane)
690 internal::TensorProductMatrixSymmetricSum::spectral_assembly<
691 Number>(mass_cbegin + nm * lane,
692 deriv_cbegin + nm * lane,
695 eigenval_begin + n_rows * lane,
696 eigenvec_begin + nm * lane);
700 for (
unsigned int vv = 0; vv < macro_size; ++vv)
701 offsets_n[vv] = n_rows * vv;
703 eigenvalues_flat.data(),
707 eigenvectors_flat.data(),
715 template <std::
size_t dim,
typename Number>
716 inline std::array<Table<2, Number>, dim>
724 template <std::
size_t dim,
typename Number>
725 inline std::array<Table<2, Number>, dim>
728 std::array<Table<2, Number>, dim> mass_copy;
742 template <std::
size_t dim,
typename Number>
743 inline std::array<Table<2, Number>, dim>
746 std::array<Table<2, Number>, dim> matrices;
748 std::fill(matrices.begin(), matrices.end(),
matrix);
755 template <
int n_rows_1d_templated, std::
size_t dim,
typename Number>
760 const unsigned int n_rows_1d_non_templated,
761 const std::array<const Number *, dim> &
mass_matrix,
762 const std::array<const Number *, dim> &derivative_matrix)
764 const unsigned int n_rows_1d = n_rows_1d_templated == 0 ?
765 n_rows_1d_non_templated :
767 const unsigned int n = Utilities::fixed_power<dim>(n_rows_1d);
770 Number *t = tmp.
begin();
777 eval({}, {}, {}, n_rows_1d, n_rows_1d);
781 const Number *
A = derivative_matrix[0];
782 eval.template apply<0, false, false>(
A, src, dst);
787 const Number *A0 = derivative_matrix[0];
789 const Number *A1 = derivative_matrix[1];
791 eval.template apply<0, false, false>(M0, src, t);
792 eval.template apply<1, false, false>(A1, t, dst);
793 eval.template apply<0, false, false>(A0, src, t);
794 eval.template apply<1, false, true>(M1, t, dst);
799 const Number *A0 = derivative_matrix[0];
801 const Number *A1 = derivative_matrix[1];
803 const Number *A2 = derivative_matrix[2];
805 eval.template apply<0, false, false>(M0, src, t + n);
806 eval.template apply<1, false, false>(M1, t + n, t);
807 eval.template apply<2, false, false>(A2, t, dst);
808 eval.template apply<1, false, false>(A1, t + n, t);
809 eval.template apply<0, false, false>(A0, src, t + n);
810 eval.template apply<1, false, true>(M1, t + n, t);
811 eval.template apply<2, false, true>(M2, t, dst);
820 template <
int n_rows_1d_templated, std::
size_t dim,
typename Number>
822 apply_inverse(Number * dst,
825 const unsigned int n_rows_1d_non_templated,
827 const std::array<const Number *, dim> &
eigenvalues,
828 const Number *inverted_eigenvalues =
nullptr)
830 const unsigned int n_rows_1d = n_rows_1d_templated == 0 ?
831 n_rows_1d_non_templated :
833 const unsigned int n = Utilities::fixed_power<dim>(n_rows_1d);
836 Number *t = tmp.
begin();
843 eval({}, {}, {}, n_rows_1d, n_rows_1d);
853 eval.template apply<0, true, false>(S, src, t);
855 for (
unsigned int i = 0; i < n_rows_1d; ++i)
856 if (inverted_eigenvalues)
857 t[i] *= inverted_eigenvalues[i];
861 eval.template apply<0, false, false>(S, t, dst);
868 eval.template apply<0, true, false>(S0, src, t);
869 eval.template apply<1, true, false>(S1, t, dst);
871 for (
unsigned int i1 = 0, c = 0; i1 < n_rows_1d; ++i1)
872 for (
unsigned int i0 = 0; i0 < n_rows_1d; ++i0, ++c)
873 if (inverted_eigenvalues)
874 dst[c] *= inverted_eigenvalues[c];
878 eval.template apply<0, false, false>(S0, dst, t);
879 eval.template apply<1, false, false>(S1, t, dst);
887 eval.template apply<0, true, false>(S0, src, t);
888 eval.template apply<1, true, false>(S1, t, dst);
889 eval.template apply<2, true, false>(S2, dst, t);
891 for (
unsigned int i2 = 0, c = 0; i2 < n_rows_1d; ++i2)
892 for (
unsigned int i1 = 0; i1 < n_rows_1d; ++i1)
893 for (
unsigned int i0 = 0; i0 < n_rows_1d; ++i0, ++c)
894 if (inverted_eigenvalues)
895 t[c] *= inverted_eigenvalues[c];
900 eval.template apply<0, false, false>(S0, t, dst);
901 eval.template apply<1, false, false>(S1, dst, t);
902 eval.template apply<2, false, false>(S2, t, dst);
911 template <
int n_rows_1d_templated, std::
size_t dim,
typename Number>
913 select_vmult(Number * dst,
916 const unsigned int n_rows_1d,
917 const std::array<const Number *, dim> &
mass_matrix,
918 const std::array<const Number *, dim> &derivative_matrix);
922 template <
int n_rows_1d_templated, std::
size_t dim,
typename Number>
924 select_apply_inverse(Number * dst,
927 const unsigned int n_rows_1d,
929 const std::array<const Number *, dim> &
eigenvalues,
930 const Number *inverted_eigenvalues =
nullptr);
935 template <
int dim,
typename Number,
int n_rows_1d>
940 for (
unsigned int d = 1;
d < dim; ++
d)
947 template <
int dim,
typename Number,
int n_rows_1d>
952 for (
unsigned int d = 1;
d < dim; ++
d)
959 template <
int dim,
typename Number,
int n_rows_1d>
965 std::lock_guard<std::mutex> lock(this->mutex);
966 this->vmult(dst_view, src_view, this->tmp_array);
971 template <
int dim,
typename Number,
int n_rows_1d>
981 Number * dst = dst_view.
begin();
982 const Number *src = src_view.
begin();
984 std::array<const Number *, dim>
mass_matrix, derivative_matrix;
986 for (
unsigned int d = 0;
d < dim; ++
d)
989 derivative_matrix[
d] = &this->derivative_matrix[
d](0, 0);
992 const unsigned int n_rows_1d_non_templated = this->mass_matrix[0].n_rows();
995 internal::TensorProductMatrixSymmetricSum::vmult<
996 n_rows_1d == -1 ? 0 : n_rows_1d>(dst,
999 n_rows_1d_non_templated,
1003 internal::TensorProductMatrixSymmetricSum::select_vmult<1>(
1007 n_rows_1d_non_templated,
1014 template <
int dim,
typename Number,
int n_rows_1d>
1020 std::lock_guard<std::mutex> lock(this->mutex);
1021 this->apply_inverse(dst_view, src_view, this->tmp_array);
1026 template <
int dim,
typename Number,
int n_rows_1d>
1036 Number * dst = dst_view.
begin();
1037 const Number *src = src_view.
begin();
1041 for (
unsigned int d = 0;
d < dim; ++
d)
1047 const unsigned int n_rows_1d_non_templated = this->mass_matrix[0].n_rows();
1049 if (n_rows_1d != -1)
1050 internal::TensorProductMatrixSymmetricSum::apply_inverse<
1051 n_rows_1d == -1 ? 0 : n_rows_1d>(
1054 internal::TensorProductMatrixSymmetricSum::select_apply_inverse<1>(
1060 template <
int dim,
typename Number,
int n_rows_1d>
1074 template <
int dim,
typename Number,
int n_rows_1d>
1075 template <
typename T>
1078 const T &derivative_matrix)
1085 template <
int dim,
typename Number,
int n_rows_1d>
1086 template <
typename T>
1090 const T &derivative_matrix)
1093 internal::TensorProductMatrixSymmetricSum::convert<dim>(
mass_matrix);
1094 this->derivative_matrix =
1095 internal::TensorProductMatrixSymmetricSum::convert<dim>(derivative_matrix);
1097 internal::TensorProductMatrixSymmetricSum::setup(this->mass_matrix,
1098 this->derivative_matrix,
1105 template <
int dim,
typename Number,
int n_rows_1d>
1108 const bool precompute_inverse_diagonal)
1109 : compress_matrices(compress_matrices)
1110 , precompute_inverse_diagonal(precompute_inverse_diagonal)
1115 template <
int dim,
typename Number,
int n_rows_1d>
1118 const AdditionalData &additional_data)
1119 : compress_matrices(additional_data.compress_matrices)
1120 , precompute_inverse_diagonal(additional_data.precompute_inverse_diagonal)
1125 template <
int dim,
typename Number,
int n_rows_1d>
1128 const unsigned int size)
1130 if (compress_matrices ==
false)
1131 mass_and_derivative_matrices.resize(size * dim);
1138 template <
int dim,
typename Number,
int n_rows_1d>
1139 template <
typename T>
1142 const unsigned int index,
1147 internal::TensorProductMatrixSymmetricSum::convert<dim>(Ms_in);
1149 internal::TensorProductMatrixSymmetricSum::convert<dim>(Ks_in);
1151 for (
unsigned int d = 0;
d < dim; ++
d)
1153 if (compress_matrices ==
false)
1155 const MatrixPairType
matrix(Ms[
d], Ks[
d]);
1156 mass_and_derivative_matrices[
index * dim +
d] =
matrix;
1160 using VectorizedArrayTrait =
1163 std::bitset<VectorizedArrayTrait::width()>
mask;
1165 for (
unsigned int v = 0; v < VectorizedArrayTrait::width(); ++v)
1167 typename VectorizedArrayTrait::value_type a = 0.0;
1169 for (
unsigned int i = 0; i < Ms[
d].size(0); ++i)
1170 for (
unsigned int j = 0; j < Ms[
d].size(1); ++j)
1172 a +=
std::abs(VectorizedArrayTrait::get(Ms[
d][i][j], v));
1173 a +=
std::abs(VectorizedArrayTrait::get(Ks[
d][i][j], v));
1176 mask[v] = (a != 0.0);
1179 const MatrixPairTypeWithMask
matrix{
mask, {Ms[
d], Ks[
d]}};
1181 const auto ptr = cache.find(
matrix);
1183 if (ptr != cache.end())
1185 const auto ptr_index = ptr->second;
1186 indices[
index * dim +
d] = ptr_index;
1189 for (
unsigned int v = 0; v < VectorizedArrayTrait::width();
1191 if ((mask[v] ==
true) && (ptr->first.first[v] ==
false))
1201 auto mask_new = ptr->first.first;
1202 auto Ms_new = ptr->first.second.first;
1203 auto Ks_new = ptr->first.second.second;
1205 for (
unsigned int v = 0; v < VectorizedArrayTrait::width();
1207 if (mask_new[v] ==
false && mask[v] ==
true)
1211 for (
unsigned int i = 0; i < Ms_new.size(0); ++i)
1212 for (
unsigned int j = 0; j < Ms_new.size(1); ++j)
1214 VectorizedArrayTrait::get(Ms_new[i][j], v) =
1215 VectorizedArrayTrait::get(Ms[
d][i][j], v);
1216 VectorizedArrayTrait::get(Ks_new[i][j], v) =
1217 VectorizedArrayTrait::get(Ks[
d][i][j], v);
1223 const MatrixPairTypeWithMask entry_new{mask_new,
1226 const auto ptr_ = cache.find(entry_new);
1229 cache[entry_new] = ptr_index;
1234 const auto size = cache.size();
1235 indices[
index * dim +
d] = size;
1244 template <
int dim,
typename Number,
int n_rows_1d>
1248 const auto store = [&](
const unsigned int index,
1249 const MatrixPairType &M_and_K) {
1253 std::array<Table<2, Number>, 1> derivative_matrix;
1254 derivative_matrix[0] = M_and_K.second;
1259 internal::TensorProductMatrixSymmetricSum::setup(
mass_matrix,
1264 for (
unsigned int i = 0, m = matrix_ptr[index], v = vector_ptr[index];
1268 for (
unsigned int j = 0; j <
mass_matrix[0].n_cols(); ++j, ++m)
1271 this->derivative_matrices[m] = derivative_matrix[0][i][j];
1279 if (compress_matrices ==
false)
1286 this->vector_ptr.resize(mass_and_derivative_matrices.size() + 1);
1287 this->matrix_ptr.resize(mass_and_derivative_matrices.size() + 1);
1289 for (
unsigned int i = 0; i < mass_and_derivative_matrices.size(); ++i)
1291 const auto &M = mass_and_derivative_matrices[i].first;
1293 this->vector_ptr[i + 1] = M.n_rows();
1294 this->matrix_ptr[i + 1] = M.n_rows() * M.n_cols();
1297 for (
unsigned int i = 0; i < mass_and_derivative_matrices.size(); ++i)
1299 this->vector_ptr[i + 1] += this->vector_ptr[i];
1300 this->matrix_ptr[i + 1] += this->matrix_ptr[i];
1303 this->mass_matrices.resize_fast(matrix_ptr.back());
1304 this->derivative_matrices.resize_fast(matrix_ptr.back());
1305 this->eigenvectors.resize_fast(matrix_ptr.back());
1306 this->eigenvalues.resize_fast(vector_ptr.back());
1308 for (
unsigned int i = 0; i < mass_and_derivative_matrices.size(); ++i)
1309 store(i, mass_and_derivative_matrices[i]);
1311 mass_and_derivative_matrices.clear();
1313 else if (cache.size() == indices.size())
1317 this->vector_ptr.resize(cache.size() + 1);
1318 this->matrix_ptr.resize(cache.size() + 1);
1320 std::map<unsigned int, MatrixPairType> inverted_cache;
1322 for (
const auto &i : cache)
1323 inverted_cache[i.second] = i.first.second;
1325 for (
unsigned int i = 0; i < indices.size(); ++i)
1327 const auto &M = inverted_cache[indices[i]].first;
1329 this->vector_ptr[i + 1] = M.n_rows();
1330 this->matrix_ptr[i + 1] = M.n_rows() * M.n_cols();
1333 for (
unsigned int i = 0; i < cache.size(); ++i)
1335 this->vector_ptr[i + 1] += this->vector_ptr[i];
1336 this->matrix_ptr[i + 1] += this->matrix_ptr[i];
1339 this->mass_matrices.resize_fast(matrix_ptr.back());
1340 this->derivative_matrices.resize_fast(matrix_ptr.back());
1341 this->eigenvectors.resize_fast(matrix_ptr.back());
1342 this->eigenvalues.resize_fast(vector_ptr.back());
1344 for (
unsigned int i = 0; i < indices.size(); ++i)
1345 store(i, inverted_cache[indices[i]]);
1354 this->vector_ptr.resize(cache.size() + 1);
1355 this->matrix_ptr.resize(cache.size() + 1);
1357 for (
const auto &i : cache)
1359 const auto &M = i.first.second.first;
1361 this->vector_ptr[i.second + 1] = M.n_rows();
1362 this->matrix_ptr[i.second + 1] = M.n_rows() * M.n_cols();
1365 for (
unsigned int i = 0; i < cache.size(); ++i)
1367 this->vector_ptr[i + 1] += this->vector_ptr[i];
1368 this->matrix_ptr[i + 1] += this->matrix_ptr[i];
1371 this->mass_matrices.resize_fast(matrix_ptr.back());
1372 this->derivative_matrices.resize_fast(matrix_ptr.back());
1373 this->eigenvectors.resize_fast(matrix_ptr.back());
1374 this->eigenvalues.resize_fast(vector_ptr.back());
1376 for (
const auto &i : cache)
1377 store(i.second, i.first.second);
1382 if (precompute_inverse_diagonal)
1387 for (
unsigned int i = 0; i < this->eigenvalues.size(); ++i)
1388 this->eigenvalues[i] = Number(1.0) / this->eigenvalues[i];
1399 std::vector<unsigned int> indices_ev;
1401 if (indices.size() > 0)
1404 const unsigned int n_cells = indices.size() / dim;
1405 std::map<std::array<unsigned int, dim>,
unsigned int> cache_ev;
1406 std::vector<unsigned int> cache_ev_idx(
n_cells);
1408 for (
unsigned int i = 0, c = 0; i <
n_cells; ++i)
1410 std::array<unsigned int, dim> id;
1412 for (
unsigned int d = 0;
d < dim; ++
d, ++c)
1415 const auto id_ptr = cache_ev.find(
id);
1417 if (id_ptr == cache_ev.end())
1419 const auto size = cache_ev.size();
1420 cache_ev_idx[i] = size;
1421 cache_ev[id] = size;
1425 cache_ev_idx[i] = id_ptr->second;
1430 std::vector<unsigned int> new_indices;
1431 new_indices.reserve(indices.size() / dim * (dim + 1));
1433 for (
unsigned int i = 0, c = 0; i <
n_cells; ++i)
1435 for (
unsigned int d = 0;
d < dim; ++
d, ++c)
1436 new_indices.push_back(indices[c]);
1437 new_indices.push_back(cache_ev_idx[i]);
1441 indices_ev.resize(cache_ev.size() * dim);
1442 for (
const auto &entry : cache_ev)
1443 for (
unsigned int d = 0;
d < dim; ++
d)
1444 indices_ev[entry.second * dim +
d] = entry.first[
d];
1450 const unsigned int n_diag =
1451 ((indices_ev.size() > 0) ? indices_ev.size() :
1452 (matrix_ptr.size() - 1)) /
1455 std::vector<unsigned int> new_vector_ptr(n_diag + 1, 0);
1456 std::vector<unsigned int> new_vector_n_rows_1d(n_diag, 0);
1458 for (
unsigned int i = 0; i < n_diag; ++i)
1460 const unsigned int c = (indices_ev.size() > 0) ?
1461 indices_ev[dim * i + 0] :
1464 const unsigned int n_rows = vector_ptr[c + 1] - vector_ptr[c];
1466 new_vector_n_rows_1d[i] = n_rows;
1470 for (
unsigned int i = 0; i < n_diag; ++i)
1471 new_vector_ptr[i + 1] += new_vector_ptr[i];
1473 this->inverted_eigenvalues.resize(new_vector_ptr.back());
1476 for (
unsigned int i = 0; i < n_diag; ++i)
1478 std::array<Number *, dim> evs;
1480 for (
unsigned int d = 0;
d < dim; ++
d)
1483 ->
eigenvalues[this->vector_ptr[(indices_ev.size() > 0) ?
1484 indices_ev[dim * i +
d] :
1487 const unsigned int mm = new_vector_n_rows_1d[i];
1490 for (
unsigned int i1 = 0, c = 0; i1 < mm; ++i1)
1491 for (
unsigned int i0 = 0; i0 < mm; ++i0, ++c)
1492 this->inverted_eigenvalues[new_vector_ptr[i] + c] =
1493 Number(1.0) / (evs[1][i1] + evs[0][i0]);
1497 for (
unsigned int i2 = 0, c = 0; i2 < mm; ++i2)
1498 for (
unsigned int i1 = 0; i1 < mm; ++i1)
1499 for (
unsigned int i0 = 0; i0 < mm; ++i0, ++c)
1500 this->inverted_eigenvalues[new_vector_ptr[i] + c] =
1501 Number(1.0) / (evs[2][i2] + evs[1][i1] + evs[0][i0]);
1506 std::swap(this->vector_ptr, new_vector_ptr);
1507 std::swap(this->vector_n_rows_1d, new_vector_n_rows_1d);
1510 this->eigenvalues.clear();
1516 template <
int dim,
typename Number,
int n_rows_1d>
1524 Number * dst = dst_in.
begin();
1525 const Number *src = src_in.
begin();
1527 if (this->eigenvalues.empty() ==
false)
1531 unsigned int n_rows_1d_non_templated = 0;
1533 for (
unsigned int d = 0;
d < dim; ++
d)
1535 const unsigned int translated_index =
1536 (indices.size() > 0) ? indices[dim * index +
d] : (dim * index +
d);
1539 this->eigenvectors.data() + matrix_ptr[translated_index];
1541 this->eigenvalues.data() + vector_ptr[translated_index];
1542 n_rows_1d_non_templated =
1543 vector_ptr[translated_index + 1] - vector_ptr[translated_index];
1546 if (n_rows_1d != -1)
1547 internal::TensorProductMatrixSymmetricSum::apply_inverse<
1548 n_rows_1d == -1 ? 0 : n_rows_1d>(dst,
1551 n_rows_1d_non_templated,
1555 internal::TensorProductMatrixSymmetricSum::select_apply_inverse<1>(
1559 n_rows_1d_non_templated,
1566 const Number * inverted_eigenvalues =
nullptr;
1567 unsigned int n_rows_1d_non_templated = 0;
1569 for (
unsigned int d = 0;
d < dim; ++
d)
1571 const unsigned int translated_index =
1572 (indices.size() > 0) ?
1573 indices[((dim == 1) ? 1 : (dim + 1)) *
index +
d] :
1577 this->eigenvectors.data() + matrix_ptr[translated_index];
1581 const unsigned int translated_index =
1582 ((indices.size() > 0) && (dim != 1)) ?
1583 indices[(dim + 1) *
index + dim] :
1586 inverted_eigenvalues =
1587 this->inverted_eigenvalues.data() + vector_ptr[translated_index];
1588 n_rows_1d_non_templated =
1590 (vector_ptr[translated_index + 1] - vector_ptr[translated_index]) :
1591 vector_n_rows_1d[translated_index];
1594 if (n_rows_1d != -1)
1595 internal::TensorProductMatrixSymmetricSum::apply_inverse<
1596 n_rows_1d == -1 ? 0 : n_rows_1d>(dst,
1599 n_rows_1d_non_templated,
1602 inverted_eigenvalues);
1604 internal::TensorProductMatrixSymmetricSum::select_apply_inverse<1>(
1608 n_rows_1d_non_templated,
1611 inverted_eigenvalues);
1617 template <
int dim,
typename Number,
int n_rows_1d>
1633 template <
int dim,
typename Number,
int n_rows_1d>
1638 if (matrix_ptr.size() == 0)
1641 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::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)
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
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
void apply_inverse(const unsigned int index, const ArrayView< Number > &dst_in, const ArrayView< const Number > &src_in, AlignedVector< Number > &tmp_array) const
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
void apply_inverse(const ArrayView< Number > &dst, const ArrayView< const Number > &src, AlignedVector< Number > &tmp) const
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 & ExcInternalError()
#define Assert(cond, exc)
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
#define AssertThrow(cond, exc)
@ matrix
Contents is actually a matrix.
@ eigenvalues
Eigenvalue vector is filled.
void mass_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const double factor=1.)
std::enable_if_t< std::is_fundamental< T >::value, std::size_t > memory_consumption(const T &t)
void swap(MemorySpaceData< T, MemorySpace > &u, MemorySpaceData< T, MemorySpace > &v)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > epsilon(const Tensor< 2, dim, Number > &Grad_u)
VectorType::value_type * begin(VectorType &V)
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
void transform(const InputIterator &begin_in, const InputIterator &end_in, OutputIterator out, const Function &function, const unsigned int grainsize)
::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< 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)