# Multiplying a Matrix with its Transpose

Currently, I am in the process of optimizing a MIPS assembly program that takes a n x n matrix and multiplies it with its transpose. I am trying to optimize my matrix calculation algorithm so that it completes in as few clock cycles as possible. I have been given the A matrix, with its values stored in RAM. I must then calculate B=A*transpose(A).

There are a few caveats:

1. The matrix multiplication must be the dot product of the ith row of A and the jth column of B. It is not meant to be the element-wise multiplication. See the Wikipedia article.
2. I am not to make my algorithm more mathematically efficient then the un-modified example I will show below. I.e. I cannot exploit the symmetric nature resulting when you multiply a matrix with its transpose.

Here is the pseudo code example I have been given:

// Given array A which is unsigned int A[n*n] (ie word or 32 bit form)
// Reset array B which is unsigned int B[n*n] (ie word or 32 bit form)
for(int i = 0; i < (n * n); i++)
{
B[i] = 0;
}

// Matrix Multiplicaiton B = A*A'
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
{
for (int k = 0; k < n; k++)
{
B[i + n * j] = B[i + n * j] + A[i + n * k] * A[j + n * k];
}
}
}


Here is my attempt at optimizing the above example:

// Given array A which is unsigned int A[n*n] (ie word or 32-bit form)

// Matrix Multiplicaiton B = A*A'
for(int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
{
temp = 0;
n_times_i = n * i;

for (int k = 0; k < (n*n); k+=n)
{
temp += A[j + k] * A[i + k];
}

B[j + n_times_i] = temp;
}
}


As you can see, I have swapped things around to avoid unnecessary calculations where possible.

However, I was wondering whether anyone can see any other way of speeding things up? I.e. Cleverly swapping the order of the loops, etc.

Any help would be greatly appreciated!

Barring any compiler heroics, you are computing n*n a total of $$\n^3\$$ times. You might want to cache that result.

const int nn = n*n;


B[j + n_times_i] is a linearly increasing address location, given that j increases by 1 for each middle loop, and and i increases once for each outer loop, which is n increases of j. Taking advantage of that, you can skip the j + n*i calculation, and B[ ] indexing.

int *pB  = &B;
// ... loops & calculation of temp omitted for brevity.
*pB++ = temp;


Result:

const int nn = n*n;
int *pB = B;

for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
temp = 0;
for (int k = 0; k < nn; k += n) {
temp += A[i+k] * A[j+k];
}
*pB++ = temp;
}
}


You may find that you can get additional speed by using pointer arithmetic for A[i+k] and A[j+k].

int *pAi = A + i;
int *pAj = A + j;
for (k=0; k<n; k++) {   // Note: n.  The nn variable is no longer needed.
temp += *pAi * *pAj;
pAi += n;
pAj += n;
}


But, you’ll need to profile to find out. It depends on the number of free registers you have ... and the compiler/optimizers these days are pretty darn good.

• "compiler heroics" is quite a strong term for a very simple optimisation (even if we don't declare n as constant, all optimisers will spot it's unmodified and re-use the expression result). – Toby Speight May 31 '19 at 13:23
• @TobySpeight Shh! They might hear you. They are subtle and quick to anger, and I appreciate all heroics efforts they make on my behalf, big and small, and dare not anger them in any way, for without them I would be forced to wrestle with the dreaded realm and the demons of assembly, ARM, Thumb, MIPS and the brethren of x86, and the pony ... he comes! – AJNeufeld May 31 '19 at 14:26
• @TobySpeight, thanks so much for your help on the solution above. Quick question for you - What might the pointer arithmetic for *pAi and *pAj look like in MIPS assembly? – Sam Talbot Jun 1 '19 at 6:06