C 矩阵简单运算(+ - * 转置 求逆 行列式 炒鸡完整版)

前一段时间用C(++)编程，由于需要进行矩阵运算，当时见识浅薄，还不知道有 eigen 这样强大的矩阵库，所以在网上拼拼凑凑，捣鼓出了一套简单的基于二维数组的 c 矩阵运算功能，现在完善了一下，希望对大家有所帮助：

1. myMatrix.h 头文件

#include <cstdio>

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

Matrix CreateMat(uint16_t row, uint16_t 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 = (DType **)malloc(row * sizeof(DType *));//先指针的指针
		if (mat.data == NULL)
		{
			printf("error, in CreateMat: mat.data==NULL");
			exit(1);
		}
		uint16_t i;
		for (i = 0; i < row; i++)
		{
			*(mat.data + i) = (DType *)malloc(column * sizeof(DType));//再分配每行的指针
			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
	uint16_t 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
	uint16_t i;
	for (i = 0; i < mat->row; i++)
		free(mat->data[i]);/*释放行*/
	free(mat->data);/*释放头指针*/
}
void SetMatData(Matrix* mat, DType(*data))//用指针变量作为形参
{   //set datas to the matrix
	uint16_t 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;
	uint16_t 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;
	uint16_t 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 TransposeMat(const Matrix mat)
{  //transpose the matrix, mat=mat'转置
	Matrix mat_T;
	mat_T = CreateMat(mat.column, mat.row);
	uint16_t 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 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);
		uint16_t i, j;
		for (i = 0; i < mat1.row; i++)
		{
			for (j = 0; j < mat2.column; j++)
			{
				uint16_t m;
				for (m = 0; m < mat1.column; m++)
				{
					mat.data[i][j] += mat1.data[i][m] * mat2.data[m][j];
				}
			}
		}
	}
	return mat;
}
Matrix CopyMat(const Matrix mat)
{
	Matrix mat_C;
	mat_C = CreateMat(mat.row, mat.column);
	uint16_t 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;
}
DType DetMat(const Matrix m)
{  //get matrix's derterminent value行列式
	char swap_f;
	DType 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 (uint16_t i = 0; i < m.column; i++)//这是运用LU分解求矩阵行列式
	{
		det = det * mat.data[i][i];
		for (uint16_t j = i + 1; j < m.column; j++)//更新一次矩阵
		{
			mat.data[j][i] = mat.data[j][i] / mat.data[i][i];
			for (uint16_t 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 uint16_t a)
{  //制造一个单位阵
	Matrix mat = CreateMat(a, a);
	for (uint16_t i = 0; i < a; i++)
	{
		for (uint16_t 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);
	}
	uint16_t a = m.column;
	Matrix mat = CopyMat(m);
	Matrix inv_mat = Eye(a);
	//step1:从上向下
	for (uint16_t i = 0; i < a; i++)
	{
		for (uint16_t j = i + 1; j < a; j++)//更新一次矩阵
		{
			DType p = mat.data[j][i] / mat.data[i][i];
			for (uint16_t 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--)//更新一次矩阵
		{
			DType p = mat.data[j][i] / mat.data[i][i];
			for (uint16_t 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 (uint16_t i = 0; i < a; i++)
	{
		DType p = mat.data[i][i];
		for (uint16_t 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 (uint16_t i = 0; i < mat.row; i++)
	{
		for (uint16_t 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

2. 使用示例

#include <cstdio>
#include <cstdlib>

#include "myMatrix.h"

void main()
{
	Matrix A=CreateMat(3,3);
	double A_data[3][3]={1,2,8,4,5,6,7,8,9};
	SetMatData(&A,*A_data);
	Matrix B=CreateMat(3,3);
	double B_data[3][3]={1,2,3,4,5,6,7,8,9};
	SetMatData(&B,*B_data);
	printf("A：\n");
	ShowMat(A);
	printf("B：\n");
	ShowMat(B);
	printf("A+B：\n");
	ShowMat(A+B);
	printf("A-B：\n");
	ShowMat(A-B);
	printf("A*B：\n");
	ShowMat(A*B);
	printf("A的行列式：\n");
	printf("%lf\n",DetMat(A));
	printf("A的转置：\n");
	ShowMat(TransposeMat(A));
	printf("A的逆：\n");
	ShowMat(~A);
	system("pause");
}

结果：

A：
     1.000        2.000        8.000
     4.000        5.000        6.000
     7.000        8.000        9.000
B：
     1.000        2.000        3.000
     4.000        5.000        6.000
     7.000        8.000        9.000
A+B：
     2.000        4.000       11.000
     8.000       10.000       12.000
    14.000       16.000       18.000
A-B：
     0.000        0.000        5.000
     0.000        0.000        0.000
     0.000        0.000        0.000
A*B：
    65.000       76.000       87.000
    66.000       81.000       96.000
   102.000      126.000      150.000
A的行列式：
-15.000000
A的转置：
     1.000        4.000        7.000
     2.000        5.000        8.000
     8.000        6.000        9.000
A的逆：
     0.200       -3.067        1.867
    -0.400        3.133       -1.733
     0.200       -0.400        0.200
请按任意键继续. . .

特别感谢此篇博客的博主：
https://blog.csdn.net/shuoyueqishilove/article/details/80427501
积分多的朋友可给点积分：
https://download.csdn.net/download/gou_hailong/11467477