Skip to content

Commit 261a627

Browse files
jvdp1jvdp1
committed
Deploying to stdlib-fpm from @ 2fdfab4 🚀
1 parent c542299 commit 261a627

File tree

3 files changed

+226
-0
lines changed

3 files changed

+226
-0
lines changed

‎example/example_kronecker_product.f90

Lines changed: 26 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,26 @@
1+
program example_kronecker_product
2+
use stdlib_linalg, only: kronecker_product
3+
implicit none
4+
integer, parameter :: m1 = 1, n1 = 2, m2 = 2, n2 = 3
5+
integer :: i, j
6+
real :: A(m1, n1), B(m2,n2)
7+
real, allocatable :: C(:,:)
8+
9+
do j = 1, n1
10+
do i = 1, m1
11+
A(i,j) = i*j ! A = [1, 2]
12+
end do
13+
end do
14+
15+
do j = 1, n2
16+
do i = 1, m2 ! B = [ 1, 2, 3 ]
17+
B(i,j) = i*j ! [ 2, 4, 6 ]
18+
end do
19+
end do
20+
21+
C = kronecker_product(A, B)
22+
! C = [ a(1,1) * B(:,:) | a(1,2) * B(:,:) ]
23+
! or in other words,
24+
! C = [ 1.00 2.00 3.00 2.00 4.00 6.00 ]
25+
! [ 2.00 4.00 6.00 4.00 8.00 12.00 ]
26+
end program example_kronecker_product

‎src/stdlib_linalg.f90

Lines changed: 41 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -12,6 +12,7 @@ module stdlib_linalg
1212
public :: eye
1313
public :: trace
1414
public :: outer_product
15+
public :: kronecker_product
1516
public :: cross_product
1617
public :: is_square
1718
public :: is_diagonal
@@ -240,6 +241,46 @@ pure module function outer_product_iint64(u, v) result(res)
240241
end function outer_product_iint64
241242
end interface outer_product
242243

244+
interface kronecker_product
245+
!! version: experimental
246+
!!
247+
!! Computes the Kronecker product of two arrays of size M1xN1, and of M2xN2, returning an (M1*M2)x(N1*N2) array
248+
!! ([Specification](../page/specs/stdlib_linalg.html#
249+
!! kronecker_product-computes-the-kronecker-product-of-two-matrices))
250+
pure module function kronecker_product_rsp(A, B) result(C)
251+
real(sp), intent(in) :: A(:,:), B(:,:)
252+
real(sp) :: C(size(A,dim=1)*size(B,dim=1),size(A,dim=2)*size(B,dim=2))
253+
end function kronecker_product_rsp
254+
pure module function kronecker_product_rdp(A, B) result(C)
255+
real(dp), intent(in) :: A(:,:), B(:,:)
256+
real(dp) :: C(size(A,dim=1)*size(B,dim=1),size(A,dim=2)*size(B,dim=2))
257+
end function kronecker_product_rdp
258+
pure module function kronecker_product_csp(A, B) result(C)
259+
complex(sp), intent(in) :: A(:,:), B(:,:)
260+
complex(sp) :: C(size(A,dim=1)*size(B,dim=1),size(A,dim=2)*size(B,dim=2))
261+
end function kronecker_product_csp
262+
pure module function kronecker_product_cdp(A, B) result(C)
263+
complex(dp), intent(in) :: A(:,:), B(:,:)
264+
complex(dp) :: C(size(A,dim=1)*size(B,dim=1),size(A,dim=2)*size(B,dim=2))
265+
end function kronecker_product_cdp
266+
pure module function kronecker_product_iint8(A, B) result(C)
267+
integer(int8), intent(in) :: A(:,:), B(:,:)
268+
integer(int8) :: C(size(A,dim=1)*size(B,dim=1),size(A,dim=2)*size(B,dim=2))
269+
end function kronecker_product_iint8
270+
pure module function kronecker_product_iint16(A, B) result(C)
271+
integer(int16), intent(in) :: A(:,:), B(:,:)
272+
integer(int16) :: C(size(A,dim=1)*size(B,dim=1),size(A,dim=2)*size(B,dim=2))
273+
end function kronecker_product_iint16
274+
pure module function kronecker_product_iint32(A, B) result(C)
275+
integer(int32), intent(in) :: A(:,:), B(:,:)
276+
integer(int32) :: C(size(A,dim=1)*size(B,dim=1),size(A,dim=2)*size(B,dim=2))
277+
end function kronecker_product_iint32
278+
pure module function kronecker_product_iint64(A, B) result(C)
279+
integer(int64), intent(in) :: A(:,:), B(:,:)
280+
integer(int64) :: C(size(A,dim=1)*size(B,dim=1),size(A,dim=2)*size(B,dim=2))
281+
end function kronecker_product_iint64
282+
end interface kronecker_product
283+
243284

