# Compute weighted average around a point in a matrix quickly

I have a short snippet of code that simply computes a weighted average for surrounding elements in a square matrix. The actual implementation that I'm working on is not an average (more complex equation), but I am using this one to figure out how to handle various performance hinderances.

A few notes about my system:

• Windows 7, x64
• Visual Studio 2010, with the Intel C++ Compiler
• Profiling it using Intel VTune Amplifier
• Running in Release, with optimizations on (Intel compiler doesn't specify which level, but with comparison I did I think it's O2 or O3).
• The timing values I am mentioning I got from running it here with -O2
• DIM=512 and ITERATIONS=1000

Complete code is provided at the end, but for the next few explanations I will just focus on this loop of interest (it's the only loop, so...). I started off with a pretty straightforward, basic implementation:

for (int iter = 0; iter < ITERATIONS; iter++)
{
for (int x = 1; x < DIM-1; x++) // avoid boundary cases for this example
{
for (int y = 1; y < DIM-1; y++)
{
f0 = d_matrix[x][y];
f1 = d_matrix[x-1][y];
f2 = d_matrix[x+1][y];
f3 = d_matrix[x][y-1];
f4 = d_matrix[x][y+1];

d_res_matrix[x][y] = f0*0.6 + f1*0.1 + f2*0.1 + f3*0.1 + f4*0.1;
}
}
for (int x = 0; x < DIM; x++)
{
for (int y = 0; y < (DIM); y++)
{
d_matrix[x][y] = d_res_matrix[x][y];
}
}
}


This attempt took ~1.9s to execute. VTune suggested that I had problems with 4K Aliasing (read-before-write memory situations) specifically on lines for the y-loop and the memory writes that preceed them (in both loops). It also identified back-end bound, core-bound problems. I figured that the 4K aliasing problems might be causing the core-bound issues as well. To address this, I decided to get rid of the constant need to fetch x and y and rewrote the code using pointers:

for (int iter = 0; iter < ITERATIONS; iter++)
{
for (int x = 1; x < DIM-1; x++) // avoid boundary cases for this example
{
pf3 = &d_matrix[x][0];
pf0 = &d_matrix[x][1];
pf1 = &d_matrix[x-1][1];
pf2 = &d_matrix[x+1][1];
pf4 = &d_matrix[x][2];
save_write_loc = &d_res_matrix[x][0];
for (int y = 1; y < DIM-1; y++)
{
f0 = *pf0; pf0++;
f1 = *pf1; pf1++;
f2 = *pf2; pf2++;
f3 = *pf3; pf3++;
f4 = *pf4; pf4++;

*save_write_loc++ = f0*0.6+f1*0.1+f2*0.1+f3*0.1+f4*0.1;
}
}
for (int x = 0; x < DIM; x++)
{
s_m = &d_matrix[x][0];
s_r_m = &d_res_matrix[x][0];
for (int y = 0; y < (DIM; y++)
{
*s_m = *s_r_m;
s_m++; s_r_m++;
}
}
}


With this the execution time became ~1.3s. I ran VTune again. It once again complained about 4K Aliasing, but now it was about all the pointer dereferencing lines (such as f1 = *pf1). I figured it was probably caused by the the incrementing of the pointer I'm doing right after. It also suggested avoiding storing intermediate values, so I collapsed the loop into two lines as below. I also unrolled the 2nd loop 8 times to become:

for (int iter = 0; iter < ITERATIONS; iter++)
{
for (int x = 1; x < DIM-1; x++) // avoid boundary cases for this example
{
pf3 = &d_matrix[x][0];
pf0 = &d_matrix[x][1];
pf1 = &d_matrix[x-1][1];
pf2 = &d_matrix[x+1][1];
pf4 = &d_matrix[x][2];
save_write_loc = &d_res_matrix[x][0];
for (int y = 1; y < DIM-1; y++)
{
*save_write_loc++ = (*pf0)*0.6 + (*pf1)*0.1 + (*pf2)*0.1 + (*pf3)*0.1 + (*pf4)*0.1;
pf0++; pf1++; pf2++; pf3++; pf4++;
}
}
for (int x = 0; x < DIM; x++)
{
s_m = &d_matrix[x][0];
s_r_m = &d_res_matrix[x][0];
for (int y = 0; y < (DIM/8); y=y+8)
{
*s_m = *s_r_m;
s_m++; s_r_m++;
*s_m = *s_r_m;
s_m++; s_r_m++;
*s_m = *s_r_m;
s_m++; s_r_m++;
*s_m = *s_r_m;
s_m++; s_r_m++;
*s_m = *s_r_m;
s_m++; s_r_m++;
*s_m = *s_r_m;
s_m++; s_r_m++;
*s_m = *s_r_m;
s_m++; s_r_m++;
*s_m = *s_r_m;
s_m++; s_r_m++;
}
}
}


This cut the execution time further to ~0.7s. (I also tried unrolling it 4 and 16 times, but those were slower, so I settled on 8).

I am now stuck. I don't really know what else I can do to make it go any faster (if anything, but I'm sure there is something).

VTune is still complaining about 4K Aliasing in the *save_write_loc++ =... line. Maybe that is caused by the pointer increments happening right after since it is such a tight loop? That same line is still triggering back-end bound, core-bound port utilization problems. Since there is so much going on (multiplications, additions, fetches, stores), I don't really know which part is causing the problem exactly.

The complete code that can be compiled is here.

I'm thinking of having a 1D array instead of a 2D matrix. In that case, the locations will be next to each other, and perhaps they can be cached more efficiently. I am going to try this and report back. But I would appreciate any sort of suggestions on how to make this code faster.

• You could just remove your previous attempts and just keep the latest one. Apr 10, 2015 at 21:01
• I've been messing with your code to try to make it faster but haven't had much luck so far. I do have a suggestion (not code related) if what you're trying to do is numerical relaxation (Jacobi method in this case), and that is to look at some other methods like Gauss–Seidel or conjugate gradient which are usually faster. Apr 10, 2015 at 22:33
• Yes, make it a 1D array, and prefer a linear access pattern (touch array indices next to each other, in linear order).
– Juho
Apr 11, 2015 at 12:52

Since I do not know what the actual implementation will be I can only comment on the code as it is presented. From the algorithm it appears that you are applying a filter multiple times in order to monitor some propagation effect - as opposed to applying the filter multiple times to stabilize time profiling. From here on out I will assume that running ITERATIONS times is part of the core algorithm.

# General C++ concepts

• Initializing memory

Are you sure that it is legal to memset double types to zero to obtain logical $0.0$? Me neither. Use std::fill instead or better yet if you have C++03 or newer use value initialization so that you can initialize the arrays to zero with double *array = new double[n]().

If T is an array type, each element of the array is value-initialized

...

double f = double(); // non-class value-init, value is 0.0

• Use std::vector so that you do not have to do the memory allocation yourself.

Just be sure to use the -O2 or higher compiler flag alongside or else vector is slower than raw arrays (in my experience).

# Bugs

Your optimized routine has a couple major bugs. I always recommend checking the output between original code and optimized code should you need to perform optimization.

• Shift bug

save_write_loc = &d_res_matrix[x][0];


You are writing output values from column 1 to column 0. This would be fine if you ran the algorithm once but since you are running this algorithm ITERATIONS times, it results in shifting the output out multiple times. In general, on iteration $i$ you are writing the $i^{th}$ column of the original d_matrix to column 0 of d_res_matrix. To fix this simply keep the matrices aligned:

save_write_loc = &d_res_matrix[x][1];

• Skipped output indices bug

for (int y = 0; y < (DIM/8); y=y+8)


When you unrolled the loop you both divided the bounds by $8$ and increased the step to $8$. You only wanted to do one or the other. You should prefer changing the bounds since DIM/8 can be computed at compile time:

for (int y = 0; y < (DIM/8); y++)


The code runs much slower with these bug fixes.

# Quick Optimization

The algorithm cannot be performed in-place so you are using d_res_matrix to store the output of applying your filter to d_matrix. But then you want the output to be placed back in d_matrix so you perform a deep copy.

However a swap of the pointers would suffice so you could use this instead:

std::swap(d_matrix, d_res_matrix);


This has a major impact on the performance of your original code, but as you will see we can do better.

# Problematic Optimization

Flattening the 2D array to 1D sounds like a good idea until you factor in the cost of post-processing. When the array is processed as 1D the border pixels change values and, in fact, exhibit border-wrap effects.

We often fix problems like this by wrapping the matrix inside a false border - i.e. we add a $1\times1$ border all around the matrix. Unfortunately this would not work here. Instead, you can fix this by applying the filter and then going back and zeroing out border pixels each iteration.

# More Optimizations

• Use the -O3 flag before attempting generic optimizations.

Readability often goes out the window when applying optimizations to source code. The built-in compiler flags can be used instead. As you will see below, applying the -O3 flag gives me comparable results to unrolling loops, using pointer-based indexing, and flattening the array to 1D.

• Consider a generic matrix element $\textrm{matrix}[a][b]$. Let's determine how many times $\textrm{matrix}[a][b]$ is multiplied by $0.1$.

$\textrm{matrix}[a][b]$ has four 4-connected neighbors: $$\textrm{matrix}[a-1][b],\enspace \textrm{matrix}[a+1][b],\enspace \textrm{matrix}[a][b-1],\enspace \textrm{matrix}[a][b+1]$$ Each of these neighbors multiplies $\textrm{matrix}[a][b]$ by $0.1$. Therefore, $\textrm{matrix}[a][b]$ is multiplied by $0.1$ four times. Instead we could cache this computation once and reuse it when we need it. This requires $\mathcal{O}(n)$ space, however, since the overall algorithm cannot be done in-place, we can use some space we already were going to need.

(Technically the original algorithm only needs one row of extra space and the cached multiply algorithm only needs three rows of extra space.)

• The real optimization here is parallelization. If we had $\mathrm{DIM} \cdot \mathrm{DIM}$ processors then each processor could be charged with applying the filter to its pixel in any order we choose (or rather don't choose). Hence there is zero serial code per iteration. Of course we would not actually want that many processors/threads due to communication costs. If you are interested in running this algorithm in parallel you could use OpenMP, MPI, boost::thread, etc.

# Code and Timings

Test machine - CentOS 6.5, libstdc++-4.4.7, -O3, Intel Core i7-2670QM (2.20 GHz)

• Your original code fixed, std::swap applied, unrolled interior loop, other minor changes
Runtime: 1.54 seconds

#include <time.h>
#include <stdio.h>
#include <cstdlib>
#include <stdlib.h>
#include <vector>
#include <cstring>

#define DIM 512
#define ITERATIONS 1000

#define START_TIMING_ND clock_t t2; t2=clock();
#define STOP_TIMING_ND {long int final_nd=clock()-t2; printf("NEW/DELETE took %li ticks (%f seconds) \n", final_nd, ((float)final_nd)/CLOCKS_PER_SEC);}

int main(void)
{
// new/delete
double ** d_matrix, ** d_res_matrix;

d_res_matrix = new double * [DIM];
d_matrix = new double * [DIM];
for (int i = 0; i < DIM; i++)
{
d_matrix[i] = new double [DIM]();
d_res_matrix[i] = new double[DIM]();
}
d_matrix[20][45] = 1; // start somewhere

// vector calculations
double * save_write_loc;
double * pf0, *pf1, *pf2, *pf3, *pf4;
double * s_m, * s_r_m;
START_TIMING_ND;

for (int iter = 0; iter < ITERATIONS; iter++)
{
for (int x = 1; x < DIM-1; x++) // avoid boundary cases for this example
{
pf3 = &d_matrix[x][0];
pf0 = &d_matrix[x][1];
pf1 = &d_matrix[x-1][1];
pf2 = &d_matrix[x+1][1];
pf4 = &d_matrix[x][2];
save_write_loc = &d_res_matrix[x][1];
for (int y = 0; y < (DIM-2)/2; y++)
{
*save_write_loc++ = (*pf0++)*0.6 + (*pf1++)*0.1 + (*pf2++)*0.1 + (*pf3++)*0.1 + (*pf4++)*0.1;
*save_write_loc++ = (*pf0++)*0.6 + (*pf1++)*0.1 + (*pf2++)*0.1 + (*pf3++)*0.1 + (*pf4++)*0.1;
}
}
std::swap(d_matrix, d_res_matrix);
}
STOP_TIMING_ND;

for(int i = 0; i < DIM; i++)
{
for(int j = 0; j < DIM; j++)
{
//printf("%lf   ", d_matrix[i][j]);
}
//printf("\n");
}

// delete dynamic stuff
for (int i = 0; i < DIM; i++)
{
delete [] d_matrix[i];
delete [] d_res_matrix[i];

}
delete [] d_matrix;
delete [] d_res_matrix;

return 0;
}




• Cached multiplication code
Runtime: 1.10 seconds

#include <time.h>
#include <stdio.h>
#include <cstdlib>
#include <stdlib.h>
#include <vector>
#include <cstring>

#define DIM 512
#define ITERATIONS 1000

#define START_TIMING_ND clock_t t2; t2=clock();
#define STOP_TIMING_ND {long int final_nd=clock()-t2; printf("NEW/DELETE took %li ticks (%f seconds) \n", final_nd, ((float)final_nd)/CLOCKS_PER_SEC);}

int main(void)
{
// new/delete
double ** d_matrix, ** cmatrix;

cmatrix = new double * [DIM];
d_matrix = new double * [DIM];
for (int i = 0; i < DIM; i++)
{
d_matrix[i] = new double [DIM]();
cmatrix[i] = new double[DIM]();
}
d_matrix[20][45] = 1; // start somewhere

// vector calculations
double * save_write_loc;
double * pf0, *pf1, *pf2, *pf3, *pf4;
START_TIMING_ND;

for (int iter = 0; iter < ITERATIONS; iter++)
{
// Store 0.1 * d_matrix
for(int i = 0; i < DIM; i++)
{
for(int j = 0; j < DIM; j++)
{
cmatrix[i][j] = 0.1 * d_matrix[i][j];
}
}

for(int i = 1; i < DIM-1; i++)
{
for(int j = 1; j < DIM-1; j++)
{
d_matrix[i][j] = 0.6 * d_matrix[i][j]
+ cmatrix[i-1][j]
+ cmatrix[i+1][j]
+ cmatrix[i][j-1]
+ cmatrix[i][j+1];
}
}
}
STOP_TIMING_ND;

for(int i = 0; i < DIM; i++)
{
for(int j = 0; j < DIM; j++)
{
//printf("%lf   ", d_matrix[i][j]);
}
//printf("\n");
}

// delete dynamic stuff
for (int i = 0; i < DIM; i++)
{
delete [] d_matrix[i];
delete [] cmatrix[i];

}
delete [] d_matrix;
delete [] cmatrix;

return 0;
}




• Cached multiplication, flattened std::vector, iterator indexing
Runtime: 1.08 seconds (the unrolled code is slightly slower for me)

#include <time.h>
#include <cstdlib>
#include <vector>
#include <iterator>
#include <stdio.h>

#define DIM 512
#define ITERATIONS 1000

#define START_TIMING_ND clock_t t2; t2=clock();
#define STOP_TIMING_ND {long int final_nd=clock()-t2; printf("NEW/DELETE took %li ticks (%f seconds) \n", final_nd, ((float)final_nd)/CLOCKS_PER_SEC);}

int main(void)
{
std::vector<double> cmatrix(DIM*DIM, 0.0);
std::vector<double> d_matrix (DIM*DIM, 0.0);

d_matrix[20*DIM+45] = 1; // start somewhere

// vector calculations
START_TIMING_ND;

for (int iter = 0; iter < ITERATIONS; iter++)
{
std::vector<double>::iterator dit = d_matrix.begin();
std::vector<double>::iterator cit = cmatrix.begin();

// Store 0.1 * d_matrix;
while(dit != d_matrix.end())
{
*cit++ = 0.1 * *dit++;
}

std::vector<double>::iterator pf0 = d_matrix.begin() + DIM+1;  // [1][1]
std::vector<double>::iterator pf1 = cmatrix.begin() + DIM;     // [1][0]
std::vector<double>::iterator pf2 = cmatrix.begin() + DIM+2;   // [1][2]
std::vector<double>::iterator pf3 = cmatrix.begin() + 1;       // [0][1]
std::vector<double>::iterator pf4 = cmatrix.begin() + 2*DIM+1; // [2][1]

std::vector<double>::iterator stop = d_matrix.begin() + DIM*DIM - DIM - 1; // [DIM-2][DIM-1]

while(pf0 != stop)
{
*pf0 = (*pf0++)*0.6 + (*pf1++) + (*pf2++) + (*pf3++) + (*pf4++);
}

dit = d_matrix.begin();
std::vector<double>::iterator deit = d_matrix.begin() + DIM - 1; // [0][DIM-1]

// Post-processing - Zero out the border pixels
while(dit != d_matrix.end())
{
*dit = 0.0;
*deit = 0.0;
dit = deit+1;
deit += DIM;
}
}
STOP_TIMING_ND;

for(int i = 0; i < DIM*DIM; i++)
{
//printf("%lf   ", d_matrix[i]);

if(not ((i+1) % DIM))
{
//printf("\n");
}
}

return 0;
}


• Separable filter
The filter in this algorithm is:

$f = \begin{bmatrix} 0 && 0.1 && 0 \\ 0.1 && 0.6 && 0.1 \\ 0 && 0.1 && 0\end{bmatrix}$

It is not separable, however, the cached multiplication above is even better than a separable filter. If a separable filter is large enough then applying the filter as two 1D filters is several times faster than applying it as a 2D filter.

• The timings above were verified with CLOCK_PROCESS_CPUTIME_ID Apr 17, 2015 at 21:30
• Wow, thank you so much! I knew I'd balls up the unrolling lol, the results did seem a bit too good to be true. You are right, I am simulating propagation. I will look at the separable filter, as it seems like it would work quite well for my application. I am rather surprised at the performance you got using <vector>. I did try it as well (I also tried malloc), but it never worked as well as new/delete or even malloc, which I found strange.
– Mewa
Apr 17, 2015 at 21:49
• I ended up trying OpenMP yesterday as well (on my original, fetch-x-y-everytime-nothing-unrolled implementation) and it was quite quick, as I expected. I am contemplating putting it on a GPU, but I think the memory overhead to copy everything over and back might be too much. Though I think if I do lots of iterations on a single thread (100?) it might be worth it. Theoretically DIM is >> 512. Oh, also the compiler I use doesn't explicitly specify O3/O2, but it is O3 as I found in documentation.
– Mewa
Apr 17, 2015 at 21:50
• @Mewa Yes that performance was what I got on my i7. My personal development machine only has a Intel T7400 Core 2 Duo in which your code (with fixes) was by far the fastest. As for parallel computing, I'm primarily an MPI guy but that is only because I have access to a computing cluster. Regarding separable filters, unfortunately most filters are not so I'm not sure how promising that will be. If you want help determining if your filter is separable or not just let me know. Apr 20, 2015 at 20:27
• IEEE 754 guarantees that a +0.0 double is represented in memory as 64 bits of zeroes. Jun 8, 2015 at 22:39