MKL普通矩阵运算示例及函数封装
本示例将介绍MKL中的矩阵乘法和求逆,使用MKL进行此类大型矩阵运算可大量节省计算时间和空间,但由于MKL中的原生API接口繁杂,因此将常用函数封装,便于后续使用,最后在实际例子中调用接口执行想要的矩阵运算。
1 MKL矩阵乘法案例
所用示例如下,矩阵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
输出为:
![](https://img2023.cnblogs.com/blog/1602810/202304/1602810-20230423171238573-313239310.png)
(2)MKL矩阵乘法
矩阵乘法接口
/*
输入:
MatrixA 矩阵A数据,类型为float型的二维数组
rowA 矩阵A的行数
colA 矩阵A的列数
MatrixB 矩阵B数据,类型为float型的二维数组
colC 矩阵C的列数
allocType 二维矩阵定义方式,
allocType=0时,定义为A=alloc(row,col), A[col][row](第一维表示列,第二维表示行)
allocType=1时,定义为A=alloc(col,row), A[row][col](第一维表示行,第二维表示列)
输出:
MatrixC 矩阵A*B数据,类型为float型的二维数组,为矩阵乘法计算结果
*/
bool MKL_MatrixMul(float **MatrixA, int rowsA, int colsA, float **MatrixB, int colsC, float **MatrixC, int allocType);
函数代码
函数将使用MKL库中的矩阵乘法接口cblas_?gemm
实现,具体用法及参数详解见MKL库矩阵乘法(cblas_?gemm) - GeoFXR - 博客园 (cnblogs.com)
#include "MKL_Matrix_Methods.h"
//矩阵乘法
bool MKL_MatrixMul(float **MatrixA, int rowsA, int colsA, float **MatrixB, int colsC, float **MatrixC, int allocType) {
if (MatrixA == NULL) {
return false;
}
if (MatrixB == NULL) {
return false;
}
if (MatrixC == NULL) {
return false;
}
int M, N, K;
int lda, ldb, ldc;
M = rowsA;
N = colsA;
K = colsC;
float *A = NULL;
float *B = NULL;
float *C = NULL;
//由于mkl的矩阵乘法函数仅支持一维数组,需对输入进行转换
A = (float*)mkl_malloc(M*N * sizeof(float), 64);
B = (float*)mkl_malloc(N*K * sizeof(float), 64);
C = (float*)mkl_malloc(M*K * sizeof(float), 64);
if (A == NULL || B == NULL || C == NULL) {
mkl_free(A);
mkl_free(B);
mkl_free(C);
return false;
}
//赋值
int i = 0;
int j = 0;
//如果alloc是行在前列在后
if (allocType == 0) {
lda = M;
ldb = N;
ldc = M;
for (i = 0; i < N; i++) {
memcpy(A + i * M, MatrixA[i], M * sizeof(float));
}
for (i = 0; i < K; i++) {
memcpy(B + i * N, MatrixB[i], N * sizeof(float));
}
cblas_sgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, M, K, N, 1, A,
lda, B, ldb, 0, C, ldc);
for (i = 0; i < K; i++) {
memcpy(MatrixC[i], C + i * M, M * sizeof(float));
}
}
//如果alloc是列在前行在后
else {
lda = N;
ldb = K;
ldc = K;
for (i = 0; i < M; i++) {
memcpy(A + i * N, MatrixA[i], N * sizeof(float));
}
for (i = 0; i < N; i++) {
memcpy(B + i * K, MatrixB[i], K * sizeof(float));
}
cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, M, K, N, 1, A, lda, B, ldb,
0, C, ldc);
for (i = 0; i < M; i++) {
memcpy(MatrixC[i], C + i * K, K * sizeof(float));
}
}
}
在执行main.cpp中的MKL_MatrixMul_Demo()
之后,
![](https://img2023.cnblogs.com/blog/1602810/202304/1602810-20230423171343001-1413034742.png)
2 MKL矩阵求逆案例
(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)
输出为:
![](https://img2023.cnblogs.com/blog/1602810/202304/1602810-20230423171430586-1250122176.png)
(2)MKL求逆
矩阵求逆接口
/*
函数功能:基于MKL库的矩阵求逆
输入:
Matrix 输入矩阵Matrix,n行n列
n 矩阵的行、列数
输出:
Matrix Matrix 的逆,n行n列
*/
bool MKL_MatrixInverse(float**Matrix, int n)
函数代码
使用MKL中的LAPACK库计算线性方程组\(AX=B\)的解,当\(B\)为单位阵时,\(X\)即为\(A\)的逆矩阵\(A^{-1}\)。函数使用的接口为LAPACKE_sgesv,具体参数详解见MKL库线性方程组求解(LAPACKE_?gesv) - GeoFXR - 博客园 (cnblogs.com)
bool MKL_MatrixInverse(float**Matrix, int n) {
MKL_INT nrhs = n, lda = n, ldb = n;
// 创建置换矩阵,长度为n的数组
int *ipiv = (int*)mkl_malloc(n * sizeof(int), 64);
if (ipiv == NULL) {
return false;
}
// 创建MKL矩阵
float *matA = (float *)mkl_malloc(n * n * sizeof(float), 64);
if (matA == NULL) {
return false;
}
//将二维数组转换为一维MKL数组
for (int i = 0; i < n; i++) {
memcpy(matA + i * n, Matrix[i], n * sizeof(float));
}
// 创建一个单位阵B
float *matEye = (float *)mkl_malloc(n * n * sizeof(float), 64);
if (matEye == NULL) {
return false;
}
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
matEye[i * n + j] = (i == j) ? 1.0 : 0.0;
}
}
// 调用求解AX=B函数
LAPACKE_sgesv(LAPACK_ROW_MAJOR, n, nrhs, matA, lda, ipiv, matEye, ldb);
// 将MKL矩阵转换回普通二维数组
for (int i = 0; i < n; i++) {
memcpy(Matrix[i], matEye + i * n, n * sizeof(float));
}
// 释放内存
mkl_free(matA);
mkl_free(ipiv);
mkl_free(matEye);
return true;
}
在执行main.cpp中的MKL_MatrixInverse_Demo()
之后,
![](https://img2023.cnblogs.com/blog/1602810/202304/1602810-20230423171505745-414617491.png)
结果与matlab一致。
完整代码
Ⅰ MKL_Matrix_Methods.h
#pragma once
#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include"mkl.h"
#include "mkl_types.h"
#include"mkl_lapacke.h"
bool MKL_MatrixMul(float **MatrixA, int rowsA, int colsA, float **MatrixB, int colsC, float **MatrixC, int allocType);
bool MKL_MatrixInverse(float**Matrix, int n);
Ⅱ MKL_Matrix_Methods.cpp
#include "MKL_Matrix_Methods.h"
//矩阵乘法
bool MKL_MatrixMul(float **MatrixA, int rowsA, int colsA, float **MatrixB, int colsC, float **MatrixC, int allocType) {
if (MatrixA == NULL) {
return false;
}
if (MatrixB == NULL) {
return false;
}
if (MatrixC == NULL) {
return false;
}
int M, N, K;
int lda, ldb, ldc;
M = rowsA;
N = colsA;
K = colsC;
float *A = NULL;
float *B = NULL;
float *C = NULL;
//由于mkl的矩阵乘法函数仅支持一维数组,需对输入进行转换
A = (float*)mkl_malloc(M*N * sizeof(float), 64);
B = (float*)mkl_malloc(N*K * sizeof(float), 64);
C = (float*)mkl_malloc(M*K * sizeof(float), 64);
if (A == NULL || B == NULL || C == NULL) {
mkl_free(A);
mkl_free(B);
mkl_free(C);
return false;
}
//赋值
int i = 0;
int j = 0;
//如果alloc是行在前列在后
if (allocType == 0) {
lda = M;
ldb = N;
ldc = M;
for (i = 0; i < N; i++) {
memcpy(A + i * M, MatrixA[i], M * sizeof(float));
}
for (i = 0; i < K; i++) {
memcpy(B + i * N, MatrixB[i], N * sizeof(float));
}
cblas_sgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, M, K, N, 1, A,
lda, B, ldb, 0, C, ldc);
for (i = 0; i < K; i++) {
memcpy(MatrixC[i], C + i * M, M * sizeof(float));
}
}
//如果alloc是列在前行在后
else {
lda = N;
ldb = K;
ldc = K;
for (i = 0; i < M; i++) {
memcpy(A + i * N, MatrixA[i], N * sizeof(float));
}
for (i = 0; i < N; i++) {
memcpy(B + i * K, MatrixB[i], K * sizeof(float));
}
cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, M, K, N, 1, A, lda, B, ldb,
0, C, ldc);
for (i = 0; i < M; i++) {
memcpy(MatrixC[i], C + i * K, K * sizeof(float));
}
}
}
//矩阵求逆
bool MKL_MatrixInverse(float**Matrix, int n) {
MKL_INT nrhs = n, lda = n, ldb = n;
// 创建置换矩阵,长度为n的数组
int *ipiv = (int*)mkl_malloc(n * sizeof(int), 64);
if (ipiv == NULL) {
return false;
}
// 创建MKL矩阵
float *matA = (float *)mkl_malloc(n * n * sizeof(float), 64);
if (matA == NULL) {
return false;
}
//将二维数组转换为一维MKL数组
for (int i = 0; i < n; i++) {
memcpy(matA + i * n, Matrix[i], n * sizeof(float));
}
// 创建一个单位阵B
float *matEye = (float *)mkl_malloc(n * n * sizeof(float), 64);
if (matEye == NULL) {
return false;
}
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
matEye[i * n + j] = (i == j) ? 1.0 : 0.0;
}
}
// 调用求解AX=B函数
LAPACKE_sgesv(LAPACK_ROW_MAJOR, n, nrhs, matA, lda, ipiv, matEye, ldb);
// 将MKL矩阵转换回普通二维数组
for (int i = 0; i < n; i++) {
memcpy(Matrix[i], matEye + i * n, n * sizeof(float));
}
// 释放内存
mkl_free(matA);
mkl_free(ipiv);
mkl_free(matEye);
return true;
}
Ⅲ main.cpp
#include "MKL_Matrix_Methods.h"
#include "alloc.h"
#define M 6
#define N 5
#define K 4
void MKL_MatrixMul_Demo();
void MKL_MatrixInverse_Demo();
int main()
{
MKL_MatrixMul_Demo(); //矩阵乘法示例
MKL_MatrixInverse_Demo(); //矩阵求逆示例
}
//矩阵乘法
void MKL_MatrixMul_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}
};
int allocType = 1; //alloc2(col,row)列在前,行在后
//将一般二维数组转换为alloc表示
float **matrixA = alloc2float(colsA, rowsA);
memset(matrixA[0], 0, rowsA*colsA * sizeof(float));
memcpy(matrixA[0], Atemp, rowsA*colsA * sizeof(float));
float **matrixB = alloc2float(colsB, rowsB);
memset(matrixB[0], 0, rowsB*colsB * sizeof(float));
memcpy(matrixB[0], Btemp, rowsB*colsB * sizeof(float));
float **matrixC = alloc2float(colsC, rowsC);
memset(matrixC[0], 0, rowsC*colsC * sizeof(float));
MKL_MatrixMul(matrixA, rowsA, colsA, matrixB, colsC, matrixC,1); //调用MKL矩阵乘法接口
/* 输出结果 */
printf("*************** MKL Matrix Multiply ***************\n ");
for (int i = 0; i < rowsC; i++) {
for (int j = 0; j < colsC; j++) {
printf("%f ", matrixC[i][j]);
}
printf("\n");
}
free(matrixA);
free(matrixB);
free(matrixC);
}
// 矩阵求逆
void MKL_MatrixInverse_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));
//复制二维数组到二级指针
memcpy(matrixA[0], Atemp, rowsA*colsA * sizeof(float));
//求逆
MKL_MatrixInverse(matrixA, rowsA);
/* 输出结果 */
printf("*************** MKL Matrix Inverse ***************\n ");
for (int i = 0; i < rowsA; i++) {
for (int j = 0; j < colsA; j++) {
printf("%f ", matrixA[j][i]);
}
printf("\n");
}
free(matrixA);
}