244285
! Cross product (of two vectors)
245286
interface cross_product

‎src/stdlib_linalg_kronecker.f90

Lines changed: 159 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,159 @@
1+
submodule (stdlib_linalg) stdlib_linalg_kronecker
2+
3+
implicit none
4+
5+
contains
6+
7+
pure module function kronecker_product_rsp(A, B) result(C)
8+
real(sp), intent(in) :: A(:,:), B(:,:)
9+
real(sp) :: C(size(A,dim=1)*size(B,dim=1),size(A,dim=2)*size(B,dim=2))
10+
integer :: m1, n1, maxM1, maxN1, maxM2, maxN2
11+
12+
maxM1 = size(A, dim=1)
13+
maxN1 = size(A, dim=2)
14+
maxM2 = size(B, dim=1)
15+
maxN2 = size(B, dim=2)
16+
17+
18+
do n1 = 1, maxN1
19+
do m1 = 1, maxM1
20+
! We use the Wikipedia convention for ordering of the matrix elements
21+
! https://en.wikipedia.org/wiki/Kronecker_product
22+
C((m1-1)*maxM2+1:m1*maxM2, (n1-1)*maxN2+1:n1*maxN2) = A(m1, n1) * B(:,:)
23+
end do
24+
end do
25+
end function kronecker_product_rsp
26+
pure module function kronecker_product_rdp(A, B) result(C)
27+
real(dp), intent(in) :: A(:,:), B(:,:)
28+
real(dp) :: C(size(A,dim=1)*size(B,dim=1),size(A,dim=2)*size(B,dim=2))
29+
integer :: m1, n1, maxM1, maxN1, maxM2, maxN2
30+
31+
maxM1 = size(A, dim=1)
32+
maxN1 = size(A, dim=2)
33+
maxM2 = size(B, dim=1)
34+
maxN2 = size(B, dim=2)
35+
36+
37+
do n1 = 1, maxN1
38+
do m1 = 1, maxM1
39+
! We use the Wikipedia convention for ordering of the matrix elements
40+
! https://en.wikipedia.org/wiki/Kronecker_product
41+
C((m1-1)*maxM2+1:m1*maxM2, (n1-1)*maxN2+1:n1*maxN2) = A(m1, n1) * B(:,:)
42+
end do
43+
end do
44+
end function kronecker_product_rdp
45+
pure module function kronecker_product_csp(A, B) result(C)
46+
complex(sp), intent(in) :: A(:,:), B(:,:)
47+
complex(sp) :: C(size(A,dim=1)*size(B,dim=1),size(A,dim=2)*size(B,dim=2))
48+
integer :: m1, n1, maxM1, maxN1, maxM2, maxN2
49+
50+
maxM1 = size(A, dim=1)
51+
maxN1 = size(A, dim=2)
52+
maxM2 = size(B, dim=1)
53+
maxN2 = size(B, dim=2)
54+
55+
56+
do n1 = 1, maxN1
57+
do m1 = 1, maxM1
58+
! We use the Wikipedia convention for ordering of the matrix elements
59+
! https://en.wikipedia.org/wiki/Kronecker_product
60+
C((m1-1)*maxM2+1:m1*maxM2, (n1-1)*maxN2+1:n1*maxN2) = A(m1, n1) * B(:,:)
61+
end do
62+
end do
63+
end function kronecker_product_csp
64+
pure module function kronecker_product_cdp(A, B) result(C)
65+
complex(dp), intent(in) :: A(:,:), B(:,:)
66+
complex(dp) :: C(size(A,dim=1)*size(B,dim=1),size(A,dim=2)*size(B,dim=2))
67+
integer :: m1, n1, maxM1, maxN1, maxM2, maxN2
68+
69+
maxM1 = size(A, dim=1)
70+
maxN1 = size(A, dim=2)
71+
maxM2 = size(B, dim=1)
72+
maxN2 = size(B, dim=2)
73+
74+
75+
do n1 = 1, maxN1
76+
do m1 = 1, maxM1
77+
! We use the Wikipedia convention for ordering of the matrix elements
78+
! https://en.wikipedia.org/wiki/Kronecker_product
79+
C((m1-1)*maxM2+1:m1*maxM2, (n1-1)*maxN2+1:n1*maxN2) = A(m1, n1) * B(:,:)
80+
end do
81+
end do
82+
end function kronecker_product_cdp
83+
pure module function kronecker_product_iint8(A, B) result(C)
84+
integer(int8), intent(in) :: A(:,:), B(:,:)
85+
integer(int8) :: C(size(A,dim=1)*size(B,dim=1),size(A,dim=2)*size(B,dim=2))
86+
integer :: m1, n1, maxM1, maxN1, maxM2, maxN2
87+
88+
maxM1 = size(A, dim=1)
89+
maxN1 = size(A, dim=2)
90+
maxM2 = size(B, dim=1)
91+
maxN2 = size(B, dim=2)
92+
93+
94+
do n1 = 1, maxN1
95+
do m1 = 1, maxM1
96+
! We use the Wikipedia convention for ordering of the matrix elements
97+
! https://en.wikipedia.org/wiki/Kronecker_product
98+
C((m1-1)*maxM2+1:m1*maxM2, (n1-1)*maxN2+1:n1*maxN2) = A(m1, n1) * B(:,:)
99+
end do
100+
end do
101+
end function kronecker_product_iint8
102+
pure module function kronecker_product_iint16(A, B) result(C)
103+
integer(int16), intent(in) :: A(:,:), B(:,:)
104+
integer(int16) :: C(size(A,dim=1)*size(B,dim=1),size(A,dim=2)*size(B,dim=2))
105+
integer :: m1, n1, maxM1, maxN1, maxM2, maxN2
106+
107+
maxM1 = size(A, dim=1)
108+
maxN1 = size(A, dim=2)
109+
maxM2 = size(B, dim=1)
110+
maxN2 = size(B, dim=2)
111+
112+
113+
do n1 = 1, maxN1
114+
do m1 = 1, maxM1
115+
! We use the Wikipedia convention for ordering of the matrix elements
116+
! https://en.wikipedia.org/wiki/Kronecker_product
117+
C((m1-1)*maxM2+1:m1*maxM2, (n1-1)*maxN2+1:n1*maxN2) = A(m1, n1) * B(:,:)
118+
end do
119+
end do
120+
end function kronecker_product_iint16
121+
pure module function kronecker_product_iint32(A, B) result(C)
122+
integer(int32), intent(in) :: A(:,:), B(:,:)
123+
integer(int32) :: C(size(A,dim=1)*size(B,dim=1),size(A,dim=2)*size(B,dim=2))
124+
integer :: m1, n1, maxM1, maxN1, maxM2, maxN2
125+
126+
maxM1 = size(A, dim=1)
127+
maxN1 = size(A, dim=2)
128+
maxM2 = size(B, dim=1)
129+
maxN2 = size(B, dim=2)
130+
131+
132+
do n1 = 1, maxN1
133+
do m1 = 1, maxM1
134+
! We use the Wikipedia convention for ordering of the matrix elements
135+
! https://en.wikipedia.org/wiki/Kronecker_product
136+
C((m1-1)*maxM2+1:m1*maxM2, (n1-1)*maxN2+1:n1*maxN2) = A(m1, n1) * B(:,:)
137+
end do
138+
end do
139+
end function kronecker_product_iint32
140+
pure module function kronecker_product_iint64(A, B) result(C)
141+
integer(int64), intent(in) :: A(:,:), B(:,:)
142+
integer(int64) :: C(size(A,dim=1)*size(B,dim=1),size(A,dim=2)*size(B,dim=2))
143+
integer :: m1, n1, maxM1, maxN1, maxM2, maxN2
144+
145+
maxM1 = size(A, dim=1)
146+
maxN1 = size(A, dim=2)
147+
maxM2 = size(B, dim=1)
148+
maxN2 = size(B, dim=2)
149+
150+
151+
do n1 = 1, maxN1
152+
do m1 = 1, maxM1
153+
! We use the Wikipedia convention for ordering of the matrix elements
154+
! https://en.wikipedia.org/wiki/Kronecker_product
155+
C((m1-1)*maxM2+1:m1*maxM2, (n1-1)*maxN2+1:n1*maxN2) = A(m1, n1) * B(:,:)
156+
end do
157+
end do
158+
end function kronecker_product_iint64
159+
end submodule stdlib_linalg_kronecker

0 commit comments

Comments
 (0)