Ho una funzione per il calcolo del prodotto matriciale:
codice:
double *gemm(double *mat1,int N, int M, double *mat2, int K)
{
    register double temp;
    double *result;
    result=(double*)malloc(N*K*sizeof(double));
    for(int i=0;i<N;i++)
    {
        for(int j=0;j<K;j++)
        {
            temp=0.0;
            for(int k=0;k<N;k++)
            {
                temp+=mat1[k+i*M]*mat2[k*K+j];
            }
            result[i*K+j]=temp;
        }
    }
    return result;
}
I parametri sono: mat1 e mat2 le matrici, N il numero di righe della matrice mat1, M il suo numero di colonne (e quindi anche il numero di righe della mat2), K il numero di colonne di mat2.
Ho messo temp come register double appunto per velocizzare i conti.
Ma il problema è che se le matrici hanno dimensioni grandissime, sicuramente una riga di mat1 non c' entra tutta nella cache.
Per cui dovrei moltiplicarla a blocchi, facendo preferibilmente i conti su elementi adiacenti in modo da "beccare" quelli che ci sono in cache.
Come posso fare per ottenere sempre il prodotto matriciale, ma moltiplicando prima gli elementi adiacenti?

PS:Le matrici sono memorizzate per righe in un array.