MKL普通矩阵运算示例及函数封装

本示例将介绍MKL中的矩阵乘法和求逆，使用MKL进行此类大型矩阵运算可大量节省计算时间和空间，但由于MKL中的原生API接口繁杂，因此将常用函数封装，便于后续使用，最后在实际例子中调用接口执行想要的矩阵运算。

1 MKL矩阵乘法案例

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

(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

输出为：

矩阵乘法接口

/*
输入：
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);

}

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

	}
	// 创建一个单位阵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函数

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

	}
	// 创建一个单位阵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}
	};

	//将一般二维数组转换为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);
}

本栏目推荐文章
• 函数的参数类型
• 【C语言】函数的递归调用
• AntDesignBlazor示例——暗黑模式
• NUS CS1101S：SICP JavaScript 描述：一、使用函数构建抽象
• asp.net mvc4 controller构造函数
• 函数返回值
• 10-函数进阶-作用域
• 11-函数进阶-参数类型
• 9-函数基础
• vim中函数的跳转与回退

矩阵 示例 函数 MKL矩阵 示例 函数mkl mkl 数组 矩阵 函数sortrows 矩阵 数值 函数 学期 矩阵 函数 语言 方法 示例 函数 常见sql requestanimationframe示例 函数amp 钩子 示例 函数golang 函数 全网 示例 源代码 示例 函数golang pipe