Added the Kronecker product functionality to stdlib_linalg, and added appropriate unit tests for all supported types.

Author: adenchfi

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

test/linalg/test_linalg.fypp

@@ -3,7 +3,7 @@
 module test_linalg
     use testdrive, only : new_unittest, unittest_type, error_type, check, skip_test
     use stdlib_kinds, only: sp, dp, xdp, qp, int8, int16, int32, int64
-    use stdlib_linalg, only: diag, eye, trace, outer_product, cross_product
+    use stdlib_linalg, only: diag, eye, trace, outer_product, cross_product, kronecker_product
 
     implicit none
 
@@ -48,6 +48,16 @@ contains
             new_unittest("trace_int16", test_trace_int16), &
             new_unittest("trace_int32", test_trace_int32), &
             new_unittest("trace_int64", test_trace_int64), &
+            new_unittest("kronecker_product_rsp", test_kronecker_product_rsp), &
+            new_unittest("kronecker_product_rdp", test_kronecker_product_rdp), &
+            new_unittest("kronecker_product_rqp", test_kronecker_product_rqp), &
+            new_unittest("kronecker_product_csp", test_kronecker_product_csp), &
+            new_unittest("kronecker_product_cdp", test_kronecker_product_cdp), &
+            new_unittest("kronecker_product_cqp", test_kronecker_product_cqp), &
+            new_unittest("kronecker_product_int8", test_kronecker_product_iint8), &
+            new_unittest("kronecker_product_int16", test_kronecker_product_iint16), &
+            new_unittest("kronecker_product_int32", test_kronecker_product_iint32), &
+            new_unittest("kronecker_product_int64", test_kronecker_product_iint64), &
             new_unittest("outer_product_rsp", test_outer_product_rsp), &
             new_unittest("outer_product_rdp", test_outer_product_rdp), &
             new_unittest("outer_product_rqp", test_outer_product_rqp), &
@@ -552,6 +562,336 @@ contains
             "trace(h) == sum(c(0:nd:2)) failed.")
 
     end subroutine test_trace_int64
