]> https://gitweb.dealii.org/ - dealii.git/commitdiff
add kronecker product to fullmatrix
authorMichał Wichrowski <mtwichrowski@gmail.com>
Thu, 1 May 2025 16:24:45 +0000 (18:24 +0200)
committerMichał Wichrowski <mtwichrowski@gmail.com>
Thu, 1 May 2025 16:24:45 +0000 (18:24 +0200)
include/deal.II/lac/full_matrix.h
include/deal.II/lac/full_matrix.templates.h

index 1be35fe6185c8a19632cf2104ae4e3bb0782ef8a..2f9ae851be178d2c081693e245c91a10e9e15e88 100644 (file)
@@ -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<number> &A,
+                    const FullMatrix<number> &B,
+                    const bool                adding = false);
+
   /**
    * Matrix-vector-multiplication.
    *
index bf163eb8e4a1630688ded798f74d695fb6d0d507..e9e160cd247cb7f577243ff255c69b05338c0c24 100644 (file)
@@ -942,6 +942,32 @@ FullMatrix<number>::triple_product(const FullMatrix<number> &A,
 }
 
 
+template <typename number>
+void
+FullMatrix<number>::kronecker_product(const FullMatrix<number> &A,
+                                      const FullMatrix<number> &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 <typename number>
 template <typename number2>
 number2

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.