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

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

(资料图片)

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

MKL中主要用到的稀疏矩阵格式有COO，CSR及CSC（与CSR类似）三种，以下将简要介绍COO格式与CSR格式：

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

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

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

例：

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

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

例：

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

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

A=[1,-1,-3,0,0;  -2,5,0,0,0;   0,0,4,6,4;  -4,0,2,7,0;   0,8,0,0,-5;   1,0,0,0,0];B=[1,-1,-3,0;  -2,5,0,0;   0,0,4,6;  -4,0,2,7;   0,8,0,0];A*B

输出为：

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

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

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

此流程简述如下：

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

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

函数代码

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) {//生成csr格式稀疏矩阵sparse_matrix_t csrA, cooA;sparse_status_t status = mkl_sparse_s_create_coo(&cooA,SPARSE_INDEX_BASE_ZERO,rowsA,    // number of rowscolsA,    // number of colsnnz,  // number of nonzerosia,ja,a);if (status != SPARSE_STATUS_SUCCESS) {printf("Error creating COO sparse matrix.\n");}mkl_sparse_convert_csr(cooA, SPARSE_OPERATION_NON_TRANSPOSE, &csrA);//调用mkl稀疏矩阵与稠密矩阵乘法  C=alpha*op(A)*B+beta*Cdouble alpha = 1.0;double beta = 0.0;int M, N, K;int ncols, ldx, ldy;if (flag == 1 || flag == 2) {//转置或共轭转置，AT或AHM = colsA;N = rowsA;K = colsC;ncols = K;ldx = K;ldy = K;}else {//默认的情况下，AM = rowsA;N = colsA;K = colsC;ncols = N;ldx = K;ldy = K;}//将二维稠密矩阵B转为一维float *denseB1D = (float*)mkl_malloc(N*K * sizeof(float), 64);for (int i = 0; i < N; i++) {memcpy(denseB1D + i * K, denseB[i], K * sizeof(float));}struct matrix_descr descr = { SPARSE_MATRIX_TYPE_GENERAL, SPARSE_FILL_MODE_FULL, SPARSE_DIAG_NON_UNIT };float *denseC1D = (float*)mkl_malloc(M * K * sizeof(float), 64);sparse_operation_t operation = SPARSE_OPERATION_NON_TRANSPOSE;  //默认 op(A)=A;if (flag == 0) {   //稀疏矩阵 op(A)=A operation = SPARSE_OPERATION_NON_TRANSPOSE;}if (flag == 1) {//稀疏矩阵 op(A)=AToperation = SPARSE_OPERATION_TRANSPOSE;}else if (flag == 2){    //稀疏矩阵 op(A)=AHoperation = SPARSE_OPERATION_CONJUGATE_TRANSPOSE;}else {//稀疏矩阵 op(A)=A operation = SPARSE_OPERATION_NON_TRANSPOSE;}mkl_sparse_s_mm(operation, alpha, csrA, descr, SPARSE_LAYOUT_ROW_MAJOR, denseB1D, ncols, ldx, beta, denseC1D, ldy);//将计算结果转为2维for (int i = 0; i < M; i++) {memcpy(denseC[i], denseC1D + i * K, K * sizeof(float));}//释放mkl_sparse_destroy(csrA);mkl_sparse_destroy(cooA);mkl_free(denseB1D);mkl_free(denseC1D);return true;}

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

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

mkl_sparse_spmm(SPARSE_OPERATION_NON_TRANSPOSE,//是否对A矩阵进行操作                csrA, //A矩阵，CSR格式                csrB,//B矩阵，CSR格式                &csrC//C矩阵，CSR格式               );

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

/*函数功能：根据已知坐标、元素值，创建CSR稀疏矩阵输入：float *a    稀疏矩阵的值  ----参照稀疏矩阵的coo格式MKL_INT *ia  稀疏矩阵的行指针MKL_INT *ja  稀疏矩阵的列索引int     nnz  稀疏矩阵的数量int     nrows稀疏矩阵的行数int     ncols稀疏矩阵的列数输出：sparse_matrix_t   CSR格式稀疏矩阵*/sparse_matrix_t MKL_Coo2Csr(int* ia, int *ja, float *a, int nrows, int ncols, int nnz) {//建立coo矩阵A 与 csrA矩阵sparse_matrix_t cooA, csrA;mkl_sparse_s_create_coo(&cooA,SPARSE_INDEX_BASE_ZERO,nrows,    // number of rowsncols,    // number of colsnnz,  // number of nonzerosia,ja,a);//coo转csrmkl_sparse_convert_csr(cooA, SPARSE_OPERATION_NON_TRANSPOSE, &csrA);//释放cooA矩阵mkl_sparse_destroy(cooA);//返回csrA矩阵return csrA;}/*函数功能：打印CSR稀疏矩阵的前n行n列元素输入：sparse_matrix_t   CSR格式稀疏矩阵int m   前m行int n前n列*/void Print_Sparse_Csr_Matrix(sparse_matrix_t csrA,int m,int n) {sparse_index_base_t indexing;int nrows;int ncols;MKL_INT* csr_row_start;MKL_INT* csr_row_end;MKL_INT* csr_col_indx;float* csr_values;mkl_sparse_s_export_csr(csrA, &indexing, &nrows, &ncols, &csr_row_start, &csr_row_end, &csr_col_indx, &csr_values);float **A_dense = alloc2float(ncols, nrows);memset(A_dense[0], 0, nrows*ncols * sizeof(float));//将value转换为普通二维数组for (int i = 0; i < nrows; i++) {for (int j = csr_row_start[i]; j < csr_row_start[i + 1]; j++) {A_dense[i][csr_col_indx[j]] = csr_values[j];}}//输出稠密矩阵for (int i = 0; i < m; i++) {for (int j = 0; j < n; j++) {printf("%f ", A_dense[i][j]);}printf("\n");}free2float(A_dense);}

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

/*输入：float *a    稀疏矩阵A的属性  ----参照稀疏矩阵的coo格式MKL_INT *ia  MKL_INT *ja  int     nnzA int    rowsA int    colsA float *b    稀疏矩阵B的属性  ----参照稀疏矩阵的coo格式MKL_INT *ib  MKL_INT *jb  int     nnzB int    rowsB int    colsB 输出：sparse_matrix_t  稀疏矩阵C(A*B)*/sparse_matrix_t MKL_Sparse_CooXCoo(int* ia, int *ja, float *a, int rowsA, int colsA, int nnzA,int* ib, int *jb, float *b, int rowsB, int colsB, int nnzB);

函数代码

sparse_matrix_t MKL_Sparse_CooXCoo(int* ia, int *ja, float *a, int rowsA, int colsA, int nnzA,int* ib, int *jb, float *b, int rowsB, int colsB, int nnzB) {//根据坐标创建csrA和csrBsparse_matrix_t csrA, csrB, csrC;csrA = MKL_Coo2Csr(ia, ja, a, rowsA, colsA, nnzA);csrB = MKL_Coo2Csr(ib, jb, b, rowsB, colsB, nnzB);//csrC创建mkl_sparse_d_create_csr(&csrC,SPARSE_INDEX_BASE_ZERO,rowsA,    // number of rowscolsB,    // number of colsNULL,  // number of nonzerosNULL,NULL,NULL);        //MKL乘法mkl_sparse_spmm(SPARSE_OPERATION_NON_TRANSPOSE, csrA, csrB, &csrC);//释放矩阵A和Bmkl_sparse_destroy(csrA);mkl_sparse_destroy(csrB);return csrC;}

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

计算结果与matlab结果一致。

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

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

输出为：

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

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

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

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

PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &n, a, ia, ja, &perm, &nrhs, iparm, &msglvl, b, x, &error);

pt：指向PARDISO内部数据结构的指针。数组长度为64，必须用零初始化，并且不再改动。

maxfct：最大因子数，通常设置为1。

mnum：与maxfct一起使用，用于区分不同的矩阵。

mtype：矩阵类型。具体取值如下：

phase：指定PARDISO的阶段。具体取值如：

n：\(AX=B\)的方程个数，简记为矩阵\(A\)的行。

a：稀疏矩阵\(A\)的非零元素（CSR格式中的values）。

ia：CSR格式中的行索引。

ja：CSR格式中的列索引。

perm：保存大小为 n 的置换向量。

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

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

msglvl：Message level information，0-不生成输出，1-打印计算信息。

b：\(B\)矩阵。

x：\(X\)矩阵。

error：错误代码。

稀疏矩阵求逆接口

/*输入：float *a    稀疏矩阵的值  ----参照稀疏矩阵的csr格式MKL_INT *ia  稀疏矩阵的行指针MKL_INT *ja  稀疏矩阵的列索引int     nnz  稀疏矩阵的数量int     n    稀疏矩阵的维度 n*nMKL_INT mtype 稀疏矩阵类型输出：float **Ainv   [稠密矩阵]稀疏矩阵的逆 n*n   */bool MKL_Sparse_Inverse(float **Ainv, float *a, MKL_INT *ia, MKL_INT *ja, int nnz, int n,MKL_INT mtype);

函数代码

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

bool MKL_Sparse_Inverse(float **Ainv, float *a, MKL_INT *ia, MKL_INT *ja, int nnz, int n, MKL_INT mtype) {/*STEP1 根据输入数组创建COO格式稀疏矩阵*/    //由于CSR格式不易表示，所以采取的路线为通过坐标创建COO格式矩阵//再通过mkl接口将COO矩阵转为CSR矩阵sparse_matrix_t csrA, cooA;sparse_status_t status = mkl_sparse_s_create_coo(&cooA,SPARSE_INDEX_BASE_ZERO,n,    // 稀疏矩阵的行、列n,nnz,  // 非零元素个数ia,//行索引                                      ja,//列索引a);//矩阵元素值if (status != SPARSE_STATUS_SUCCESS) {printf("Error creating COO sparse matrix.\n");}//coo转csr格式mkl_sparse_convert_csr(cooA, SPARSE_OPERATION_NON_TRANSPOSE, &csrA);/*STEP2 根据CSR格式稀疏矩阵得到其ia,ja,a三组数据*/sparse_index_base_t indexing;int nrows;int ncols;MKL_INT* ia_csr = (MKL_INT*)mkl_malloc(nnz * sizeof(MKL_INT), 64);//csr格式的行指针MKL_INT* csr_row_end = (MKL_INT*)mkl_malloc(nnz * sizeof(MKL_INT), 64);MKL_INT* ja_csr = (MKL_INT*)mkl_malloc(nnz * sizeof(MKL_INT), 64);//csr格式的列索引float* a_csr = (float*)mkl_malloc(nnz * sizeof(float), 64);//csr格式的矩阵值//利用mkl_sparse_s_export_csr接口实现mkl_sparse_s_export_csr(csrA, &indexing, &nrows, &ncols, &ia_csr, &csr_row_end, &ja_csr, &a_csr);/*Step3 设置稀疏矩阵参数*///初始化B矩阵和X矩阵float *b = NULL;   //保存单位矩阵用于求逆float *x = NULL;   //解矩阵b = (float*)mkl_malloc(n*n * sizeof(float), 64);x = (float*)mkl_malloc(n*n * sizeof(float), 64);//将B矩阵初始化为单位阵for (int i = 0; i < n; i++) {for (int j = 0; j < n; j++) {if (i == j) {b[i*n + j] = 1;}else {b[i*n + j] = 0;}}}//初始化ptvoid *pt[64];for (int i = 0; i < 64; i++){pt[i] = 0;}//初始化矩阵相关控制参数MKL_INT maxfct = 1;MKL_INT mnum = 1;MKL_INT phase = 1;MKL_INT perm = 0;MKL_INT nrhs = n;MKL_INT error = 0;MKL_INT msglvl = 0;//设置iparm参数MKL_INT iparm[64];//批量初始化for (int i = 0; i < 64; i++){iparm[i] = 0;}iparm[0] = 1;//开启自定义iparm[12] = 1;  //提高准确性iparm[27] = 1;  //为float型iparm[34] = 1;  //初始索引为0/*Step4 分析矩阵、数值分解、求解*/phase = 13;//phase设置为13，表示分析、数值分解、求解阶段PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &n, a_csr, ia_csr, ja_csr, &perm, &nrhs, iparm, &msglvl, b, x, &error);//将解矩阵copy给输出,即A的逆矩阵for (int i = 0; i < n; i++) {memcpy(Ainv[i], x + i * n, n * sizeof(float));}//释放矩阵内存phase = -1;           //phase设置为-1PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &n, a_csr, ia_csr, ja_csr, &perm, &nrhs, iparm, &msglvl, b, x, &error);//释放内存mkl_sparse_destroy(cooA);mkl_sparse_destroy(csrA);mkl_free(x);mkl_free(b);return true;}

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

#pragma once#include#include#include "alloc.h"#include"mkl.h"#include "mkl_types.h"#include"mkl_lapacke.h"#include "mkl_spblas.h"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);sparse_matrix_t MKL_Sparse_CooXCoo(int* ia, int *ja, float *a, int rowsA, int colsA, int nnzA,int* ib, int *jb, float *b, int rowsB, int colsB, int nnzB);sparse_matrix_t MKL_Coo2Csr(int* ia, int *ja, float *a, int nrows, int ncols, int nnz);void Print_Sparse_Csr_Matrix(sparse_matrix_t csrA,int m,int n);bool MKL_Sparse_Inverse(float **Ainv, float *a, MKL_INT *ia, MKL_INT *ja, int nnz, int n, MKL_INT mtype);

#include "MKL_Sparse_Methods.h"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) {//生成csr格式稀疏矩阵sparse_matrix_t csrA, cooA;sparse_status_t status = mkl_sparse_s_create_coo(&cooA,SPARSE_INDEX_BASE_ZERO,rowsA,    // number of rowscolsA,    // number of colsnnz,  // number of nonzerosia,ja,a);if (status != SPARSE_STATUS_SUCCESS) {printf("Error creating COO sparse matrix.\n");}mkl_sparse_convert_csr(cooA, SPARSE_OPERATION_NON_TRANSPOSE, &csrA);//调用mkl稀疏矩阵与稠密矩阵乘法  C=alpha*op(A)*B+beta*Cdouble alpha = 1.0;double beta = 0.0;int M, N, K;int ncols, ldx, ldy;if (flag == 1 || flag == 2) {M = colsA;N = rowsA;K = colsC;ncols = K;ldx = K;ldy = K;}else {//默认的情况下M = rowsA;N = colsA;K = colsC;ncols = N;ldx = K;ldy = K;}//将二维稠密矩阵B转为一维float *denseB1D = (float*)mkl_malloc(N*K * sizeof(float), 64);for (int i = 0; i < N; i++) {memcpy(denseB1D + i * K, denseB[i], K * sizeof(float));}struct matrix_descr descr = { SPARSE_MATRIX_TYPE_GENERAL, SPARSE_FILL_MODE_FULL, SPARSE_DIAG_NON_UNIT };float *denseC1D = (float*)mkl_malloc(M * K * sizeof(float), 64);sparse_operation_t operation = SPARSE_OPERATION_NON_TRANSPOSE;  //默认 op(A)=A;if (flag == 0) {   //稀疏矩阵 op(A)=A operation = SPARSE_OPERATION_NON_TRANSPOSE;}if (flag == 1) {//稀疏矩阵 op(A)=AToperation = SPARSE_OPERATION_TRANSPOSE;}else if (flag == 2){    //稀疏矩阵 op(A)=AHoperation = SPARSE_OPERATION_CONJUGATE_TRANSPOSE;}else {//稀疏矩阵 op(A)=A operation = SPARSE_OPERATION_NON_TRANSPOSE;}mkl_sparse_s_mm(operation, alpha, csrA, descr, SPARSE_LAYOUT_ROW_MAJOR, denseB1D, ncols, ldx, beta, denseC1D, ldy);//将计算结果转为2维for (int i = 0; i < M; i++) {memcpy(denseC[i], denseC1D + i * K, K * sizeof(float));}//释放mkl_sparse_destroy(csrA);mkl_sparse_destroy(cooA);mkl_free(denseB1D);mkl_free(denseC1D);return true;}//根据已知坐标、元素值，创建CSR稀疏矩阵sparse_matrix_t MKL_Coo2Csr(int* ia, int *ja, float *a, int nrows, int ncols, int nnz) {//建立coo矩阵A 与 csrA矩阵sparse_matrix_t cooA, csrA;mkl_sparse_s_create_coo(&cooA,SPARSE_INDEX_BASE_ZERO,nrows,    // number of rowsncols,    // number of colsnnz,  // number of nonzerosia,ja,a);//coo转csrmkl_sparse_convert_csr(cooA, SPARSE_OPERATION_NON_TRANSPOSE, &csrA);//释放cooA矩阵mkl_sparse_destroy(cooA);//返回csrA矩阵return csrA;}void Print_Sparse_Csr_Matrix(sparse_matrix_t csrA,int m,int n) {sparse_index_base_t indexing;int nrows;int ncols;MKL_INT* csr_row_start;MKL_INT* csr_row_end;MKL_INT* csr_col_indx;float* csr_values;mkl_sparse_s_export_csr(csrA, &indexing, &nrows, &ncols, &csr_row_start, &csr_row_end, &csr_col_indx, &csr_values);float **A_dense = alloc2float(ncols, nrows);memset(A_dense[0], 0, nrows*ncols * sizeof(float));//将value转换为普通二维数组for (int i = 0; i < nrows; i++) {for (int j = csr_row_start[i]; j < csr_row_start[i + 1]; j++) {A_dense[i][csr_col_indx[j]] = csr_values[j];}}//输出稠密矩阵for (int i = 0; i < m; i++) {for (int j = 0; j < n; j++) {printf("%f ", A_dense[i][j]);}printf("\n");}free2float(A_dense);}//Coo格式×Coo格式矩阵乘法sparse_matrix_t MKL_Sparse_CooXCoo(int* ia, int *ja, float *a, int rowsA, int colsA, int nnzA,int* ib, int *jb, float *b, int rowsB, int colsB, int nnzB) {//根据坐标创建csrA和csrBsparse_matrix_t csrA, csrB, csrC;csrA = MKL_Coo2Csr(ia, ja, a, rowsA, colsA, nnzA);csrB = MKL_Coo2Csr(ib, jb, b, rowsB, colsB, nnzB);//csrC创建mkl_sparse_d_create_csr(&csrC,SPARSE_INDEX_BASE_ZERO,rowsA,    // number of rowscolsB,    // number of colsNULL,  // number of nonzerosNULL,NULL,NULL);mkl_sparse_spmm(SPARSE_OPERATION_NON_TRANSPOSE, csrA, csrB, &csrC);//释放矩阵A和Bmkl_sparse_destroy(csrA);mkl_sparse_destroy(csrB);return csrC;}//稀疏矩阵求逆bool MKL_Sparse_Inverse(float **Ainv, float *a, MKL_INT *ia, MKL_INT *ja, int nnz, int n, MKL_INT mtype) {/*STEP1 根据输入数组创建COO格式稀疏矩阵*///由于CSR格式不易表示，所以采取的路线为通过坐标创建COO格式矩阵//再通过mkl接口将COO矩阵转为CSR矩阵sparse_matrix_t csrA, cooA;sparse_status_t status = mkl_sparse_s_create_coo(&cooA,SPARSE_INDEX_BASE_ZERO,n,    // 稀疏矩阵的行、列n,nnz,  // 非零元素个数ia,//行索引                                      ja,//列索引a);//矩阵元素值if (status != SPARSE_STATUS_SUCCESS) {printf("Error creating COO sparse matrix.\n");}//coo转csr格式mkl_sparse_convert_csr(cooA, SPARSE_OPERATION_NON_TRANSPOSE, &csrA);/*STEP2 根据CSR格式稀疏矩阵得到其ia,ja,a三组数据*/sparse_index_base_t indexing;int nrows;int ncols;MKL_INT* ia_csr = (MKL_INT*)mkl_malloc(nnz * sizeof(MKL_INT), 64);//csr格式的行指针MKL_INT* csr_row_end = (MKL_INT*)mkl_malloc(nnz * sizeof(MKL_INT), 64);MKL_INT* ja_csr = (MKL_INT*)mkl_malloc(nnz * sizeof(MKL_INT), 64);//csr格式的列索引float* a_csr = (float*)mkl_malloc(nnz * sizeof(float), 64);//csr格式的矩阵值//利用mkl_sparse_s_export_csr接口实现mkl_sparse_s_export_csr(csrA, &indexing, &nrows, &ncols, &ia_csr, &csr_row_end, &ja_csr, &a_csr);/*Step3 设置稀疏矩阵参数*///初始化B矩阵和X矩阵float *b = NULL;   //保存单位矩阵用于求逆float *x = NULL;   //解矩阵b = (float*)mkl_malloc(n*n * sizeof(float), 64);x = (float*)mkl_malloc(n*n * sizeof(float), 64);//将B矩阵初始化为单位阵for (int i = 0; i < n; i++) {for (int j = 0; j < n; j++) {if (i == j) {b[i*n + j] = 1;}else {b[i*n + j] = 0;}}}//初始化ptvoid *pt[64];for (int i = 0; i < 64; i++){pt[i] = 0;}//初始化矩阵相关控制参数MKL_INT maxfct = 1;MKL_INT mnum = 1;MKL_INT phase = 1;MKL_INT perm = 0;MKL_INT nrhs = n;MKL_INT error = 0;MKL_INT msglvl = 0;//设置iparm参数MKL_INT iparm[64];//批量初始化for (int i = 0; i < 64; i++){iparm[i] = 0;}iparm[0] = 1;//开启自定义iparm[12] = 1;  //提高准确性iparm[27] = 1;  //为float型iparm[34] = 1;  //初始索引为0/*Step4 分析矩阵、数值分解、求解*/phase = 13;//phase设置为13，表示分析、数值分解、求解阶段PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &n, a_csr, ia_csr, ja_csr, &perm, &nrhs, iparm, &msglvl, b, x, &error);//将解矩阵copy给输出,即A的逆矩阵for (int i = 0; i < n; i++) {memcpy(Ainv[i], x + i * n, n * sizeof(float));}//释放矩阵内存phase = -1;           //phase设置为-1PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &n, a_csr, ia_csr, ja_csr, &perm, &nrhs, iparm, &msglvl, b, x, &error);//释放内存mkl_sparse_destroy(cooA);mkl_sparse_destroy(csrA);mkl_free(x);mkl_free(b);return true;}

#include "MKL_Sparse_Methods.h"#include "alloc.h"#define M 6#define N 5#define K 4void MKL_Sparse_CooXDense_Demo();void MKL_Sparse_CooXCoo_Demo();void MKL_Sparse_Inverse_Demo();int main() {    MKL_Sparse_CooXDense_Demo();//稀疏乘稠密    MKL_Sparse_CooXCoo_Demo();//稀疏×稀疏    MKL_Sparse_Inverse_Demo();//稀疏矩阵求逆    return 0;}//稀疏乘稠密void MKL_Sparse_CooXDense_Demo() {int flag = 0;/*flag=0时表示A(COO)*B(Dense)flag=1时表示AT(COO)*B(Dense)*/int rowsA, colsA;if (flag == 0) {rowsA = M, colsA = N;}else if (flag == 1) {rowsA = N, colsA = M;}int rowsB = N, colsB = K;int rowsC = M, colsC = K;float Atemp[M][N] = {{1,-1,-3,0,0},{-2,5,0,0,0},{0,0,4,6,4},{-4,0,2,7,0},{0,8,0,0,-5},{1,0,0,0,0},};float ATtemp[N][M] = {{1,-2,0,-4,0,1},{-1,5,0,0,8,0},{-3,0,4,2,0,0},{0,0,6,7,0,0},{0,0,4,0,-5,0},};float Btemp[N][K] = {{1,-1,-3,0},{-2,5,0,0},{0,0,4,6},{-4,0,2,7},{0,8,0,0}};//将一般二维数组转换为alloc表示float **matrixA = alloc2float(colsA, rowsA);memset(matrixA[0], 0, rowsA*colsA * sizeof(float));float **matrixB = alloc2float(colsB, rowsB);memset(matrixB[0], 0, rowsB*colsB * sizeof(float));float **matrixC = alloc2float(colsC, rowsC);memset(matrixC[0], 0, rowsC*colsC * sizeof(float));//复制二维数组到二级指针if (flag == 0) {memcpy(matrixA[0], Atemp, rowsA*colsA * sizeof(float));}else if (flag == 1) {memcpy(matrixA[0], ATtemp, rowsA*colsA * sizeof(float));}memcpy(matrixB[0], Btemp, rowsB*colsB * sizeof(float));//统计二维数组的非零元素个数int nnz = 0;for (int i = 0; i < rowsA; i++) {for (int j = 0; j < colsA; j++) {//统计nnzif (matrixA[i][j] != 0) {nnz++;}}}//获取稠密矩阵的coo形式MKL_INT *ia = (MKL_INT*)mkl_malloc(nnz * sizeof(MKL_INT), 64);MKL_INT *ja = (MKL_INT*)mkl_malloc(nnz * sizeof(MKL_INT), 64);float *a = (float*)mkl_malloc(nnz * sizeof(float), 64);//获取coo数据int k = 0;for (int i = 0; i < rowsA; i++) {for (int j = 0; j < colsA; j++) {if (matrixA[i][j] != 0.0) {a[k] = matrixA[i][j];ia[k] = i;ja[k] = j;k++;}}}//执行矩阵乘法MKL_Sparse_CooXDense(ia, ja, a, nnz, matrixB, matrixC, rowsA, colsA, colsC, flag);/* 输出结果 */printf("*************** MKL Sparse X Dense ***************\n ");printf("===============  flag=%d  ================\n", flag);for (int i = 0; i < rowsC; i++) {for (int j = 0; j < colsC; j++) {printf("%f ", matrixC[i][j]);}printf("\n");}free2float(matrixA);free2float(matrixB);free2float(matrixC);mkl_free(ia);mkl_free(ja);mkl_free(a);}//稀疏×稀疏void MKL_Sparse_CooXCoo_Demo() {int rowsA = M, colsA = N;int rowsB = N, colsB = K;int rowsC = M, colsC = K;float Atemp[M][N] = {{1,-1,-3,0,0},{-2,5,0,0,0},{0,0,4,6,4},{-4,0,2,7,0},{0,8,0,0,-5},{1,0,0,0,0},};float Btemp[N][K] = {{1,-1,-3,0},{-2,5,0,0},{0,0,4,6},{-4,0,2,7},{0,8,0,0}};//将一般二维数组转换为alloc表示float **matrixA = alloc2float(colsA, rowsA);memset(matrixA[0], 0, rowsA*colsA * sizeof(float));float **matrixB = alloc2float(colsB, rowsB);memset(matrixB[0], 0, rowsB*colsB * sizeof(float));//复制二维数组到二级指针memcpy(matrixA[0], Atemp, rowsA*colsA * sizeof(float));memcpy(matrixB[0], Btemp, rowsB*colsB * sizeof(float));//***********A矩阵稀疏表示//统计二维数组的非零元素个数int nnzA = 0;//A矩阵for (int i = 0; i < rowsA; i++) {for (int j = 0; j < colsA; j++) {//统计nnzif (matrixA[i][j] != 0) {nnzA++;}}}//获取稠密矩阵的coo形式MKL_INT *ia = (MKL_INT*)mkl_malloc(nnzA * sizeof(MKL_INT), 64);MKL_INT *ja = (MKL_INT*)mkl_malloc(nnzA * sizeof(MKL_INT), 64);float *a = (float*)mkl_malloc(nnzA * sizeof(float), 64);//获取coo数据int k = 0;for (int i = 0; i < rowsA; i++) {for (int j = 0; j < colsA; j++) {if (matrixA[i][j] != 0.0) {a[k] = matrixA[i][j];ia[k] = i;ja[k] = j;k++;}}}//***********B矩阵稀疏表示int nnzB = 0;//B矩阵for (int i = 0; i < rowsB; i++) {for (int j = 0; j < colsB; j++) {//统计nnzif (matrixB[i][j] != 0) {nnzB++;}}}//获取稠密矩阵的coo形式MKL_INT *ib = (MKL_INT*)mkl_malloc(nnzB * sizeof(MKL_INT), 64);MKL_INT *jb = (MKL_INT*)mkl_malloc(nnzB * sizeof(MKL_INT), 64);float *b = (float*)mkl_malloc(nnzB * sizeof(float), 64);//获取coo数据k = 0;for (int i = 0; i < rowsB; i++) {for (int j = 0; j < colsB; j++) {if (matrixB[i][j] != 0.0) {b[k] = matrixB[i][j];ib[k] = i;jb[k] = j;k++;}}}sparse_matrix_t csrC = MKL_Sparse_CooXCoo(ia, ja, a, rowsA, colsA, nnzA, ib, jb, b, rowsB, colsB, nnzB);printf("*************** MKL Sparse X Sparse ***************\n ");Print_Sparse_Csr_Matrix(csrC, rowsC, colsC);//打印矩阵C结果mkl_sparse_destroy(csrC);mkl_free(ia);mkl_free(ja);mkl_free(a);mkl_free(ib);mkl_free(jb);mkl_free(b);free2float(matrixA);free2float(matrixB);}//稀疏矩阵求逆void MKL_Sparse_Inverse_Demo() {int rowsA = N, colsA = N;float Atemp[N][N] = {{1,2,4,0,0},{2,2,0,0,0},{4,0,3,0,0},{0,0,0,4,0},{0,0,0,0,5},};//将一般二维数组转换为alloc表示float **matrixA = alloc2float(colsA, rowsA);memset(matrixA[0], 0, rowsA*colsA * sizeof(float));float **matrixA_inv = alloc2float(colsA, rowsA);memset(matrixA[0], 0, rowsA*colsA * sizeof(float));//复制二维数组到二级指针memcpy(matrixA[0], Atemp, rowsA*colsA * sizeof(float));//统计二维数组的非零元素个数int nnz = 0;for (int i = 0; i < rowsA; i++) {for (int j = 0; j < colsA; j++) {//统计nnzif (matrixA[i][j] != 0) {nnz++;}}}//获取稠密矩阵的coo形式MKL_INT *ia = (MKL_INT*)mkl_malloc(nnz * sizeof(MKL_INT), 64);MKL_INT *ja = (MKL_INT*)mkl_malloc(nnz * sizeof(MKL_INT), 64);float *a = (float*)mkl_malloc(nnz * sizeof(float), 64);//获取coo数据int k = 0;for (int i = 0; i < rowsA; i++) {for (int j = 0; j < colsA; j++) {if (matrixA[i][j] != 0.0) {a[k] = matrixA[i][j];ia[k] = i;ja[k] = j;k++;}}}MKL_INT mtype = 11;//设置矩阵类型为一般实矩阵MKL_Sparse_Inverse(matrixA, a, ia, ja, nnz, rowsA, mtype);printf("*************** MKL Sparse Matrix Inverse ***************\n ");for (int i = 0; i < N; i++) {for (int j = 0; j < N; j++) {printf("%f ", matrixA[i][j]);}printf("\n");}free2float(matrixA);free2float(matrixA_inv);mkl_free(ia);mkl_free(ja);mkl_free(a);}