Fortran嵌套Do循环计算Kronecker乘积结果异常，求问题排查

问题描述 #

我尝试用Fortran的do循环计算Kronecker乘积（张量积），参考Wolfram Mathworld的符号规范，但得到的结果矩阵C与Python中numpy.kron生成的12x6矩阵不一致，以下是我的Fortran代码和Python示例结果。

Fortran代码 #

program test

implicit none
real, dimension (3,3) :: A
real, dimension (4,2) :: B
real, dimension (:,:), allocatable :: C
integer :: i, j, k, l, m, n, p, q, alpha, beta
integer :: Ccols, Crows

A = reshape( (/ 1, 0, 2, 1, 9, 5, 1, 4, -1 /), shape(A) )
B = reshape( (/ 1, -2, 9, 5, 4, -1, 3, 9 /), shape( B ) )

m = size(A,1) ! no. A rows
n = size(A,2) ! no. A cols
p = size(B,1) ! no. B rows
q = size(B,2) ! no. B cols
allocate(C(m*p, n*q)) ! m*p columns, n*q rows
C = 0

do i = 1, n  ! iterate over A cols
    do j = 1, m  ! iterate over A rows
        do k = 1, q  ! iterate over B cols
            do l = 1, p  ! iterate over B rows
                alpha = p*(i-1) + k  ! C row index
                beta = q*(j-1) + l  ! C col index
                C(beta, alpha) = A(i,j) * B(k,l)
            end do
        end do
    end do
end do

print *, C

end program test

Python示例结果 #

>>> import numpy as np
>>> A = np.array([[1,0,2],[1,9,5],[1,4,-1]])
>>> B = np.array([[1,-2],[9,5],[4,-1],[3,9]])
>>> np.kron(A,B)

   [[  1,  -2,   0,   0,   2,  -4],
    [  9,   5,   0,   0,  18,  10],
    [  4,  -1,   0,   0,   8,  -2],
    [  3,   9,   0,   0,   6,  18],
    [  1,  -2,   9, -18,   5, -10],
    [  9,   5,  81,  45,  45,  25],
    [  4,  -1,  36,  -9,  20,  -5],
    [  3,   9,  27,  81,  15,  45],
    [  1,  -2,   4,  -8,  -1,   2],
    [  9,   5,  36,  20,  -9,  -5],
    [  4,  -1,  16,  -4,  -4,   1],
    [  3,   9,  12,  36,  -3,  -9]]

错误分析 #

代码主要存在三个核心问题：

1. 矩阵元素索引颠倒
   • Fortran中矩阵采用(行号, 列号)的索引方式，你在调用A(i,j)时，i是A的列序号、j是A的行序号，正确写法应为A(j,i)；同理，B(k,l)中k是B的列序号、l是B的行序号，正确写法应为B(l,k)。
2. C矩阵的索引计算错误
   • Kronecker乘积中，C矩阵的行索引应由A的行序号和B的行序号计算：(j-1)*p + l；列索引应由A的列序号和B的列序号计算：(i-1)*q + k。你当前的alpha和beta计算逻辑完全颠倒，且公式错误。
3. allocate语句注释错误
   • 注释标注“mp columns, nq rows”，但Fortran数组声明是(行数, 列数)，实际C(m*p, n*q)是正确的维度（12行6列），只是注释描述错误，不影响运行但易造成混淆。

修正后的Fortran代码 #

program test

implicit none
real, dimension (3,3) :: A
real, dimension (4,2) :: B
real, dimension (:,:), allocatable :: C
integer :: i, j, k, l, m, n, p, q, c_row, c_col

A = reshape( (/ 1, 0, 2, 1, 9, 5, 1, 4, -1 /), shape(A) )
B = reshape( (/ 1, -2, 9, 5, 4, -1, 3, 9 /), shape( B ) )

m = size(A,1) ! A的行数
n = size(A,2) ! A的列数
p = size(B,1) ! B的行数
q = size(B,2) ! B的列数
allocate(C(m*p, n*q)) ! C的维度：m*p行，n*q列
C = 0.0

do i = 1, n  ! 遍历A的列
    do j = 1, m  ! 遍历A的行
        do k = 1, q  ! 遍历B的列
            do l = 1, p  ! 遍历B的行
                c_row = (j-1)*p + l  ! C的行索引
                c_col = (i-1)*q + k  ! C的列索引
                C(c_row, c_col) = A(j,i) * B(l,k)
            end do
        end do
    end do
end do

! 按矩阵格式打印结果，方便对比
do j = 1, size(C,1)
    write(*, '(6F6.0)') C(j,:)
end do

end program test

内容的提问来源于stack exchange，提问作者Medulla Oblongata