I'm trying to calculate a tensor/Kronecker product in Fortran for two vectors x=[1,4] and complex y=[2+i,3+2i,3-i]. The solution z should be an array [2+i, 3+2i, 3-i, 8+4i, 12+8i, 12-4i]. In my MWE,
program myfun
implicit none
real, dimension (:), allocatable :: x
complex, dimension (:), allocatable :: y
complex, dimension (:), allocatable :: z
integer :: i, j ! indexing
integer :: k=1 ! indexing
x = (/ 1,4 /)
y = (/ 2,3,3 /) + (/ 1,2,-1 /) * cmplx(0,1)
allocate(z(size(x)*size(y)))
do i = 1,size(x)
do j = 1,size(y)
z(k) = x(i) * y(j)
print *, x(i) * y(j)
k = k+1
end do
end do
print *, z
print *, KronProd(x,y)
contains
function KronProd(A,B) result(C)
IMPLICIT NONE
real, dimension (:,:), intent(in) :: A, B
real, dimension (:,:), allocatable :: C
integer :: i = 0, j = 0, k = 0, l = 0
integer :: m = 0, n = 0, p = 0, q = 0
allocate(C(size(A,1)*size(B,1),size(A,2)*size(B,2)))
C = 0
do i = 1,size(A,1)
do j = 1,size(A,2)
n=(i-1)*size(B,1) + 1
m=n+size(B,1) - 1
p=(j-1)*size(B,2) + 1
q=p+size(B,2) - 1
C(n:m,p:q) = A(i,j)*B
enddo
enddo
end function KronProd
end program myfun
Edit: After editing the code using the comments below, the z vector gives the correct result.
I don't want a do loop every time I calculate the Kronecker product, so I'm trying the KronProd function from ECoulter's Kronecker.f95.
In my code, KronProd is taken unchanged from Kronecker.f95. But I get Error: Rank mismatch in argument 'a' (rank-2 and rank-1). As pointed out by other users, there could be an issue with the dimension (:,:), intent(in). How do I make this statement more general so that the A and B variables can take either matrices or arrays?
The function KronProd cannot be directly used because
It is intended for 2D (rank 2)arrays, not 1D (rank 1) arrays.
It is intended for real arrays, not complex ones.
There are multiple approaches how to enable the function.
First would be pointer remapping of the 1D arrays to a 2D arrays. The result would have to be remapped back. Because only one of.the arrays is complex, you could do the real part separately and the imaginary part separately. Or you can create a complex version of the subroutine.
Another approach would be creating multiple specific versions of the function (different array ranks and real/complex versions) and a generic interface.
For the remapping consider
real, dimension (:), allocatable, target :: x
complex, dimension (:), allocatable, target :: y
complex, dimension (:), allocatable, target :: z
real, pointer, dimension(:,:) :: xp
complex, pointer, dimension(:,:) :: yp, zp
...
xp(1:size(x),1:1) => x
yp(1:size(y),1:1) => y
zp(1:size(z),1:1) => z
zp = KronProd_real_complex(xp,yp)
where the function is now
function KronProd_real_complex(A,B) result(C)
IMPLICIT NONE
real, dimension (:,:), intent(in) :: A
complex, dimension (:,:), intent(in) :: B
complex, dimension (:,:), allocatable :: C
integer :: i = 0, j = 0, k = 0, l = 0
integer :: m = 0, n = 0, p = 0, q = 0
allocate(C(size(A,1)*size(B,1),size(A,2)*size(B,2)))
C = 0
do i = 1,size(A,1)
do j = 1,size(A,2)
n=(i-1)*size(B,1) + 1
m=n+size(B,1) - 1
p=(j-1)*size(B,2) + 1
q=p+size(B,2) - 1
C(n:m,p:q) = A(i,j)*B
enddo
enddo
end function KronProd_real_complex
Another option, and quite a good one for small arrays where performance is not an issue, is to just copy the 1D arrays into 2D arrays, separate the real and imaginary components in this way and call the original function.
real, allocatable, dimension(:,:) :: x2d, y2d, z2d
allocate(x2d(1:size(x),1))
allocate(y2d(1:size(y),1))
allocate(z2d(1:size(z),1))
x2d(:,1) = x
y2d(:,1) = y%re
z2d = KronProd(x2d ,y2d)
z%re = z2d(:,1)
y2d(:,1) = y%im
z2d = KronProd(x2d ,y2d)
z%im = z2d(:,1)
print *, z
This can also be hidden inside a function, potentially with automatic arrays instead of the allocatable ones.