MKL稀疏矩陣運算範例及函數封裝
2023-04-23 21:01:05
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 < 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)=AT
operation = SPARSE_OPERATION_TRANSPOSE;
}
else if (flag == 2)
{ //稀疏矩陣 op(A)=AH
operation = 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()後，
(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 < 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和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 < n; i++) {
for (int j = 0; j < n; j++) {
if (i == j) {
b[i*n + j] = 1;
}
else {
b[i*n + j] = 0;
}
}
}
//初始化pt
void *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設定為-1
PARDISO(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結果一致：
完整程式碼
Ⅰ MKL_Sparse_Methods.h
#pragma once
#include<stdio.h>
#include<stdlib.h>
#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 < 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)=AT
operation = SPARSE_OPERATION_TRANSPOSE;
}
else if (flag == 2)
{ //稀疏矩陣 op(A)=AH
operation = 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 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;
}
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和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_sparse_spmm(SPARSE_OPERATION_NON_TRANSPOSE, csrA, csrB, &csrC);
//釋放矩陣A和B
mkl_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;
}
}
}
//初始化pt
void *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設定為-1
PARDISO(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
#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 < rowsA; i++) {
for (int j = 0; j < colsA; j++) {
//統計nnz
if (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++) {
//統計nnz
if (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++) {
//統計nnz
if (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++) {
//統計nnz
if (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);
}