我有一个代码，它正在对 3D 矩阵中的每个点进行一些处理。数组 input_vec_1D 由一个不寻常的索引 ind_prime 访问，该索引取决于循环变量（对于上下文，索引由我在 this paper 的等式 42e 中使用的算法确定，我的完整代码是here）。通过首先将矩阵转换为一维数组，计算正确的索引，进行处理，然后重新整形回 3D，我设法使其正常工作：

Nx = 8; Ny = 6; Nz = 4; Ntot = Nx*Ny*Nz;                    % Number of points
xvals = rand(1,Nx); yvals = rand(1,Ny); zvals = rand(1,Nz); % Grid vectors

input_vec_3D = rand(Ny,Nx,Nz); % Dummy 3D array

factor1 = 3.6*xvals; % some constant times xvals
factor2 = 1.2*yvals;
factor3 = 8.5*zvals;

input_vec_1D = reshape( permute(input_vec_3D,[3,1,2]),[Ntot 1]); % Reshape to 1D for loop
output_vec = zeros(Ntot,1);
for ind = 1:Ntot
    j1 = floor( floor( (ind-1)/Nz ) /Ny ) + 1;   
    j2 = mod( floor( (ind-1)/Nz ),Ny ) + 1; 
    j3 = mod( (ind-1),Nz ) + 1;
    n1 = mod( 5*(j1-1),Nx);
    n2 = mod( 3*(j2-1),Ny);
    n3 = mod( 2*(j3-1),Nz);
    ind_prime = mod( ( n3 + Nz*(n2 + Ny*n1) ),Ntot ) + 1; % a different index for input_vec
    output_vec(ind) = output_vec(ind) + input_vec_1D(ind_prime) * factor1(j1)*factor2(j2)*factor3(j3);
end
output_vec = permute( reshape( output_vec,[Nz,Ny,Nx] ),[2,3,1] );  % Reshape back to 3D

这个对所有元素的循环是我代码中最慢的部分，所以我想加快速度 - 通过矢量化或其他方式。

我的数组通常是 512x512x1024 的复数双精度数，因此对于我的应用程序至关重要的是，由于 RAM 有限（大约 6 GB），我没有存储任何临时的超大矩阵，这排除了使用 meshgrid() 生成因子（请注意，factor1、factor2、factor3 只是一维向量，因此它们的内存使用量很小）。

我得到了一个非常相似的循环 here 的帮助，在这种情况下使用 Matlab 的隐式扩展解决了这个问题。然而，这更复杂，因为在处理线中使用了不同的索引 ind_prime、ind 和 j。