MKL稀疏矩阵运算示例及函数封装

2023-04-27

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

0 稀疏矩阵

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

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

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

MKL中主要用到的稀疏矩阵格式有COO，CSR及CSC（与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计算结果

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

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

输出为：

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

用于计算稀疏矩阵(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的第一维
               );

此流程简述如下：

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

稀疏矩阵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 rows
        colsA,    // number of cols
        nnz,  // number of nonzeros
        ia,
        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*C
    double alpha = 1.0;
    double beta = 0.0;
    int M, N, K;
    int ncols, ldx, ldy;

    if (flag == 1 || flag == 2) {   //转置或共轭转置，AT或AH
        M = colsA;
        N = rowsA;
        K = colsC;
        ncols = K;
        ldx = K;
        ldy = K;
    }
    else {  //默认的情况下，A
        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

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

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

两个稀疏矩阵的乘法，使用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 rows
        ncols,    // number of cols
        nnz,  // number of nonzeros
        ia,
        ja,
        a);

    //coo转csr
    mkl_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

稀疏矩阵(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和csrB
    sparse_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 rows
        colsB,    // number of cols
        NULL,  // number of nonzeros
        NULL,
        NULL,
        NULL);
        //MKL乘法
    mkl_sparse_spmm(SPARSE_OPERATION_NON_TRANSPOSE, csrA, csrB, &csrC);
    //释放矩阵A和B
    mkl_sparse_destroy(csrA);
    mkl_sparse_destroy(csrB);
    return csrC;
}

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

计算结果与matlab结果一致。

2 稀疏矩阵求逆

(1)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)

输出为：

(2)MKL计算

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

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

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

在我们使用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：矩阵类型。具体取值如下：
- 1 - 实对称矩阵；
- 2 - 实对称正定矩阵；
- -2 - 实对称不定矩阵；
- 3 - 复对称矩阵；
- 11 - 实数、非对称矩阵；
phase：指定PARDISO的阶段。具体取值如：
- 11-分析阶段；
- 12-分析、数值分解阶段；
- 13-分析、数值分解、求解阶段；
- 22-数值分解阶段；
- 23-数值分解、求解阶段；
- 33-求解、迭代阶段；
- -1-释放所有矩阵内存；
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)，以下仅列出一些常用且便于理解的参数：
- 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开始索引；
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*n
MKL_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

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

完整代码

Ⅰ MKL_Sparse_Methods.h

#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);

Ⅱ MKL_Sparse_Methods.cpp

#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 rows
        colsA,    // number of cols
        nnz,  // number of nonzeros
        ia,
        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*C
    double 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

Ⅲ main.cpp

#include "MKL_Sparse_Methods.h"
#include "alloc.h"

#define M 6
#define N 5
#define K 4

void 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

文章来源于互联网：MKL稀疏矩阵运算示例及函数封装

THE END

二维码