如何将C中的两个稀疏矩阵相乘

Question

我有一个稀疏矩阵D，我想把D_transpose和D相乘得到L，如下所示：

L=D'*D；

我使用稀疏LAS来处理稀疏矩阵，但文档中说没有什么可以乘以两个稀疏矩阵。

我完全被卡住了，不知道该怎么办。

D的尺寸通常在500000乘250000之间。我根本无法分配那么多内存，所以只能使用稀疏矩阵来完成。

我是用MATLAB做这件事的，我不明白如果MATLAB在接口下面也使用稀疏LAS，它是怎么做的——或者是这样？如果没有，它用什么？我也可以用它。

感谢您的阅读！

编辑：已解决。我需要L矩阵与一个向量相乘。因此，我没有首先计算L，而是简单地计算D'*（D*x），从而避免了两个稀疏矩阵相乘的需要。我现在只做稀疏矩阵和稠密向量乘法，这是由稀疏LAS支持的。

Answer 1

英特尔的库MKL支持两个稀疏矩阵的乘积。您可以从MKL的参考中检查例程mkl_sparse_spmm。（旧版本的MKL使用了程序mkl_?csrmultcsr，现在已弃用。）

此外，还有一个古老但广泛使用的开源线程安全稀疏矩阵库SPARSKIT（及其更新版本SPARSKIT2）。它们都是由F77编写的。有一个程序amub可以用来得到两个稀疏矩阵的乘积。您可以检查是否有人将此库重写为C.

Answer 2

它实际上在发布的文档中有说明。

第11页

5.2使用稀疏BLAS矩阵

一旦稀疏BLAS矩阵句柄已经完全构建完成（可以通过检查属性blas_valid_handle进行测试）使用矩阵句柄执行操作此时显示的四个操作支持表3.2和3.3中的除了使用稀疏BLAS矩阵执行操作之外通过其句柄查询其属性。表5.5列出了可以通过调用get属性例程获得。

表3.3第4页

USMM稀疏矩阵矩阵乘法

因此，支持似乎是存在的。我就是找不到BLAS_usmm函数的签名。也许你可以查一下标题。

编辑：如果您从NIST获得sparseLas，您可以检查blas_sparse_proto.h文件中的BLAS_*usmm函数的签名和参数。

 /* 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 3

据我所知，您的问题主要是在内存中存储巨大的矩阵。您可以将值存储在（行、列）对中。例如，

1 0 0
0 0 2
0 4 0

该矩阵可以存储在std::map<pair<int, int>, int>中，作为：

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

现在是计算部分。假设第一个矩阵存储在map1中，第二个矩阵存储于map2中，并且答案存储于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

如果你想在C.中使用类似地图的自定义数据结构，你需要使用类似地图