如何在 for 循环 c++ 犰狳中加快对巨大矢量的调用

Question

>我正在使用犰狳库。我的程序太慢并且我需要加快它的部分如下

for(int q(0); q < Nk*Nk; q++){
    for(int k(0); k < Nk*Nk; k++){
        int kq = (k+q) % (Nk*Nk);
        cx_mat Gx = ((Eigveck.slice(k)).t())*(Vxk.slice(k)-Vxk.slice(kq))*Eigveck.slice(kq);
        cx_mat Gy = ((Eigveck.slice(k)).t())*(Vyk.slice(k)-Vyk.slice(kq))*Eigveck.slice(kq);
        vec ek = Eigvalk.col(k);
        vec ekq = Eigvalk.col(kq);
        for(int i(0); i < Ltot; i++){
            for(int j(0); j < Ltot; j++){
                chi = chi + (abs(Gx(i,j))*abs(Gx(i,j))+abs(Gy(i,j))*abs(Gy(i,j)))*(1.0/(1.0+exp(ekq(j)/T))-1.0/(1.0+exp(ek(i)/T)))*((ekq(j)-ek(i))/((ekq(j)-ek(i))*(ekq(j)-ek(i))+eta*eta))/(Nk*Nk);
            }
        }
    }
    double qx = (G1(0)*floor(q/Nk)/Nk+G2(0)*(q % Nk)/Nk);
    double qy = (G1(1)*floor(q/Nk)/Nk+G2(1)*(q % Nk)/Nk);
    lindhard << qx << "     " << qy << "     " << -chi << "    " << endl;
}

在这部分之前，我定义了一个巨大的矩阵，Eigvalk和巨大的立方体，Eigveck，Vxk，Vyk。

现在，在 for 循环中对它们的值的调用非常慢，需要很长时间。 立方体包含特征向量和给定问题的其他数量。问题是，对于 Nk=10(测试代码的非常小的 Nk(，计算 Nk*Nk=100 乘以 47 个特征向量需要 0.1 秒，执行使用它们的循环需要 4.5 秒。我已经检查了花费时间的部分是呼叫

cx_mat Gx = .....

我还尝试定义向量或巨大的cx_mat(通过矢量化矩阵(而不是cx_cube，但没有任何变化。

有没有更好的方法解决这个问题？

Answer 1

我没有看到市长错误。矩阵的遍历顺序是可以的。

我认为您的代码可以使用这样的 openMP 缩减有效地并行计算

for(int q(0); q < Nk*Nk; q++){
    #pragma omp parallel for default(shared) reduction(+:chi)
    for(int k(0); k < Nk*Nk; k++){
        int kq = (k+q) % (Nk*Nk);
        cx_mat Gx = ((Eigveck.slice(k)).t())*(Vxk.slice(k)-Vxk.slice(kq))*Eigveck.slice(kq);
        cx_mat Gy = ((Eigveck.slice(k)).t())*(Vyk.slice(k)-Vyk.slice(kq))*Eigveck.slice(kq);
        vec ek = Eigvalk.col(k);
        vec ekq = Eigvalk.col(kq);
        for(int i(0); i < Ltot; i++){
            for(int j(0); j < Ltot; j++){
                chi = chi + (abs(Gx(i,j))*abs(Gx(i,j))+abs(Gy(i,j))*abs(Gy(i,j)))*(1.0/(1.0+exp(ekq(j)/T))-1.0/(1.0+exp(ek(i)/T)))*((ekq(j)-ek(i))/((ekq(j)-ek(i))*(ekq(j)-ek(i))+eta*eta))/(Nk*Nk);
            }
        }
    }
    double qx = (G1(0)*floor(q/Nk)/Nk+G2(0)*(q % Nk)/Nk);
    double qy = (G1(1)*floor(q/Nk)/Nk+G2(1)*(q % Nk)/Nk);
    lindhard << qx << "     " << qy << "     " << -chi << "    " << endl;
}

附言。

也许您定义了一些常量局部变量，例如

const auto delta = ekq(j)-ek(i);

您是如何测量热点的？

您使用哪些编译器选项？我假设您打开了适当的优化级别，对吧？