经验首页 前端设计 程序设计 Java相关 移动开发 数据库/运维 软件/图像 大数据/云计算 其他经验
当前位置:技术经验 » 大数据/云/AI » 人工智能基础 » 查看文章
MKL稀疏矩阵运算示例及函数封装
来源:cnblogs  作者:GeoFXR  时间:2023/4/24 8:59:00  对本文有异议

Intel MKL库提供了大量优化程度高、效率快的稀疏矩阵算法,使用MKL库的将大型矩阵进行稀疏表示后,利用稀疏矩阵运算可大量节省计算时间和空间,但由于MKL中的原生API接口繁杂,因此将常用函数封装,便于后续使用,最后在实际例子中调用接口执行想要的矩阵运算。

0 稀疏矩阵

稀疏矩阵是指矩阵中大部分元素为0的矩阵。存储这些0值数据会耗费大量的存储空间,并且计算时也会产生不必要的时间浪费。为了更有效地存储和处理这种类型的矩阵,有几种不同的稀疏矩阵格式。

下面是几种常见的稀疏矩阵格式:

  • COO格式:COO格式(坐标格式)用三个数组存储非零元素的行、列索引以及值。
  • CSR格式:CSR格式(压缩行格式)用三个数组存储矩阵的非零元素值、列索引和行指针。行指针数组指示每行中第一个非零元素的位置。
  • CSC格式:CSC格式(压缩列格式)与CSR格式类似,但是是按列存储非零元素。
  • DIA格式:DIA格式(对角线格式)使用一个二维数组存储非零元素。数组中的每一行表示矩阵的一个对角线,并且只有矩阵中存在的对角线上的元素才被存储。
  • BSR格式:BSR格式(块压缩行格式)用四个数组存储矩阵的非零元素。其中三个数组与CSR格式相同,第四个数组存储块的大小。
  • ELL格式:ELL格式(行程格式)使用两个数组存储矩阵的非零元素。其中一个数组存储元素的值,另一个数组存储元素在每行中的位置。每行中最大非零元素数量相同。

MKL中主要用到的稀疏矩阵格式有COOCSRCSC(与CSR类似)三种,以下将简要介绍COO格式与CSR格式:

(1)COO(Coordinate,坐标格式)

也被称为三元组格式,在 COO 格式中,每一个非零元素都用一个三元组 (row, column, value) 来表示,其中 row 和 column 分别代表该元素所在的行和列的索引,value 则代表该元素的值。由于 COO 格式中的非零元素的存储是无序的,因此在进行矩阵向量乘法等操作时,需要对 COO 格式进行排序。

COO 格式的优点:非常简单直观,易于理解和实现,同时可以处理任意稀疏度的矩阵。

缺点:存储开销较大,需要存储每个非零元素的行列索引,同时由于无序存储的缘故,在进行一些稀疏矩阵的计算时会需要排序,因此在效率上可能不如其他稀疏矩阵格式。

例:

(2)CSR(Compressed Sparse Row,行压缩格式)

常用于稀疏矩阵的存储和计算,CSR格式通过将矩阵的非零元素存储在一个一维数组中,并用两个一维数组存储行指针和列指针,来有效地压缩稀疏矩阵的存储空间。

CSR格式的一维数组包含三个部分:数据、列索引和行指针。假设稀疏矩阵的大小为m × n,其中非零元素个数为nnz。分别介绍这三个数组的含义:

  • 数据数组(values array):存储非零元素的值,大小为nnz。
  • 列索引数组(column indices array):存储每个非零元素所在的列数,大小为nnz。
  • 行指针数组(row pointer array):存储每一行的第一个非零元素在数据数组中的下标,大小为m+1。

例:

1 稀疏矩阵乘法

所用示例如下,矩阵A、B分别为

\[A = {\left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} 1&{ - 1} \end{array}}&{ - 3}&0&0\{\begin{array}{*{20}{c}} { - 2}&5 \end{array}}&0&0&0\{\begin{array}{*{20}{c}} 0&0 \end{array}}&4&6&4\{\begin{array}{*{20}{c}} {\begin{array}{*{20}{l}} { - 4}\0\1 \end{array}}&{\begin{array}{*{20}{l}} 0\8\0 \end{array}} \end{array}}&{\begin{array}{*{20}{l}} 2\0\0 \end{array}}&{\begin{array}{*{20}{l}} 7\0\0 \end{array}}&{\begin{array}{*{20}{l}} 0\{ - 5}\0 \end{array}} \end{array}} \right]_{6 \times 5}}{\begin{array}{*{20}{c}} {}&{B = \left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} 1\{ - 2} \end{array}}&{\begin{array}{*{20}{c}} { - 1}\5 \end{array}}&{\begin{array}{*{20}{c}} { - 3}\0 \end{array}}&{\begin{array}{*{20}{c}} 0\0 \end{array}}\0&0&4&6\{ - 4}&0&2&7\0&8&0&0 \end{array}} \right]} \end{array}_{5 \times 4}} \]

(1)matlab计算结果

作为标准答案,验证后续调用的正确性。

  1. A=[1,-1,-3,0,0;
  2. -2,5,0,0,0;
  3. 0,0,4,6,4;
  4. -4,0,2,7,0;
  5. 0,8,0,0,-5;
  6. 1,0,0,0,0];
  7. B=[1,-1,-3,0;
  8. -2,5,0,0;
  9. 0,0,4,6;
  10. -4,0,2,7;
  11. 0,8,0,0];
  12. A*B

输出为:

(2)稀疏矩阵X稠密矩阵

用于计算稀疏矩阵(COO格式表示)的矩阵\(A\),与稠密矩阵\(B\)的乘积。

函数将使用MKL库中的稀疏矩阵乘法接口mkl_sparse_s_mm实现\(y = alpha * A * x + beta * y\),具体用法及参数详解如下:

  1. mkl_sparse_s_mm(operation, //表示矩阵乘法的操作类型,可以是普通/转置/共轭转置
  2. alpha, //乘法系数
  3. A, //稀疏矩阵A
  4. descr, //结构体,描述矩阵的属性,包括存储格式、存储顺序等
  5. SPARSE_LAYOUT_ROW_MAJOR,//矩阵存储顺序
  6. x, //X矩阵,稠密
  7. columns, // 矩阵x的列数
  8. ldx, //矩阵x的第一维
  9. beta, //加法后系数
  10. y, //y矩阵,即输出矩阵
  11. ldy //矩阵y的第一维
  12. );

