How can I optimize a large number of copy and access operations in my numerical solver?

Question

I am working on a small pice of code which transports a density from one timestep to another iteratively until the destination time is reachted

The code I am working on has the following form:

1. initialize

2. ping-pong between two buffers now and next

   1. compute new refinement level for next
   2. for each value of the now vector some quantity gets addet to next.

My question: is there a general pattern how to speed up such a piece of code?

Side-question: can I implement this in a nicer/better way using stl algorithms?

---

std::vector<double> ping;
std::vector<double> pong;
ping.reserve(1000000);
pong.reserve(1000000);


std::vector* now= &ping;
std::vector* next = &pong;

Initialize(now);  // Fill first timestep
for(size_t t = 0; t < 1000; t++)  // timesteps
{
    size_t M = now->size();
    size_t N = calcNewRefinement(t, now);
    next->resize(N);
    for(size_t i = 0; i < N; i++) // all elements
    {
        for(size_t j = 0; j < now->size(); j++)
        {
            if (j > 0 && j < N)
                (*next)[i] += ExpensiveFunction((*now)[j-1], (*now)[j], (*now)[j+1], (*next)[i])
            else if (j == 0)
                (*next)[i] += ExpensiveFunction2((*now)[j], (*now)[j+1], (*next)[i])
            else if (j == M-1)
                (*next)[i] += ExpensiveFunction3((*now)[j-1], (*now)[j], (*next)[i])
        }
    }
    vector<double> *intermediate = now;
    now = next;
    next = intermediate;
}

Answer 1

The general advice for optimizations is to do profiling first. I assume you already did it, and found that your "copy and access operations" (as described in the question) have to be optimized.

In this case, let me note that the name ExpensiveFunction is misleading, because it cannot possibly be expensive when a few copy and access operations are so significant in your code.

The "general pattern" for optimizing is: look at your inner loop, and try to remove unnecessary operations.

In your case, you have the following there:

• for (...; j < now->size(); ...) - try replacing now->size() by M - there is a good chance your compiler already did it, but you never know...
• if (j > 0 && < j < N) - you can remove these checks completely if you separate your loop into 3 parts (first iteration; middle iterations; last iteration)
• now[j-1], now[j], now[j+1] - some c++ implementations insist on array bound checking for each access (this is not required by c++); if yours is like this, try disabling the check, or replace your std::vector by std::array or (if it doesn't help) by a C-style array
• next[i] = ... - as above
• You can try optimizing the code in your expensive functions instead...