c - 递归 block 矩阵乘法

Question

尝试递归地实现分块矩阵乘法。它适用于 2x2 的矩阵，但增加到 4x4 等尺寸时答案差异很大

3个for循环的结果

1.53 0.89 0.53 1.33 
1.75 1.09 0.72 1.17 
1.78 1.43 0.57 1.69 
1.73 1.04 0.62 1.51 

递归结果

1.34 1.49 0.30 1.45 
2.02 1.93 0.79 1.30 
2.70 2.75 0.87 2.21 
1.81 1.84 0.59 1.47

如果矩阵中的 block 数量大于 4，我将 block 分成四个较大的 block 并取平方根以获得新的维度，然后进行 8 次递归调用。

void myRecMat(float** MatrixA, float** MatrixB, float** MatrixC, int srA, int scA, int srB, int scB, int srC, int scC, int blocks,int dim){
 if(blocks > 4) 
{ blocks=blocks/4;
      int newDim = dim/2;

        myRecMat(MatrixA,MatrixB,MatrixC, srA,scA,srB,scB,srC,scC,blocks,newDim);
        myRecMat(MatrixA,MatrixB,MatrixC, srA,scA+newDim,srB+newDim,scB,srC,scC,blocks,newDim);
        myRecMat(MatrixA,MatrixB,MatrixC, srA,scA,srB,scB+newDim,srC,scC+newDim,blocks,newDim);
        myRecMat(MatrixA,MatrixB,MatrixC, srA,scA+newDim,srB+newDim,scB,srC+newDim,scC,blocks,newDim);
        myRecMat(MatrixA,MatrixB,MatrixC, srA+newDim,scA,srB,scB,srC+newDim,scC,blocks,newDim);
        myRecMat(MatrixA,MatrixB,MatrixC, srA+newDim,scA+newDim,srB+newDim,scB,srC+newDim,scC,blocks,newDim);
        myRecMat(MatrixA,MatrixB,MatrixC, srA+newDim,scA+newDim,srB,scB+newDim,srC+newDim,scC+newDim,blocks,newDim);
        myRecMat(MatrixA,MatrixB,MatrixC, srA+newDim,scA+newDim,srB+newDim,scB+newDim,srC+newDim,scC+newDim,blocks,newDim); }                   
else
{
 int i,j,k,endR,endC;
 endR=srC+dim;
 endC=scC+dim;


 for(i=srC; i< endR; i++)
        for(j=scC;j< endC;j++)
            for(k=0; k<newDim; k++) 
                    c[i][j] += a[i][k]*b[k][j];

}
}

sr 和 sc 用于起始行和列。间距应该是正确的，所以老实说我在这里没有线索。提前致谢。

Answer 1

我已经编译并仔细调试了您的代码。如果您只打算在 2^k*2^k 的矩阵 上使用此函数，这两个修改会有所帮助。

首先:

for(i=srC; i< endR; i++) {
    for(j=scC;j< endC;j++) {
        for(k=0; k<newDim; k++) 
            /*c[i][j] += a[i][k]*b[k][j];*/
            c[i][j] += a[i][scA+k] * b[srB+k][j];
    }
}

第二个:

myRecMat(MatrixA,MatrixB,MatrixC,srA,scA,srB,scB,srC,scC,blocks,newDim);
myRecMat(MatrixA,MatrixB,MatrixC,srA,scA+newDim,srB+newDim,scB,srC,scC,blocks,newDim);
myRecMat(MatrixA,MatrixB,MatrixC,srA,scA,srB,scB+newDim,srC,scC+newDim,blocks,newDim); 
/*myRecMat(MatrixA,MatrixB,MatrixC,srA,scA+newDim,srB+newDim,scB,srC+newDim, scC,blocks,newDim);*/
myRecMat(MatrixA,MatrixB,MatrixC,srA,scA+newDim,srB+newDim,scB+newDim,srC, scC+newDim,blocks,newDim);
myRecMat(MatrixA,MatrixB,MatrixC,srA+newDim,scA,srB,scB,srC+newDim,scC,blocks,newDim);
myRecMat(MatrixA,MatrixB,MatrixC,srA+newDim,scA+newDim,srB+newDim,scB,srC+newDim,scC,blocks,newDim);
/*myRecMat(MatrixA,MatrixB,MatrixC,srA+newDim,scA+newDim,srB,scB+newDim,srC+newDim,scC+newDim,blocks,newDim);*/
myRecMat(MatrixA,MatrixB,MatrixC,srA+newDim,scA,srB,scB+newDim,srC+newDim,scC+newDim,blocks,newDim);
myRecMat(MatrixA,MatrixB,MatrixC,srA+newDim,scA+newDim,srB+newDim,scB+newDim,srC+newDim,scC+newDim,blocks,newDim);