此流程简述如下:

  1. 获取待稀疏表示矩阵\(A\)的COO格式(ia,ja,value),以及非零元素个数nnz;
  2. 根据三组数据创建COO格式稀疏矩阵,并通过MKL转换接口将其转为CSR格式;
  3. 执行稀疏矩阵csrA与稠密矩阵denseB的乘积,使用mkl_sparse_s_mm接口计算矩阵乘法,结果为稠密矩阵C;
  4. 将计算结果转为需要的尺寸(此例为二维数组)返回。

稀疏矩阵coo乘稠密矩阵接口

  1. /*
  2. 输入:
  3. ia 稀疏矩阵A的行索引,一维MKL整型
  4. ja 稀疏矩阵A的列索引,一维MKL整型
  5. a 稀疏矩阵A的数据值,一维浮点型
  6. nnz 非零元素个数
  7. denseB 稠密矩阵B数据,类型为float型的二维数组
  8. rowsA 稀疏矩阵A的行数
  9. colsA 稀疏矩阵A的列数
  10. colsC A、B两矩阵相乘结果C的列数
  11. flag 对稀疏矩阵A的操作,选项为0、1、2。0-A 1-AT(A矩阵的转置) 2-AH(A矩阵的共轭转置) 默认为0
  12. 输出:
  13. denseC 稠密矩阵C数据,为A与B相乘后的结果,类型为float型的二维数组
  14. */
  15. bool MKL_Sparse_CooXDense(int *ia, int *ja, float *a, int nnz, float **denseB, float **denseC, int rowsA, int colsA, int colsC, int flag);

函数代码

  1. bool MKL_Sparse_CooXDense(MKL_INT *ia, MKL_INT *ja, float *a, int nnz, float **denseB, float **denseC, int rowsA, int colsA, int colsC, int flag) {
  2. //生成csr格式稀疏矩阵
  3. sparse_matrix_t csrA, cooA;
  4. sparse_status_t status = mkl_sparse_s_create_coo(&cooA,
  5. SPARSE_INDEX_BASE_ZERO,
  6. rowsA, // number of rows
  7. colsA, // number of cols
  8. nnz, // number of nonzeros
  9. ia,
  10. ja,
  11. a);
  12. if (status != SPARSE_STATUS_SUCCESS) {
  13. printf("Error creating COO sparse matrix.\n");
  14. }
  15. mkl_sparse_convert_csr(cooA, SPARSE_OPERATION_NON_TRANSPOSE, &csrA);
  16. //调用mkl稀疏矩阵与稠密矩阵乘法 C=alpha*op(A)*B+beta*C
  17. double alpha = 1.0;
  18. double beta = 0.0;
  19. int M, N, K;
  20. int ncols, ldx, ldy;
  21. if (flag == 1 || flag == 2) { //转置或共轭转置,AT或AH
  22. M = colsA;
  23. N = rowsA;
  24. K = colsC;
  25. ncols = K;
  26. ldx = K;
  27. ldy = K;
  28. }
  29. else { //默认的情况下,A
  30. M = rowsA;
  31. N = colsA;
  32. K = colsC;
  33. ncols = N;
  34. ldx = K;
  35. ldy = K;
  36. }
  37. //将二维稠密矩阵B转为一维
  38. float *denseB1D = (float*)mkl_malloc(N*K * sizeof(float), 64);
  39. for (int i = 0; i < N; i++) {
  40. memcpy(denseB1D + i * K, denseB[i], K * sizeof(float));
  41. }
  42. struct matrix_descr descr = { SPARSE_MATRIX_TYPE_GENERAL, SPARSE_FILL_MODE_FULL, SPARSE_DIAG_NON_UNIT };
  43. float *denseC1D = (float*)mkl_malloc(M * K * sizeof(float), 64);
  44. sparse_operation_t operation = SPARSE_OPERATION_NON_TRANSPOSE; //默认 op(A)=A;
  45. if (flag == 0) { //稀疏矩阵 op(A)=A
  46. operation = SPARSE_OPERATION_NON_TRANSPOSE;
  47. }
  48. if (flag == 1) {//稀疏矩阵 op(A)=AT
  49. operation = SPARSE_OPERATION_TRANSPOSE;
  50. }
  51. else if (flag == 2)
  52. { //稀疏矩阵 op(A)=AH
  53. operation = SPARSE_OPERATION_CONJUGATE_TRANSPOSE;
  54. }
  55. else {//稀疏矩阵 op(A)=A
  56. operation = SPARSE_OPERATION_NON_TRANSPOSE;
  57. }
  58. mkl_sparse_s_mm(operation, alpha, csrA, descr, SPARSE_LAYOUT_ROW_MAJOR, denseB1D, ncols, ldx, beta, denseC1D, ldy);
  59. //将计算结果转为2维
  60. for (int i = 0; i < M; i++) {
  61. memcpy(denseC[i], denseC1D + i * K, K * sizeof(float));
  62. }
  63. //释放
  64. mkl_sparse_destroy(csrA);
  65. mkl_sparse_destroy(cooA);
  66. mkl_free(denseB1D);
  67. mkl_free(denseC1D);
  68. return true;
  69. }

执行main.cpp中的MKL_Sparse_CooXDense_Demo()后,

(3)稀疏矩阵X稀疏矩阵

两个稀疏矩阵的乘法,使用mkl_sparse_spmm接口实现\(C = op(A) * B\),该接口的用法相对简单,

  1. mkl_sparse_spmm(SPARSE_OPERATION_NON_TRANSPOSE,//是否对A矩阵进行操作
  2. csrA, //A矩阵,CSR格式
  3. csrB, //B矩阵,CSR格式
  4. &csrC //C矩阵,CSR格式
  5. );

