Kronecker Product addition to stdlib_linalg

src/CMakeLists.txt

@@ -22,6 +22,7 @@ set(fppFiles
     stdlib_linalg.fypp
     stdlib_linalg_diag.fypp
     stdlib_linalg_outer_product.fypp
+    stdlib_linalg_kronecker.fypp
     stdlib_linalg_cross_product.fypp
     stdlib_optval.fypp
     stdlib_selection.fypp

src/stdlib_linalg.fypp

@@ -14,6 +14,7 @@ module stdlib_linalg
   public :: eye
   public :: trace
   public :: outer_product
+  public :: kronecker_product
   public :: cross_product
   public :: is_square
   public :: is_diagonal
@@ -93,6 +94,20 @@ module stdlib_linalg
     #:endfor
   end interface outer_product
 
+  interface kronecker_product
+    !! version: experimental
+    !!
+    !! Computes the Kronecker product of two arrays size M1xN1, M2xN2, returning an (M1*M2)x(N1*N2) array
+    !! ([Specification](../page/specs/stdlib_linalg.html#
+    !! kronecker_product-computes-the-kronecker-product-of-two-matrices))
+    #:for k1, t1 in RCI_KINDS_TYPES
+      pure module function kronecker_product_${t1[0]}$${k1}$(A, B) result(C)
+        ${t1}$, intent(in) :: A(:,:), B(:,:)
+        ${t1}$ :: C(size(A,dim=1)*size(B,dim=1),size(A,dim=2)*size(B,dim=2))
+      end function kronecker_product_${t1[0]}$${k1}$
+    #:endfor
+  end interface kronecker_product
+
 
   ! Cross product (of two vectors)
   interface cross_product

src/stdlib_linalg_kronecker.fypp

@@ -0,0 +1,28 @@
+#:include "common.fypp"
+#:set RCI_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES + INT_KINDS_TYPES
+submodule (stdlib_linalg) stdlib_linalg_kronecker
+
+  implicit none
+
+contains
+
+  #:for k1, t1 in RCI_KINDS_TYPES
+    pure module function kronecker_product_${t1[0]}$${k1}$(A, B) result(C)
+      ${t1}$, intent(in) :: A(:,:), B(:,:)
+      ${t1}$ :: C(size(A,dim=1)*size(B,dim=1),size(A,dim=2)*size(B,dim=2))
+      integer :: m1, n1, m2, n2, maxM1, maxN1, maxM2, maxN2
+
+      maxM1 = size(A, dim=1)
+      maxN1 = size(A, dim=2)
+      maxM2 = size(B, dim=1)
+      maxN2 = size(B, dim=2)
+
+      do n1=1, maxN1
+         do m1=1, maxM1
+            ! We use the numpy.kron convention for ordering of the matrix elements
+            C((m1-1)*maxM2+1:m1*maxM2, (n1-1)*maxN2+1:n1*maxN2) = A(m1, n1) * B(:,:)
+         end do
+      end do
+    end function kronecker_product_${t1[0]}$${k1}$
+  #:endfor
+end submodule stdlib_linalg_kronecker