+    
+    subroutine test_kronecker_product_rsp(error) 
+      !> Error handling
+      type(error_type), allocatable, intent(out) :: error
+      integer, parameter :: m1=1, n1=2, m2=2, n2=3
+      real(sp), parameter :: tol = 1.e-6
+
+      real(sp) :: A(m1,n1), B(m2,n2)
+      real(sp) :: C(m1*m2,n1*n2), expected(m1*m2,n1*n2), diff(m1*m2,n1*n2)
+
+      integer :: i,j
+
+      do j=1, n1
+         do i=1, m1
+            A(i,j) = i*j ! A = [1, 2]
+         end do
+      end do
+
+      do j=1, n2
+         do i=1, m2
+            B(i,j) = i*j ! B = [[1, 2, 3], [2, 4, 6]]
+         end do
+      end do
+
+      C = kronecker_product(A,B)
+      expected = transpose(reshape([1,2,3, 2,4,6, 2,4,6, 4, 8, 12], [m2*n2, m1*n1]))
+      diff = C - expected
+
+      call check(error, all(abs(diff) .le. abs(tol)), "all(abs(diff) .le. abs(tol)) failed")
+      ! Expected: C = [1*B, 2*B] = [[1,2,3, 2,4,6], [2,4,6, 4, 8, 12]]
+
+    end subroutine test_kronecker_product_rsp
+    
+    subroutine test_kronecker_product_rdp(error) 
+      !> Error handling
+      type(error_type), allocatable, intent(out) :: error
+      integer, parameter :: m1=1, n1=2, m2=2, n2=3
+      real(dp), parameter :: tol = 1.e-6
+
+      real(dp) :: A(m1,n1), B(m2,n2)
+      real(dp) :: C(m1*m2,n1*n2), expected(m1*m2,n1*n2), diff(m1*m2,n1*n2)
+
+      integer :: i,j
+
+      do j=1, n1
+         do i=1, m1
+            A(i,j) = i*j ! A = [1, 2]
+         end do
+      end do
+
+      do j=1, n2
+         do i=1, m2
+            B(i,j) = i*j ! B = [[1, 2, 3], [2, 4, 6]]
+         end do
+      end do
+
+      C = kronecker_product(A,B)
+      expected = transpose(reshape([1,2,3, 2,4,6, 2,4,6, 4, 8, 12], [m2*n2, m1*n1]))
+      diff = C - expected
+
+      call check(error, all(abs(diff) .le. abs(tol)), "all(abs(diff) .le. abs(tol)) failed")
+      ! Expected: C = [1*B, 2*B] = [[1,2,3, 2,4,6], [2,4,6, 4, 8, 12]]
+
+    end subroutine test_kronecker_product_rdp
+
+    subroutine test_kronecker_product_rqp(error) 
+      !> Error handling
+      type(error_type), allocatable, intent(out) :: error
+#:if WITH_QP
+      integer, parameter :: m1=1, n1=2, m2=2, n2=3
+      real(qp), parameter :: tol = 1.e-6
+
+      real(qp) :: A(m1,n1), B(m2,n2)
+      real(qp) :: C(m1*m2,n1*n2), expected(m1*m2,n1*n2), diff(m1*m2,n1*n2)
+
+      integer :: i,j
+
+      do j=1, n1
+         do i=1, m1
+            A(i,j) = i*j ! A = [1, 2]
+         end do
+      end do
+
+      do j=1, n2
+         do i=1, m2
+            B(i,j) = i*j ! B = [[1, 2, 3], [2, 4, 6]]
+         end do
+      end do
+
+      C = kronecker_product(A,B)
+      expected = transpose(reshape([1,2,3, 2,4,6, 2,4,6, 4, 8, 12], [m2*n2, m1*n1]))
+      diff = C - expected
+
+      call check(error, all(abs(diff) .le. abs(tol)), "all(abs(diff) .le. abs(tol)) failed")
+      ! Expected: C = [1*B, 2*B] = [[1,2,3, 2,4,6], [2,4,6, 4, 8, 12]]
+#:else
+        call skip_test(error, "Quadruple precision is not enabled")
+#:endif
+
+    end subroutine test_kronecker_product_rqp
+    
+    subroutine test_kronecker_product_csp(error) 
+      !> Error handling
+      type(error_type), allocatable, intent(out) :: error
+      integer, parameter :: m1=1, n1=2, m2=2, n2=3
+      complex(sp), parameter :: tol = 1.e-6
+
+      complex(sp) :: A(m1,n1), B(m2,n2)
+      complex(sp) :: C(m1*m2,n1*n2), expected(m1*m2,n1*n2), diff(m1*m2,n1*n2)
+
+      integer :: i,j
+
+      do j=1, n1
+         do i=1, m1
+            A(i,j) = i*j ! A = [1, 2]
+         end do
+      end do
+
+      do j=1, n2
+         do i=1, m2
+            B(i,j) = i*j ! B = [[1, 2, 3], [2, 4, 6]]
+         end do
+      end do
+
+      C = kronecker_product(A,B)
+      expected = transpose(reshape([1,2,3, 2,4,6, 2,4,6, 4, 8, 12], [m2*n2, m1*n1]))
+      diff = C - expected
+
+      call check(error, all(abs(diff) .le. abs(tol)), "all(abs(diff) .le. abs(tol)) failed")
+      ! Expected: C = [1*B, 2*B] = [[1,2,3, 2,4,6], [2,4,6, 4, 8, 12]]
+
+    end subroutine test_kronecker_product_csp
+
+    subroutine test_kronecker_product_cdp(error) 
+      !> Error handling
+      type(error_type), allocatable, intent(out) :: error
+      integer, parameter :: m1=1, n1=2, m2=2, n2=3
+      complex(dp), parameter :: tol = 1.e-6
+
+      complex(dp) :: A(m1,n1), B(m2,n2)
+      complex(dp) :: C(m1*m2,n1*n2), expected(m1*m2,n1*n2), diff(m1*m2,n1*n2)
+
+      integer :: i,j
+
+      do j=1, n1
+         do i=1, m1
+            A(i,j) = i*j ! A = [1, 2]
+         end do
+      end do
+
+      do j=1, n2
+         do i=1, m2
+            B(i,j) = i*j ! B = [[1, 2, 3], [2, 4, 6]]
+         end do
+      end do
+
+      C = kronecker_product(A,B)
+      expected = transpose(reshape([1,2,3, 2,4,6, 2,4,6, 4, 8, 12], [m2*n2, m1*n1]))
+      diff = C - expected
+
+      call check(error, all(abs(diff) .le. abs(tol)), "all(abs(diff) .le. abs(tol)) failed")
+      ! Expected: C = [1*B, 2*B] = [[1,2,3, 2,4,6], [2,4,6, 4, 8, 12]]
+
+    end subroutine test_kronecker_product_cdp
+
+    subroutine test_kronecker_product_cqp(error) 
+      !> Error handling
+      type(error_type), allocatable, intent(out) :: error
+#:if WITH_QP
+      integer, parameter :: m1=1, n1=2, m2=2, n2=3
+      complex(qp), parameter :: tol = 1.e-6
+
+      complex(qp) :: A(m1,n1), B(m2,n2)
+      complex(qp) :: C(m1*m2,n1*n2), expected(m1*m2,n1*n2), diff(m1*m2,n1*n2)
+
+      integer :: i,j
+
+      do j=1, n1
+         do i=1, m1
+            A(i,j) = i*j ! A = [1, 2]
+         end do
+      end do
+
+      do j=1, n2
+         do i=1, m2
+            B(i,j) = i*j ! B = [[1, 2, 3], [2, 4, 6]]
+         end do
+      end do
+
+      C = kronecker_product(A,B)
+      expected = transpose(reshape([1,2,3, 2,4,6, 2,4,6, 4, 8, 12], [m2*n2, m1*n1]))
+      diff = C - expected
+
+      call check(error, all(abs(diff) .le. abs(tol)), "all(abs(diff) .le. abs(tol)) failed")
+      ! Expected: C = [1*B, 2*B] = [[1,2,3, 2,4,6], [2,4,6, 4, 8, 12]]
+#:else
+        call skip_test(error, "Quadruple precision is not enabled")
+#:endif
+
+    end subroutine test_kronecker_product_cqp
+    
+    subroutine test_kronecker_product_iint8(error) 
+      !> Error handling
+      type(error_type), allocatable, intent(out) :: error
+      integer, parameter :: m1=1, n1=2, m2=2, n2=3
+      integer(int8), parameter :: tol = 1.e-6
+
+      integer(int8) :: A(m1,n1), B(m2,n2)
+      integer(int8) :: C(m1*m2,n1*n2), expected(m1*m2,n1*n2), diff(m1*m2,n1*n2)
+
+      integer :: i,j
+
+      do j=1, n1
+         do i=1, m1
+            A(i,j) = i*j ! A = [1, 2]
+         end do
+      end do
+
+      do j=1, n2
+         do i=1, m2
+            B(i,j) = i*j ! B = [[1, 2, 3], [2, 4, 6]]
+         end do
+      end do
+
+      C = kronecker_product(A,B)
+      expected = transpose(reshape([1,2,3, 2,4,6, 2,4,6, 4, 8, 12], [m2*n2, m1*n1]))
+      diff = C - expected
+
+      call check(error, all(abs(diff) .le. abs(tol)), "all(abs(diff) .le. abs(tol)) failed")
+      ! Expected: C = [1*B, 2*B] = [[1,2,3, 2,4,6], [2,4,6, 4, 8, 12]]
+
+    end subroutine test_kronecker_product_iint8
+    
+    subroutine test_kronecker_product_iint16(error) 
+      !> Error handling
+      type(error_type), allocatable, intent(out) :: error
+      integer, parameter :: m1=1, n1=2, m2=2, n2=3
+      integer(int16), parameter :: tol = 1.e-6
+
+      integer(int16) :: A(m1,n1), B(m2,n2)
+      integer(int16) :: C(m1*m2,n1*n2), expected(m1*m2,n1*n2), diff(m1*m2,n1*n2)
+
+      integer :: i,j
+
+      do j=1, n1
+         do i=1, m1
+            A(i,j) = i*j ! A = [1, 2]
+         end do
+      end do
+
+      do j=1, n2
+         do i=1, m2
+            B(i,j) = i*j ! B = [[1, 2, 3], [2, 4, 6]]
+         end do
+      end do
+
+      C = kronecker_product(A,B)
+      expected = transpose(reshape([1,2,3, 2,4,6, 2,4,6, 4, 8, 12], [m2*n2, m1*n1]))
+      diff = C - expected
+
+      call check(error, all(abs(diff) .le. abs(tol)), "all(abs(diff) .le. abs(tol)) failed")
+      ! Expected: C = [1*B, 2*B] = [[1,2,3, 2,4,6], [2,4,6, 4, 8, 12]]
+
+    end subroutine test_kronecker_product_iint16
+    
+    subroutine test_kronecker_product_iint32(error) 
+      !> Error handling
+      type(error_type), allocatable, intent(out) :: error
+      integer, parameter :: m1=1, n1=2, m2=2, n2=3
+      integer(int32), parameter :: tol = 1.e-6
+
+      integer(int32) :: A(m1,n1), B(m2,n2)
+      integer(int32) :: C(m1*m2,n1*n2), expected(m1*m2,n1*n2), diff(m1*m2,n1*n2)
+
+      integer :: i,j
+
+      do j=1, n1
+         do i=1, m1
+            A(i,j) = i*j ! A = [1, 2]
+         end do
+      end do
+
+      do j=1, n2
+         do i=1, m2
+            B(i,j) = i*j ! B = [[1, 2, 3], [2, 4, 6]]
+         end do
+      end do
+
+      C = kronecker_product(A,B)
+      expected = transpose(reshape([1,2,3, 2,4,6, 2,4,6, 4, 8, 12], [m2*n2, m1*n1]))
+      diff = C - expected
+
+      call check(error, all(abs(diff) .le. abs(tol)), "all(abs(diff) .le. abs(tol)) failed")
+      ! Expected: C = [1*B, 2*B] = [[1,2,3, 2,4,6], [2,4,6, 4, 8, 12]]
+
+    end subroutine test_kronecker_product_iint32
+    
+    subroutine test_kronecker_product_iint64(error) 
+      !> Error handling
+      type(error_type), allocatable, intent(out) :: error
+      integer, parameter :: m1=1, n1=2, m2=2, n2=3
+      integer(int64), parameter :: tol = 1.e-6
+
+      integer(int64) :: A(m1,n1), B(m2,n2)
+      integer(int64) :: C(m1*m2,n1*n2), expected(m1*m2,n1*n2), diff(m1*m2,n1*n2)
+
+      integer :: i,j
+
+      do j=1, n1
+         do i=1, m1
+            A(i,j) = i*j ! A = [1, 2]
+         end do
+      end do
+
+      do j=1, n2
+         do i=1, m2
+            B(i,j) = i*j ! B = [[1, 2, 3], [2, 4, 6]]
+         end do
+      end do
+
+      C = kronecker_product(A,B)
+
+      expected = transpose(reshape([1,2,3, 2,4,6, 2,4,6, 4, 8, 12], [m2*n2, m1*n1]))
+      
+      diff = C - expected
+
+      call check(error, all(abs(diff) .le. abs(tol)), "all(abs(diff) .le. abs(tol)) failed")
+      ! Expected: C = [1*B, 2*B] = [[1,2,3, 2,4,6], [2,4,6, 4, 8, 12]]
+
+    end subroutine test_kronecker_product_iint64
 
 
     subroutine test_outer_product_rsp(error)