Deploying to stdlib-fpm from @ 2fdfab4 🚀

example/example_kronecker_product.f90

@@ -0,0 +1,26 @@
+program example_kronecker_product
+  use stdlib_linalg, only: kronecker_product
+  implicit none
+  integer, parameter :: m1 = 1, n1 = 2, m2 = 2, n2 = 3
+  integer :: i, j
+  real :: A(m1, n1), B(m2,n2)
+  real, allocatable :: C(:,:)
+
+  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 = [ 1, 2, 3 ]
+        B(i,j) = i*j !     [ 2, 4, 6 ]
+     end do
+  end do
+
+  C = kronecker_product(A, B)
+  ! C =     [ a(1,1) * B(:,:) | a(1,2) * B(:,:) ]
+  ! or in other words, 
+  ! C =     [  1.00      2.00      3.00      2.00      4.00      6.00  ]
+  !         [  2.00      4.00      6.00      4.00      8.00     12.00  ]
+end program example_kronecker_product

src/stdlib_linalg.f90

@@ -12,6 +12,7 @@ module stdlib_linalg
   public :: eye
   public :: trace
   public :: outer_product
+  public :: kronecker_product
   public :: cross_product
   public :: is_square
   public :: is_diagonal
@@ -240,6 +241,46 @@ pure module function outer_product_iint64(u, v) result(res)
       end function outer_product_iint64
   end interface outer_product
 
+  interface kronecker_product
+    !! version: experimental
+    !!
+    !! Computes the Kronecker product of two arrays of size M1xN1, and of M2xN2, returning an (M1*M2)x(N1*N2) array
+    !! ([Specification](../page/specs/stdlib_linalg.html#
+    !! kronecker_product-computes-the-kronecker-product-of-two-matrices))
+      pure module function kronecker_product_rsp(A, B) result(C)
+        real(sp), intent(in) :: A(:,:), B(:,:)
+        real(sp) :: C(size(A,dim=1)*size(B,dim=1),size(A,dim=2)*size(B,dim=2))
+      end function kronecker_product_rsp
+      pure module function kronecker_product_rdp(A, B) result(C)
+        real(dp), intent(in) :: A(:,:), B(:,:)
+        real(dp) :: C(size(A,dim=1)*size(B,dim=1),size(A,dim=2)*size(B,dim=2))
+      end function kronecker_product_rdp
+      pure module function kronecker_product_csp(A, B) result(C)
+        complex(sp), intent(in) :: A(:,:), B(:,:)
+        complex(sp) :: C(size(A,dim=1)*size(B,dim=1),size(A,dim=2)*size(B,dim=2))
+      end function kronecker_product_csp
+      pure module function kronecker_product_cdp(A, B) result(C)
+        complex(dp), intent(in) :: A(:,:), B(:,:)
+        complex(dp) :: C(size(A,dim=1)*size(B,dim=1),size(A,dim=2)*size(B,dim=2))
+      end function kronecker_product_cdp
+      pure module function kronecker_product_iint8(A, B) result(C)
+        integer(int8), intent(in) :: A(:,:), B(:,:)
+        integer(int8) :: C(size(A,dim=1)*size(B,dim=1),size(A,dim=2)*size(B,dim=2))
+      end function kronecker_product_iint8
+      pure module function kronecker_product_iint16(A, B) result(C)
+        integer(int16), intent(in) :: A(:,:), B(:,:)
+        integer(int16) :: C(size(A,dim=1)*size(B,dim=1),size(A,dim=2)*size(B,dim=2))
+      end function kronecker_product_iint16
+      pure module function kronecker_product_iint32(A, B) result(C)
+        integer(int32), intent(in) :: A(:,:), B(:,:)
+        integer(int32) :: C(size(A,dim=1)*size(B,dim=1),size(A,dim=2)*size(B,dim=2))
+      end function kronecker_product_iint32
+      pure module function kronecker_product_iint64(A, B) result(C)
+        integer(int64), intent(in) :: A(:,:), B(:,:)
+        integer(int64) :: C(size(A,dim=1)*size(B,dim=1),size(A,dim=2)*size(B,dim=2))
+      end function kronecker_product_iint64
+  end interface kronecker_product
+
 
   ! Cross product (of two vectors)
   interface cross_product

