Sparse matrix multiplication for billion by billion matrices

I have to do multiplication of two billion by billion sparse matrices (on a CPU), hence any help or hints in optimizing the below given code would be extremely useful.

Note: I am only showing the code to compute number of non-zeros in the output matrix. Buffer based adjustment of memory is not an option for such a large matrix. Matrix multiplication will be performed in parts. The line with the problem is workarray[matrixb->colInd[k]] = 1 as it is repeating some calculations but I have not been able to come up with a method to reduce the computations.

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <omp.h>
struct sparsemat
{
int nzmax;
int rows;
int cols;
int *rowPtr;
int *colInd;
double *values;
};
struct darray
{
double * array;
int rows;
int cols;
};
struct sparsemat *create(int rows, int cols, int nzmax)
{
struct sparsemat *matrix;
matrix = (struct sparsemat*) calloc(1, sizeof (struct sparsemat));
matrix->rows = rows;
matrix->nzmax = nzmax;
matrix->cols = cols;
matrix->colInd = (int*) calloc(nzmax, sizeof (int));
matrix->rowPtr = (int*) calloc(rows + 1, sizeof (int));
matrix->values = (double*) calloc(nzmax, sizeof (double));
return matrix;
}
void _destroy2(struct sparsemat *matrix)
{
if (matrix != 0)
{
if (matrix->rowPtr != 0)
{
free(matrix->rowPtr);
matrix->rowPtr = NULL;
}
if (matrix->colInd != 0)
{
free(matrix->colInd);
matrix->colInd = NULL;
}
if (matrix->values != 0)
{
free(matrix->values);
matrix->values = NULL;
}
matrix->nzmax = 0; matrix->rows = 0; matrix->cols = 0; free(matrix); matrix = NULL;
}
}
void sparse_dual2(const struct sparsemat * const matrixa, const struct sparsemat * const matrixb, struct sparsemat * matrixc)
/* loop counters and scratch variables*/
{
int i, j, k, index;
double scalar;
/* Workarray that contains indicator of non zero entries*/
unsigned char *workarray = 0;
double *multiply = 0;
/* Initilize different fields of output sparse matrix without a function */
matrixc->rows = matrixa->rows;
matrixc->cols = matrixb->cols;
matrixc->rowPtr = (int*)calloc(matrixa->rows + 1, sizeof(int));
/**************Parallel Part *******************************************/
/* Compute number of non-zero entries in each row of the output matrix
and also number of non-zero entries in the output matrix that is C */
# pragma omp parallel private (i, j, k, index, scalar) firstprivate(workarray, multiply)
{
workarray = (unsigned char*)calloc(matrixb->cols, sizeof(unsigned char));
# pragma omp for
/*matrixa->rows i.e., number of rows in the output matrix */
for (i = 0; i<matrixa->rows; ++i)
{
index = i+1;
for (j = matrixa->rowPtr[i]; j <= matrixa->rowPtr[index] - 1; ++j)
{
# pragma GCC ivdep
for (k = matrixb->rowPtr[matrixa->colInd[j]]; k <= matrixb->rowPtr[matrixa->colInd[j] + 1] - 1; ++k)
workarray[matrixb->colInd[k]] = 1;
//workarray[matrixb->colInd[k]]^= (0 ^ 1) & 1; bit based
}
# pragma GCC ivdep
/* Vectorize Loop; Compute NNZ in each row of the output matrix*/
for (j = 0; j <matrixb->cols; ++j)
{
matrixc->rowPtr[index] += workarray[j];
workarray[j] = 0;
}
}
/* Set WorkArray to 0 in each thread*/
free(workarray); workarray = NULL;
}
}
int main()
{
struct sparsemat *matrixa = 0, *matrixb = 0, *matrixc = 0;
int nrows_a = 5;
int ncols_a = 10;
int nzmax_a = 19;
int nrows_b = 10;
int ncols_b = 5;
int nzmax_b = 12;

int i;
double start, finish;

int rowPtrA = { 0, 3, 5, 8, 14, 19 };
int colIndA = { 1, 7, 9, 5, 9, 7, 8, 9, 1, 4, 6, 7, 8, 9, 0, 1, 5, 6, 7 };
double valuesA = { 0.68, 0.99, 0.10, 0.65, 0.84, 0.08, 0.31, 0.97, 0.53, 0.47, 0.50, 0.07, 0.04, 0.45, 0.25, 0.93, 0.81, 0.64, 0.87 };

int rowPtrB = { 0, 1, 2, 2, 3, 5, 5, 7, 9, 10, 12 };
int colIndB = { 3, 1, 2, 0, 3, 3, 4, 0, 2, 0, 0, 1 };
double valuesB = { 0.78, 0.006, 0.79, 0.01, 0.92, 0.95, 0.56, 0.99, 0.82, 0.89, 0.17, 0.43 };

matrixa = create(nrows_a, ncols_a, nzmax_a);
for (i = 0; i < nzmax_a; ++i)
{
matrixa->colInd[i] = colIndA[i];
matrixa->values[i] = valuesA[i];
if (i <= nrows_a)
matrixa->rowPtr[i] = rowPtrA[i];
}
matrixb = create(nrows_b, ncols_b, nzmax_b);
for (i = 0; i < nzmax_b; ++i)
{
matrixb->colInd[i] = colIndB[i];
matrixb->values[i] = valuesB[i];
if (i <= nrows_b)
matrixb->rowPtr[i] = rowPtrB[i];
}

start = omp_get_wtime();
matrixc = (struct sparsemat*) calloc(1, sizeof (struct sparsemat));
sparse_dual2(matrixa, matrixb, matrixc);
finish = omp_get_wtime();

printf("Rows in C = %i \n", matrixc->rows);
printf("Cols in C = %i \n", matrixc->cols);
printf("Non Zeros in C = %i  \n", matrixc->nzmax);
printf("Time Taken = %f\n", finish - start);

_destroy2(matrixa);
_destroy2(matrixb);
_destroy2(matrixc);
getchar();
return 0;
}
• Are you sure you mean c++ and not c? – nwp Dec 31 '14 at 2:49
• It is C or C++ does not matter. Anyhow the code given above can be compiled by both a C and C++ compiler. Calloc rather than new is used for a reason as this code is part of a larger code that uses realloc – Vineet Yadav Dec 31 '14 at 17:12
• It matters for the code review. I started writing a review but then thought writing that C is bad C++ will not tell you anything useful. It makes more sense to re-tag this as C to avoid useless reviews. – nwp Dec 31 '14 at 19:48

