c++ - 使用 MKL 稀疏矩阵 vector 相乘

Question

有没有人有使用 MKL 稀疏矩阵 vector 乘法例程的简单 C++ 代码示例？我需要使用“mkl_zcsrsymv”将一个复数对称矩阵(存储在下三角矩阵中)与一个复数 vector 相乘，但我找不到这方面的单一示范性示例。无法为复杂情况编译我的代码。

Answer 1

在 Intel's homepage 上阅读 mkl_zcsrsymv 的文档.这里对称是按字面意思理解的! (不是Hermitian的意思)

使用 icpc -mkl test.c 编译以获得最大的便利。

#include <stdio.h>
#include "mkl_spblas.h"

int main()
{
  /* Matrix in CRS format
   *
   * { {  0,  0,  i } 
   *   {  0, -1,  2 } 
   *   {  i,  2,  0 } } 
   */
  int m = 3;
  MKL_Complex16 a[] = { {0,1}, {-1,0}, {2,0}, {0,1}, {2,0} };
  int ia[] = { 0, 1, 3, 5 };
  int ja[] = { 2, 1, 2, 0, 1 };

  MKL_Complex16 x[] = { {1,0}, {2,0}, {3,0} };
  MKL_Complex16 y[] = { {0,0}, {0,0}, {0,0} };

  char uplo = 'L';
  // Use MKL to compute
  // y = A*x
  mkl_cspblas_zcsrsymv(&uplo, &m, a, ia, ja, x, y);

  printf("y = { (%g,%g), (%g,%g), (%g,%g) }\n",
         y[0].real, y[0].imag,
         y[1].real, y[1].imag,
         y[2].real, y[2].imag
    );
}

输出为 y = { (0,3), (4,0), (4,1) }。查看 WolframAlpha .

---

这也是 mkl_dcsrmv 的示例.

#include <stdio.h>
#include "mkl_spblas.h"

int main()
{
  /* Matrix in CRS format
   *
   * { { 0,  0,  1 } 
   *   { 0, -1,  2 } 
   *   { 1,  0,  0 } } 
   */
  int m = 3;
  int k = 3;
  double val[] = { 1, -1, 2, 1 };
  int indx[]   = { 2, 1, 2, 0 };
  int pntrb[]  = { 0, 1, 3 };
  int pntre[]  = { 1, 3, 4 };

  double x[] = { 1, 2, 3 };
  double y[] = { 0, 0, 0 };

  double alpha = 1;
  double beta  = 0;
  char transa = 'N';
  char matdescra[] = {
    'G', // type of matrix
    ' ', // triangular indicator (ignored in multiplication)
    ' ', // diagonal indicator (ignored in multiplication)
    'C'  // type of indexing
  };

  // Use MKL to compute
  // y = alpha*A*x + beta*y
  mkl_dcsrmv(&transa, &m, &k, &alpha, matdescra, val, indx, pntrb, pntre, x, &beta, y);

  printf("y = { %g, %g, %g }\n", y[0], y[1], y[2]);
}

输出是 y = { 3, 4, 1 }。查看 WolframAlpha .

---

在玩这个的时候我发现它直接兼容 Armadillo .这使得在 C++ 中使用起来非常方便。这里我先用 Armadillo 生成一个随机对称矩阵并将其转换为稀疏矩阵。我将其与随机 vector 相乘。最后，我将结果与 Armadillo 的稀疏矩阵 vector 积进行比较。精度差异很大。

#include <iostream>
#include <armadillo>
#define MKL_Complex16 arma::cx_double
#include "mkl_spblas.h"

int main()
{
  /* Matrix in CRS format
   *
   * { {  0,  0,  i } 
   *   {  0, -1,  2 } 
   *   {  i,  2,  0 } } 
   */
  int dim = 1000;
  arma::sp_cx_mat a(arma::randu<arma::cx_mat>(dim,dim));
  a += a.st();
  arma::cx_vec x = arma::randu<arma::cx_vec>(dim);
  arma::cx_vec y(dim);

  char uplo = 'L';
  // Use MKL to compute
  // y = A*x
  mkl_cspblas_zcsrsymv(&uplo, &dim,
                       a.values, (int*)a.col_ptrs, (int*)a.row_indices,
                       x.memptr(), y.memptr());

  std::cout << std::boolalpha
            << arma::approx_equal(y, a*x, "absdiff", 1e-10)
            << '\n';
}