首页 > 其他分享 >MKL普通矩阵运算示例及函数封装

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

时间：2023-04-23 17:47:43浏览次数：50  
标签：MKL 封装 示例 int mkl float 矩阵 sizeof

本示例将介绍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));

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

标签：MKL,封装,示例,int,mkl,float,矩阵,sizeof
From： https://www.cnblogs.com/GeophysicsWorker/p/17347228.html

相关文章

  • ai问答：使用 Vue3 组合式API 和 TypeScript 封装 echarts 折线图
    <template><divref="chart"style="height:500px;"></div></template><scriptlang="ts">import{ref,onMounted,watch}from'vue'import*asechartsfrom'echarts'e......
  • 授权实现 封装权限信息
    限制访问资源所需权限​SpringSecurity为我们提供了基于注解的权限控制方案，这也是我们项目中主要采用的方式。我们可以使用注解去指定访问对应的资源所需的权限。​但是要使用它我们需要先开启相关配置我们前面在写UserDetailsServiceImpl的时候说过，在查询出用户后还要获取......
  • php 无限极分类 封装
    <?phpnamespaceApp\Services;useIlluminate\Http\Request;/***ClassPendingService*@packageApp\Service*无限分类公共类*/classLimitlessService{protected$_request;//publicfunction__construct(Request$request)//{//......
  • python matplotlib 散点图的拟合直线的简单示例
     #samplepointsX=[0,5,10,15,20]Y=[0,7,10,13,20]#solveforaandbdefbest_fit(X,Y):xbar=sum(X)/len(X)ybar=sum(Y)/len(Y)n=len(X)#orlen(Y)numer=sum([xi*yiforxi,yiinzip(X,Y)])-n*xbar*y......
  • 设置权限所需资源、封装权限信息
    设置权限所需资源SpringSecurity为我们提供了基于注解的权限控制方案，这也是我们项目中主要采用的方式。我们可以使用注解去指定访问对应的资源所需的权限。但是要使用它我们需要先开启相关配置。@EnableGlobalMethodSecurity(prePostEnabled=true)在SecurityConfig配置类......
  • SPI机制的简单示例？
    我们现在需要使用一个内容搜索接口，搜索的实现可能是基于文件系统的搜索，也可能是基于数据库的搜索。    可以看到输出结果：文件搜索helloworld如果在com.cainiao.ys.spi.learn.Search文件里写上两个实现类，那最后的输出结果就是两行了。这就是因为ServiceLoader.lo......
  • 设计模式-模板模式在Java中的使用示例-悍马模型制造示例
    场景设计模式-模板模式在Java中的使用示例：设计模式-模板模式在Java中的使用示上面整理了模板模式的使用示例，为加强理解特记录另一个使用示例，以下示例摘自设计模式之禅第二版。模板方法模式定义一个操作中的算法的框架，而将一些步骤延迟到子类中。使得子类可以不改变一个算法的结构即......
  • 使用mybatisPlus修改数据-示例
    mapperimportcom.atguigu.yygh.model.hosp.HospitalSet;importcom.baomidou.mybatisplus.core.mapper.BaseMapper;publicinterfaceHospitalSetMapperextendsBaseMapper<HospitalSet>{} serviceimportcom.atguigu.yygh.model.hosp.HospitalSet;impor......
  • 2023-04-21：用go语言重写ffmpeg的metadata.c示例。
    2023-04-21：用go语言重写ffmpeg的metadata.c示例。答案2023-04-21：这段Go代码演示了如何使用ffmpeg-go库中的函数来读取多媒体文件元数据，包括视频、音频等信息。它的大体过程如下：设置环境变量以加载FFmpeg动态链接库这里将FFmpeg库中的各个动态链接库路径添加到环境变......
  • 高频RF电路里面都见不到大封装的电容？
    俗话说：低频选大电容，高频选小电容，所以在高频电路中都见不到大封装电容电容的封装尺寸越小，其所包含的等效电感越小，因为这样电容里面的金属板或导线就越小。电感是阻碍电流的变化，也就是频率越高，阻抗越大，所以在高频的世界里，如果需要电容，就要使用尺寸小的，要不然大尺寸的电容可能会变成......

赞助商

阅读排行

网站主页关于我们联系我们网站地图

本网站内容转载自其他媒体，侵权联系[admin##ips99.com]。

Copyright © 2020-2023 IPS99 版权所有 IPS99