在进行稀疏矩阵示例之前,先补充两个封装函数:MKL_Coo2Csr()Print_Sparse_Csr_Matrix()

  1. /*
  2. 函数功能:根据已知坐标、元素值,创建CSR稀疏矩阵
  3. 输入:
  4. float *a 稀疏矩阵的值 ----参照稀疏矩阵的coo格式
  5. MKL_INT *ia 稀疏矩阵的行指针
  6. MKL_INT *ja 稀疏矩阵的列索引
  7. int nnz 稀疏矩阵的数量
  8. int nrows稀疏矩阵的行数
  9. int ncols稀疏矩阵的列数
  10. 输出:
  11. sparse_matrix_t CSR格式稀疏矩阵
  12. */
  13. sparse_matrix_t MKL_Coo2Csr(int* ia, int *ja, float *a, int nrows, int ncols, int nnz) {
  14. //建立coo矩阵A 与 csrA矩阵
  15. sparse_matrix_t cooA, csrA;
  16. mkl_sparse_s_create_coo(&cooA,
  17. SPARSE_INDEX_BASE_ZERO,
  18. nrows, // number of rows
  19. ncols, // number of cols
  20. nnz, // number of nonzeros
  21. ia,
  22. ja,
  23. a);
  24. //coo转csr
  25. mkl_sparse_convert_csr(cooA, SPARSE_OPERATION_NON_TRANSPOSE, &csrA);
  26. //释放cooA矩阵
  27. mkl_sparse_destroy(cooA);
  28. //返回csrA矩阵
  29. return csrA;
  30. }
  31. /*
  32. 函数功能:打印CSR稀疏矩阵的前n行n列元素
  33. 输入:
  34. sparse_matrix_t CSR格式稀疏矩阵
  35. int m 前m行
  36. int n 前n列
  37. */
  38. void Print_Sparse_Csr_Matrix(sparse_matrix_t csrA,int m,int n) {
  39. sparse_index_base_t indexing;
  40. int nrows;
  41. int ncols;
  42. MKL_INT* csr_row_start;
  43. MKL_INT* csr_row_end;
  44. MKL_INT* csr_col_indx;
  45. float* csr_values;
  46. mkl_sparse_s_export_csr(csrA, &indexing, &nrows, &ncols, &csr_row_start, &csr_row_end, &csr_col_indx, &csr_values);
  47. float **A_dense = alloc2float(ncols, nrows);
  48. memset(A_dense[0], 0, nrows*ncols * sizeof(float));
  49. //将value转换为普通二维数组
  50. for (int i = 0; i < nrows; i++) {
  51. for (int j = csr_row_start[i]; j < csr_row_start[i + 1]; j++) {
  52. A_dense[i][csr_col_indx[j]] = csr_values[j];
  53. }
  54. }
  55. //输出稠密矩阵
  56. for (int i = 0; i < m; i++) {
  57. for (int j = 0; j < n; j++) {
  58. printf("%f ", A_dense[i][j]);
  59. }
  60. printf("\n");
  61. }
  62. free2float(A_dense);
  63. }

稀疏矩阵(csr)乘稀疏矩阵(csr)接口

  1. /*
  2. 输入:
  3. float *a 稀疏矩阵A的属性 ----参照稀疏矩阵的coo格式
  4. MKL_INT *ia
  5. MKL_INT *ja
  6. int nnzA
  7. int rowsA
  8. int colsA
  9. float *b 稀疏矩阵B的属性 ----参照稀疏矩阵的coo格式
  10. MKL_INT *ib
  11. MKL_INT *jb
  12. int nnzB
  13. int rowsB
  14. int colsB
  15. 输出:
  16. sparse_matrix_t 稀疏矩阵C(A*B)
  17. */
  18. sparse_matrix_t MKL_Sparse_CooXCoo(
  19. int* ia, int *ja, float *a, int rowsA, int colsA, int nnzA,
  20. int* ib, int *jb, float *b, int rowsB, int colsB, int nnzB);

函数代码

  1. sparse_matrix_t MKL_Sparse_CooXCoo(
  2. int* ia, int *ja, float *a, int rowsA, int colsA, int nnzA,
  3. int* ib, int *jb, float *b, int rowsB, int colsB, int nnzB) {
  4. //根据坐标创建csrA和csrB
  5. sparse_matrix_t csrA, csrB, csrC;
  6. csrA = MKL_Coo2Csr(ia, ja, a, rowsA, colsA, nnzA);
  7. csrB = MKL_Coo2Csr(ib, jb, b, rowsB, colsB, nnzB);
  8. //csrC创建
  9. mkl_sparse_d_create_csr(&csrC,
  10. SPARSE_INDEX_BASE_ZERO,
  11. rowsA, // number of rows
  12. colsB, // number of cols
  13. NULL, // number of nonzeros
  14. NULL,
  15. NULL,
  16. NULL);
  17. //MKL乘法
  18. mkl_sparse_spmm(SPARSE_OPERATION_NON_TRANSPOSE, csrA, csrB, &csrC);
  19. //释放矩阵A和B
  20. mkl_sparse_destroy(csrA);
  21. mkl_sparse_destroy(csrB);
  22. return csrC;
  23. }

执行main.cpp中的MKL_Sparse_CooXCoo_Demo()后,

计算结果与matlab结果一致。

2 稀疏矩阵求逆

(1)matlab计算结果

作为标准答案,验证后续调用的正确性。

  1. A = [1 2 4 0 0;
  2. 2 2 0 0 0;
  3. 4 0 3 0 0;
  4. 0 0 0 4 0;
  5. 0 0 0 0 5];
  6. A_inv = inv(A)

输出为:

(2)MKL计算

MKL的求逆计算相对复杂,以下将介绍MKL中的高性能稀疏线性方程组求解器PARDISO(Parallel Direct Sparse Solver Interface),PARDISO 实现了高效的并行算法和内存管理技术,可以处理大规模、高度稀疏的线性方程组求解,并具有高性能和可扩展性。

PARDISO 的求解过程包括以下几个步骤:

  1. 输入矩阵:提供线性方程组的稀疏矩阵\(A\),稀疏CSR格式。
  2. 分析矩阵:对矩阵进行预处理和分解,生成求解器所需的数据结构和信息,包括消元树、消元顺序、LU 分解等。
  3. 求解线性方程组:使用 LU 分解求解线性方程组,可以直接求解$ A X = B \(或者\) AX = λX $问题。
  4. 输出解向量:输出求解得到的解向量 \(X\)

在我们使用PARDISO接口求逆时,思路为将\(B\)设为单位阵,此时求解\(X\)即为矩阵\(A\)的逆\(A^{-1}\)

