问题

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?

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

答案1

得分: 0

以下是您要翻译的内容：

"The function KronProd cannot be directly used because

1. It is intended for 2D (rank 2) arrays, not 1D (rank 1) arrays.

2. 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 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."

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)

  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

  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