src/stdlib_linalg_kronecker.f90

@@ -0,0 +1,159 @@
+submodule (stdlib_linalg) stdlib_linalg_kronecker
+
+  implicit none
+
+contains
+
+    pure module function kronecker_product_rsp(A, B) result(C)
+      real(sp), intent(in) :: A(:,:), B(:,:)
+      real(sp) :: C(size(A,dim=1)*size(B,dim=1),size(A,dim=2)*size(B,dim=2))
+      integer :: m1, n1, 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 Wikipedia convention for ordering of the matrix elements
+	    ! https://en.wikipedia.org/wiki/Kronecker_product
+            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_rsp
+    pure module function kronecker_product_rdp(A, B) result(C)
+      real(dp), intent(in) :: A(:,:), B(:,:)
+      real(dp) :: C(size(A,dim=1)*size(B,dim=1),size(A,dim=2)*size(B,dim=2))
+      integer :: m1, n1, 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 Wikipedia convention for ordering of the matrix elements
+	    ! https://en.wikipedia.org/wiki/Kronecker_product
+            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_rdp
+    pure module function kronecker_product_csp(A, B) result(C)
+      complex(sp), intent(in) :: A(:,:), B(:,:)
+      complex(sp) :: C(size(A,dim=1)*size(B,dim=1),size(A,dim=2)*size(B,dim=2))
+      integer :: m1, n1, 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 Wikipedia convention for ordering of the matrix elements
+	    ! https://en.wikipedia.org/wiki/Kronecker_product
+            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_csp
+    pure module function kronecker_product_cdp(A, B) result(C)
+      complex(dp), intent(in) :: A(:,:), B(:,:)
+      complex(dp) :: C(size(A,dim=1)*size(B,dim=1),size(A,dim=2)*size(B,dim=2))
+      integer :: m1, n1, 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 Wikipedia convention for ordering of the matrix elements
+	    ! https://en.wikipedia.org/wiki/Kronecker_product
+            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_cdp
+    pure module function kronecker_product_iint8(A, B) result(C)
+      integer(int8), intent(in) :: A(:,:), B(:,:)
+      integer(int8) :: C(size(A,dim=1)*size(B,dim=1),size(A,dim=2)*size(B,dim=2))
+      integer :: m1, n1, 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 Wikipedia convention for ordering of the matrix elements
+	    ! https://en.wikipedia.org/wiki/Kronecker_product
+            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_iint8
+    pure module function kronecker_product_iint16(A, B) result(C)
+      integer(int16), intent(in) :: A(:,:), B(:,:)
+      integer(int16) :: C(size(A,dim=1)*size(B,dim=1),size(A,dim=2)*size(B,dim=2))
+      integer :: m1, n1, 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 Wikipedia convention for ordering of the matrix elements
+	    ! https://en.wikipedia.org/wiki/Kronecker_product
+            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_iint16
+    pure module function kronecker_product_iint32(A, B) result(C)
+      integer(int32), intent(in) :: A(:,:), B(:,:)
+      integer(int32) :: C(size(A,dim=1)*size(B,dim=1),size(A,dim=2)*size(B,dim=2))
+      integer :: m1, n1, 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 Wikipedia convention for ordering of the matrix elements
+	    ! https://en.wikipedia.org/wiki/Kronecker_product
+            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_iint32
+    pure module function kronecker_product_iint64(A, B) result(C)
+      integer(int64), intent(in) :: A(:,:), B(:,:)
+      integer(int64) :: C(size(A,dim=1)*size(B,dim=1),size(A,dim=2)*size(B,dim=2))
+      integer :: m1, n1, 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 Wikipedia convention for ordering of the matrix elements
+	    ! https://en.wikipedia.org/wiki/Kronecker_product
+            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_iint64
+end submodule stdlib_linalg_kronecker