How do I multiply two sparse matrices in C?

Question

I have a sparse matrix D, and I want to multiply D_transpose and D to get L as follows:

L = D'*D;

I am using sparseBLAS to deal with sparse matrices, but the documentation says there's nothing to multiply two sparse matrices.

I am completely stuck and have no idea how to proceed.

The dimensions of D are typically around 500,000 by 250,000. I cannot assign that much memory at all, so it just has to be done using sparse matrices.

I did this using MATLAB, and I don't understand how MATLAB does it if it also uses sparseBLAS underneath the interface - or does it? If not, what does it use? I could use that too.

Thank you for reading!

EDIT: Solved. I needed the L matrix to multiply with a vector. So instead of first calculating L, I simply did D'*(D*x), thus avoiding the need for multiplication of two sparse matrices. I now only do sparse matrix and dense vector multiplicatin, which is supported by sparseBLAS.

Answer 1

It actually is stated in the documentation posted.

Page 11

5.2 Using Sparse BLAS matrices

Once a Sparse BLAS matrix handle has been completely constructed (something that can be tested by checking the property blas_valid_handle), it is possible to use the matrix handle to perform operations. At this time the four operations shown in Tables 3.2 and 3.3 are supported. In addition to performing operations with a Sparse BLAS matrix, it is possible to query its properties through its handle. Table 5.5 lists the properties that can be obtained by calling the get properties routine.

Table 3.3 Page 4

USMM sparse matrix-matrix multiplication

So support seems to be there. I just can't find the signature for the BLAS_usmm function. Maybe you can check into the headers.

Edit: If you got your sparseBLas from NIST you can check the blas_sparse_proto.h file for the BLAS_*usmm functions for the signatures and parameters.

 /* Level 3 Computational Routines */

int BLAS_susmm( enum blas_order_type order, enum blas_trans_type transa,
    int nrhs, float alpha, blas_sparse_matrix A, const float *b, int ldb,
        float *c, int ldc );
int BLAS_dusmm( enum blas_order_type order, enum blas_trans_type transa,
        int nrhs, double alpha, blas_sparse_matrix A, const double *b,
        int ldb, double *c, int ldc );
int BLAS_cusmm( enum blas_order_type order, enum blas_trans_type transa,
         int nrhs, const void *alpha, blas_sparse_matrix A, const void *b, 
     int ldb, void *c, int ldc );
int BLAS_zusmm( enum blas_order_type order, enum blas_trans_type transa,
         int nrhs, const void *alpha, blas_sparse_matrix A, const void *b, 
     int ldb, void *c, int ldc );

Answer 2

As far as I can understand, your problem is mainly storing huge matrix in the memory. You can store the values in (row, column) pair. For example,

1 0 0
0 0 2
0 4 0

This matrix can be stored in a std::map<pair<int, int>, int> as:

map[make_pair(1, 1)] = 1
map[make_pair(2, 3)] = 2
map[make_pair(3, 2)] = 4

Now the computation part. Assuming that first matrix is stored in map1, second matrix is stored in map2 and the answer is stored in mapAns,

for each element x in map1:
    for each element y in map2:
        if x.column == y.row:
            mapAns[x.row, y.column] += x.value * y.value

You need to use map like custom data structure if you want to do the same in C.