关于PARDISO接口参数的详细描述参考oneMKL PARDISO Parameters in Tabular Form (intel.com),以下简单介绍:

  1. PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &n, a, ia, ja, &perm, &nrhs, iparm, &msglvl, b, x, &error);
  1. pt:指向PARDISO内部数据结构的指针。数组长度为64,必须用零初始化,并且不再改动。

  2. maxfct:最大因子数,通常设置为1。

  3. mnum:与maxfct一起使用,用于区分不同的矩阵。

  4. mtype:矩阵类型。具体取值如下:

    • 1 - 实对称矩阵;
    • 2 - 实对称正定矩阵;
    • -2 - 实对称不定矩阵;
    • 3 - 复对称矩阵;
    • 11 - 实数、非对称矩阵;
  5. phase:指定PARDISO的阶段。具体取值如:

    • 11-分析阶段;
    • 12-分析、数值分解阶段;
    • 13-分析、数值分解、求解阶段;
    • 22-数值分解阶段;
    • 23-数值分解、求解阶段;
    • 33-求解、迭代阶段;
    • -1-释放所有矩阵内存;
  6. n:\(AX=B\)的方程个数,简记为矩阵\(A\)的行。

  7. a:稀疏矩阵\(A\)的非零元素(CSR格式中的values)。

  8. ia:CSR格式中的行索引。

  9. ja:CSR格式中的列索引。

  10. perm:保存大小为 n 的置换向量。

  11. nrhs:官方解释为:需要求解的右侧数(Number of right-hand sides that need to be solved for),一般为1。

  12. iparm:iparm是PARDISO中的一个长度为64的整数数组,用于控制PARDISO求解器的行为。iparm中每个参数的详细说明参见pardiso iparm Parameter (intel.com),以下仅列出一些常用且便于理解的参数:

    • iparm[0]:0-使用默认值,非0-使用自定义参数;
    • iparm[11]:对稀疏矩阵A进行操作后求解。0-求解\(AX=B\),1-求解\(A^HX=B\),2-求解\(A^TX=B\)
    • iparm[12]: 使用(非)对称加权匹配提高准确性,0-禁用,1-开启;
    • iparm[27]: 单精度/双精度,0-double,1-float;
    • iparm[34]: 以0或1作为初始索引,0-从1开始索引,1-从0开始索引;
  13. msglvl:Message level information,0-不生成输出,1-打印计算信息。

  14. b:\(B\)矩阵。

  15. x:\(X\)矩阵。

  16. error:错误代码。

稀疏矩阵求逆接口

  1. /*
  2. 输入:
  3. float *a 稀疏矩阵的值 ----参照稀疏矩阵的csr格式
  4. MKL_INT *ia 稀疏矩阵的行指针
  5. MKL_INT *ja 稀疏矩阵的列索引
  6. int nnz 稀疏矩阵的数量
  7. int n 稀疏矩阵的维度 n*n
  8. MKL_INT mtype 稀疏矩阵类型
  9. 输出:
  10. float **Ainv [稠密矩阵]稀疏矩阵的逆 n*n
  11. */
  12. bool MKL_Sparse_Inverse(float **Ainv, float *a, MKL_INT *ia, MKL_INT *ja, int nnz, int n,MKL_INT mtype);

函数代码

在代码中对求逆步骤进行解释

  1. bool MKL_Sparse_Inverse(float **Ainv, float *a, MKL_INT *ia, MKL_INT *ja, int nnz, int n, MKL_INT mtype) {
  2. /*STEP1 根据输入数组创建COO格式稀疏矩阵*/
  3. //由于CSR格式不易表示,所以采取的路线为通过坐标创建COO格式矩阵
  4. //再通过mkl接口将COO矩阵转为CSR矩阵
  5. sparse_matrix_t csrA, cooA;
  6. sparse_status_t status = mkl_sparse_s_create_coo(&cooA,
  7. SPARSE_INDEX_BASE_ZERO,
  8. n, // 稀疏矩阵的行、列
  9. n,
  10. nnz, // 非零元素个数
  11. ia,//行索引
  12. ja,//列索引
  13. a);//矩阵元素值
  14. if (status != SPARSE_STATUS_SUCCESS) {
  15. printf("Error creating COO sparse matrix.\n");
  16. }
  17. //coo转csr格式
  18. mkl_sparse_convert_csr(cooA, SPARSE_OPERATION_NON_TRANSPOSE, &csrA);
  19. /*STEP2 根据CSR格式稀疏矩阵得到其ia,ja,a三组数据*/
  20. sparse_index_base_t indexing;
  21. int nrows;
  22. int ncols;
  23. MKL_INT* ia_csr = (MKL_INT*)mkl_malloc(nnz * sizeof(MKL_INT), 64);//csr格式的行指针
  24. MKL_INT* csr_row_end = (MKL_INT*)mkl_malloc(nnz * sizeof(MKL_INT), 64);
  25. MKL_INT* ja_csr = (MKL_INT*)mkl_malloc(nnz * sizeof(MKL_INT), 64);//csr格式的列索引
  26. float* a_csr = (float*)mkl_malloc(nnz * sizeof(float), 64);//csr格式的矩阵值
  27. //利用mkl_sparse_s_export_csr接口实现
  28. mkl_sparse_s_export_csr(csrA, &indexing, &nrows, &ncols, &ia_csr, &csr_row_end, &ja_csr, &a_csr);
  29. /*Step3 设置稀疏矩阵参数*/
  30. //初始化B矩阵和X矩阵
  31. float *b = NULL; //保存单位矩阵用于求逆
  32. float *x = NULL; //解矩阵
  33. b = (float*)mkl_malloc(n*n * sizeof(float), 64);
  34. x = (float*)mkl_malloc(n*n * sizeof(float), 64);
  35. //将B矩阵初始化为单位阵
  36. for (int i = 0; i < n; i++) {
  37. for (int j = 0; j < n; j++) {
  38. if (i == j) {
  39. b[i*n + j] = 1;
  40. }
  41. else {
  42. b[i*n + j] = 0;
  43. }
  44. }
  45. }
  46. //初始化pt
  47. void *pt[64];
  48. for (int i = 0; i < 64; i++)
  49. {
  50. pt[i] = 0;
  51. }
  52. //初始化矩阵相关控制参数
  53. MKL_INT maxfct = 1;
  54. MKL_INT mnum = 1;
  55. MKL_INT phase = 1;
  56. MKL_INT perm = 0;
  57. MKL_INT nrhs = n;
  58. MKL_INT error = 0;
  59. MKL_INT msglvl = 0;
  60. //设置iparm参数
  61. MKL_INT iparm[64];
  62. //批量初始化
  63. for (int i = 0; i < 64; i++)
  64. {
  65. iparm[i] = 0;
  66. }
  67. iparm[0] = 1; //开启自定义
  68. iparm[12] = 1; //提高准确性
  69. iparm[27] = 1; //为float型
  70. iparm[34] = 1; //初始索引为0
  71. /*Step4 分析矩阵、数值分解、求解*/
  72. phase = 13; //phase设置为13,表示分析、数值分解、求解阶段
  73. PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &n, a_csr, ia_csr, ja_csr, &perm, &nrhs, iparm, &msglvl, b, x, &error);
  74. //将解矩阵copy给输出,即A的逆矩阵
  75. for (int i = 0; i < n; i++) {
  76. memcpy(Ainv[i], x + i * n, n * sizeof(float));
  77. }
  78. //释放矩阵内存
  79. phase = -1; //phase设置为-1
  80. PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &n, a_csr, ia_csr, ja_csr, &perm, &nrhs, iparm, &msglvl, b, x, &error);
  81. //释放内存
  82. mkl_sparse_destroy(cooA);
  83. mkl_sparse_destroy(csrA);
  84. mkl_free(x);
  85. mkl_free(b);
  86. return true;
  87. }

