An assignment algorithm in C

Question

As a post processing step of a similarity matrix computation that leads to the non-negative matrix mt4_ in gpu, I am performing an assignment step to determine which element in rows/columns is most similar to which element in columns/rows. So this is more like a first-come-first-serve algorithm in which the largest similarity is picked up and then that particular row and column is excluded (by flagging the first element of that row and column --using -1.0f) and the process continues as such until all rows/columns are exhausted. The results are exported to two contiguous sub-arrays of sms one for each row-to-column and column-to-row assignment. Here is the function in which sZ is a const unsigned int equal to 5000 (mt4_, and by extension mt2, are of size sZ x sZ):

void asn(float *mt2, float *mt4_, unsigned int *sms) {

  unsigned int i, j;
  unsigned int a, b;
  register unsigned long k;
  register float mx;
  register unsigned int mm = 0, nn = 0;

  for (j=0; j<2; j++) {
    k = j*sZ;
    cudaMemcpy(mt2, mt4_, mtb, cudaMemcpyDeviceToHost);
    for (i=0; i<sZ; i++) {
      mx = -1.0f;
      for (a=0; a<sZ; a++) {
        if (mt2[a*sZ] != -1.0f)
          for (b=0; b<sZ; b++) {
            if (mt2[a*sZ+b] > mx && mt2[b] != -1.0f) {
              mx = mt2[a*sZ+b];
              mm = a;
              nn = b;
            }
          }
      }
      if (fabsf(mt2[mm*sZ+mm] - mt2[mm*sZ+nn]) <= thd)
        nn = mm;
      sms[k+mm] = nn+1;
      mt2[mm*sZ] = mt2[nn] = -1.0f;
    }
  }
}

Here thd is the const float threshold that introduces a bias towards assigning a column index to its identical row index if their similarity is a within a tolerance thd of the best available score (this is to compensate for numerical errors).

The code runs fine but my problem is that it takes too long (about 100 seconds on my compute node), and this is something I can't afford due to the large number of calls to this function.

I am new to C and it is very well probable that I am missing straightforward performance optimizations, so I would really appreciate if you would comment in case you see a window for improvement.

Background: Blondel et al., A measure of similarity between graph vertices

Answer 1

You may try to prepare an array of

struct {
    float value;
    int row;
    int col;
} data[sZ * sZ];

and two arrays

bool valid_column[sZ];
bool valid_row[sZ];

Sort data in the descending order, and scan it once, marking rows and columns invalid as you pick up the values. That shall drive the performance from \$O(sZ^3)\$ to \$O(sZ^2 \ln{sZ})\$.

Answer 2

All of what I'm suggesting here will be done by a good optimizing compiler. Some of what I'm suggesting is computer architecture dependent so you should depend on the compiler rather than doing it if you can. Generate assembly code and look at it to better optimize.

First, register assignment isn't guaranteed, it is only a recommendation to the compiler. You want to move your register variable declarations up. The compiler ignores your register allocations when it runs out of registers.

void asn(float *mt2, float *mt4_, unsigned int *sms) {
  register float mx;
  register unsigned int mm = 0, nn = 0;
  register unsigned int a, b;
  unsigned int i, j;
  unsigned long k;

You want to base your register allocations on what is being used most, since k is only used twice it shouldn't be a register. The variables mm, nn, a and b are used the most (sZ * sZ * sZ) so they should be registers.

On many computers counting down (auto decrement) is faster than counting up because decrement and test for zero is one opcode:

    for (i=sZ; --i; ) {
      mx = -1.0f;
      for (a=sZ; --a; ) {
        if (mt2[a*sZ] != -1.0f)
          for (b=sZ; --b; ) {
            if (mt2[a*sZ+b] > mx && mt2[b] != -1.0f) {
              mx = mt2[a*sZ+b];
              mm = a;
              nn = b;
            }
          }
      }

This does impact readability and maintainability so it should not be used unless the optimizing compiler doesn't do it for you.

Testing for equivalency of floating points is to be discouraged because of rounding errors (floating point errors).

Variable names and constants should be longer and clearer, sZ => MatrixSize.

Answer 3

"There's your problem!"

You say that sZ=5000, this means that:

for (j=0; j<2; j++) {
  ...
  for (i=0; i<sZ; i++) {
    mx = -1.0f;
    for (a=0; a<sZ; a++) {
      if (mt2[a*sZ] != -1.0f) // (1)
        for (b=0; b<sZ; b++) {
          if (mt2[a*sZ+b] > mx && mt2[b] != -1.0f) {
            mx = mt2[a*sZ+b];
            mm = a;
            nn = b;
          }

the inner most 3 lines of this code may be executed at most 2*5000^3=2.5E11 times. Provided that the inner loop will complete in say 10 instructions and you have 3GHz CPU with 5 IPC, the loop will take about: 2.5E11 * 10 / (5*3E9) = 5/3E2 = 166s to complete so I'm guessing the if statement (1) reduces that by about 40-ish%, give or take.

I don't believe that you can make this code faster by optimizing, and there is nothing that's inefficient in your C code.

You have to approach this by changing to a better algorithm. I unfortunately don't have time to help you with that right now but I hope that I have at least pointed you in the right direction.