简易C语言矩阵运算库

参考网址：

异想家纯C语言矩阵运算库 - Sandeepin - 博客园

这次比opencv快⑥倍！！！

参考上述网址，整理了一下代码：

//main.c#include <stdio.h>
#include <stdlib.h>
#include <string.h>#include "matlib.h"typedef struct 
{int x;int y;
} Point_t;typedef struct 
{Point_t stSrc;Point_t stDsc;
} Perspective_map_t;//解非齐次线性方程组Ax=b求x
int get_perspective_transform(Perspective_map_t stPointMap[4] ,Matrix_t *pstPerspectiveMat)
{int ret = 0;Matrix_t * pstMatrix_A = NULL;Matrix_t * pstMatrix_InvA = NULL;Matrix_t * pstMatrix_x = NULL;Matrix_t * pstMatrix_b = NULL;#define GET_2D_ARRAY_VAL(buf, i, j, column) *((buf) + (i)*(column) + (j))if(!pstPerspectiveMat || pstPerspectiveMat->row != 3 || pstPerspectiveMat->column!= 3){printf("Param err\n");return -1;}pstMatrix_A = MatCreate(8, 8);pstMatrix_InvA = MatCreate(8, 8);pstMatrix_x = MatCreate(8, 1);pstMatrix_b = MatCreate(8, 1);for(int i = 0; i < 8; i += 2){GET_2D_ARRAY_VAL(pstMatrix_A->matrix_data, i, 0, 8) = stPointMap[i /2 ].stSrc.x;GET_2D_ARRAY_VAL(pstMatrix_A->matrix_data, i, 1, 8) = stPointMap[i /2 ].stSrc.y;GET_2D_ARRAY_VAL(pstMatrix_A->matrix_data, i, 2, 8) = 1;GET_2D_ARRAY_VAL(pstMatrix_A->matrix_data, i, 6, 8) = -stPointMap[i /2 ].stSrc.x * stPointMap[i /2 ].stDsc.x;GET_2D_ARRAY_VAL(pstMatrix_A->matrix_data, i, 7, 8) = -stPointMap[i /2 ].stDsc.x * stPointMap[i /2 ].stSrc.y;GET_2D_ARRAY_VAL(pstMatrix_A->matrix_data, i+1, 3, 8) = stPointMap[i /2 ].stSrc.x;GET_2D_ARRAY_VAL(pstMatrix_A->matrix_data, i+1, 4, 8) = stPointMap[i /2 ].stSrc.y;GET_2D_ARRAY_VAL(pstMatrix_A->matrix_data, i+1, 5, 8) = 1;GET_2D_ARRAY_VAL(pstMatrix_A->matrix_data, i+1, 6, 8) = -stPointMap[i /2 ].stSrc.x * stPointMap[i /2 ].stDsc.y;GET_2D_ARRAY_VAL(pstMatrix_A->matrix_data, i+1, 7, 8) = -stPointMap[i /2 ].stSrc.y * stPointMap[i /2 ].stDsc.y;  pstMatrix_b->matrix_data[i] = stPointMap[i /2 ].stDsc.x;pstMatrix_b->matrix_data[i + 1] = stPointMap[i /2 ].stDsc.y;}ret = MatInv(pstMatrix_A, pstMatrix_InvA);ret |= MatMul(pstMatrix_InvA, pstMatrix_b, pstMatrix_x);memcpy(pstPerspectiveMat->matrix_data, pstMatrix_x->matrix_data, pstMatrix_x->row * pstMatrix_x->column * sizeof(double));GET_2D_ARRAY_VAL(pstPerspectiveMat->matrix_data, 2, 2, 3) = 1;MatRelease(pstMatrix_A);MatRelease(pstMatrix_InvA);MatRelease(pstMatrix_x);MatRelease(pstMatrix_b);#undef GET_2D_ARRAY_VALreturn ret;
}int main(int argc, char** argv)
{Matrix_t Matrix_A;double A[] = { -3, 2, -5, -1, 0, -2, 3, -4, 1 };Matrix_A.row = 3;Matrix_A.column = 3;Matrix_A.matrix_data = A;Matrix_t Matrix_B;double B[] = { 1, 4, 7, 3, 0, 5, -1, 9, 11 };Matrix_B.row = 3;Matrix_B.column = 3;Matrix_B.matrix_data = B;Matrix_t Matrix_C;double C[] = { 1, 2, 3, 4, 5, 6, 7, 8, 9 };Matrix_C.row = 3;Matrix_C.column = 3;Matrix_C.matrix_data = C;Matrix_t Matrix_Adj;double Adj[9] = {0};Matrix_Adj.row = 3;Matrix_Adj.column = 3;Matrix_Adj.matrix_data = Adj;Matrix_t Matrix_Out;double Out[9] = {0};Matrix_Out.row = 3;Matrix_Out.column = 3;Matrix_Out.matrix_data = Out;printf("A+B：\n");MatAdd(&Matrix_A, &Matrix_B, &Matrix_Out); MatShow(&Matrix_Out);printf("A-B：\n");MatSub(&Matrix_A, &Matrix_B, &Matrix_Out); MatShow(&Matrix_Out);printf("A*B：\n");MatMul(&Matrix_A, &Matrix_B, &Matrix_Out);MatShow(&Matrix_Out);printf("2*C：\n");MatMulk(&Matrix_C, 2, &Matrix_Out);MatShow(&Matrix_Out);printf("C的转置：\n");MatT(&Matrix_C, &Matrix_Out); MatShow(&Matrix_Out);double det = 0.0;printf("B的行列式值：\n");MatDet(&Matrix_B, &det);printf("%16lf\n", det);printf("C的行列式值：\n");MatDet(&Matrix_C, &det);printf("%16lf\n", det);//矩阵的逆printf("A的逆：\n");MatInv(&Matrix_A, &Matrix_Out); MatShow(&Matrix_Out);printf("C的逆：\n");MatInv(&Matrix_C, &Matrix_Out);MatShow(&Matrix_Out);//矩阵代数余子式printf("C的(0,0)元素的代数余子式：\n");printf("%16lf\n", MatACof(C, 3, 0, 0));//矩阵伴随矩阵printf("A的伴随矩阵：\n");MatAdj(&Matrix_A, &Matrix_Adj); MatShow(&Matrix_Adj);//蒋方正矩阵库应用：多元线性回归//多元线性回归公式：β=(X'X)^(-1)X'Ydouble X[15][5] = {1, 316, 1536, 874, 981,//第一列要补11, 385, 1771, 777, 1386,1, 299, 1565, 678, 1672,1, 326, 1970, 785, 1864,1, 441, 1890, 785, 2143,1, 460, 2050, 709, 2176,1, 470, 1873, 673, 1769,1, 504, 1955, 793, 2207,1, 348, 2016, 968, 2251,1, 400, 2199, 944, 2390,1, 496, 1328, 749, 2287,1, 497, 1920, 952, 2388,1, 533, 1400, 1452, 2093,1, 506, 1612, 1587, 2083,1, 458, 1613, 1485, 2390};double Y[15][1] = {3894,4628,4569,5340,5449,5599,5010,5694,5792,6126,5025,5924,5657,6019,6141};Matrix_t *pstMatrix_X = NULL;Matrix_t *pstMatrix_Y = NULL;Matrix_t *pstMatrix_XT = NULL;Matrix_t *pstMatrix_XTX = NULL;Matrix_t *pstMatrix_InvXTX = NULL;Matrix_t *pstMatrix_InvXTXXT = NULL;Matrix_t *pstMatrix_InvXTXXTY = NULL;pstMatrix_X = MatCreateInitFromArray(15, 5, (double *)X);pstMatrix_Y = MatCreateInitFromArray(15, 1, (double *)Y);pstMatrix_XT = MatCreate(5, 15);pstMatrix_XTX = MatCreate(5, 5);pstMatrix_InvXTX = MatCreate(5, 5);pstMatrix_InvXTXXT = MatCreate(5, 15);pstMatrix_InvXTXXTY = MatCreate(5, 1);//多元线性回归公式：β=(X'X)^(-1)X'YMatT(pstMatrix_X, pstMatrix_XT);MatMul(pstMatrix_XT, pstMatrix_X, pstMatrix_XTX);MatInv(pstMatrix_XTX, pstMatrix_InvXTX);MatMul(pstMatrix_InvXTX, pstMatrix_XT, pstMatrix_InvXTXXT);MatMul(pstMatrix_InvXTXXT, pstMatrix_Y, pstMatrix_InvXTXXTY);printf("多元线性回归β系数：\n");MatShow(pstMatrix_InvXTXXTY);//保存矩阵到csvMatWrite("MulLinearBetaCoef.csv", pstMatrix_InvXTXXTY);MatRelease(pstMatrix_X);MatRelease(pstMatrix_Y);MatRelease(pstMatrix_XT);MatRelease(pstMatrix_XTX);MatRelease(pstMatrix_InvXTX);MatRelease(pstMatrix_InvXTXXT);MatRelease(pstMatrix_InvXTXXTY);printf("透视变换：\n");Perspective_map_t stPointMap[4] = {{{ 50 , 50 }  ,{ 0  , 0  }},{{ 200, 50 }  ,{ 300, 0  }},{{ 50 , 200}  ,{ 0  , 400}},{{ 200, 200}  ,{ 300, 400}},};Matrix_t *pstPerspectiveMat = NULL;Matrix_t *pstPerspectiveInvMat = NULL;pstPerspectiveMat = MatCreate(3, 3);pstPerspectiveInvMat = MatCreate(3, 3);	//获取透视变换矩阵get_perspective_transform(stPointMap , pstPerspectiveMat);MatShow(pstPerspectiveMat);//反透视变换矩阵MatInv(pstPerspectiveMat, pstPerspectiveInvMat);MatShow(pstPerspectiveInvMat);//保存矩阵到csvMatWrite("PerspectiveMat.csv", pstPerspectiveInvMat);MatRelease(pstPerspectiveMat);MatRelease(pstPerspectiveInvMat);//读取csv文件数据到矩阵printf("读取csv文件：\n");Matrix_t *pstReadMat = NULL;pstReadMat = MatCreate(3, 3);MatRead("PerspectiveMat.csv", pstReadMat);MatShow(pstReadMat);return 0;
}

