RTKLIB——matmul(矩阵乘法函数）

笔者个人喜欢使用malloc开辟二维矩阵进行计算，在C语言线性代数专栏中对C语言中实验矩阵乘积、转置、求行列式、求逆等方法进行了详细的介绍。
RTKLIB中矩阵乘积函数matmul提供了另一种非常棒的思路，更加高效、便捷，且能够保证相对程度上的精度，以下是对matmul函数的介绍：

matmul是用来进行矩阵乘法的函数，其中传入参数如下：

判断是否需要转置标志：tr，“NN”、“NT”、“TN”、“TT”,其中，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;ifor (j=0;jd=0.0;switch (f) {case 1: for (x=0;xd+=A[i+x*n]*B[x+j*m]; }break;case 2: for (x=0;xd+=A[i+x*n]*B[j+x*k]; }break;case 3: for (x=0;xd+=A[x+i*m]*B[x+j*m]; }break;case 4: for (x=0;xd+=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的缩放因子。
  上述代码需要注意几点：

• 矩阵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的同学都知道，列优先存储时逻辑上会更加清晰）

[123456789]\left[ \begin{matrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \end{matrix} \right] ​147​258​369​​

搞懂列优先后，再看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为空，需要将其初始化。