c++ armadillo 如何在for循环中加速调用一个巨大的向量

How to speed up the call to a huge vector in a for loop c++ armadillo

我正在使用犰狳库。我的程序太慢了,我需要加快速度的部分如下

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,但没有任何变化。

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

我没有看到市长错误。矩阵的遍历顺序没问题

我认为您的代码可以像这样使用 openMP reduction 高效并行计算

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;
}

Ps。

也许你定义了一些const局部变量,比如

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

您是如何测量热点的?

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