在执行main.cpp中的MKL_Sparse_Inverse_Demo()之后,输出如下,与matlab结果一致:

完整代码

Ⅰ MKL_Sparse_Methods.h

  1. #pragma once
  2. #include<stdio.h>
  3. #include<stdlib.h>
  4. #include "alloc.h"
  5. #include"mkl.h"
  6. #include "mkl_types.h"
  7. #include"mkl_lapacke.h"
  8. #include "mkl_spblas.h"
  9. bool MKL_Sparse_CooXDense(MKL_INT *ia, MKL_INT *ja, float *a, int nnz, float **denseB, float **denseC, int rowsA, int colsA, int colsC, int flag);
  10. sparse_matrix_t MKL_Sparse_CooXCoo(
  11. int* ia, int *ja, float *a, int rowsA, int colsA, int nnzA,
  12. int* ib, int *jb, float *b, int rowsB, int colsB, int nnzB);
  13. sparse_matrix_t MKL_Coo2Csr(int* ia, int *ja, float *a, int nrows, int ncols, int nnz);
  14. void Print_Sparse_Csr_Matrix(sparse_matrix_t csrA,int m,int n);
  15. bool MKL_Sparse_Inverse(float **Ainv, float *a, MKL_INT *ia, MKL_INT *ja, int nnz, int n, MKL_INT mtype);

