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>
297 const auto &M_0 = left.first;
298 const auto &K_0 = left.second;
299 const auto &M_1 = right.first;
300 const auto &K_1 = right.second;
302 std::bitset<width>
mask;
304 for (
unsigned int v = 0; v <
width; ++v)
309 for (
unsigned int i = 0; i < M_0.size(0); ++i)
310 for (
unsigned int j = 0; j < M_0.size(1); ++j)
318 mask[v] = (a != 0.0) && (
b != 0.0);
324 if (comparator(M_0, M_1))
326 else if (comparator(M_1, M_0))
328 else if (comparator(K_0, K_1))
365 template <
int dim,
typename Number,
int n_rows_1d = -1>
405 template <
typename T>
407 insert(
const unsigned int index,
const T &Ms,
const T &Ks);
521 template <
typename Number>
524 const Number * derivative_matrix,
525 const unsigned int n_rows,
526 const unsigned int n_cols,
532 std::vector<bool> constrained_dofs(n_rows,
false);
534 for (
unsigned int i = 0; i < n_rows; ++i)
538 Assert(derivative_matrix[i + i * n_rows] == 0.0,
541 for (
unsigned int j = 0; j < n_rows; ++j)
543 Assert(derivative_matrix[i + j * n_rows] == 0,
545 Assert(derivative_matrix[j + i * n_rows] == 0,
549 constrained_dofs[i] =
true;
553 const auto transpose_fill_nm = [&constrained_dofs](Number * out,
555 const unsigned int n,
556 const unsigned int m) {
557 for (
unsigned int mm = 0, c = 0; mm < m; ++mm)
558 for (
unsigned int nn = 0; nn < n; ++nn, ++c)
560 (mm == nn && constrained_dofs[mm]) ? Number(1.0) : in[c];
563 std::vector<::Vector<Number>> eigenvecs(n_rows);
567 transpose_fill_nm(&(mass_copy(0, 0)),
mass_matrix, n_rows, n_cols);
568 transpose_fill_nm(&(deriv_copy(0, 0)), derivative_matrix, n_rows, n_cols);
573 for (
unsigned int i = 0, c = 0; i < n_rows; ++i)
574 for (
unsigned int j = 0; j < n_cols; ++j, ++c)
575 if (constrained_dofs[i] ==
false)
578 for (
unsigned int i = 0; i < n_rows; ++i, ++
eigenvalues)
584 template <std::
size_t dim,
typename Number>
591 const unsigned int n_rows_1d =
mass_matrix[0].n_cols();
594 for (
unsigned int dir = 0; dir < dim; ++dir)
599 derivative_matrix[dir].n_rows());
601 derivative_matrix[dir].n_cols());
606 internal::TensorProductMatrixSymmetricSum::spectral_assembly<Number>(
608 &(derivative_matrix[dir](0, 0)),
618 template <std::
size_t dim,
typename Number, std::
size_t n_lanes>
629 const unsigned int n_rows_1d =
mass_matrix[0].n_cols();
630 constexpr
unsigned int macro_size =
632 const std::size_t nm_flat_size_max = n_rows_1d * n_rows_1d * macro_size;
633 const std::size_t n_flat_size_max = n_rows_1d * macro_size;
635 std::vector<Number> mass_matrix_flat;
636 std::vector<Number> deriv_matrix_flat;
637 std::vector<Number> eigenvalues_flat;
638 std::vector<Number> eigenvectors_flat;
639 mass_matrix_flat.resize(nm_flat_size_max);
640 deriv_matrix_flat.resize(nm_flat_size_max);
641 eigenvalues_flat.resize(n_flat_size_max);
642 eigenvectors_flat.resize(nm_flat_size_max);
643 std::array<unsigned int, macro_size> offsets_nm;
644 std::array<unsigned int, macro_size> offsets_n;
645 for (
unsigned int dir = 0; dir < dim; ++dir)
650 derivative_matrix[dir].n_rows());
652 derivative_matrix[dir].n_cols());
654 const unsigned int n_rows =
mass_matrix[dir].n_rows();
655 const unsigned int n_cols =
mass_matrix[dir].n_cols();
656 const unsigned int nm = n_rows * n_cols;
657 for (
unsigned int vv = 0; vv < macro_size; ++vv)
658 offsets_nm[vv] = nm * vv;
664 mass_matrix_flat.data());
667 &(derivative_matrix[dir](0, 0)),
669 deriv_matrix_flat.data());
671 const Number *mass_cbegin = mass_matrix_flat.data();
672 const Number *deriv_cbegin = deriv_matrix_flat.data();
673 Number * eigenvec_begin = eigenvectors_flat.data();
674 Number * eigenval_begin = eigenvalues_flat.data();
675 for (
unsigned int lane = 0; lane < macro_size; ++lane)
676 internal::TensorProductMatrixSymmetricSum::spectral_assembly<
677 Number>(mass_cbegin + nm * lane,
678 deriv_cbegin + nm * lane,
681 eigenval_begin + n_rows * lane,
682 eigenvec_begin + nm * lane);
686 for (
unsigned int vv = 0; vv < macro_size; ++vv)
687 offsets_n[vv] = n_rows * vv;
689 eigenvalues_flat.data(),
693 eigenvectors_flat.data(),
701 template <std::
size_t dim,
typename Number>
702 inline std::array<Table<2, Number>, dim>
710 template <std::
size_t dim,
typename Number>
711 inline std::array<Table<2, Number>, dim>
714 std::array<Table<2, Number>, dim> mass_copy;
728 template <std::
size_t dim,
typename Number>
729 inline std::array<Table<2, Number>, dim>
732 std::array<Table<2, Number>, dim> matrices;
734 std::fill(matrices.begin(), matrices.end(),
matrix);
741 template <
int n_rows_1d_templated, std::
size_t dim,
typename Number>
746 const unsigned int n_rows_1d_non_templated,
747 const std::array<const Number *, dim> &
mass_matrix,
748 const std::array<const Number *, dim> &derivative_matrix)
750 const unsigned int n_rows_1d = n_rows_1d_templated == 0 ?
751 n_rows_1d_non_templated :
753 const unsigned int n = Utilities::fixed_power<dim>(n_rows_1d);
756 Number *t = tmp.
begin();
763 eval({}, {}, {}, n_rows_1d, n_rows_1d);
767 const Number *
A = derivative_matrix[0];
768 eval.template apply<0, false, false>(
A, src, dst);
773 const Number *A0 = derivative_matrix[0];
775 const Number *A1 = derivative_matrix[1];
777 eval.template apply<0, false, false>(M0, src, t);
778 eval.template apply<1, false, false>(A1, t, dst);
779 eval.template apply<0, false, false>(A0, src, t);
780 eval.template apply<1, false, true>(M1, t, dst);
785 const Number *A0 = derivative_matrix[0];
787 const Number *A1 = derivative_matrix[1];
789 const Number *A2 = derivative_matrix[2];
791 eval.template apply<0, false, false>(M0, src, t + n);
792 eval.template apply<1, false, false>(M1, t + n, t);
793 eval.template apply<2, false, false>(A2, t, dst);
794 eval.template apply<1, false, false>(A1, t + n, t);
795 eval.template apply<0, false, false>(A0, src, t + n);
796 eval.template apply<1, false, true>(M1, t + n, t);
797 eval.template apply<2, false, true>(M2, t, dst);
806 template <
int n_rows_1d_templated, std::
size_t dim,
typename Number>
808 apply_inverse(Number * dst,
811 const unsigned int n_rows_1d_non_templated,
813 const std::array<const Number *, dim> &
eigenvalues)
815 const unsigned int n_rows_1d = n_rows_1d_templated == 0 ?
816 n_rows_1d_non_templated :
818 const unsigned int n = Utilities::fixed_power<dim>(n_rows_1d);
821 Number *t = tmp.
begin();
828 eval({}, {}, {}, n_rows_1d, n_rows_1d);
838 eval.template apply<0, true, false>(S, src, t);
839 for (
unsigned int i = 0; i < n_rows_1d; ++i)
841 eval.template apply<0, false, false>(S, t, dst);
848 eval.template apply<0, true, false>(S0, src, t);
849 eval.template apply<1, true, false>(S1, t, dst);
850 for (
unsigned int i1 = 0, c = 0; i1 < n_rows_1d; ++i1)
851 for (
unsigned int i0 = 0; i0 < n_rows_1d; ++i0, ++c)
853 eval.template apply<0, false, false>(S0, dst, t);
854 eval.template apply<1, false, false>(S1, t, dst);
862 eval.template apply<0, true, false>(S0, src, t);
863 eval.template apply<1, true, false>(S1, t, dst);
864 eval.template apply<2, true, false>(S2, dst, t);
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)
870 eval.template apply<0, false, false>(S0, t, dst);
871 eval.template apply<1, false, false>(S1, dst, t);
872 eval.template apply<2, false, false>(S2, t, dst);
881 template <
int n_rows_1d_templated, std::
size_t dim,
typename Number>
883 select_vmult(Number * dst,
886 const unsigned int n_rows_1d,
887 const std::array<const Number *, dim> &
mass_matrix,
888 const std::array<const Number *, dim> &derivative_matrix);
892 template <
int n_rows_1d_templated, std::
size_t dim,
typename Number>
894 select_apply_inverse(Number * dst,
897 const unsigned int n_rows_1d,
899 const std::array<const Number *, dim> &
eigenvalues);
904 template <
int dim,
typename Number,
int n_rows_1d>
909 for (
unsigned int d = 1;
d < dim; ++
d)
916 template <
int dim,
typename Number,
int n_rows_1d>
921 for (
unsigned int d = 1;
d < dim; ++
d)
928 template <
int dim,
typename Number,
int n_rows_1d>
934 std::lock_guard<std::mutex> lock(this->mutex);
935 this->vmult(dst_view, src_view, this->tmp_array);
940 template <
int dim,
typename Number,
int n_rows_1d>
950 Number * dst = dst_view.
begin();
951 const Number *src = src_view.
begin();
953 std::array<const Number *, dim>
mass_matrix, derivative_matrix;
955 for (
unsigned int d = 0;
d < dim; ++
d)
958 derivative_matrix[
d] = &this->derivative_matrix[
d](0, 0);
961 const unsigned int n_rows_1d_non_templated = this->mass_matrix[0].n_rows();
964 internal::TensorProductMatrixSymmetricSum::vmult<
965 n_rows_1d == -1 ? 0 : n_rows_1d>(dst,
968 n_rows_1d_non_templated,
972 internal::TensorProductMatrixSymmetricSum::select_vmult<1>(
976 n_rows_1d_non_templated,
983 template <
int dim,
typename Number,
int n_rows_1d>
989 std::lock_guard<std::mutex> lock(this->mutex);
990 this->apply_inverse(dst_view, src_view, this->tmp_array);
995 template <
int dim,
typename Number,
int n_rows_1d>
1005 Number * dst = dst_view.
begin();
1006 const Number *src = src_view.
begin();
1010 for (
unsigned int d = 0;
d < dim; ++
d)
1016 const unsigned int n_rows_1d_non_templated = this->mass_matrix[0].n_rows();
1018 if (n_rows_1d != -1)
1019 internal::TensorProductMatrixSymmetricSum::apply_inverse<
1020 n_rows_1d == -1 ? 0 : n_rows_1d>(
1023 internal::TensorProductMatrixSymmetricSum::select_apply_inverse<1>(
1029 template <
int dim,
typename Number,
int n_rows_1d>
1043 template <
int dim,
typename Number,
int n_rows_1d>
1044 template <
typename T>
1047 const T &derivative_matrix)
1054 template <
int dim,
typename Number,
int n_rows_1d>
1055 template <
typename T>
1059 const T &derivative_matrix)
1062 internal::TensorProductMatrixSymmetricSum::convert<dim>(
mass_matrix);
1063 this->derivative_matrix =
1064 internal::TensorProductMatrixSymmetricSum::convert<dim>(derivative_matrix);
1066 internal::TensorProductMatrixSymmetricSum::setup(this->mass_matrix,
1067 this->derivative_matrix,
1074 template <
int dim,
typename Number,
int n_rows_1d>
1077 : compress_matrices(compress_matrices)
1082 template <
int dim,
typename Number,
int n_rows_1d>
1085 const AdditionalData &additional_data)
1086 : compress_matrices(additional_data.compress_matrices)
1091 template <
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);
1104 template <
int dim,
typename Number,
int n_rows_1d>
1105 template <
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 const MatrixPairType
matrix(Ms[
d], Ks[
d]);
1121 if (compress_matrices ==
false)
1123 mass_and_derivative_matrices[
index * dim +
d] =
matrix;
1127 const auto ptr = cache.find(
matrix);
1129 if (ptr != cache.end())
1130 indices[
index * dim +
d] = ptr->second;
1133 const auto size = cache.size();
1134 indices[
index * dim +
d] = size;
1143 template <
int dim,
typename Number,
int n_rows_1d>
1147 const auto store = [&](
const unsigned int index,
1148 const MatrixPairType &M_and_K) {
1152 std::array<Table<2, Number>, 1> derivative_matrix;
1153 derivative_matrix[0] = M_and_K.second;
1158 internal::TensorProductMatrixSymmetricSum::setup(
mass_matrix,
1163 for (
unsigned int i = 0, m = matrix_ptr[index], v = vector_ptr[index];
1167 for (
unsigned int j = 0; j <
mass_matrix[0].n_cols(); ++j, ++m)
1170 this->derivative_matrices[m] = derivative_matrix[0][i][j];
1178 if (compress_matrices ==
false)
1185 this->vector_ptr.resize(mass_and_derivative_matrices.size() + 1);
1186 this->matrix_ptr.resize(mass_and_derivative_matrices.size() + 1);
1188 for (
unsigned int i = 0; i < mass_and_derivative_matrices.size(); ++i)
1190 const auto &M = mass_and_derivative_matrices[i].first;
1192 this->vector_ptr[i + 1] = M.n_rows();
1193 this->matrix_ptr[i + 1] = M.n_rows() * M.n_cols();
1196 for (
unsigned int i = 0; i < mass_and_derivative_matrices.size(); ++i)
1198 this->vector_ptr[i + 1] += this->vector_ptr[i];
1199 this->matrix_ptr[i + 1] += this->matrix_ptr[i];
1202 this->mass_matrices.resize_fast(matrix_ptr.back());
1203 this->derivative_matrices.resize_fast(matrix_ptr.back());
1204 this->eigenvectors.resize_fast(matrix_ptr.back());
1205 this->eigenvalues.resize_fast(vector_ptr.back());
1207 for (
unsigned int i = 0; i < mass_and_derivative_matrices.size(); ++i)
1208 store(i, mass_and_derivative_matrices[i]);
1210 mass_and_derivative_matrices.clear();
1212 else if (cache.size() == indices.size())
1216 this->vector_ptr.resize(cache.size() + 1);
1217 this->matrix_ptr.resize(cache.size() + 1);
1219 std::map<unsigned int, MatrixPairType> inverted_cache;
1221 for (
const auto &i : cache)
1222 inverted_cache[i.second] = i.first;
1224 for (
unsigned int i = 0; i < indices.size(); ++i)
1226 const auto &M = inverted_cache[indices[i]].first;
1228 this->vector_ptr[i + 1] = M.n_rows();
1229 this->matrix_ptr[i + 1] = M.n_rows() * M.n_cols();
1232 for (
unsigned int i = 0; i < cache.size(); ++i)
1234 this->vector_ptr[i + 1] += this->vector_ptr[i];
1235 this->matrix_ptr[i + 1] += this->matrix_ptr[i];
1238 this->mass_matrices.resize_fast(matrix_ptr.back());
1239 this->derivative_matrices.resize_fast(matrix_ptr.back());
1240 this->eigenvectors.resize_fast(matrix_ptr.back());
1241 this->eigenvalues.resize_fast(vector_ptr.back());
1243 for (
unsigned int i = 0; i < indices.size(); ++i)
1244 store(i, inverted_cache[indices[i]]);
1253 this->vector_ptr.resize(cache.size() + 1);
1254 this->matrix_ptr.resize(cache.size() + 1);
1256 for (
const auto &i : cache)
1258 const auto &M = i.first.first;
1260 this->vector_ptr[i.second + 1] = M.n_rows();
1261 this->matrix_ptr[i.second + 1] = M.n_rows() * M.n_cols();
1264 for (
unsigned int i = 0; i < cache.size(); ++i)
1266 this->vector_ptr[i + 1] += this->vector_ptr[i];
1267 this->matrix_ptr[i + 1] += this->matrix_ptr[i];
1270 this->mass_matrices.resize_fast(matrix_ptr.back());
1271 this->derivative_matrices.resize_fast(matrix_ptr.back());
1272 this->eigenvectors.resize_fast(matrix_ptr.back());
1273 this->eigenvalues.resize_fast(vector_ptr.back());
1275 for (
const auto &i : cache)
1276 store(i.second, i.first);
1284 template <
int dim,
typename Number,
int n_rows_1d>
1292 Number * dst = dst_in.
begin();
1293 const Number *src = src_in.
begin();
1296 unsigned int n_rows_1d_non_templated = 0;
1298 for (
unsigned int d = 0;
d < dim; ++
d)
1300 const unsigned int translated_index =
1301 (indices.size() > 0) ? indices[dim * index +
d] : (dim * index +
d);
1304 this->eigenvectors.data() + matrix_ptr[translated_index];
1305 eigenvalues[
d] = this->eigenvalues.data() + vector_ptr[translated_index];
1306 n_rows_1d_non_templated =
1307 vector_ptr[translated_index + 1] - vector_ptr[translated_index];
1310 if (n_rows_1d != -1)
1311 internal::TensorProductMatrixSymmetricSum::apply_inverse<
1312 n_rows_1d == -1 ? 0 : n_rows_1d>(
1315 internal::TensorProductMatrixSymmetricSum::select_apply_inverse<1>(
1321 template <
int dim,
typename Number,
int n_rows_1d>
1337 template <
int dim,
typename Number,
int n_rows_1d>
1342 if (matrix_ptr.size() == 0)
1345 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
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
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
std::map< MatrixPairType, unsigned int, internal::TensorProductMatrixSymmetricSum::MatrixPairComparator< Number > > cache
void reserve(const unsigned int size)
void insert(const unsigned int index, const T &Ms, const T &Ks)
std::vector< unsigned int > indices
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)
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)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
VectorType::value_type * begin(VectorType &V)
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 Predicate &predicate, const unsigned int grainsize)
AdditionalData(const bool compress_matrices=true)
std::pair< Table< 2, Number >, Table< 2, Number > > MatrixPairType
typename VectorizedArrayTrait::value_type ScalarNumber
bool operator()(const MatrixPairType &left, const MatrixPairType &right) const
static constexpr std::size_t width
static constexpr std::size_t width
static Number & get(Number &value, unsigned int c)
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)