Let's start from scratch, please.
The naïve algorithm for multiplying an n-row, k-column matrix A and a k-row, m-column matrix B is
Code:
Cr,c = ∑i=1..k Ar,i Bi,c
and the result is an n-row, m-column matrix C.
For details and background, please look at the Wikipedia article on matrix multiplication. The Strassen algorithm is much more complicated, but only beats the naïve algorithm for larger matrices. In other words, the naïve approach is quite acceptable for this.
Given 16 eight-by-eight matrices M1 to M16, their product is
Code:
M1 × M2 × M3 × ... × M[sub]15[sub] × M16
= (((((((((((((((M1 × M2) × M3) × ... × M15) × M16
Or, in pseudocode:
Code:
T1 and T2 are temporary matrices, R is the result matrix
T1 = M1 × M2
T2 = T1 × M3
T1 = T2 × M4
T2 = T1 × M5
T1 = T2 × M6
T2 = T1 × M7
T1 = T2 × M8
T2 = T1 × M9
T1 = T2 × M10
T2 = T1 × M11
T1 = T2 × M12
T2 = T1 × M13
T1 = T2 × M14
T2 = T1 × M15
R = T2 × M16
So, you only need a function which multiplies two matrices.
Assume you have
Code:
struct matrix8x8 {
double element[8][8];
};
you then only need
Code:
void multiply8x8(struct matrix8x8 *const result,
const struct matrix8x8 *const left,
const struct matrix8x8 *const right)
{
double s;
int r, c, i;
for (r = 0; r < 8; r++) {
for (c = 0; c < 8; c++) {
s = left->element[r][0] * right->element[0][c];
for (i = 1; i < 8; i++)
s += left->element[r][i] * right->element[i][c];
result->element[r][c] = s;
}
}
}
To implement the pseudocode above for any number of 8x8 matrices, you could use
Code:
void product8x8(struct matrix8x8 *const result,
const struct matrix8x8 *const *const matrixarray,
const int count)
{
struct matrix8x8 temp1, temp2;
int n;
if (count < 1)
return;
if (count == 1) {
*result = *matrixarray[0];
return;
} else
if (count == 2) {
multiply8x8(result, matrixarray[0], matrixarray[1]);
return;
}
multiply8x8(temp1, matrixarray[0], matrixarray[1]);
for (n = 2; n + 2 < count; n += 2) {
multiply8x8(temp2, temp1, matrixarray[n]);
multiply8x8(temp1, temp2, matrixarray[n + 1]);
}
if (count & 1) {
multiply8x8(result, temp1, matrixarray[count - 1]);
} else {
multiply8x8(temp2, temp1, matrixarray[count - 2]);
multiply8x8(result, temp2, matrixarray[count - 1]);
}
return;
}
The loop does all even numbered multiplications except for the last one, and all odd numbered multiplications except the first and last. (I know that sounds odd, but if you look at the pseudocode, that's just the way it repeats.)
The final if statement takes care of the final multiplication (if odd number of matrices) or final two multiplications (if even number of matrices).
If the matrices were of different sizes, then it would be useful to pick the multiplication order in such a fashion to keep the total number of multiplications as small as necessary -- but I guess that is a subject for some other thread.