Ⅱ MKL_Sparse_Methods.cpp

  1. #include "MKL_Sparse_Methods.h"
  2. bool MKL_Sparse_CooXDense(MKL_INT *ia, MKL_INT *ja, float *a, int nnz, float **denseB, float **denseC, int rowsA, int colsA, int colsC, int flag) {
  3. //生成csr格式稀疏矩阵
  4. sparse_matrix_t csrA, cooA;
  5. sparse_status_t status = mkl_sparse_s_create_coo(&cooA,
  6. SPARSE_INDEX_BASE_ZERO,
  7. rowsA, // number of rows
  8. colsA, // number of cols
  9. nnz, // number of nonzeros
  10. ia,
  11. ja,
  12. a);
  13. if (status != SPARSE_STATUS_SUCCESS) {
  14. printf("Error creating COO sparse matrix.\n");
  15. }
  16. mkl_sparse_convert_csr(cooA, SPARSE_OPERATION_NON_TRANSPOSE, &csrA);
  17. //调用mkl稀疏矩阵与稠密矩阵乘法 C=alpha*op(A)*B+beta*C
  18. double alpha = 1.0;
  19. double beta = 0.0;
  20. int M, N, K;
  21. int ncols, ldx, ldy;
  22. if (flag == 1 || flag == 2) {
  23. M = colsA;
  24. N = rowsA;
  25. K = colsC;
  26. ncols = K;
  27. ldx = K;
  28. ldy = K;
  29. }
  30. else { //默认的情况下
  31. M = rowsA;
  32. N = colsA;
  33. K = colsC;
  34. ncols = N;
  35. ldx = K;
  36. ldy = K;
  37. }
  38. //将二维稠密矩阵B转为一维
  39. float *denseB1D = (float*)mkl_malloc(N*K * sizeof(float), 64);
  40. for (int i = 0; i < N; i++) {
  41. memcpy(denseB1D + i * K, denseB[i], K * sizeof(float));
  42. }
  43. struct matrix_descr descr = { SPARSE_MATRIX_TYPE_GENERAL, SPARSE_FILL_MODE_FULL, SPARSE_DIAG_NON_UNIT };
  44. float *denseC1D = (float*)mkl_malloc(M * K * sizeof(float), 64);
  45. sparse_operation_t operation = SPARSE_OPERATION_NON_TRANSPOSE; //默认 op(A)=A;
  46. if (flag == 0) { //稀疏矩阵 op(A)=A
  47. operation = SPARSE_OPERATION_NON_TRANSPOSE;
  48. }
  49. if (flag == 1) {//稀疏矩阵 op(A)=AT
  50. operation = SPARSE_OPERATION_TRANSPOSE;
  51. }
  52. else if (flag == 2)
  53. { //稀疏矩阵 op(A)=AH
  54. operation = SPARSE_OPERATION_CONJUGATE_TRANSPOSE;
  55. }
  56. else {//稀疏矩阵 op(A)=A
  57. operation = SPARSE_OPERATION_NON_TRANSPOSE;
  58. }
  59. mkl_sparse_s_mm(operation, alpha, csrA, descr, SPARSE_LAYOUT_ROW_MAJOR, denseB1D, ncols, ldx, beta, denseC1D, ldy);
  60. //将计算结果转为2维
  61. for (int i = 0; i < M; i++) {
  62. memcpy(denseC[i], denseC1D + i * K, K * sizeof(float));
  63. }
  64. //释放
  65. mkl_sparse_destroy(csrA);
  66. mkl_sparse_destroy(cooA);
  67. mkl_free(denseB1D);
  68. mkl_free(denseC1D);
  69. return true;
  70. }
  71. //根据已知坐标、元素值,创建CSR稀疏矩阵
  72. sparse_matrix_t MKL_Coo2Csr(int* ia, int *ja, float *a, int nrows, int ncols, int nnz) {
  73. //建立coo矩阵A 与 csrA矩阵
  74. sparse_matrix_t cooA, csrA;
  75. mkl_sparse_s_create_coo(&cooA,
  76. SPARSE_INDEX_BASE_ZERO,
  77. nrows, // number of rows
  78. ncols, // number of cols
  79. nnz, // number of nonzeros
  80. ia,
  81. ja,
  82. a);
  83. //coo转csr
  84. mkl_sparse_convert_csr(cooA, SPARSE_OPERATION_NON_TRANSPOSE, &csrA);
  85. //释放cooA矩阵
  86. mkl_sparse_destroy(cooA);
  87. //返回csrA矩阵
  88. return csrA;
  89. }
  90. void Print_Sparse_Csr_Matrix(sparse_matrix_t csrA,int m,int n) {
  91. sparse_index_base_t indexing;
  92. int nrows;
  93. int ncols;
  94. MKL_INT* csr_row_start;
  95. MKL_INT* csr_row_end;
  96. MKL_INT* csr_col_indx;
  97. float* csr_values;
  98. mkl_sparse_s_export_csr(csrA, &indexing, &nrows, &ncols, &csr_row_start, &csr_row_end, &csr_col_indx, &csr_values);
  99. float **A_dense = alloc2float(ncols, nrows);
  100. memset(A_dense[0], 0, nrows*ncols * sizeof(float));
  101. //将value转换为普通二维数组
  102. for (int i = 0; i < nrows; i++) {
  103. for (int j = csr_row_start[i]; j < csr_row_start[i + 1]; j++) {
  104. A_dense[i][csr_col_indx[j]] = csr_values[j];
  105. }
  106. }
  107. //输出稠密矩阵
  108. for (int i = 0; i < m; i++) {
  109. for (int j = 0; j < n; j++) {
  110. printf("%f ", A_dense[i][j]);
  111. }
  112. printf("\n");
  113. }
  114. free2float(A_dense);
  115. }
  116. //Coo格式×Coo格式矩阵乘法
  117. sparse_matrix_t MKL_Sparse_CooXCoo(
  118. int* ia, int *ja, float *a, int rowsA, int colsA, int nnzA,
  119. int* ib, int *jb, float *b, int rowsB, int colsB, int nnzB) {
  120. //根据坐标创建csrA和csrB
  121. sparse_matrix_t csrA, csrB, csrC;
  122. csrA = MKL_Coo2Csr(ia, ja, a, rowsA, colsA, nnzA);
  123. csrB = MKL_Coo2Csr(ib, jb, b, rowsB, colsB, nnzB);
  124. //csrC创建
  125. mkl_sparse_d_create_csr(&csrC,
  126. SPARSE_INDEX_BASE_ZERO,
  127. rowsA, // number of rows
  128. colsB, // number of cols
  129. NULL, // number of nonzeros
  130. NULL,
  131. NULL,
  132. NULL);
  133. mkl_sparse_spmm(SPARSE_OPERATION_NON_TRANSPOSE, csrA, csrB, &csrC);
  134. //释放矩阵A和B
  135. mkl_sparse_destroy(csrA);
  136. mkl_sparse_destroy(csrB);
  137. return csrC;
  138. }
  139. //稀疏矩阵求逆
  140. bool MKL_Sparse_Inverse(float **Ainv, float *a, MKL_INT *ia, MKL_INT *ja, int nnz, int n, MKL_INT mtype) {
  141. /*STEP1 根据输入数组创建COO格式稀疏矩阵*/
  142. //由于CSR格式不易表示,所以采取的路线为通过坐标创建COO格式矩阵
  143. //再通过mkl接口将COO矩阵转为CSR矩阵
  144. sparse_matrix_t csrA, cooA;
  145. sparse_status_t status = mkl_sparse_s_create_coo(&cooA,
  146. SPARSE_INDEX_BASE_ZERO,
  147. n, // 稀疏矩阵的行、列
  148. n,
  149. nnz, // 非零元素个数
  150. ia,//行索引
  151. ja,//列索引
  152. a);//矩阵元素值
  153. if (status != SPARSE_STATUS_SUCCESS) {
  154. printf("Error creating COO sparse matrix.\n");
  155. }
  156. //coo转csr格式
  157. mkl_sparse_convert_csr(cooA, SPARSE_OPERATION_NON_TRANSPOSE, &csrA);
  158. /*STEP2 根据CSR格式稀疏矩阵得到其ia,ja,a三组数据*/
  159. sparse_index_base_t indexing;
  160. int nrows;
  161. int ncols;
  162. MKL_INT* ia_csr = (MKL_INT*)mkl_malloc(nnz * sizeof(MKL_INT), 64);//csr格式的行指针
  163. MKL_INT* csr_row_end = (MKL_INT*)mkl_malloc(nnz * sizeof(MKL_INT), 64);
  164. MKL_INT* ja_csr = (MKL_INT*)mkl_malloc(nnz * sizeof(MKL_INT), 64);//csr格式的列索引
  165. float* a_csr = (float*)mkl_malloc(nnz * sizeof(float), 64);//csr格式的矩阵值
  166. //利用mkl_sparse_s_export_csr接口实现
  167. mkl_sparse_s_export_csr(csrA, &indexing, &nrows, &ncols, &ia_csr, &csr_row_end, &ja_csr, &a_csr);
  168. /*Step3 设置稀疏矩阵参数*/
  169. //初始化B矩阵和X矩阵
  170. float *b = NULL; //保存单位矩阵用于求逆
  171. float *x = NULL; //解矩阵
  172. b = (float*)mkl_malloc(n*n * sizeof(float), 64);
  173. x = (float*)mkl_malloc(n*n * sizeof(float), 64);
  174. //将B矩阵初始化为单位阵
  175. for (int i = 0; i < n; i++) {
  176. for (int j = 0; j < n; j++) {
  177. if (i == j) {
  178. b[i*n + j] = 1;
  179. }
  180. else {
  181. b[i*n + j] = 0;
  182. }
  183. }
  184. }
  185. //初始化pt
  186. void *pt[64];
  187. for (int i = 0; i < 64; i++)
  188. {
  189. pt[i] = 0;
  190. }
  191. //初始化矩阵相关控制参数
  192. MKL_INT maxfct = 1;
  193. MKL_INT mnum = 1;
  194. MKL_INT phase = 1;
  195. MKL_INT perm = 0;
  196. MKL_INT nrhs = n;
  197. MKL_INT error = 0;
  198. MKL_INT msglvl = 0;
  199. //设置iparm参数
  200. MKL_INT iparm[64];
  201. //批量初始化
  202. for (int i = 0; i < 64; i++)
  203. {
  204. iparm[i] = 0;
  205. }
  206. iparm[0] = 1; //开启自定义
  207. iparm[12] = 1; //提高准确性
  208. iparm[27] = 1; //为float型
  209. iparm[34] = 1; //初始索引为0
  210. /*Step4 分析矩阵、数值分解、求解*/
  211. phase = 13; //phase设置为13,表示分析、数值分解、求解阶段
  212. PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &n, a_csr, ia_csr, ja_csr, &perm, &nrhs, iparm, &msglvl, b, x, &error);
  213. //将解矩阵copy给输出,即A的逆矩阵
  214. for (int i = 0; i < n; i++) {
  215. memcpy(Ainv[i], x + i * n, n * sizeof(float));
  216. }
  217. //释放矩阵内存
  218. phase = -1; //phase设置为-1
  219. PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &n, a_csr, ia_csr, ja_csr, &perm, &nrhs, iparm, &msglvl, b, x, &error);
  220. //释放内存
  221. mkl_sparse_destroy(cooA);
  222. mkl_sparse_destroy(csrA);
  223. mkl_free(x);
  224. mkl_free(b);
  225. return true;
  226. }

