异想家纯C语言矩阵运算库

异想家纯C语言矩阵运算库

　　Sandeepin最近做的项目中需要在嵌入式芯片里跑一些算法，而这些单片机性能不上不下，它能跑些简单的程序，但又还没到上Linux系统的地步。所以只好用C语言写一些在高级语言里一个函数就解决的算法了，由于算法需要运用矩阵运算，自己就先用纯C语言写了个简单的矩阵运算库。

　　代码只实现了矩阵最基本的运算，包括矩阵的加、减、乘、数乘、转置、行列式、逆矩阵、代数余子式、伴随矩阵等运算。此外增加了一些实用函数，如显示矩阵、从csv文件读取保存矩阵等函数。具体的例子在主函数中体现，其中还用自己这个矩阵运算库做了一个简单的应用，利用公式β=(X'X)^(-1)X'Y来进行多元线性回归系数计算，大家可以参考参考，欢迎批评。

　　JfzMatLib.c文件代码：

#include "JfzMatLib.h"int main(int argc, char** argv)
{//矩阵的基本运算：加、减、乘、数乘、转置、行列式、逆矩阵、代数余子式、伴随矩阵//初始实验矩阵double A[] = { -3, 2, -5, -1, 0, -2, 3, -4, 1 };double B[] = { 1, 4, 7, 3, 0, 5, -1, 9, 11 };double C[] = { 1, 2, 3, 4, 5, 6, 7, 8, 9 };//计算结果矩阵double *Add = (double*)malloc(sizeof(double) * 9);double *Sub = (double*)malloc(sizeof(double) * 9);double *Mul = (double*)malloc(sizeof(double) * 9);double *kMat = (double*)malloc(sizeof(double) * 9);double *CT = (double*)malloc(sizeof(double) * 9);double *AInv = (double*)malloc(sizeof(double) * 9);double *Adj = (double*)malloc(sizeof(double) * 9);//显示矩阵A、B、Cprintf("A：\n"); MatShow(A, 3, 3);printf("B：\n"); MatShow(B, 3, 3);printf("C：\n"); MatShow(C, 3, 3);//矩阵相加printf("A+B：\n");Add = MatAdd(A, B, 3, 3); MatShow(Add, 3, 3);//矩阵相减printf("A-B：\n");Sub = MatSub(A, B, 3, 3); MatShow(Sub, 3, 3);//矩阵相乘printf("A*B：\n");Mul = MatMul(A, 3, 3, B, 3, 3); MatShow(Mul, 3, 3);//矩阵数乘printf("2*C：\n");kMat = MatMulk(C, 3, 3, 2); MatShow(kMat, 3, 3);//矩阵转置printf("C的转置：\n");CT = MatT(C, 3, 3); MatShow(CT, 3, 3);//矩阵行列式值printf("B的行列式值：\n");printf("%16lf\n", MatDet(B, 3));printf("C的行列式值：\n");printf("%16lf\n", MatDet(C, 3));//矩阵的逆printf("A的逆：\n");AInv = MatInv(A, 3, 3); MatShow(AInv, 3, 3);printf("C的逆：\n");MatInv(C, 3, 3);//矩阵代数余子式printf("C的(0,0)元素的代数余子式：\n");printf("%16lf\n", MatACof(C, 3, 0, 0));//矩阵伴随矩阵printf("A的伴随矩阵：\n");Adj = MatAdj(A, 3, 3); MatShow(Adj, 3, 3);//蒋方正矩阵库应用：多元线性回归//多元线性回归公式：β=(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};//多元线性回归公式：β=(X'X)^(-1)X'Ydouble *XT = (double*)malloc(sizeof(double) * 5 * 15);double *XTX = (double*)malloc(sizeof(double) * 5 * 5);double *InvXTX = (double*)malloc(sizeof(double) * 5 * 5);double *InvXTXXT = (double*)malloc(sizeof(double) * 5 * 15);double *InvXTXXTY = (double*)malloc(sizeof(double) * 5 * 1);XT = MatT((double*)X, 15, 5);XTX = MatMul(XT, 5, 15, (double*)X, 15, 5);InvXTX = MatInv(XTX, 5, 5);InvXTXXT = MatMul(InvXTX, 5, 5, XT, 5, 15);InvXTXXTY = MatMul(InvXTXXT, 5, 15, (double*)Y, 15, 1);printf("多元线性回归β系数：\n");MatShow(InvXTXXTY, 5, 1);//保存矩阵到csvMatWrite("XTX.csv", XTX, 5, 5);MatWrite("InvXTXXTY.csv", InvXTXXTY, 5, 1);MatWrite("Fuck.csv", A, 3, 3);MatWrite("Fuck2.csv", A, 9, 1);//从csv读取矩阵double *Fuck = (double*)malloc(sizeof(double) * 3 * 3);Fuck = MatRead("Fuck.csv");MatShow(Fuck, 3, 3);double *Fuck2 = (double*)malloc(sizeof(double) * 9 * 1);Fuck2 = MatRead("Fuck2.csv");MatShow(Fuck2, 9, 1);double *InvXTXXTYread = (double*)malloc(sizeof(double) * 5 * 1);InvXTXXTYread = MatRead("InvXTXXTY.csv");MatShow(InvXTXXTYread, 5, 1);double *XTXread = (double*)malloc(sizeof(double) * 5 * 5);XTXread = MatRead("XTX.csv");MatShow(XTXread, 5, 5);system("pause");
}

　　JfzMatLib.h头文件代码：

 

#include <stdio.h>
#include <stdlib.h>
#include <string.h>void    MatShow(double* Mat, int row, int col);
double* MatAdd(double* A, double* B, int row, int col);
double* MatSub(double* A, double* B, int row, int col);
double* MatMul(double* A, int Arow, int Acol, double* B, int Brow, int Bcol);
double* MatMulk(double *A, int row, int col, double k);
double* MatT(double* A, int row, int col);
double  MatDet(double *A, int row);
double* MatInv(double *A, int row, int col);
double  MatACof(double *A, int row, int m, int n);
double* MatAdj(double *A, int row, int col);
double *MatRead(char *csvFileName, int row, int col);
void MatWrite(char *A, int row, int col);// (det用)功能：求逆序对的个数
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
int sgn(int order)
{return order % 2 ? -1 : 1;
}// (det用)功能：交换两整数a、b的值
void swap(int *a, int *b)
{int m;m = *a;*a = *b;*b = m;
}// 功能：求矩阵行列式的核心函数
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;
}// 功能：矩阵显示
// 形参：(输入)矩阵首地址指针Mat，矩阵行数row和列数col。
// 返回：无
void MatShow(double* Mat, int row, int col)
{for (int i = 0; i < row*col; i++){printf("%16lf ", Mat[i]);if (0 == (i + 1) % col) printf("\n");}
}// 功能：矩阵相加
// 形参：(输入)矩阵A首地址指针A，矩阵B首地址指针B，矩阵A(也是矩阵B)行数row和列数col
// 返回：A+B
double* MatAdd(double* A, double* B, int row, int col)
{double *Out = (double*)malloc(sizeof(double) * row * col);for (int i = 0; i < row; i++)for (int j = 0; j < col; j++)Out[col*i + j] = A[col*i + j] + B[col*i + j];return Out;
}// 功能：矩阵相减
// 形参：(输入)矩阵A，矩阵B，矩阵A(也是矩阵B)行数row和列数col
// 返回：A-B
double* MatSub(double* A, double* B, int row, int col)
{double *Out = (double*)malloc(sizeof(double) * row * col);for (int i = 0; i < row; i++)for (int j = 0; j < col; j++)Out[col*i + j] = A[col*i + j] - B[col*i + j];return Out;
}// 功能：矩阵相乘
// 形参：(输入)矩阵A，矩阵A行数row和列数col，矩阵B，矩阵B行数row和列数col
// 返回：A*B
double* MatMul(double* A, int Arow, int Acol, double* B, int Brow, int Bcol)
{double *Out = (double*)malloc(sizeof(double) * Arow * Bcol);if (Acol != Brow){printf("        Shit!矩阵不可乘!\n");return NULL;}if (Acol == Brow){int i, j, k;for (i = 0; i < Arow; i++)for (j = 0; j < Bcol; j++){Out[Bcol*i + j] = 0;for (k = 0; k < Acol; k++)Out[Bcol*i + j] = Out[Bcol*i + j] + A[Acol*i + k] * B[Bcol*k + j];}return Out;}
}// 功能：矩阵数乘(实数k乘以矩阵A)
// 形参：(输入)矩阵A首地址指针，矩阵行数row和列数col，实数k
// 返回：kA
double* MatMulk(double *A, int row, int col, double k)
{double *Out = (double*)malloc(sizeof(double) * row * col);for (int i = 0; i < row * col; i++){*Out = *A * k;Out++;A++;}Out = Out - row * col;return Out;
}// 功能：矩阵转置
// 形参：(输入)矩阵A首地址指针A，行数row和列数col
// 返回：A的转置
double* MatT(double* A, int row, int col)
{double *Out = (double*)malloc(sizeof(double) * row * col);for (int i = 0; i < row; i++)for (int j = 0; j < col; j++)Out[row*j + i] = A[col*i + j];return Out;
}// 功能：求行列式值
// 形参：(输入)矩阵A首地址指针A，行数row
// 返回：A的行列式值
double MatDet(double *A, int row)
{int *list = (int*)malloc(sizeof(int) * row);for (int i = 0; i < row; i++)list[i] = i;double Out = det(A, row, 0, list, 0.0);free(list);return Out;
}// 功能：矩阵的逆
// 形参：(输入)矩阵A首地址指针A，行数row和列数col
// 返回：A的逆
double *MatInv(double *A, int row, int col)
{double *Out = (double*)malloc(sizeof(double) * row * col);double det = MatDet(A, row); //求行列式if (det == 0){printf("        Fuck!矩阵不可逆!\n");return NULL;}if (det != 0){Out = MatAdj(A, row, col); //求伴随矩阵int len = row * row;for (int i = 0; i < len; i++)*(Out + i) /= det;return Out;}
}// 功能：求代数余子式
// 形参：(输入)矩阵A首地址指针A，矩阵行数row, 元素a的下标m，n(从0开始)，
// 返回：NxN 矩阵中元素A(mn)的代数余子式
double MatACof(double *A, int row, int m, int n)
{int len = (row - 1) * (row - 1);double *cofactor = (double*)malloc(sizeof(double) * len);int count = 0;int raw_len = row * row;for (int i = 0; i < raw_len; i++)if (i / row != m && i % row != n)*(cofactor + count++) = *(A + i);double ret = MatDet(cofactor, row - 1);if ((m + n) % 2)ret = -ret;free(cofactor);return ret;
}// 功能：求伴随矩阵
// 形参：(输入)矩阵A首地址指针A，行数row和列数col
// 返回：A的伴随矩阵
double *MatAdj(double *A, int row, int col)
{double *Out = (double*)malloc(sizeof(double) * row * col);int len = row * row;int count = 0;for (int i = 0; i < len; i++){*(Out + count++) = MatACof(A, row, i % row, i / row);}return Out;
}// 读取文件行数
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)
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;
}// 逗号间隔数据提取
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
double *MatRead(char *csvFileName)
{int row = FileReadRow(csvFileName); int col = FileReadCol(csvFileName);double *Out = (double*)malloc(sizeof(double) * row * col);FILE *f = fopen(csvFileName, "r");char buffer[4096];while (fgets(buffer, sizeof(buffer), f)){//printf("buffer[%s]\n",buffer);double double_Out[128] = { 0 };GetCommaData(buffer, double_Out);for (int i = 0; i < col; i++){//printf("double_Out:%lf\n", double_Out[i]);*Out = double_Out[i];Out++;}}Out = Out - row * col;//指针移回数据开头fclose(f);return Out;
}// 功能：将矩阵A存入csv文件中
// 形参：(输入)保存的csv文件名，矩阵A首地址指针A，行数row和列数col
// 返回：无
void MatWrite(const char *csvFileName, double *A, int row, int col)
{FILE *DateFile;double *Ap = A;DateFile = fopen(csvFileName, "w");//追加的方式保存生成的时间戳for (int i = 0; i < row*col; i++){if ((i+1) % col == 0) fprintf(DateFile, "%lf\n", *Ap);//保存到文件，到列数换行else fprintf(DateFile, "%lf,", *Ap);//加逗号Ap++;}fclose(DateFile);
}

 

　　运行结果如图：

转载于:.html