使用numpy memmap时出现Bus Error: Core dumped问题求助

首先得戳破一个误区：虽然你用了np.memmap，但中间的FFT和fftshift操作并没有真正做到全程核外计算——很多步骤会偷偷把整个超大数组加载到物理内存里，这才是触发Bus Error的核心原因，哪怕你有128GB内存也顶不住。

问题成因拆解 #

1. fftshift直接触发全数组加载 #

np.fft.fftshift(X)接收的是memmap对象，但这个操作需要把数组的上下/左右块交换位置，没法直接通过内存映射的偏移量实现（原文件是连续存储的，fftshift后的数组是不连续的）。所以numpy会把整个65536×65536的float64数组（32GB）完全加载到物理内存，这已经占了你内存的四分之一。

2. FFT中间临时数组的叠加开销 #

接下来的rfft2、A*B、irfft2每一步都会生成超大的临时数组：

• rfft2处理后会得到65536×32769的复数数组（约16GB）；
• A*B会生成同样大小的复数数组（再占16GB）；
• irfft2又会把结果还原回65536×65536的float64数组（32GB）；
  再加上FFT算法本身需要的工作空间（多线程情况下会翻倍），这些加起来已经接近甚至超过128GB的物理内存上限。当内存耗尽时，系统会触发Bus Error（而非普通OOM），因为numpy尝试访问无法映射的内存区域。

3. 合并写法的雪上加霜 #

你把所有步骤合并成一行的写法，反而会让numpy无法及时释放中间临时数组的内存——多个超大数组同时驻留内存，直接压垮内存资源。

解决办法 #

1. 替代fftshift：频域相位调整，避免全数组加载 #

fftshift的本质是给原数组每个元素乘以(-1)^(i+j)（i,j是行列索引），这个操作可以移到频域完成，不用提前对原数组做移位。这样就能完全避免加载整个X数组到内存：

import numpy as np

rows, cols = 65536, 65536
X = np.memmap('X.bin', dtype='float64', mode='r', shape=(rows, cols))
Y = np.memmap('Y.bin', dtype='float64', mode='r', shape=(rows, cols))

# 生成(-1)^(i+j)的相位因子，用memmap存储避免加载到内存
phase_mmap = np.memmap('phase.bin', dtype='float64', mode='w+', shape=(rows, cols))
i = np.arange(rows)[:, np.newaxis]
j = np.arange(cols)
phase_mmap[:] = (-1) ** (i + j)

# 直接对memmap做FFT，再通过频域乘法替代fftshift
X_fft = np.fft.rfft2(X)
phase_fft = np.fft.rfft2(phase_mmap)
A = X_fft * phase_fft  # 等价于原代码的np.fft.rfft2(np.fft.fftshift(X))

B = np.fft.rfft2(Y)
C = np.fft.irfft2(A * B)
alpha_x, alpha_y = np.gradient(C, edge_order=2)

# 清理临时文件
import os
os.remove('phase.bin')

2. 分块处理：把大数组拆成小块计算 #

如果频域相位调整的方式不适用，你可以把65536×65536的数组拆成多个8192×8192的小块（你已经验证过这个大小可行），逐块处理后再合并结果：

import numpy as np

rows, cols = 65536, 65536
block_size = 8192

X = np.memmap('X.bin', dtype='float64', mode='r', shape=(rows, cols))
Y = np.memmap('Y.bin', dtype='float64', mode='r', shape=(rows, cols))

# 预先创建结果的memmap文件
alpha_x = np.memmap('alpha_x.bin', dtype='float64', mode='w+', shape=(rows, cols))
alpha_y = np.memmap('alpha_y.bin', dtype='float64', mode='w+', shape=(rows, cols))

for i in range(0, rows, block_size):
    for j in range(0, cols, block_size):
        # 读取小块数据
        X_block = X[i:i+block_size, j:j+block_size]
        Y_block = Y[i:i+block_size, j:j+block_size]
        # 执行单块计算
        A_block = np.fft.rfft2(np.fft.fftshift(X_block))
        B_block = np.fft.rfft2(Y_block)
        C_block = np.fft.irfft2(A_block * B_block)
        dx_block, dy_block = np.gradient(C_block, edge_order=2)
        # 写回结果到memmap
        alpha_x[i:i+block_size, j:j+block_size] = dx_block
        alpha_y[i:i+block_size, j:j+block_size] = dy_block

这种方式每块的内存占用和你之前成功的8192×8192场景一致，完全不会触发内存溢出。

3. 优化内存与FFT配置 #

• 关闭numpy多线程：多线程FFT会额外占用大量工作内存，你可以通过环境变量关闭：

  export OMP_NUM_THREADS=1
  

  或者在代码里设置：

  import os
  os.environ['OMP_NUM_THREADS'] = '1'
  

• 检查memmap文件完整性：确保X.bin和Y.bin的大小正好是65536*65536*8字节（32GB），文件损坏或大小不足也会触发Bus Error。

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