Ⅲ main.cpp

  1. #include "MKL_Sparse_Methods.h"
  2. #include "alloc.h"
  3. #define M 6
  4. #define N 5
  5. #define K 4
  6. void MKL_Sparse_CooXDense_Demo();
  7. void MKL_Sparse_CooXCoo_Demo();
  8. void MKL_Sparse_Inverse_Demo();
  9. int main() {
  10. MKL_Sparse_CooXDense_Demo();//稀疏乘稠密
  11. MKL_Sparse_CooXCoo_Demo(); //稀疏×稀疏
  12. MKL_Sparse_Inverse_Demo();//稀疏矩阵求逆
  13. return 0;
  14. }
  15. //稀疏乘稠密
  16. void MKL_Sparse_CooXDense_Demo() {
  17. int flag = 0; /* flag=0时表示A(COO)*B(Dense)
  18. flag=1时表示AT(COO)*B(Dense)*/
  19. int rowsA, colsA;
  20. if (flag == 0) {
  21. rowsA = M, colsA = N;
  22. }
  23. else if (flag == 1) {
  24. rowsA = N, colsA = M;
  25. }
  26. int rowsB = N, colsB = K;
  27. int rowsC = M, colsC = K;
  28. float Atemp[M][N] = {
  29. {1,-1,-3,0,0},
  30. {-2,5,0,0,0},
  31. {0,0,4,6,4},
  32. {-4,0,2,7,0},
  33. {0,8,0,0,-5},
  34. {1,0,0,0,0},
  35. };
  36. float ATtemp[N][M] = {
  37. {1,-2,0,-4,0,1},
  38. {-1,5,0,0,8,0},
  39. {-3,0,4,2,0,0},
  40. {0,0,6,7,0,0},
  41. {0,0,4,0,-5,0},
  42. };
  43. float Btemp[N][K] = {
  44. {1,-1,-3,0},
  45. {-2,5,0,0},
  46. {0,0,4,6},
  47. {-4,0,2,7},
  48. {0,8,0,0}
  49. };
  50. //将一般二维数组转换为alloc表示
  51. float **matrixA = alloc2float(colsA, rowsA);
  52. memset(matrixA[0], 0, rowsA*colsA * sizeof(float));
  53. float **matrixB = alloc2float(colsB, rowsB);
  54. memset(matrixB[0], 0, rowsB*colsB * sizeof(float));
  55. float **matrixC = alloc2float(colsC, rowsC);
  56. memset(matrixC[0], 0, rowsC*colsC * sizeof(float));
  57. //复制二维数组到二级指针
  58. if (flag == 0) {
  59. memcpy(matrixA[0], Atemp, rowsA*colsA * sizeof(float));
  60. }
  61. else if (flag == 1) {
  62. memcpy(matrixA[0], ATtemp, rowsA*colsA * sizeof(float));
  63. }
  64. memcpy(matrixB[0], Btemp, rowsB*colsB * sizeof(float));
  65. //统计二维数组的非零元素个数
  66. int nnz = 0;
  67. for (int i = 0; i < rowsA; i++) {
  68. for (int j = 0; j < colsA; j++) {
  69. //统计nnz
  70. if (matrixA[i][j] != 0) {
  71. nnz++;
  72. }
  73. }
  74. }
  75. //获取稠密矩阵的coo形式
  76. MKL_INT *ia = (MKL_INT*)mkl_malloc(nnz * sizeof(MKL_INT), 64);
  77. MKL_INT *ja = (MKL_INT*)mkl_malloc(nnz * sizeof(MKL_INT), 64);
  78. float *a = (float*)mkl_malloc(nnz * sizeof(float), 64);
  79. //获取coo数据
  80. int k = 0;
  81. for (int i = 0; i < rowsA; i++) {
  82. for (int j = 0; j < colsA; j++) {
  83. if (matrixA[i][j] != 0.0) {
  84. a[k] = matrixA[i][j];
  85. ia[k] = i;
  86. ja[k] = j;
  87. k++;
  88. }
  89. }
  90. }
  91. //执行矩阵乘法
  92. MKL_Sparse_CooXDense(ia, ja, a, nnz, matrixB, matrixC, rowsA, colsA, colsC, flag);
  93. /* 输出结果 */
  94. printf("*************** MKL Sparse X Dense ***************\n ");
  95. printf("=============== flag=%d ================\n", flag);
  96. for (int i = 0; i < rowsC; i++) {
  97. for (int j = 0; j < colsC; j++) {
  98. printf("%f ", matrixC[i][j]);
  99. }
  100. printf("\n");
  101. }
  102. free2float(matrixA);
  103. free2float(matrixB);
  104. free2float(matrixC);
  105. mkl_free(ia);
  106. mkl_free(ja);
  107. mkl_free(a);
  108. }
  109. //稀疏×稀疏
  110. void MKL_Sparse_CooXCoo_Demo() {
  111. int rowsA = M, colsA = N;
  112. int rowsB = N, colsB = K;
  113. int rowsC = M, colsC = K;
  114. float Atemp[M][N] = {
  115. {1,-1,-3,0,0},
  116. {-2,5,0,0,0},
  117. {0,0,4,6,4},
  118. {-4,0,2,7,0},
  119. {0,8,0,0,-5},
  120. {1,0,0,0,0},
  121. };
  122. float Btemp[N][K] = {
  123. {1,-1,-3,0},
  124. {-2,5,0,0},
  125. {0,0,4,6},
  126. {-4,0,2,7},
  127. {0,8,0,0}
  128. };
  129. //将一般二维数组转换为alloc表示
  130. float **matrixA = alloc2float(colsA, rowsA);
  131. memset(matrixA[0], 0, rowsA*colsA * sizeof(float));
  132. float **matrixB = alloc2float(colsB, rowsB);
  133. memset(matrixB[0], 0, rowsB*colsB * sizeof(float));
  134. //复制二维数组到二级指针
  135. memcpy(matrixA[0], Atemp, rowsA*colsA * sizeof(float));
  136. memcpy(matrixB[0], Btemp, rowsB*colsB * sizeof(float));
  137. //***********A矩阵稀疏表示
  138. //统计二维数组的非零元素个数
  139. int nnzA = 0; //A矩阵
  140. for (int i = 0; i < rowsA; i++) {
  141. for (int j = 0; j < colsA; j++) {
  142. //统计nnz
  143. if (matrixA[i][j] != 0) {
  144. nnzA++;
  145. }
  146. }
  147. }
  148. //获取稠密矩阵的coo形式
  149. MKL_INT *ia = (MKL_INT*)mkl_malloc(nnzA * sizeof(MKL_INT), 64);
  150. MKL_INT *ja = (MKL_INT*)mkl_malloc(nnzA * sizeof(MKL_INT), 64);
  151. float *a = (float*)mkl_malloc(nnzA * sizeof(float), 64);
  152. //获取coo数据
  153. int k = 0;
  154. for (int i = 0; i < rowsA; i++) {
  155. for (int j = 0; j < colsA; j++) {
  156. if (matrixA[i][j] != 0.0) {
  157. a[k] = matrixA[i][j];
  158. ia[k] = i;
  159. ja[k] = j;
  160. k++;
  161. }
  162. }
  163. }
  164. //***********B矩阵稀疏表示
  165. int nnzB = 0; //B矩阵
  166. for (int i = 0; i < rowsB; i++) {
  167. for (int j = 0; j < colsB; j++) {
  168. //统计nnz
  169. if (matrixB[i][j] != 0) {
  170. nnzB++;
  171. }
  172. }
  173. }
  174. //获取稠密矩阵的coo形式
  175. MKL_INT *ib = (MKL_INT*)mkl_malloc(nnzB * sizeof(MKL_INT), 64);
  176. MKL_INT *jb = (MKL_INT*)mkl_malloc(nnzB * sizeof(MKL_INT), 64);
  177. float *b = (float*)mkl_malloc(nnzB * sizeof(float), 64);
  178. //获取coo数据
  179. k = 0;
  180. for (int i = 0; i < rowsB; i++) {
  181. for (int j = 0; j < colsB; j++) {
  182. if (matrixB[i][j] != 0.0) {
  183. b[k] = matrixB[i][j];
  184. ib[k] = i;
  185. jb[k] = j;
  186. k++;
  187. }
  188. }
  189. }
  190. sparse_matrix_t csrC = MKL_Sparse_CooXCoo(ia, ja, a, rowsA, colsA, nnzA, ib, jb, b, rowsB, colsB, nnzB);
  191. printf("*************** MKL Sparse X Sparse ***************\n ");
  192. Print_Sparse_Csr_Matrix(csrC, rowsC, colsC); //打印矩阵C结果
  193. mkl_sparse_destroy(csrC);
  194. mkl_free(ia);
  195. mkl_free(ja);
  196. mkl_free(a);
  197. mkl_free(ib);
  198. mkl_free(jb);
  199. mkl_free(b);
  200. free2float(matrixA);
  201. free2float(matrixB);
  202. }
  203. //稀疏矩阵求逆
  204. void MKL_Sparse_Inverse_Demo() {
  205. int rowsA = N, colsA = N;
  206. float Atemp[N][N] = {
  207. {1,2,4,0,0},
  208. {2,2,0,0,0},
  209. {4,0,3,0,0},
  210. {0,0,0,4,0},
  211. {0,0,0,0,5},
  212. };
  213. //将一般二维数组转换为alloc表示
  214. float **matrixA = alloc2float(colsA, rowsA);
  215. memset(matrixA[0], 0, rowsA*colsA * sizeof(float));
  216. float **matrixA_inv = alloc2float(colsA, rowsA);
  217. memset(matrixA[0], 0, rowsA*colsA * sizeof(float));
  218. //复制二维数组到二级指针
  219. memcpy(matrixA[0], Atemp, rowsA*colsA * sizeof(float));
  220. //统计二维数组的非零元素个数
  221. int nnz = 0;
  222. for (int i = 0; i < rowsA; i++) {
  223. for (int j = 0; j < colsA; j++) {
  224. //统计nnz
  225. if (matrixA[i][j] != 0) {
  226. nnz++;
  227. }
  228. }
  229. }
  230. //获取稠密矩阵的coo形式
  231. MKL_INT *ia = (MKL_INT*)mkl_malloc(nnz * sizeof(MKL_INT), 64);
  232. MKL_INT *ja = (MKL_INT*)mkl_malloc(nnz * sizeof(MKL_INT), 64);
  233. float *a = (float*)mkl_malloc(nnz * sizeof(float), 64);
  234. //获取coo数据
  235. int k = 0;
  236. for (int i = 0; i < rowsA; i++) {
  237. for (int j = 0; j < colsA; j++) {
  238. if (matrixA[i][j] != 0.0) {
  239. a[k] = matrixA[i][j];
  240. ia[k] = i;
  241. ja[k] = j;
  242. k++;
  243. }
  244. }
  245. }
  246. MKL_INT mtype = 11; //设置矩阵类型为一般实矩阵
  247. MKL_Sparse_Inverse(matrixA, a, ia, ja, nnz, rowsA, mtype);
  248. printf("*************** MKL Sparse Matrix Inverse ***************\n ");
  249. for (int i = 0; i < N; i++) {
  250. for (int j = 0; j < N; j++) {
  251. printf("%f ", matrixA[i][j]);
  252. }
  253. printf("\n");
  254. }
  255. free2float(matrixA);
  256. free2float(matrixA_inv);
  257. mkl_free(ia);
  258. mkl_free(ja);
  259. mkl_free(a);
  260. }

原文链接:https://www.cnblogs.com/GeophysicsWorker/p/17347266.html

 友情链接:直通硅谷  点职佳  北美留学生论坛

本站QQ群:前端 618073944 | Java 606181507 | Python 626812652 | C/C++ 612253063 | 微信 634508462 | 苹果 692586424 | C#/.net 182808419 | PHP 305140648 | 运维 608723728

W3xue 的所有内容仅供测试,对任何法律问题及风险不承担任何责任。通过使用本站内容随之而来的风险与本站无关。
关于我们  |  意见建议  |  捐助我们  |  报错有奖  |  广告合作、友情链接(目前9元/月)请联系QQ:27243702 沸活量
皖ICP备17017327号-2 皖公网安备34020702000426号