c语言实现 基础矩阵运算

参考链接：

https://blog.csdn.net/qq_41649694/article/details/81432663

https://blog.csdn.net/Gou_Hailong/article/details/98329113?utm_medium=distribute.wap_relevant.none-task-blog-BlogCommendFromMachineLearnPai2-3.nonecase&depth_1-utm_source=distribute.wap_relevant.none-task-blog-BlogCommendFromMachineLearnPai2-3.nonecase

math.h

#include <stdio.h>
#include <stdbool.h> 
#include <math.h>
#include <stdlib.h>

//数据类型定义 
typedef struct
{
	int row;
	int column;
	double **data;
}Matrix;

Matrix CreateMat(int row, int column)
{  //create a matrix with 0
	Matrix mat;
	if (row <= 0 || column <= 0)
	{
		printf("error, in CreateMat: row <= 0||column<=0\n");
		exit(1);
	}
	if (row > 0 && column > 0)
	{
		mat.row = row;
		mat.column = column;
		mat.data = (double **)malloc(row * sizeof(double *));//先指针的指针
		if (mat.data == NULL)
		{
			printf("error, in CreateMat: mat.data==NULL");
			exit(1);
		}
		int i;
		for (i = 0; i < row; i++)
		{
			*(mat.data + i) = (double *)malloc(column * sizeof(double));//再分配每行的指针
			if (mat.data[i] == NULL)
			{
				printf("error, in CreateMat: mat.data==NULL");
				exit(1);
			}
		}
	}
	return mat;
}

void ClearMat(Matrix* mat)
{  //set all datas in matrix to zero
	int i, j;
	for (i = 0; i < mat->row; i++)
	{
		for (j = 0; j < mat->column; j++)
		{
			mat->data[i][j] = 0;
		}
	}
	return;
}

void FreeMat(Matrix *mat)
{  //free a matrix
	int i;
	for (i = 0; i < mat->row; i++)
		free(mat->data[i]);/*释放行*/
	free(mat->data);/*释放头指针*/
}

void SetMatData(Matrix* mat, double(*data))//用指针变量作为形参
{   //set datas to the matrix
	int i, j;
	for (i = 0; i < mat->row; i++)
	{
		for (j = 0; j < mat->column; j++)
		{
			mat->data[i][j] = data[i*mat->column + j];
		}
	}
}

Matrix AddMat(const Matrix mat1, const Matrix mat2)
{  //mat=mat1+mat2矩阵相加
	if (mat1.row != mat2.row)
	{
		printf("error, in AddMat: mat1->row != mat2->row\n");
		exit(1);
	}
	if (mat1.column != mat2.column)
	{
		printf("error, in AddMat: mat1->column != mat2->column\n");
		exit(1);
	}
	Matrix mat;
	int i, j;
	mat = CreateMat(mat1.row, mat1.column);
	for (i = 0; i < mat1.row; i++)
	{
		for (j = 0; j < mat1.column; j++)
			mat.data[i][j] = mat1.data[i][j] + mat2.data[i][j];
	}
	return mat;
}

Matrix SubMat(const Matrix mat1, const Matrix mat2)
{  //mat=mat1-mat2矩阵相减
	if (mat1.row != mat2.row)
	{
		printf("error, in AddMat: mat1->row != mat2->row\n");
		exit(1);
	}
	if (mat1.column != mat2.column)
	{
		printf("error, in AddMat: mat1->column != mat2->column\n");
		exit(1);
	}
	Matrix mat;
	int i, j;
	mat = CreateMat(mat1.row, mat1.column);
	for (i = 0; i < mat1.row; i++)
	{
		for (j = 0; j < mat1.column; j++)
			mat.data[i][j] = mat1.data[i][j] - mat2.data[i][j];
	}
	return mat;
}

Matrix MultMat(const Matrix mat1, const Matrix mat2)
{  //mat=mat1*mat2矩阵相乘
	Matrix mat;
	if (mat1.column != mat2.row)
	{
		printf("error,In MultMat: mat1.column != mat2.row\n");
		exit(1);
	}
	else
	{
		mat = CreateMat(mat1.row, mat2.column);
		ClearMat(&mat);
		int i, j;
		for (i = 0; i < mat1.row; i++)
		{
			for (j = 0; j < mat2.column; j++)
			{
				int m;
				for (m = 0; m < mat1.column; m++)
				{
					mat.data[i][j] += mat1.data[i][m] * mat2.data[m][j];
				}
			}
		}
	}
	return mat;
}

