異想家純C語言矩陣運算庫

時間  2020-01-28
標籤 異想 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、C
	printf("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'Y
	double X[15][5] = {
		1, 316, 1536, 874, 981,//第一列要補1
		1, 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'Y
	double *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);

	//保存矩陣到csv
	MatWrite("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++;// 數據數=逗號數+1
	printf("數組列數：%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);
}

 

　　運行結果如圖：spa

相關文章
相關標籤/搜索
每日一句
    每一个你不满意的现在，都有一个你没有努力的曾经。
本站公眾號
   歡迎關注本站公眾號,獲取更多信息
聯繫我們 最近搜索 最新文章 沪ICP备13005482号-10 MyBatis教程 SQL 教程 MySQL教程 Java 教程 Thymeleaf 教程 Hibernate教程 Spring教程 Redis教程