Fortran中能否实现三维矩阵卷积？具体实现方法是什么

Fortran完全可以实现三维矩阵卷积，不管是直接计算还是FFT频域加速的方案都有成熟落地路径，作为原生面向高性能科学计算的语言，Fortran做这类矩阵运算的性能通常优于多数通用编程语言，非常适合大尺寸三维数据的卷积计算场景。

方案1：直接计算法（适合卷积核尺寸远小于输入矩阵的场景） #

这个方案逻辑最直观，不需要依赖任何第三方库，新手最容易快速跑通。
三维卷积的核心逻辑是滑动窗口加权求和：对于尺寸为(Nx,Ny,Nz)的输入矩阵A，尺寸为(Kx,Ky,Kz)的卷积核B，full模式下输出矩阵C的尺寸为(Nx+Kx-1, Ny+Ky-1, Nz+Kz-1)，每个位置的取值为对应窗口内元素与卷积核的乘积累加和。
以下是可直接编译运行的基础实现代码：

program direct_3d_conv
    implicit none
    ! 可根据自身需求修改矩阵尺寸参数
    integer, parameter :: Nx=32, Ny=32, Nz=32, Kx=3, Ky=3, Kz=3
    real :: A(Nx,Ny,Nz), B(Kx,Ky,Kz), C(Nx+Kx-1, Ny+Ky-1, Nz+Kz-1)
    integer :: i,j,k, di,dj,dk
    real :: val

    ! 此处替换为自定义的输入矩阵、卷积核初始化逻辑
    call random_number(A)
    call random_number(B)
    C = 0.0

    ! 遍历输出矩阵的所有位置
    do i = 1, Nx+Kx-1
        do j = 1, Ny+Ky-1
            do k = 1, Nz+Kz-1
                val = 0.0
                ! 遍历卷积核完成乘加计算
                do di = 1, Kx
                    do dj = 1, Ky
                        do dk = 1, Kz
                            ! 边界判定：超出输入矩阵范围的位置按0填充
                            if (i-di+1 >=1 .and. i-di+1 <=Nx .and. &
                                j-dj+1 >=1 .and. j-dj+1 <=Ny .and. &
                                k-dk+1 >=1 .and. k-dk+1 <=Nz) then
                                val = val + A(i-di+1, j-dj+1, k-dk+1) * B(di,dj,dk)
                            end if
                        end do
                    end do
                end do
                C(i,j,k) = val
            end do
        end do
    end do

    print *, "直接法三维卷积计算完成，输出矩阵尺寸：", shape(C)
end program direct_3d_conv

使用注意点：

• 上述代码实现的是深度学习、图像处理场景常用的互相关逻辑（卷积核不翻转），如果需要数学定义上的严格卷积，只需要把卷积核的遍历顺序倒序即可。
• 如果需要same（输出与输入同尺寸）、valid（仅计算窗口完全覆盖输入的位置）模式，只需要调整输出矩阵尺寸、循环索引范围、边界截取逻辑即可。
• 编译时开启高阶优化选项可以大幅提升性能：gfortran添加-O3 -ffast-math参数，Intel Fortran添加-O3参数，编译器会自动完成循环展开、向量化等优化，3x3x3、5x5x5这类小核场景下性能足够支撑生产使用。

方案2：FFT加速法（适合大尺寸卷积核、大输入矩阵场景） #

根据卷积定理，两个矩阵空域的卷积等价于频域的逐点相乘，再通过逆傅里叶变换转回空域。当输入矩阵和卷积核尺寸都较大（比如单边长度≥64）时，FFT方法的计算效率远高于直接法。
Fortran生态有成熟的高性能FFT库FFTW可供调用，不需要手写FFT逻辑，核心实现步骤如下：

1. 计算统一的FFT计算尺寸：每个维度的长度≥输入对应维度长度+卷积核对应维度长度-1，通常取2的整数次幂可以获得最优FFT性能。
2. 将输入矩阵、卷积核分别补零到上述统一尺寸。
3. 对两个补零后的矩阵分别做三维正FFT得到频域矩阵。
4. 两个频域矩阵做逐点复数相乘。
5. 对相乘结果做三维逆FFT，归一化后取实部就是卷积结果。
   以下是基于FFTW库的参考实现代码：

program fft_3d_conv
    use, intrinsic :: iso_c_binding
    implicit none
    include 'fftw3.f03'
    ! 输入矩阵、卷积核尺寸
    integer, parameter :: Nx=128, Ny=128, Nz=128, Kx=64, Ky=64, Kz=64
    ! FFT计算统一尺寸
    integer, parameter :: Lx=192, Ly=192, Lz=192
    real :: A(Nx,Ny,Nz), B(Kx,Ky,Kz)
    complex(C_DOUBLE_COMPLEX) :: FA(Lx/2+1, Ly, Lz), FB(Lx/2+1, Ly, Lz), FC(Lx/2+1, Ly, Lz)
    real(C_DOUBLE) :: C(Lx,Ly,Lz)
    type(C_PTR) :: plan_fwdA, plan_fwdB, plan_inv
    integer :: i,j,k

    ! 此处替换为自定义的输入矩阵、卷积核初始化逻辑
    call random_number(A)
    call random_number(B)

    ! 提前创建FFT计划，自动选择最优计算路径
    plan_fwdA = fftw_plan_dft_r2c_3d(Lx, Ly, Lz, C, FA, FFTW_ESTIMATE)
    plan_fwdB = fftw_plan_dft_r2c_3d(Lx, Ly, Lz, C, FB, FFTW_ESTIMATE)
    plan_inv = fftw_plan_dft_c2r_3d(Lx, Ly, Lz, FC, C, FFTW_ESTIMATE)

    ! 补零填充输入矩阵A并做正变换
    C = 0.0
    C(1:Nx, 1:Ny, 1:Nz) = A
    call fftw_execute_dft_r2c(plan_fwdA, C, FA)

    ! 补零填充卷积核B并做正变换
    C = 0.0
    C(1:Kx, 1:Ky, 1:Kz) = B
    call fftw_execute_dft_r2c(plan_fwdB, C, FB)

    ! 频域逐点相乘
    FC = FA * FB

    ! 逆变换回空域并做归一化
    call fftw_execute_dft_c2r(plan_inv, FC, C)
    C = C / real(Lx*Ly*Lz)

    ! 销毁FFT计划释放资源
    call fftw_destroy_plan(plan_fwdA)
    call fftw_destroy_plan(plan_fwdB)
    call fftw_destroy_plan(plan_inv)

    print *, "FFT法三维卷积计算完成"
end program fft_3d_conv

使用注意点：

• 编译时需要链接FFTW3库，以gfortran为例，编译命令为gfortran fft_conv.f90 -o fft_conv -lfftw3 -lm -O3。
• 计算完成后根据自身需求截取对应位置的子矩阵，即可得到full/same/valid模式的输出结果。
• FFT方法存在微小的浮点精度损失，对精度要求极高的小核场景优先用直接法。

• 刚接触Fortran时优先跑通直接法的示例代码，确认索引逻辑、边界处理、卷积核翻转规则符合自身业务需求后，再根据性能要求选择是否切换到FFT方案。
• 如果不想手写底层逻辑，也可以直接调用成熟Fortran数值库中封装好的三维卷积接口，比如BLAS/LAPACK扩展库、OpenCV Fortran绑定中都有现成实现，传入矩阵参数即可直接调用。

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