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矩阵的行：
nB矩阵的列：
kA矩阵的行=B矩阵的列：
mA*B精度缩放因子：
alphaC矩阵缩放因子：
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=2
double arr2[]={1,1,1,1,1,1};//row=2,col=3
double* 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）手动赋值时较为麻烦文章来源：https://www.toymoban.com/news/detail-621617.html
3）内存管理不如上述方式文章来源地址https://www.toymoban.com/news/detail-621617.html
//矩阵结构体声明
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;
else
C.data[i][j] = alpha * d + beta * C.data[i][j];
}
}
}
到了这里，关于RTKLIB——matmul(矩阵乘法函数）的文章就介绍完了。如果您还想了解更多内容，请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章，希望大家以后多多支持TOY模板网！
