Is there a point in transforming Sparse Matrix Multiplication into block form?

In an assignment for a parallel computing class we have been assigned to program Sparse Binary Matrix-Matrix multiplication (SpGEMM) in C. Julia has a relatively easy to follow implementation based on Gustavson's algorithm that works great.

Thing is we also need to do the multiplication in block form, which I already did, but I don't really see any speedup in doing so. From what I understand you're supposed to use the result of A(i,k)*B(k,j), where (i,j) are coordinates in the block matrix, as a mask/filter for the next block multiplication in the sum C(i,j) = Σ( A(i,k)*B(k,j) ).

Julia's implementation though, which I followed, already has a dense boolean array when computing each row that acts as a "flag" for when not to add something again in the resulting matrix.

My question is, is there any merit in turning this into block matrix multiplication or is there something that I might be doing wrong myself.

Keep in mind my C code currently runs in half the time Matlab takes in multiplying a 5,000,000 x 5,000,000 sparse matrix. The blocked version, which I really tried to optimize and I'm also doing in the Gustavson order, gets slower and slower the smaller the block-size is set.

Here is my current code

//C=D+(A*B) (basically OR)
bool SpGEMM_dor(int  *Acol, int *Arow, int An, 
               int  *Bcol, int *Brow, int Bm,
               int **Ccol, int *Crow, int *Csize,//output
               int  *Dcol, int *Drow)//previous
{
    //printCSR(Arow,Acol,An,An,An);
    int nnzcum=0;
    bool *xb = calloc(An,sizeof(bool)); //boolean flag
    for(int i=0; i<An; i++){
        int nnzpv = nnzcum;//nnz of previous row;
        Crow[i] = nnzcum;
        if(nnzcum + An > *Csize){ //make sure theres enough space
            *Csize += MAX(An, *Csize/4);
            *Ccol = realloc(*Ccol,*Csize*sizeof(int));
        }

        //---OR---
        //add previous row items in order to exist in the next block
        for(int jj=Drow[i]; jj<Drow[i+1]; jj++){ 
            int j = Dcol[jj];
            xb[j] = true;
            (*Ccol)[nnzcum] = j;
            nnzcum++;
        }
        //--------

        //add new row items
        for(int jj=Arow[i]; jj<Arow[i+1]; jj++){
            int j = Acol[jj];

            for(int kp=Brow[j]; kp<Brow[j+1]; kp++){
                int k = Bcol[kp];
                if(!xb[k]){
                    xb[k] = true;
                    (*Ccol)[nnzcum] = k;
                    nnzcum++;
                }
            }

        }

        if(nnzcum > nnzpv){
            quickSort(*Ccol,nnzpv,nnzcum-1);
            for(int p=nnzpv; p<nnzcum; p++){
                xb[ (*Ccol)[p] ] = false;
            }
        }

    }
    Crow[An] = nnzcum;

    free(xb);
    return Crow[An];
}

The part of code that I have inside of the ----OR---- section only happens in the block version in order to add the previous block to the now-calculating one. It basically does C = D+(A*B). I've also tried calculating the next block and then merging the 2 sorted arrays of each row of the 2 CSR matrices, which seems to be slower. Also all matrices are in CSR format.