From 261a6274477ec8e23f94bd8be384158149ed6628 Mon Sep 17 00:00:00 2001 From: jvdp1 Date: Wed, 8 Mar 2023 13:49:13 +0000 Subject: [PATCH] =?UTF-8?q?Deploying=20to=20stdlib-fpm=20from=20@=20fortra?= =?UTF-8?q?n-lang/stdlib@2fdfab47e263249ab80952c79e7a7638e481ddba=20?= =?UTF-8?q?=F0=9F=9A=80?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- example/example_kronecker_product.f90 | 26 +++++ src/stdlib_linalg.f90 | 41 +++++++ src/stdlib_linalg_kronecker.f90 | 159 ++++++++++++++++++++++++++ 3 files changed, 226 insertions(+) create mode 100644 example/example_kronecker_product.f90 create mode 100644 src/stdlib_linalg_kronecker.f90 diff --git a/example/example_kronecker_product.f90 b/example/example_kronecker_product.f90 new file mode 100644 index 000000000..98bab0eca --- /dev/null +++ b/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 diff --git a/src/stdlib_linalg.f90 b/src/stdlib_linalg.f90 index cff487272..54ea03394 100644 --- a/src/stdlib_linalg.f90 +++ b/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 diff --git a/src/stdlib_linalg_kronecker.f90 b/src/stdlib_linalg_kronecker.f90 new file mode 100644 index 000000000..77ad16a9b --- /dev/null +++ b/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