我建议您查看 Eigen C++ 矩阵类库 http://eigen.tuxfamily.org/index.php?title=Main_Page。

我用你所说的大小和稀疏度的稀疏矩阵做了一个小测试，在我的中等功率 Windows 机器上，每个矩阵乘法大约需要 1 毫秒。

我的实验代码如下。如您所见，大部分代码用于设置测试矩阵。实际的矩阵乘法是一个简单的单行。（我将矩阵乘以 500 次以使时间更可靠。）

#include <iostream>
using std::cout;
using std::endl;
#include <random>

#include <Eigen/Core>
#include <Eigen/Sparse>

#include <boost/timer.hpp>

int main()
{
  const int n = 100000;
  const int nnz = static_cast<int>(.001*n);
  typedef Eigen::Triplet<double> Triplet;
  std::vector<Triplet> at(nnz), bt(nnz);
  std::random_device rd; // obtain a random number from hardware
  std::mt19937 eng(rd()); // seed the generator
  std::uniform_int_distribution<> distr(0, nnz-1); // define the range
  for (int i = 0; i < nnz; i++) {
    int ii = distr(eng), jj = distr(eng);
    at[i] = Triplet(ii, jj, 1);
    bt[i] = Triplet(jj, ii, 3);
  }
  Eigen::SparseMatrix<double> as(n, n), bs(n, n), cs(n, n);
  as.setFromTriplets(at.begin(), at.end());
  bs.setFromTriplets(bt.begin(), bt.end());
  boost::timer timer;
  const int nt = 500;
  for (int i = 0; i < nt; i++)
     cs = as.cwiseProduct(bs);
  cout << "elapsed time=" << timer.elapsed()/(double) nt  << endl;
}

如果可以选择编程语言，一种选择是使用Julia，它内置了对稀疏矩阵的支持（通过 Suitsparse）。在我的笔记本电脑上计时大约为 1.5 毫秒，您可以使用交互式动态语言，这在某些情况下可能很有用。

julia> a=sprand(1000, 1000, 0.1)
1000x1000 sparse matrix with 99749 Float64 entries:
    [5   ,    1]  =  0.725824
    [17  ,    1]  =  0.420022
    [19  ,    1]  =  0.404282
    [21  ,    1]  =  0.0307138
    [52  ,    1]  =  0.453376
    [55  ,    1]  =  0.30054
    [69  ,    1]  =  0.360203
    [74  ,    1]  =  0.346881
    [94  ,    1]  =  0.312849
    ⋮
    [932 , 1000]  =  0.978966
    [933 , 1000]  =  0.149551
    [954 , 1000]  =  0.417852
    [959 , 1000]  =  0.722707
    [964 , 1000]  =  0.519931
    [967 , 1000]  =  0.567152
    [971 , 1000]  =  0.964192
    [979 , 1000]  =  0.88494
    [987 , 1000]  =  0.286723
    [988 , 1000]  =  0.24282

julia> b=sprand(1000, 1000, 0.1)
1000x1000 sparse matrix with 99998 Float64 entries:
    [1   ,    1]  =  0.920533
    [3   ,    1]  =  0.879179
    [7   ,    1]  =  0.267203
    [25  ,    1]  =  0.522407
    [34  ,    1]  =  0.656031
    [41  ,    1]  =  0.280885
    [44  ,    1]  =  0.735824
    [68  ,    1]  =  0.433098
    [69  ,    1]  =  0.124862
    ⋮
    [932 , 1000]  =  0.505959
    [939 , 1000]  =  0.983413
    [947 , 1000]  =  0.418157
    [949 , 1000]  =  0.884657
    [963 , 1000]  =  0.412645
    [964 , 1000]  =  0.544348
    [966 , 1000]  =  0.709398
    [983 , 1000]  =  0.260483
    [989 , 1000]  =  0.1218
    [1000, 1000]  =  0.468975

julia> a.*b
1000x1000 sparse matrix with 9876 Float64 entries:
    [69  ,    1]  =  0.0449757
    [102 ,    1]  =  0.0340867
    [137 ,    1]  =  0.0794594
    [247 ,    1]  =  0.108002
    [376 ,    1]  =  0.248346
    [609 ,    1]  =  0.241789
    [633 ,    1]  =  0.224115
    [658 ,    1]  =  0.379804
    [754 ,    1]  =  0.272618
    ⋮
    [224 , 1000]  =  0.0122434
    [301 , 1000]  =  0.163899
    [309 , 1000]  =  0.0972784
    [403 , 1000]  =  0.0245688
    [659 , 1000]  =  0.0801249
    [700 , 1000]  =  0.158823
    [760 , 1000]  =  0.388442
    [926 , 1000]  =  0.193808
    [932 , 1000]  =  0.495317
    [964 , 1000]  =  0.283024

julia> @time a.*b
  0.001649 seconds (25 allocations: 3.056 MB)