Matrix TransposeMat(const Matrix mat)
{  //transpose the matrix, mat=mat'转置
	Matrix mat_T;
	mat_T = CreateMat(mat.column, mat.row);
	int i, j;
	for (i = 0; i < mat.row; i++)
	{
		for (j = 0; j < mat.column; j++)
		{
			mat_T.data[j][i] = mat.data[i][j];
		}
	}
	return mat_T;
}

Matrix CopyMat(const Matrix mat)
{
	Matrix mat_C;
	mat_C = CreateMat(mat.row, mat.column);
	int i, j;
	for (i = 0; i < mat.row; i++)
	{
		for (j = 0; j < mat.column; j++)
		{
			mat_C.data[i][j] = mat.data[i][j];
		}
	}
	return mat_C;
}

double DetMat(const Matrix m)
{  //get matrix's derterminent value行列式
	char swap_f;
	double det = 1;
	swap_f = 0;
	if (m.column != m.row)
	{
		printf("error:In DetMat (m.column != m.row)\n");
		exit(1);
	}
	Matrix mat = CopyMat(m);
	for (int i = 0; i < m.column; i++)//这是运用LU分解求矩阵行列式
	{
		det = det * mat.data[i][i];
		for (int j = i + 1; j < m.column; j++)//更新一次矩阵
		{
			mat.data[j][i] = mat.data[j][i] / mat.data[i][i];
			for (int k = i + 1; k < m.column; k++)//更新第j行
			{
				mat.data[j][k] = mat.data[j][k] - mat.data[i][k] * mat.data[j][i];
			}
		}
	}
	//最后的mat矩阵保存着LU分解的结果
	FreeMat(&mat);
	return det;
}

Matrix Eye(const int a)
{  //制造一个单位阵
	Matrix mat = CreateMat(a, a);
	for (int i = 0; i < a; i++)
	{
		for (int j = 0; j < a; j++)
		{
			if (i == j)
				mat.data[i][j] = 1;
			else
				mat.data[i][j] = 0;
		}
	}
	return mat;
}

/*
get inverse matrix
use main column element of Gauss-Jordan algrithm: A|I  --->  I|A^(-1)
*/
Matrix InverseMat(const Matrix m)
{  //求逆
	if (DetMat(m) == 0)
	{
		printf("error: In InverseMat(DetMat(mat) == 0)\n");
		exit(1);
	}
	int a = m.column;
	Matrix mat = CopyMat(m);
	Matrix inv_mat = Eye(a);
	//step1:从上向下
	for (int i = 0; i < a; i++)
	{
		for (int j = i + 1; j < a; j++)//更新一次矩阵
		{
			double p = mat.data[j][i] / mat.data[i][i];
			for (int k = 0; k < m.column; k++)//更新第j行
			{
				mat.data[j][k] = mat.data[j][k] - mat.data[i][k] * p;
				inv_mat.data[j][k] = inv_mat.data[j][k] - inv_mat.data[i][k] * p;
			}
		}
	}
	//step2:从下向上
	for (int i = a - 1; i >= 0; i--)
	{
		for (int j = i - 1; j >= 0; j--)//更新一次矩阵
		{
			double p = mat.data[j][i] / mat.data[i][i];
			for (int k = 0; k < m.column; k++)//更新第j行
			{
				mat.data[j][k] = mat.data[j][k] - mat.data[i][k] * p;
				inv_mat.data[j][k] = inv_mat.data[j][k] - inv_mat.data[i][k] * p;
			}
		}
	}
	//step3:单位化
	for (int i = 0; i < a; i++)
	{
		double p = mat.data[i][i];
		for (int j = 0; j < a; j++)
		{
			inv_mat.data[i][j] = inv_mat.data[i][j] / p;
			mat.data[i][j] = mat.data[i][j] / p;
		}
	}
	FreeMat(&mat);
	return inv_mat;
}

void ShowMat(const Matrix mat)
{   //打印出矩阵
	for (int i = 0; i < mat.row; i++)
	{
		for (int j = 0; j < mat.column; j++)
			printf("%10.3f   ", mat.data[i][j]);
		printf("\n");
	}
	return;
}


运算符重载
//Matrix operator+ (const Matrix mat1, const Matrix mat2) { return AddMat(mat1, mat2); }
//Matrix operator- (const Matrix mat1, const Matrix mat2) { return SubMat(mat1, mat2); }
//Matrix operator* (const Matrix mat1, const Matrix mat2) { return MultMat(mat1, mat2); }
//Matrix operator~ (const Matrix mat) { return InverseMat(mat); }//求逆
#pragma once