# An assignment algorithm in C

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.

• Are you using an optimization flag? If you are, which one? – pacmaninbw Mar 16 '16 at 14:52
• Have you profiled the code to see where most of the time is spent? Profiling should show you if you are spending more time in cudaMemcpy() then in your own function. – pacmaninbw Mar 16 '16 at 15:12
• @pacmaninbw : Yes, but only the O3 flag. As for profiling, cudaMemcpy() is in the milliseconds category and I am only performing two copies; so I am not worried about that particular operation. – oligilo Mar 16 '16 at 15:21
• @Alien Herb Nite : Thanks for pointing that out, I corrected the title. – oligilo Mar 16 '16 at 15:25

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})$.

• Thank you for the suggestion; a rough version of this did actually occur to me at some point but I happened to forget about it not knowing the complexity of available/reachable sort algorithms in C at the moment. I will implement this and report back. – oligilo Mar 16 '16 at 17:46
• I finally got a chance to sit down and implement your suggestion and just wanted to state the obvious that it is significantly faster. Thank you so much! – oligilo Mar 28 '16 at 22:00

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.

• Most compilers will just plain ignore the register keyword, it's pretty useless. Link: stackoverflow.com/questions/10675072/… – Emily L. Mar 16 '16 at 16:25
• @EmilyL. That's absolutely true for C++, however, this is C and you might want to look at stackoverflow.com/questions/578202/register-keyword-in-c – pacmaninbw Mar 16 '16 at 16:33
• @pacmaninbw Thank you for the suggestions, I did make the changes you suggested and it turns out counting down the for loops actually slows down the code (at least in this case). The register keywords seem to make a few-second difference (I also added to the declaration of k at a later line). I thought the restrict keyword would also make an impact but apparently it doesn't here. – oligilo Mar 16 '16 at 17:36

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.

• Thank you for the answer; yes, you're right: the algorithm itself is expensive. People seem to suggest the Hungarian algorithm in the literature as a post processing step, and although the Wikipedia page was pointing to O(n^3) on the complexity of this assignment, I think there are faster implementations available. However, I don't see any particular reason why the Hungarian algorithm should be the 'relevant' assignment here and that is why I am pushing this FCFS approach. – oligilo Mar 16 '16 at 17:42