From c14f6ebc3a5e127049fbd7bce0f4910fbb721639 Mon Sep 17 00:00:00 2001 From: Adam Denchfield Date: Fri, 3 Mar 2023 15:32:46 -0600 Subject: [PATCH 01/11] Added the Kronecker product functionality to stdlib_linalg, and added appropriate unit tests for all supported types. --- src/CMakeLists.txt | 1 + src/stdlib_linalg.fypp | 15 ++ test/linalg/test_linalg.fypp | 342 ++++++++++++++++++++++++++++++++++- 3 files changed, 357 insertions(+), 1 deletion(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 8f512af56..ceb1bd2b9 100644 --- a/src/CMakeLists.txt +++ b/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 diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index cfa43d3d9..c2c29775b 100644 --- a/src/stdlib_linalg.fypp +++ b/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 diff --git a/test/linalg/test_linalg.fypp b/test/linalg/test_linalg.fypp index 2ffd2d7de..3a416c034 100644 --- a/test/linalg/test_linalg.fypp +++ b/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) From 04caad8a80e3bff2221a44d573ce46d7ed33fa55 Mon Sep 17 00:00:00 2001 From: Adam Denchfield Date: Fri, 3 Mar 2023 15:36:49 -0600 Subject: [PATCH 02/11] Update README.md --- src/stdlib_linalg_kronecker.fypp | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) create mode 100644 src/stdlib_linalg_kronecker.fypp diff --git a/src/stdlib_linalg_kronecker.fypp b/src/stdlib_linalg_kronecker.fypp new file mode 100644 index 000000000..49114b294 --- /dev/null +++ b/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 From dc814f679e50460f7a280b7a4dd6243a2bb1a693 Mon Sep 17 00:00:00 2001 From: adenchfi Date: Sat, 4 Mar 2023 15:36:54 -0600 Subject: [PATCH 03/11] Update test/linalg/test_linalg.fypp Co-authored-by: Jeremie Vandenplas --- test/linalg/test_linalg.fypp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/linalg/test_linalg.fypp b/test/linalg/test_linalg.fypp index 3a416c034..fe6601433 100644 --- a/test/linalg/test_linalg.fypp +++ b/test/linalg/test_linalg.fypp @@ -580,8 +580,8 @@ contains end do end do - do j=1, n2 - do i=1, m2 + do j = 1, n2 + do i = 1, m2 B(i,j) = i*j ! B = [[1, 2, 3], [2, 4, 6]] end do end do From 9612c71dbfe85b7c9b7ff7fc310d0cfb341da654 Mon Sep 17 00:00:00 2001 From: adenchfi Date: Sat, 4 Mar 2023 15:41:20 -0600 Subject: [PATCH 04/11] Apply suggestions from code review Adding spaces to = signs Co-authored-by: Jeremie Vandenplas --- src/stdlib_linalg_kronecker.fypp | 4 ++-- test/linalg/test_linalg.fypp | 6 +++--- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/stdlib_linalg_kronecker.fypp b/src/stdlib_linalg_kronecker.fypp index 49114b294..be29e8b58 100644 --- a/src/stdlib_linalg_kronecker.fypp +++ b/src/stdlib_linalg_kronecker.fypp @@ -17,8 +17,8 @@ contains maxM2 = size(B, dim=1) maxN2 = size(B, dim=2) - do n1=1, maxN1 - do m1=1, maxM1 + 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 diff --git a/test/linalg/test_linalg.fypp b/test/linalg/test_linalg.fypp index fe6601433..d48f2cba0 100644 --- a/test/linalg/test_linalg.fypp +++ b/test/linalg/test_linalg.fypp @@ -566,7 +566,7 @@ contains 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 + 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) @@ -574,8 +574,8 @@ contains integer :: i,j - do j=1, n1 - do i=1, m1 + do j = 1, n1 + do i = 1, m1 A(i,j) = i*j ! A = [1, 2] end do end do From 0ce039b096cdd6dd3fa0d72cc4617d29b6b0adb8 Mon Sep 17 00:00:00 2001 From: Adam Denchfield Date: Sat, 4 Mar 2023 16:01:30 -0600 Subject: [PATCH 05/11] Added spec and example program. --- doc/specs/stdlib_linalg.md | 31 ++++++++++++++++++++ example/linalg/example_kronecker_product.f90 | 25 ++++++++++++++++ src/stdlib_linalg_kronecker.fypp | 3 +- 3 files changed, 58 insertions(+), 1 deletion(-) create mode 100644 example/linalg/example_kronecker_product.f90 diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 29eb9878d..1e0fa2ecd 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -160,6 +160,37 @@ Returns a rank-2 array equal to `u v^T` (where `u, v` are considered column vect {!example/linalg/example_outer_product.f90!} ``` +## `kronecker_product` - Computes the kronecker product of two rank-2 arrays + +### Status + +Experimental + +### Description + +Computes the Kronecker product of two rank-2 arrays + +### Syntax + +`C = [[stdlib_linalg(module):kronecker_product(interface)]](A, B)` + +### Arguments + +`A`: Shall be a rank-2 array with dimensions M1, N1 + +`v`: Shall be a rank-2 array with dimensions M2, N2 + +### Return value + +Returns a rank-2 array equal to `A \otimes B`. The shape of the returned array is `[M1*M2, N1*N2]`. + +### Example + +```fortran +{!example/linalg/example_kronecker_product.f90!} +``` + + ## `cross_product` - Computes the cross product of two vectors ### Status diff --git a/example/linalg/example_kronecker_product.f90 b/example/linalg/example_kronecker_product.f90 new file mode 100644 index 000000000..be9d72a0d --- /dev/null +++ b/example/linalg/example_kronecker_product.f90 @@ -0,0 +1,25 @@ +program example_kronecker_product + use stdlib_linalg, only: kronecker_product + implicit none + integer, parameter :: m1=1, n1=2, m2=2, n2=3 + 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_kronecker.fypp b/src/stdlib_linalg_kronecker.fypp index 49114b294..5b47ba90c 100644 --- a/src/stdlib_linalg_kronecker.fypp +++ b/src/stdlib_linalg_kronecker.fypp @@ -19,7 +19,8 @@ contains do n1=1, maxN1 do m1=1, maxM1 - ! We use the numpy.kron convention for ordering of the matrix elements + ! 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 From b3755cc8b203dbba1d0005492e94d84cfdd59997 Mon Sep 17 00:00:00 2001 From: Adam Denchfield Date: Sat, 4 Mar 2023 17:13:02 -0600 Subject: [PATCH 06/11] Condensed the tests into a single fypp loop, and made the expected output array a parameter. --- test/linalg/test_linalg.fypp | 317 ++--------------------------------- 1 file changed, 12 insertions(+), 305 deletions(-) diff --git a/test/linalg/test_linalg.fypp b/test/linalg/test_linalg.fypp index d48f2cba0..f634e69a4 100644 --- a/test/linalg/test_linalg.fypp +++ b/test/linalg/test_linalg.fypp @@ -562,15 +562,20 @@ contains "trace(h) == sum(c(0:nd:2)) failed.") end subroutine test_trace_int64 - - subroutine test_kronecker_product_rsp(error) + + +#:set RCI_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES + INT_KINDS_TYPES + #:for k1, t1 in RCI_KINDS_TYPES + subroutine test_kronecker_product_${t1[0]}$${k1}$(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 + ${t1}$, dimension(m1*m2,n1*n2), parameter :: expected & + = transpose(reshape([1,2,3, 2,4,6, 2,4,6, 4,8,12], [m2*n2, m1*n1])) + ${t1}$, 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) + ${t1}$ :: A(m1,n1), B(m2,n2) + ${t1}$ :: C(m1*m2,n1*n2), diff(m1*m2,n1*n2) integer :: i,j @@ -587,312 +592,14 @@ contains 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 - + end subroutine test_kronecker_product_${t1[0]}$${k1}$ + #:endfor subroutine test_outer_product_rsp(error) !> Error handling From 25b219da0b08bbe78033c67e8624a9d7aa78b6ae Mon Sep 17 00:00:00 2001 From: adenchfi Date: Sat, 4 Mar 2023 17:20:04 -0600 Subject: [PATCH 07/11] Apply suggestions from code review Adding spaces and capitalization. Co-authored-by: Jeremie Vandenplas --- doc/specs/stdlib_linalg.md | 2 +- example/linalg/example_kronecker_product.f90 | 10 +++++----- src/stdlib_linalg.fypp | 2 +- 3 files changed, 7 insertions(+), 7 deletions(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 1e0fa2ecd..088901cf2 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -160,7 +160,7 @@ Returns a rank-2 array equal to `u v^T` (where `u, v` are considered column vect {!example/linalg/example_outer_product.f90!} ``` -## `kronecker_product` - Computes the kronecker product of two rank-2 arrays +## `kronecker_product` - Computes the Kronecker product of two rank-2 arrays ### Status diff --git a/example/linalg/example_kronecker_product.f90 b/example/linalg/example_kronecker_product.f90 index be9d72a0d..00b5e2d81 100644 --- a/example/linalg/example_kronecker_product.f90 +++ b/example/linalg/example_kronecker_product.f90 @@ -1,18 +1,18 @@ program example_kronecker_product use stdlib_linalg, only: kronecker_product implicit none - integer, parameter :: m1=1, n1=2, m2=2, n2=3 + integer, parameter :: m1 = 1, n1 = 2, m2 = 2, n2 = 3 real :: A(m1, n1), B(m2,n2) real, allocatable :: C(:,:) - do j=1, n1 - do i=1, m1 + 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 ] + do j = 1, n2 + do i = 1, m2 ! B = [ 1, 2, 3 ] B(i,j) = i*j ! [ 2, 4, 6 ] end do end do diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index c2c29775b..3ed905c56 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -97,7 +97,7 @@ module stdlib_linalg interface kronecker_product !! version: experimental !! - !! Computes the Kronecker product of two arrays size M1xN1, M2xN2, returning an (M1*M2)x(N1*N2) array + !! 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)) #:for k1, t1 in RCI_KINDS_TYPES From 22536ec235f777fa7be799bb9c4a87de1b89bce8 Mon Sep 17 00:00:00 2001 From: adenchfi Date: Sun, 5 Mar 2023 10:18:09 -0600 Subject: [PATCH 08/11] Apply suggestions from code review Fixing minor typo v -> B Co-authored-by: Jeremie Vandenplas --- doc/specs/stdlib_linalg.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 088901cf2..671cfee2f 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -178,7 +178,7 @@ Computes the Kronecker product of two rank-2 arrays `A`: Shall be a rank-2 array with dimensions M1, N1 -`v`: Shall be a rank-2 array with dimensions M2, N2 +`B`: Shall be a rank-2 array with dimensions M2, N2 ### Return value From 3fd0defb4c569171ca29ba679180a446f2c8b961 Mon Sep 17 00:00:00 2001 From: adenchfi Date: Sun, 5 Mar 2023 10:19:39 -0600 Subject: [PATCH 09/11] Update example/linalg/example_kronecker_product.f90 Forgot to add the loop integers. Co-authored-by: Ian Giestas Pauli --- example/linalg/example_kronecker_product.f90 | 1 + 1 file changed, 1 insertion(+) diff --git a/example/linalg/example_kronecker_product.f90 b/example/linalg/example_kronecker_product.f90 index 00b5e2d81..98bab0eca 100644 --- a/example/linalg/example_kronecker_product.f90 +++ b/example/linalg/example_kronecker_product.f90 @@ -2,6 +2,7 @@ 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(:,:) From 94a555a662fb179d42eebe6249f5cdb6c78a7bc5 Mon Sep 17 00:00:00 2001 From: Adam Denchfield Date: Sun, 5 Mar 2023 10:31:54 -0600 Subject: [PATCH 10/11] Organized the unit test list into a fypp loop as well. --- test/linalg/test_linalg.fypp | 16 +++++----------- 1 file changed, 5 insertions(+), 11 deletions(-) diff --git a/test/linalg/test_linalg.fypp b/test/linalg/test_linalg.fypp index f634e69a4..6fdf7f17d 100644 --- a/test/linalg/test_linalg.fypp +++ b/test/linalg/test_linalg.fypp @@ -1,4 +1,5 @@ #:include "common.fypp" +#:set RCI_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES + INT_KINDS_TYPES module test_linalg use testdrive, only : new_unittest, unittest_type, error_type, check, skip_test @@ -48,16 +49,9 @@ 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), & + #:for k1, t1 in RCI_KINDS_TYPES + new_unittest("kronecker_product_${t1[0]}$${k1}$", test_kronecker_product_${t1[0]}$${k1}$), & + #:endfor 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), & @@ -564,7 +558,7 @@ contains end subroutine test_trace_int64 -#:set RCI_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES + INT_KINDS_TYPES + #:for k1, t1 in RCI_KINDS_TYPES subroutine test_kronecker_product_${t1[0]}$${k1}$(error) !> Error handling From a6584f75ec2af3dd9c1b580bb8d194dbe3640ec5 Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Wed, 8 Mar 2023 08:42:03 -0500 Subject: [PATCH 11/11] Update src/stdlib_linalg_kronecker.fypp Co-authored-by: Ian Giestas Pauli --- src/stdlib_linalg_kronecker.fypp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/stdlib_linalg_kronecker.fypp b/src/stdlib_linalg_kronecker.fypp index 3ab1f264f..38895f73e 100644 --- a/src/stdlib_linalg_kronecker.fypp +++ b/src/stdlib_linalg_kronecker.fypp @@ -10,7 +10,7 @@ contains 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 + integer :: m1, n1, maxM1, maxN1, maxM2, maxN2 maxM1 = size(A, dim=1) maxN1 = size(A, dim=2)