Una operación muy frecuente con matrices consiste en multiplicar una matriz por la traspuesta de otra. La implementación más sencilla de esto sería trasponer la segunda matriz y luego multiplicar la primera matriz por la traspuesta. Lamentablemente, trasponer una matriz es poco eficiente. No hay un patrón sencillo de acceso a la caché que sea bueno al mismo tiempo: hay que descomponer la matriz en trozos pequeños… en pocas palabras, un lío. Y resulta que hay una forma más sencilla de conseguirlo: invertir la posición de los índices al acceder a la segunda matriz. El código «naïve» para la multiplicación directa de dos matrices se resumía en un triple bucle:
for (int i = 0; i < m; i++)
for (int j = 0; j < p; j++)
{
double d = 0;
for (int k = 0; k < n; k++)
d += a[i, k] * b[k, j];
result[i, j] = d;
}
En la entrada sobre multiplicación de matrices, invertíamos el orden de los dos bucles interiores para mejorar el rendimiento de la caché de memoria. Además, añadíamos punteros e instrucciones AVX. Para la nueva rutina, lo que invertimos son los índices de la segunda matrix y, en consecuencia, renunciamos al truco de intercambiar los bucles internos. Perdemos algo de velocidad pero, en cualquier caso, el resultado es mejor que trasponer explícitamente y luego multiplicar.
La rutina que finalmente está incluida en Austra es la siguiente:
public unsafe Matrix MultiplyTranspose(Matrix m)
{
Contract.Requires(IsInitialized);
Contract.Requires(m.IsInitialized);
if (Cols != m.Cols)
throw new MatrixSizeException();
Contract.Ensures(Contract.Result<Matrix>().Rows == Rows);
Contract.Ensures(Contract.Result<Matrix>().Cols == m.Rows);
int r = Rows;
int n = Cols;
int c = m.Rows;
double[,] result = new double[r, c];
fixed (double* pA = values, pB = m.values, pC = result)
{
int lastBlockIndex = n & CommonMatrix.AVX_MASK;
double* pAi = pA;
double* pCi = pC;
for (int i = 0; i < r; i++)
{
double* pBj = pB;
for (int j = 0; j < c; j++)
{
int k = 0;
double acc = 0;
if (Avx.IsSupported)
{
Vector256<double> sum = Vector256<double>.Zero;
for (; k < lastBlockIndex; k += 4)
sum = Avx.Add(
left: Avx.Multiply(
left: Avx.LoadVector256(pAi + k),
right: Avx.LoadVector256(pBj + k)),
right: sum);
sum = Avx.HorizontalAdd(sum, sum);
acc = sum.ToScalar() + sum.GetElement(2);
}
for (; k < n; k++)
acc += pAi[k] * pBj[k];
pCi[j] = acc;
pBj += n;
}
pAi += n;
pCi += c;
}
}
return result;
}
Una parte importante del truco es hacer que el lenguaje que implementa Austra reconozca la multiplicación por la traspuesta y sepa usar la operación más eficiente. Pero eso es fácil de conseguir, y lo explicaré en otro momento.