user2243369
user2243369

Reputation: 63

Recursive Block Matrix Multiplcation

Trying to implement block matrix multiplication recursively. It works fine for matrices of 2x2 but increase to sizes such as 4x4 and the answers differ vastly

Result of 3 for loops

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 

Result of recursion

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

If the amount of blocks within the matrix is greater than 4 I divide blocks into four larger ones and take the square root to get the new dimension like so then make the 8 recursive calls.

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];

}
}

The sr and sc are for starting row and col. The spacing should be right so I'm honestly out of leads here. Thanks in advanced.

Upvotes: 1

Views: 1593

Answers (2)

dexjq23
dexjq23

Reputation: 386

I've compiled and carefully debugged your code. If you only intend to use this function on matrices of 2^k*2^k, these 2 modifications will help.

First:

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];
    }
}

Second:

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);

Upvotes: 1

cgledezma
cgledezma

Reputation: 622

I believe your problem here is not as much on the implementation of your method, but on the loss of precision of floating-point operations. Sometimes one may think this imprecision is neglectable, but when we do intense operations over a floating point variable, like your triple nested loop, these imprecisions become significant.

One way to go around this is to scale your floating point numbers so they "lose" their decimal part. That is, for example, if you know your matrix won't have numbers with more than two decimal digits, then multiply them all by 100 and get their integer representation. Then perform the arithmetics on integers (which are precise), and in the end get the floating point representation of the result and divide it by 100.

Hope this helps.

Upvotes: 0

Related Questions