并行化matlab算法答案

【问题标题】：Parallelize matlab algorithm并行化matlab算法
【发布时间】：2014-07-07 11:16:28
【问题描述】：

我正在使用 SART 方法研究断层扫描算法。我正在使用 2D 投影来获得 3D 体积。

我的主要问题是计算时间真的很长..

例如，对于大小为 88*75 像素的投影，我的程序需要 5 天才能完成。

通过使用 matlab 分析器，我们可以看到与氡矩阵（很重）相乘的行花费了很多时间。

这就是为什么我想并行化我的算法，但我并不真正熟悉它。

这是我要并行化的算法部分：

line which need a lot of time: Correction = mtimesx(facteur, Matrice_Radon_Transposee, 'SPEED');


for k = 1: ite % Boucle définie par le nombre d 'itérations choisi par l'
utilisateur
Soustraction = Projection2 - Projection_old;
Erreur(ligne, k) = mean(Soustraction);
Err = Erreur(ligne, k);
for z = 1: L * Nproj
for a = 1: L * L
Somme = Somme + Matrice_Radon(z, a) * Matrice_Radon(z, a);
end
facteur = landa / Somme;
Correction = mtimesx(facteur, Matrice_Radon_Transposee, 'SPEED'); % Formule de correction additive de l
  'image

            Correction = Correction*Soustraction;
            Ik1 = Ik + Correction;
            Ik=Ik1;                                                             % Stockage de l'
image corrigée dans l
  'image_itération_précédente

             Projection_old=Matrice_Radon*Ik;                          % Calcul de la projection à partir de l'
image_itération_précédente
end
pos3 = 1;
for incr2 = 1: L
for incr3 = 1: L
Ik2(incr2, incr3) = Ik1(pos3, 1); % Transformation vecteur en matrice
pos3 = pos3 + 1;
end
end
Coupe_SIRT(ligne, : , : ) = Ik2; % Stockage des coupes générées
end
fprintf('Progression : %2d/%d\n', ligne, H); % Affichage de l
  'avancement du traitement
 end

【问题讨论】：

请将您的代码格式化为更易读的格式：英文中的缩进和 cmets。
这就是人们（ASTRA、TIGRE、RTK、....）构建 GPU 断层重建算法的原因！否则它们会非常慢，没有 CPU 可以足够快地完成这项工作。

标签： matlab parallel-processing tomography-reconstruction

【解决方案1】：

我在这里看到很多嵌套的 for 循环。您应该在 MATLAB 中尽可能使用矢量化代码而不是 for 循环，以加快您的算法。也不能使用嵌套 parfor。

只要您的循环彼此独立（当然您需要安装 Parallel Computing Toolbox），就可以进行并行化。在这种情况下使用它很简单：

  result = zeros(1,n); % preallocate result storing array
  parfor i=1:n % parallel for loop
    %your code here, e.g.:
    result(i) = ...;
  end

可以手动设置一个 MATALB 工作池，但也可以只使用 parfor 并让 MATLAB 自动处理该池（至少在 r2014a 中）。

请注意，并行化并不总能加快速度。

【讨论】：

【解决方案2】：

我真的不知道如何向量化这段代码..

我不需要对所有代码进行矢量化，只需这 2 行：

Correction = mtimesx(facteur, Matrice_Radon_Transposee, 'SPEED'); 

Correction = Correction*Soustraction;

所以也许并行化 for 循环会很有趣？：

for k = 1: ite

这样写对吗：

result_Projection_old = zeros(1,ite);
result_Ik = zeros(1,ite);

parfor  k = 1 : ite                

    Somme = 0;      
    Projection_old = result_Projection_old(k);
    Ik = result_Ik(k);

    Soustraction = Projection2 - Projection_old;


    for z=1:LNproj
        for a=1:L_carre

            Somme = Somme + Matrice_Radon(z,a)*Matrice_Radon(z,a);

        end

        facteur = landa/Somme;
        Correction = mtimesx(facteur,Matrice_Radon_Transposee,'SPEED');    
        Correction = Correction*Soustraction;
        Ik1 = Ik + Correction;
        Ik=Ik1; 
        result_Ik(k) = Ik;                                                          

        Projection_old=Matrice_Radon*Ik;
        result_Projection_old(k) = Projection_old;                     

    end
end

使用此代码，重建不起作用，而不是编写 result_Ik(k) = Ik 我需要编写 Result_Ik(k+1) = Ik 以便在下一次迭代中从迭代中获得 Ik 的值k-1（上一次迭代）

【讨论】：