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

输出为：

### (2)MKL矩阵乘法

矩阵乘法接口

/*输入：MatrixA  矩阵A数据，类型为float型的二维数组rowA  矩阵A的行数colA  矩阵A的列数MatrixB  矩阵B数据，类型为float型的二维数组colC  矩阵C的列数输出：MatrixC  矩阵A*B数据，类型为float型的二维数组，为矩阵乘法计算结果*/bool MKL_MatrixMul(float **MatrixA, int rowA, int colA, float **MatrixB, int colC, float **MatrixC) ;

函数代码

函数将使用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) {if (MatrixA == NULL) {return false;}if (MatrixB == NULL) {return false;}if (MatrixC == NULL) {return false;}int M = rowsA; int N = colsA;int K = colsC;float *A = NULL;float *B = NULL;float *C = NULL;int lda = N;int ldb = K;int ldc = K;//由于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;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));}mkl_free(A);mkl_free(B);mkl_free(C);return true;}

main函数

#include "MKL_Matrix_Methods.h"#include "alloc.h"#define M 6#define N 5#define K 4// 矩阵乘法int main(){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));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);//调用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);}

结果与matlab一致。

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)

输出为：

(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));}// 创建一个单位阵Bfloat *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函数

#include "MKL_Matrix_Methods.h"#include "alloc.h"#define N 5// 矩阵乘法int main(){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);}

结果与matlab一致。

完整代码

Ⅰ MKL_Matrix_Methods.h

#pragma once#include#include#include#include"mkl.h"#include "mkl_types.h"#include"mkl_lapacke.h"bool MKL_MatrixMul(float **MatrixA, int rowA, int colA, float **MatrixB, int colC, float **MatrixC);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) {if (MatrixA == NULL) {return false;}if (MatrixB == NULL) {return false;}if (MatrixC == NULL) {return false;}int M = rowsA; int N = colsA;int K = colsC;float *A = NULL;float *B = NULL;float *C = NULL;int lda = N;int ldb = K;int ldc = K;//由于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;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));}mkl_free(A);mkl_free(B);mkl_free(C);return true;}//矩阵求逆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));}// 创建一个单位阵Bfloat *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 4void 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}};//将一般二维数组转换为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);//调用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);}

标签：