//matlib.h#ifndef __MATLIB_H__
#define __MATLIB_H__//定义csv文件读写矩阵函数，用于测试用例
#define USE_CSV_FILE_TESTtypedef struct 
{int row; //行int column; //列double *matrix_data; //缓冲区
} Matrix_t;Matrix_t *MatCreate(int row, int column);
Matrix_t *MatCreateInitFromArray(int row, int column, double *Array);
void MatRelease(Matrix_t *Mat);
int MatCopy(Matrix_t *dsc, Matrix_t *src);int MatAdd(Matrix_t* A, Matrix_t* B, Matrix_t *Out);
int MatSub(Matrix_t* A, Matrix_t* B, Matrix_t *Out);
int MatMul(Matrix_t* A, Matrix_t* B, Matrix_t *Out);
int MatMulk(Matrix_t* A, double k, Matrix_t *Out);
int MatT(Matrix_t* A, Matrix_t *Out);
int MatDet(Matrix_t* A, double *Out);
int MatInv(Matrix_t* A, Matrix_t *Out);
double MatACof(double *A, int row, int m, int n);
int MatAdj(Matrix_t* A, Matrix_t *Out);void MatShow(Matrix_t* Mat);#ifdef USE_CSV_FILE_TEST
int MatRead(char *InCsvFileName, Matrix_t *Out);
int MatWrite(const char *OutCsvFileName, Matrix_t *In);
#endif //USE_CSV_FILE_TEST#endif //__MATLIB_H__

