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:

C_{r,c} = ∑_{i=1..k} A_{r,i} B_{i,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 M_{1} to M_{16}, their product is

Code:

M_{1} × M_{2} × M_{3} × ... × M[sub]15[sub] × M_{16}
= (((((((((((((((M_{1} × M_{2}) × M_{3}) × ... × M_{15}) × M_{16}

Or, in pseudocode:

Code:

T_{1} *and* T_{2} *are temporary matrices*, R *is the result matrix*
T_{1} = M_{1} × M_{2}
T_{2} = T_{1} × M_{3}
T_{1} = T_{2} × M_{4}
T_{2} = T_{1} × M_{5}
T_{1} = T_{2} × M_{6}
T_{2} = T_{1} × M_{7}
T_{1} = T_{2} × M_{8}
T_{2} = T_{1} × M_{9}
T_{1} = T_{2} × M_{10}
T_{2} = T_{1} × M_{11}
T_{1} = T_{2} × M_{12}
T_{2} = T_{1} × M_{13}
T_{1} = T_{2} × M_{14}
T_{2} = T_{1} × M_{15}
R = T_{2} × M_{16}

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.