扰动相位振荡器积分

Question

我正在集成一个被扰动的相位振荡器系统。定义了方程组和Jacobian矩阵。我必须将系统状态的一维向量重塑为二维向量，然后进行矩阵产生。

有没有使程序更快的单线程方式？

using dim1 = std::vector<double>;
using dim2 = std::vector<std::vector<double>>;

dim2 multiply_matrices(dim2 &A, dim2 &B)
{
    int r1 = A.size();
    int c1 = A[0].size();
    int r2 = B.size();
    int c2 = B[0].size();

    assert(c1==r2 && "bad dimension for matrix multiplication");

    dim2 mult(r1, dim1(c2));

    for (int i = 0; i < r1; ++i) {
        for (int j = 0; j < c2; ++j) {
            mult[i][j] = 0;
        }
    }

    // Multiplying matrix A and B and storing in array mult.
    for (int i = 0; i < r1; ++i) {
        for (int j = 0; j < c2; ++j) {
            for (int k = 0; k < c1; ++k) {
                mult[i][j] += A[i][k] * B[k][j];
            }
        }
    }

    return mult;
}

class kuramoto_exe 
{
    vector<vector<double>> &adj;
    dim1 ω
    int kind;
    int N;
    double coupling;

public:
    kuramoto_exe (
        int N_, int kind_, double coupling_,
        dim1& omega_, dim2& adj_) : adj(adj_), omega(omega_), 
            kind(kind_), coupling(coupling_), N(N_) { }

        void operator() (const dim1 &X, dim1& dXdt, const double /*t*/)
    {


        dim1 x(X.begin(), X.begin() + N); // initial condition

        // reshaping:every colomn is a vector (reshape order "F")
        dim2 Y(N, dim1(N));
        for (int j=0; j<N; j++)
            for (int i=0; i<N; i++)
                Y[i][j] = X[N*j+i+N];  


        dim2 Jac(N, dim1(N));

        // definition of unperturbed model
        for (int i=0; i<N; i++) {
            double sumj =0.0;
            for (int j=0; j<N; j++) {
                if ((i != j) && (adj[i][j] > 1e-8))
                        sumj += adj[i][j] * sin(x[j]-x[i]);

                dXdt[i] =  omega[i]+ coupling * sumj;
            }
        }
        // calculation of the Jacobian 
        for (int i=0; i<N; i++) {
            for(int j=0; j<N; j++) {
                if (i!=j) 
                    Jac[i][j] = coupling * adj[i][j] * cos(x[j] - x[i]);
                else {
                    double sumj = 0.0;
                    for (int jj=0; jj<N; jj++) {
                        if (jj!=i)
                            sumj += adj[i][jj] * cos(x[jj] - x[i]);
                    }
                    Jac[i][j] = -1.0 * coupling * sumj;
                }
            }
        }

        dim2 M = multiply_matrices(Jac, Y);

        for (int j=0; j<N; j++)
            for (int i=0; i<N; i++)
                dXdt[N*j+i+N] = M[i][j]; // revese of reshape order "F"
    }
};

Answer 1

对于N = 1000 (如注释所建议的)，矩阵乘法在整个算法中花费的时间最多。这在很大程度上是因为它比它所需要的慢得多，但是即使有了艾根(见下文)，它仍然需要超过85%的时间。对于较低的N，这一比率也较低。

与可能的情况相比，矩阵乘法的实现非常缓慢。它是矩阵乘法基本定义的转录，没有任何优化应用。要从零开始编写一个好的实现，需要一个高级别的实现：

1. 高度优化的内环。
2. 调优贴图，以避免过多的缓存丢失。
3. 重新包装瓷砖，以避免过多的TLB失误。

第2点和第3点对于“足够大”的矩阵是必需的，1000x1000当然足够大，造成了许多额外的复杂性。

将第1点分开，这反过来意味着：

• 使用SIMD。忽略使用SIMD会立即使代码相对于可以实现的代码处于极大的劣势。
• 确保有足够的独立计算链。任何现代CPU都有很好的浮点计算吞吐量，从这个意义上讲，FP的计算并不慢。然而，这些操作单独需要花费大量的时间，因此获得高吞吐量的唯一方法是将独立操作的执行重叠起来。例如，在Haswell上，您至少需要10个独立的融合乘法-加法。
• 减少负载的数量。大多数现代CPU每个周期可以执行2次加载和2次FMA操作。这意味着，为了使其充满FMA，负载的数量不得超过FMA的数量。在实际应用中，进一步降低载荷与FMA的比值往往是很好的。
• 低级代码调优。

您可以手动完成这些操作(如果您需要，我可以向您展示一些东西)，但是不使用SIMD本质(甚至组装代码)的一种简单方法是使用现有的高效实现，例如从Eigen库或Intel MKL或其他竞争对手那里获得的实现。

这里使用特征的一种“最小更改”方法就是将Y和Jac转换为Eigen::MatrixXd对象，将它们相乘，并将结果提取为dXdt。也许最好完全避免dim2类型，避免重新打包，也因为嵌套向量并不是最有效的矩阵表示。

Answer 2

首先，确保您正在编译时启用了优化。

在multiply_matrices中，您不需要显式地用零填充矩阵，因为向量构造函数会这样做。

在kuramoto_exe构造函数中，成员初始化程序列表的成员变量顺序错误。它们将按照声明的顺序构造，因此N将在coupling之前构造，尽管它在列表后面出现。虽然这里没有问题，但如果成员变量之间存在依赖关系，则可能是这样，而且当您这样做时，一些编译器会发出警告。

在矩阵乘法中，将结果累加到局部变量，然后将结果存储在mult[i][j]中。这避免了多个索引计算的可能性。

在计算雅可比时，对角构件是水平构件之和的负值。你做的工作是你需要做的两倍。在处理行时积累和(跳过对角线元素)，然后在完成行(Jac[i][i] = -rowsum;)时存储对角线元素。

当您正在进行初始的整形(填充Y矩阵)时，最好将i循环放在外部，以便您的写入是连续的。CPU处理非顺序读取比非顺序写入更容易。

与其使用向量向量，不如考虑创建一个矩阵类，它连续地存储所有元素(在一个向量中)并计算适当的索引。这避免了当前的双重查找。