Compute weighted average around a point in a matrix quickly

Question

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.

Answer 1

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

See also

• 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.