//matlib.c#include <stdio.h>
#include <stdlib.h>
#include <string.h>#include "matlib.h"// (det用)功能：求逆序对的个数
static int inver_order(int *list, int n)
{int ret = 0;for (int i = 1; i < n; i++)for (int j = 0; j < i; j++)if (list[j] > list[i])ret++;return ret;
}// (det用)功能：符号函数，正数返回1，负数返回-1
static int sgn(int order)
{return order % 2 ? -1 : 1;
}// (det用)功能：交换两整数a、b的值
static void swap(int *a, int *b)
{int m;m = *a;*a = *b;*b = m;
}// 功能：求矩阵行列式的核心函数
static double det(double *p, int n, int k, int *list, double sum)
{if (k >= n){int order = inver_order(list, n);double item = (double)sgn(order);for (int i = 0; i < n; i++){//item *= p[i][list[i]];item *= *(p + i * n + list[i]);}return sum + item;}else{for (int i = k; i < n; i++){swap(&list[k], &list[i]);sum = det(p, n, k + 1, list, sum);swap(&list[k], &list[i]);}}return sum;
}Matrix_t *MatCreate(int row, int column)
{unsigned int len = row*column*sizeof(double);Matrix_t *pMatrix = malloc(sizeof(Matrix_t));if(!pMatrix){return NULL;}pMatrix->row = row;pMatrix->column = column;pMatrix->matrix_data = malloc(len);if(!pMatrix->matrix_data){free(pMatrix);return NULL;}memset(pMatrix->matrix_data, 0, len);return pMatrix;
}Matrix_t *MatCreateInitFromArray(int row, int column, double *Array)
{Matrix_t *pMatrix = MatCreate(row, column);memcpy(pMatrix->matrix_data, Array, row * column * sizeof(double));return pMatrix;
}void MatRelease(Matrix_t *Mat)
{if(Mat){if(Mat->matrix_data){free(Mat->matrix_data);}free(Mat);}
}// 功能：矩阵显示
// 形参：(输入)矩阵首地址指针Mat，矩阵行数row和列数col。
// 返回：无
void MatShow(Matrix_t* Mat)
{printf("\n");for (int i = 0; i < Mat->row * Mat->column; i++){printf("%16lf ", Mat->matrix_data[i]);if (0 == (i + 1) % Mat->column) {printf("\n");}}printf("\n");
}// 功能：矩阵相加
// 形参：(输入)矩阵A首地址指针A，矩阵B首地址指针B，矩阵A(也是矩阵B)行数row和列数col
// 返回：A+B
int MatAdd(Matrix_t* A, Matrix_t* B, Matrix_t* Out)
{if(!A || !B || !Out ||A->row != B->row || A->column != B->column || A->column != Out->column || A->row != Out->row ){printf("Input Matrix can not Add!\n");return -1;}for (int i = 0; i < A->row; i++)for (int j = 0; j < A->column; j++){Out->matrix_data[A->column*i + j] = A->matrix_data[A->column*i + j] + B->matrix_data[A->column*i + j];}return 0;
}// 功能：矩阵相减
// 形参：(输入)矩阵A，矩阵B，矩阵A(也是矩阵B)行数row和列数col
// 返回：A-B
int MatSub(Matrix_t* A, Matrix_t* B, Matrix_t *Out)
{if(!A || !B || !Out ||A->row != B->row || A->column != B->column || A->column != Out->column || A->row != Out->row ){printf("Input Matrix can not Sub!\n");return -1;}for (int i = 0; i < A->row; i++)for (int j = 0; j < A->column; j++){Out->matrix_data[A->column*i + j] = A->matrix_data[A->column*i + j] - B->matrix_data[A->column*i + j];}return 0;
}// 功能：矩阵相乘
// 形参：(输入)矩阵A，矩阵A行数row和列数col，矩阵B，矩阵B行数row和列数col
// 返回：A*B (Arow * Bcol)
int MatMul(Matrix_t* A, Matrix_t* B, Matrix_t *Out)
{if (!A || !B || !Out ||A->column != B->row){printf("Input Matrix can not Mul!\n");return -1;}if(Out->row != A->row ||  Out->column != B->column){printf("Output Matrix parameter err!\n");return -1;}int i, j, k;for (i = 0; i < A->row; i++)for (j = 0; j < B->column; j++){Out->matrix_data[B->column*i + j] = 0;for (k = 0; k < A->column; k++)Out->matrix_data[B->column*i + j] = Out->matrix_data[B->column*i + j] + A->matrix_data[A->column*i + k] * B->matrix_data[B->column*k + j];}return 0;
}// 功能：矩阵数乘(实数k乘以矩阵A)
// 形参：(输入)矩阵A首地址指针，矩阵行数row和列数col，实数k
// 返回：kA
int MatMulk(Matrix_t* A, double k, Matrix_t *Out)
{if(!A || !Out ||A->column != Out->column || A->row != Out->row ){printf("Input Matrix can not Mulk!\n");return -1;}for (int i = 0; i < A->row * A->column; i++){Out->matrix_data[i] = A->matrix_data[i] * k;}return 0;
}// 功能：矩阵转置
// 形参：(输入)矩阵A首地址指针A，行数row和列数col
// 返回：A的转置
int MatT(Matrix_t* A, Matrix_t *Out)
{if(!A || !Out ||Out->row != A->column || Out->column != A->row){printf("Input Matrix can not T! (%d %d)->(%d %d)\n", Out->row, Out->column, A->row, A->column);return -1;}for (int i = 0; i < A->row; i++)for (int j = 0; j < A->column; j++)Out->matrix_data[A->row*j + i] = A->matrix_data[A->column*i + j];return 0;
}// 功能：求行列式值
// 形参：(输入)矩阵A首地址指针A，行数row
// 返回：A的行列式值
int MatDet(Matrix_t* A, double *Out)
{if(!A || !Out ||A->column != A->row){printf("Input Matrix can not Det!\n");return -1;}int *list = malloc(A->row * sizeof(int));for (int i = 0; i < A->row; i++)list[i] = i;*Out = det(A->matrix_data, A->row, 0, list, 0.0);free(list);return 0;
}// 功能：矩阵的逆
// 形参：(输入)矩阵A首地址指针A，行数row和列数col
// 返回：A的逆
int MatInv(Matrix_t* A, Matrix_t *Out)
{if(!A || !Out ||A->column != A->row ||Out->row != A->column || Out->column != A->row){printf("Input Matrix can not Inv!\n");return -1;}double det;MatDet(A, &det); //求行列式if (det == 0){printf("Matrix can not Inv!\n");return -1;}MatAdj(A, Out); //求伴随矩阵int len = A->row * A->row;for (int i = 0; i < len; i++)Out->matrix_data[i] /= det;return 0;
}// 功能：求代数余子式
// 形参：(输入)矩阵A首地址指针A，矩阵行数row, 元素a的下标m，n(从0开始)，
// 返回：NxN 矩阵中元素A(mn)的代数余子式
double MatACof(double *A, int row, int m, int n)
{Matrix_t *pstMatrix_cofactor = MatCreate(row - 1, row - 1);int count = 0;int raw_len = row * row;double det = 0;for (int i = 0; i < raw_len; i++)if (i / row != m && i % row != n){*(pstMatrix_cofactor->matrix_data + count++) = *(A + i);}MatDet(pstMatrix_cofactor, &det);if ((m + n) % 2)det = -det;MatRelease(pstMatrix_cofactor);return det;
}// 功能：求伴随矩阵
// 形参：(输入)矩阵A首地址指针A，行数row和列数col
// 返回：A的伴随矩阵
int MatAdj(Matrix_t* A, Matrix_t *Out)
{if(!A || !Out ||Out->row != A->row || Out->column != A->column){printf("Input Matrix can not Adj!\n");return -1;}int len = A->row * A->row;for (int i = 0; i < len; i++){Out->matrix_data[i] = MatACof(A->matrix_data, A->row, i % A->row, i / A->row);}return 0;
}#ifdef USE_CSV_FILE_TEST
// 读取文件行数
static int FileReadRow(const char *filename)
{FILE *f = fopen(filename, "r");int i = 0;char str[4096];while (NULL != fgets(str, 4096, f))++i;printf("数组行数：%d\n", i);return i;
}// 读取文件每行数据数(逗号数+1)
static int FileReadCol(const char *filename)
{FILE *f = fopen(filename, "r");int i = 0;char str[4096];fgets(str, 4096, f);for (int j = 0; j < strlen(str); j++){if (',' == str[j]) i++;}i++;// 数据数=逗号数+1printf("数组列数：%d\n", i);return i;
}// 逗号间隔数据提取
static void GetCommaData(char str_In[4096], double double_Out[1024])
{int str_In_len = strlen(str_In); //printf("str_In_len:%d\n", str_In_len);char str_Data_temp[128];int j = 0;int double_Out_num = 0;for (int i = 0; i < str_In_len; i++){//不是逗号，则是数据，存入临时数组中if (',' != str_In[i]) str_Data_temp[j++] = str_In[i];//是逗号或\n(最后一个数据)，则数据转换为double，保存到输出数组if (',' == str_In[i] || '\n' == str_In[i]) { str_Data_temp[j] = '\0'; j = 0; /*printf("str_Data_temp:%s\n", str_Data_temp); */double_Out[double_Out_num++] = atof(str_Data_temp); memset(str_Data_temp, 0, sizeof(str_Data_temp)); }}
}// 功能：从csv文件读矩阵，保存到指针中
// 形参：(输入)csv文件名，预计行数row和列数col
// 返回：矩阵指针A
int MatRead(char *InCsvFileName, Matrix_t *Out)
{int row = FileReadRow(InCsvFileName); int col = FileReadCol(InCsvFileName);FILE *f = fopen(InCsvFileName, "r");char buffer[4096];while (fgets(buffer, sizeof(buffer), f)){//printf("buffer[%s]\n",buffer);double double_Out[1024] = { 0 };GetCommaData(buffer, double_Out);for (int i = 0; i < col; i++){//printf("double_Out:%lf\n", double_Out[i]);*Out->matrix_data++ = double_Out[i];}}Out->matrix_data = Out->matrix_data - row * col;//指针移回数据开头fclose(f);return 0;
}// 功能：将矩阵A存入csv文件中
// 形参：(输入)保存的csv文件名，矩阵A首地址指针A，行数row和列数col
// 返回：无
int MatWrite(const char *OutCsvFileName, Matrix_t *In)
{FILE *DateFile;double *Ap = In->matrix_data;DateFile = fopen(OutCsvFileName, "w");//追加的方式保存生成的时间戳for (int i = 0; i < In->row*In->column; i++){if ((i+1) % In->column == 0) fprintf(DateFile, "%lf\n", *Ap);//保存到文件，到列数换行else fprintf(DateFile, "%lf,", *Ap);//加逗号Ap++;}fclose(DateFile);return 0;
}
#endif //USE_CSV_FILE_TEST

我测试了一下解非齐次线性方程组（求透视变换矩阵），实现简单，但速度慢且内存需求较大，有兴趣可以尝试去优化。

C++矩阵运算库中，Eigen库是比较高效的

https://eigen.tuxfamily.org/index.php?title=Main_Pagehttps://eigen.tuxfamily.org/index.php?title=Main_Pagehttps://eigen.tuxfamily.org/index.php?title=Main_Page