Perturbed phase oscillator integration

I am integrating a system of perturbed phase oscillators. I define the system of equation and also the Jacobian matrix. I have to reshape the one dimensional vector of system state to a 2 dimensional vector then do the matrix production.

Is there any single-threaded way to make the program faster?

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
{
dim1 &omega;
int kind;
int N;
double coupling;

public:
kuramoto_exe (
int N_, int kind_, double coupling_,
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))

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

• About how big are the matrixes (different size classes require different kinds of optimization)? Are SIMD intrinsics within the scope of the question and if so which version? Mar 2, 2019 at 20:32
• Usually less than 1000 nodes, so I have (N*N)+N equation to solve. I don't know much about SIMD intrinsics. I think using odeint on Cuda is an option here. Mar 4, 2019 at 3:39

For N = 1000 (as suggested by the comment), the matrix multiplication takes by far the most time out of the whole algorithm. Largely because it is much slower than it needs to be, but even with Eigen (see below) it still takes over 85% of the time. For lower N that ratio is also lower.

The implementation of matrix multiplication is very slow compared to what is possible. It is a transcription of the basic definition of matrix multiplication without any optimization applied. To write a good implementation from scratch it would take, at a high level:

1. A highly optimized inner loop.
2. Tuned tiling, to avoid excessive cache misses.
3. Repacking tiles, to avoid excessive TLB misses.

Points 2 and 3 are necessary for "large enough" matrixes only, 1000x1000 certainly is large enough, creating a lot of additional complexity..

Splitting out point 1, that would in turn mean:

• Using SIMD. Neglecting to use SIMD immediately puts the code at a huge disadvantage compared to what is achievable.
• Ensuring there are enough independent chains of calculation. Any modern CPU has a good throughput for floating point calculations, in that sense FP math is not slow. However, these operations take a lot of time individually, so the only way to get high throughput is by overlapping the execution of independent operations. For example on Haswell you would need at least 10 independent fused multiply-adds.
• Reducing the number of loads. Most modern CPUs can do 2 loads and 2 FMAs per cycle. That means that in order to saturate it with FMAs, the number of loads must not exceed the number of FMAs. In practice it is often good to reduce the ratio of loads to FMAs even further.
• Low-level code tuning.

You could do these things manually (if you want I can show you some things), but a simple way to improve that without using SIMD intrinsics (or assembly code even) is by using an existing efficient implementation, for example from the Eigen library or Intel® MKL or some other competitor.

A "minimal changes" way to use Eigen here is just converting Y and Jac into Eigen::MatrixXd objects, multiplying them, and extracting the result into dXdt. It is probably better to avoid your dim2 type entirely though, to avoid some re-packing and also because nested vectors are already not the most efficient matrix representation.

First off, be sure you're compiling with optimizations enabled.

In multiply_matrices you don't need to explicitly fill your matrix with zeros, as the vector constructor will do that.

In the kuramoto_exe constructor, the member initializer list has the member variables in the wrong order. They will be constructed in the order they are declared, so N will be constructed before coupling despite its appearing later in the list. While not an issue here, it can be if there are dependencies among the member variables, and some compilers will issue a warning when you do this.

In your matrix multiply, accumulate the result to a local variable, then store the result in mult[i][j]. This avoids the potential of multiple index calculations.

When calculating the Jacobian, the diagonal members are the negative of the sum of the horizontal members. You're doing twice as much work as you need to. Accumulate the sum as you process the row (skipping the diagonal element) then store the diagonal element when you're done with the row (Jac[i][i] = -rowsum;).

When you're doing the initial reshaping (filling the Y matrix), it is probably better to put the i loop on the outside, so that your writes are sequential. It is easier for the CPU to handle nonsequential reads then nonsequential writes.

Rather than using a vector of vectors, consider creating a matrix class that stores all elements contiguously (in one vector) and computes the appropriate index. This avoids the double lookup you currently have.