A few things immediately stand out to me:

1. You aren't handling allocation failures. In create() the initial creation of the struct could fail. So could allocating any of the contained pointers. You need to deal with that case. The easiest thing is probably to check for failure of each allocation, and if any of them fail, free the already allocated ones and return NULL. This also applies to the allocation of workarray and matrixc.
2. Your naming is not very good. A name like darray may tell you it's an array of doubles, but it doesn't give you any clue as to what it's used for. (But apparently, it's not used - so why is it there?) Likewise, why call the destruction method _destroy2()? Was there a _destroy1()? If so, what's the difference between that and this? If not, why does this one end in 2? Same with sparse_dual2(). And the matrixc argument to sparse_dual2() should indicate that it's the return value. Maybe name it result or dst? (Some standard library calls, like memcpy() take the destination as the first argument, so it could be confusing.)
3. Instead of using const pointers to const data, would it make sense to use const references? Something like const struct sparsemat& foo? (Or would it need to be const struct sparsemat& const foo?) It seems like that could be simplified, though, I admit, I could be misunderstanding it.
4. I find your loop end condition on the 2 inner loops confusing. I assume it's saying, "The item that's 1 before the first item of the next row", so basically, the last item of this row? It would be clearer if you made a variable with the appropriate name and put that value into it. Or at least, change it to < matrixa->rowPtr[index] instead of <= that minus 1.
5. In terms of optimizing the particular line of interest, have you tried using pointers instead of array indices?

Something like this:

int* nextCol = matrixb->colInd[matrixb->rowPtr[matrixa->colInd[j]]];
for (k = matrixb->rowPtr[matrixa->colInd[j]]; k <= matrixb->rowPtr[matrixa->colInd[j] + 1] - 1; ++k, ++nextCol)
workarray[*nextCol] = 1;

You might even be able to keep a pointer to the next place in the work array, so something like this:

int* nextCol = &matrixb->colInd[matrixb->rowPtr[matrixa->colInd[j]]];
unsigned char* nextWorkCell = &workarray[*nextCol];
for (k = matrixb->rowPtr[matrixa->colInd[j]]; k <= matrixb->rowPtr[matrixa->colInd[j] + 1] - 1; ++k, ++nextCol)
{
nextWorkCell += nextCol;
*nextWorkCell = 1;
}

I'm not sure if that would be any better than your compiler's optimizer.