From: Michał Wichrowski Date: Thu, 1 May 2025 16:24:45 +0000 (+0200) Subject: add kronecker product to fullmatrix X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=f69a3d5c9ac85ba8b39f9115d9e685a053ebead8;p=dealii.git add kronecker product to fullmatrix --- diff --git a/include/deal.II/lac/full_matrix.h b/include/deal.II/lac/full_matrix.h index 1be35fe618..2f9ae851be 100644 --- a/include/deal.II/lac/full_matrix.h +++ b/include/deal.II/lac/full_matrix.h @@ -975,6 +975,33 @@ public: const bool transpose_D = false, const number scaling = number(1.)); + + /** + * @brief Compute the Kronecker product of two matrices. + * + * This function computes the Kronecker product of two matrices A and B, and + * stores the result in the current matrix. The Kronecker product of an + * m x n matrix A and a p x q matrix B is an (m*p) x (n*q) matrix defined as: + * + * ``` + * A ⊗ B = | a11*B a12*B ... a1n*B | + * | a21*B a22*B ... a2n*B | + * | ... ... ... ... | + * | am1*B am2*B ... amn*B | + * ``` + * + * where aij are the elements of the matrix A. + * + * @param A The first matrix (m x n). + * @param B The second matrix (p x q). + * @param adding If `true`, the result is added to the current matrix. If + * `false` (default), the current matrix is overwritten. + */ + void + kronecker_product(const FullMatrix &A, + const FullMatrix &B, + const bool adding = false); + /** * Matrix-vector-multiplication. * diff --git a/include/deal.II/lac/full_matrix.templates.h b/include/deal.II/lac/full_matrix.templates.h index bf163eb8e4..e9e160cd24 100644 --- a/include/deal.II/lac/full_matrix.templates.h +++ b/include/deal.II/lac/full_matrix.templates.h @@ -942,6 +942,32 @@ FullMatrix::triple_product(const FullMatrix &A, } +template +void +FullMatrix::kronecker_product(const FullMatrix &A, + const FullMatrix &B, + const bool adding) +{ + Assert(!A.empty(), ExcEmptyMatrix()); + Assert(!B.empty(), ExcEmptyMatrix()); + + const size_type m = A.m() * B.m(); + const size_type n = A.n() * B.n(); + + if (adding) + { + AssertDimension(m, this->m()); + AssertDimension(n, this->n()); + } + else + this->reinit(m, n); + + for (size_type i = 0; i < m; ++i) + for (size_type j = 0; j < n; ++j) + (*this)(i, j) += A(i / B.m(), j / B.n()) * B(i % B.m(), j % B.n()); +} + + template template number2