RTKLIB——matmul(矩阵乘法函数)
RTKLIB——matmul(矩阵乘法函数)
笔者个人喜欢使用malloc开辟二维矩阵进行计算,在C语言线性代数专栏中对C语言中实验矩阵乘积、转置、求行列式、求逆等方法进行了详细的介绍。
RTKLIB中矩阵乘积函数matmul提供了另一种非常棒的思路,更加高效、便捷,且能够保证相对程度
上的精度,以下是对matmul
函数的介绍:
matmul是用来进行矩阵乘法的函数,其中传入参数如下:
判断是否需要转置标志:
tr
,“NN”(case 1)、“NT”(case 2)、“TN”(case 3)、“TT”(case 4),其中,T为需要转置、N为不需要转置A矩阵的行:
n
B矩阵的列:
k
A矩阵的行=B矩阵的列:
m
A*B精度缩放因子:
alpha
C矩阵缩放因子:
beta
A
、B
、C
分别代表左矩阵A,右矩阵C,结果矩阵/原矩阵
/* multiply matrix -----------------------------------------------------------*/
extern void matmul(const char *tr, int n, int k, int m, double alpha,const double *A, const double *B, double beta, double *C)
{double d;int i,j,x,f=tr[0]=='N'?(tr[1]=='N'?1:2):(tr[1]=='N'?3:4);for (i=0;i<n;i++) {for (j=0;j<k;j++) {d=0.0;switch (f) {case 1: for (x=0;x<m;x++) {d+=A[i+x*n]*B[x+j*m]; }break;case 2: for (x=0;x<m;x++) {d+=A[i+x*n]*B[j+x*k]; }break;case 3: for (x=0;x<m;x++) {d+=A[x+i*m]*B[x+j*m]; }break;case 4: for (x=0;x<m;x++) {d+=A[x+i*m]*B[j+x*k]; }break;}if (beta==0.0) C[i+j*n]=alpha*d; else C[i+j*n]=alpha*d+beta*C[i+j*n];}}
}
上述代码需要注意几点:
- 矩阵A、B、C均是以一维数组形式来存储二维矩阵,假设A是i行j列的矩阵,
A[i][x]
的元素可用A[i+xn]
表示 - 代码中使用的是一维数组形式,因此传入参数可以是一维数组
arr[9]
,也可以是动态内存开辟的数组double* arr=(double*)malloc(sizeof(double)*9);
- 代码中判定是否需要转置是根据“N”来判定,可依据自己需求加入大小写判断等
- switch语句中务必加入break语句,若采用更严谨的写法,可加入default判断
alpha
与beta
分别为矩阵AB乘积的缩放因子和C的缩放因子- 代码中存储数据是以列优先进行存储
什么是列优先
通常来讲,C语言中存储矩阵是以行优先的形式进行存储的。
例如下面3*3的矩阵,在C语言中可以使用如下方法:arr[3][3]={1,2,3,4,5,6,7,8,9},如果采用一维数组存储下面矩阵,通常会想到arr[9]={1,2,3,4,5,6,7,8,9},
然而riklib中并不是这样的!!!rtklib中采用列主元的方式进行存储,也就是arr[9]={1,4,7,2,5,6,3,6,9},这种存储方式与Fortran中矩阵存储的形式相同(如果学过Fortran的同学都知道,列优先存储时逻辑上会更加清晰)
[ 1 2 3 4 5 6 7 8 9 ] \left[ \begin{matrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \end{matrix} \right] 147258369
搞懂列优先后,再看rktlib中矩阵运算规则就一目了然了。
为什么要对矩阵进行缩放?
在矩阵乘法中,如果乘积结果的值非常大或非常小,那么在计算机中表示这个值可能会导致精度损失或溢出。因此,在矩阵乘法中,常常需要对结果进行缩放,以保持数值的合适范围,同时尽可能地保留精度。
缩放的方法通常是对结果矩阵中的每个元素都乘以一个常数因子,这个常数因子通常被称为缩放因子或比例因子。上述代码中,如果
beta
等于 0,表示不需要对原矩阵C进行缩放,只需乘对乘积结果乘以一个常数因子alpha
;如果beta
不等于 0,则需要将原有的结果矩阵 C 乘以一个常数因子beta
,再加上新计算出的结果矩阵alpha*d
,以获得最终的结果矩阵 C。
下面是使用matmul的简单示例:
void Test()
{double arr1[]={1,1,2,2,3,3};//row=3,col=2double arr2[]={1,1,1,1,1,1};//row=2,col=3double* arr=(double*)malloc(sizeof(double)*9);memset(arr,0,sizeof(double)*9);matmul("NN",3, 3, 2, 1.0,arr1,arr2,0.0,arr);for(i=0;i<3;i++){for(j=0;j<3;j++){printf("%2.0lf\t",arr[i+j*3]);}putchar('\n');}
}
其中:matmul("NN",3, 3, 2, 1.0,arr1,arr2,0.0,arr);
表示矩阵AB不需要转置,其中A为3行2列,B为2行3列,AB乘积的缩放因子为1.0(不需要缩放),C的缩放因子因子为0,表示不需要对C进行缩放。
注意:因为有缩放因子的存在,如果原矩阵C为空,需要将其初始化。
个人改进方法
RTKLIB中的乘法函数采用一维数组表达二维矩阵,提高了开辟内存时的效率与访问速度,但也因此会有以下几个问题:
1)赋值时需要以
列优先
原则赋值,代码可读性较差2)传入参数必须手动输入矩阵的行列、容易造成人工输入bug
以下是笔者比较喜爱的一种运算方式,如果数据量不是非常庞大,下面代码更好理解与方便:
1)创建了
Matrix
结构体,用于存储矩阵信息(如果是Winodows系统,可以使用_msize函数访问内存进行维数判定,那么就需要专门创建结构体存储信息了,但是该方法Mac/Linux不适用)2)加入了几个暴力检查,主要针对输入错误
3)代码可读性高,代码中存储数值的逻辑符合正常逻辑
但是需要注意以下几点:
1)
MakeMatrix
创建矩阵的逻辑为一列一列循环创建,因此要释放内存时,需要使用free_Matrix
依次循环释放内存2)手动赋值时较为麻烦
3)内存管理不如上述方式
//矩阵结构体声明
typedef struct Matrix
{int row;int col;double** data;
}Matrix,*pMatrix;
//循环释放内存函数
void free_Matrix(Matrix arr)
{int i;for (i = 0; i < arr.row; i++){free(arr.data[i]);}
}
//创建矩阵函数
Matrix MakeMatrix(int row, int col)
{int i = 0;Matrix arr = {0};arr.row = row;arr.col = col;arr.data = (double **)malloc(sizeof(double *) * arr.row);if (arr.data == NULL)exit(-1);for (i = 0; i < arr.row; i++){arr.data[i] = (double *)malloc(sizeof(double) * arr.col);memset(arr.data[i], 0, sizeof(double) * arr.col);}return arr;
}
//打印矩阵函数
void PrintMatrix(Matrix arr)
{int i, j;for (i = 0; i < arr.row; i++){for (j = 0; j < arr.col; j++){printf("%lf\t", arr.data[i][j]);}putchar('\n');}
}
//利用Matrix结构体实现matmul矩阵乘法
extern void matmul2(const char *tr, const Matrix A, const Matrix B, const Matrix C, double alpha, double beta)
{double d;int i, j, x, f, row, col, k;f = tr[0] == 'N' ? (tr[1] == 'N' ? 1 : 2) : (tr[1] == 'N' ? 3 : 4);//这里对比一下就能看出来RTKLIB中的精妙之处,因为函数采用二维矩阵的形式,因此在计算前要根据转置情况重新判断行列情况。row = A.row, col = B.col, k = A.col;if (f == 2)col = B.row;else if (f == 3){row = A.col;k = A.row;}else if (f == 4){row = A.col;col = B.row;k = A.row;}for (i = 0; i < row; i++){for (j = 0; j < col; j++){d = 0.0;switch (f){case 1:for (x = 0; x < k; x++){d += A.data[i][x] * B.data[x][j];}break;case 2:for (x = 0; x < k; x++){d += A.data[i][x] * B.data[j][x];}break;case 3:for (x = 0; x < k; x++){d += A.data[x][i] * B.data[x][j];}break;case 4:for (x = 0; x < k; x++){d += A.data[x][i] * B.data[j][x];}break;}if (beta == 0.0)C.data[i][j] = alpha * d;elseC.data[i][j] = alpha * d + beta * C.data[i][j];}}
}
RTKLIB——matmul(矩阵乘法函数)相关推荐
- tf.matmul - 矩阵乘法
tf.matmul - 矩阵乘法 https://github.com/tensorflow/docs/tree/r1.4/site/en/api_docs/api_docs/python/tf si ...
- Python:矩阵乘法函数(高教社,《Python编程基础及应用》习题4-11)
学习记录与分享 PTA教学平台 题目 设计一个Python函数,计算两个矩阵(二维列表)的乘积. 函数接口定义 def multiply(a,b,p,q,r) a是一个p行q列的二维列表:b是一个q行 ...
- 6-2 矩阵乘法函数(高教社,《Python编程基础及应用》习题4-11)
设计一个Python函数,计算两个矩阵(二维列表)的乘积. a1.png 函数接口定义: def multiply(a,b,p,q,r) a是一个p行q列的二维列表:b是一个q行r列的二维列表: 应返 ...
- C语言 · 矩阵乘法
问题描述 给定一个N阶矩阵A,输出A的M次幂(M是非负整数) 例如: A = 1 2 3 4 A的2次幂 7 10 15 22 输入格式 第一行是一个正整数N.M(1<=N<=30, 0& ...
- 循环取矩阵的某行_1.2 震惊! 某大二本科生写的矩阵乘法吊打Mathematica-线性代数库BLAS-矩阵 (上)...
本文是 1. 线性代数库BLASzhuanlan.zhihu.com 系列的第二篇, 将讲述矩阵类的结构和矩阵基础运算的AVX2加速算法. 1. 矩阵类的结构 在讲述矩阵各种算法之前很有必要详解一下 ...
- 稀疏矩阵加法运算_1.2 震惊! 某大二本科生写的矩阵乘法吊打Mathematica-线性代数库BLAS-矩阵 (上)...
本文是 1. 线性代数库BLASzhuanlan.zhihu.com 系列的第二篇, 将讲述矩阵类的结构和矩阵基础运算的AVX2加速算法. 1. 矩阵类的结构 在讲述矩阵各种算法之前很有必要详解一下 ...
- cuBLAS矩阵乘法
cuBLAS是cuda封装好的一个数学库,头文件为<cublas_v2.h>,其中的矩阵乘法函数是我们做深度学习绕不开的函数, 下面是一个常用的函数,后面又有扩展,但是核心函数使用逻辑一样 ...
- RTKLIB 矩阵相乘函数matmul
在对RTKLIB进行二次开发时一定会用到矩阵相乘函数matmul. extern void matmul(const char *tr, int n, int k, int m, double alp ...
- tensorflow tf.matmul() (多维)矩阵相乘(多维矩阵乘法)
@tf_export("matmul") def matmul(a,b,transpose_a=False,transpose_b=False,adjoint_a=False,ad ...
最新文章
- 用eclipse创建WebService Step by Step
- ASP.NET3.5问题集
- vba 判断文本框内容是否为空_校验数据一旦失败,VBA代码自动控制焦点返回的另一备选方案...
- 向maven中央仓库提交jar
- 树莓派i2c python_树莓派2 python i2cPython中chr、unichr、ord字符函数之间的对比
- SAP Spartacus ProductService.get的几个调用场景
- php switch 函数,php switch case用法与实例教程
- 人力资源社会保障部关于公布国家职业资格目录的通知
- 贵阳中职计算机学什么区别,贵阳中职计算机专业
- 【vue】 vue中的query 路由传值的方式
- 基于nvidia的ffmpeg编解码加速
- matlab画圆的命令_matlab画圆命令资料
- 2021年T电梯修理免费试题及T电梯修理考试试卷
- 赛事相关 | 腾讯觅影×腾讯云TI平台,锁了
- 第043篇:VBA之单元格简写与引用、值与地址
- 获取今天是星期几的四种写法
- Docker安装Redis方法
- Opencv-Python-导向滤波快速导向滤波
- 热电阻 热电偶 测量电路_热电偶热电阻原理及常见故障处理
- 解释太阳能量来源《张朝阳的物理课》估算太阳寿命约百亿年
热门文章
- 1242 -- 绩点换算
- HANA tenant backup
- 绝地求生体验服显示服务器以满,绝地求生:体验服新增8排,信号枪可以在这里找到!...
- 百度云新用户组队送100元代金券,半年付100元1核1GVPS
- python数据结构考试题库_2020中国大学慕课moocPython编程基础题目答案
- 通过ROS控制真实机械臂(8)---延时时间精确控制
- 《有限元法》学习笔记2 求表面力产生的节点载荷
- Python黑魔法手册记录
- Java实现英文段落分句,python:对英文段落进行分句(对一段英语进行整句切分,切分句子)...
- C#笔记之移位运算符