# 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!

## 1 Answer

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 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 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 at 6:06