基于C语言的矩阵运算库

Alexon Xu

已于 2023-09-02 23:35:08 修改

阅读量1.9w

点赞数 13

分类专栏：算法文章标签： C语言矩阵运算

于 2018-05-23 22:46:43 首次发布

本文链接：https://blog.csdn.net/shuoyueqishilove/article/details/80427501

版权

算法专栏收录该内容

7 篇文章 1 订阅

订阅专栏

最近本着锻炼自己编程能力的目的，捣鼓了一下矩阵运算，看到网上很多的矩阵运算库都是基于C++的，很少有基于C语言的，于是自己想要写一个纯C语言的矩阵运算库。大神看到了不要喷我，我只是个小白。
个人感觉矩阵运算最重要的是数据结构的安排，我看到有些使用C语言编写的矩阵运算库里面的矩阵元素是使用以为数组存储的，这样的好处是确定了矩阵的行和列，直接将矩阵结构中的指针指向对应的数组即可；其形式如下：

typedef struct 
{
  int row;
  int column;
  float* data;
}Matrix_t;

这种一维数组有个缺点就是，运算的时候下标不好操作，需要程序员编程功底较好，对于我这种半吊子，还是使用简单点儿的二维数组吧，下标和各个元素一一对应，方便操作。这里使用动态内存分配的方式来管理矩阵的数据存储的分配和释放。

typedef unsigned int uint16_t;
typedef struct
{
	uint16_t row;
	uint16_t column;
	float **data;
}Matrix_t;

首先创建矩阵，创建矩阵时只需要矩阵的行数和列数，矩阵创建函数会自动根据矩阵的行和列分配一片动态内存给矩阵。为了方便起见，将所有的矩阵元素全部置为0；

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

矩阵元素置0函数很简单，输入指向对应矩阵的指针，将每个元素设置为0；

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

由于创建矩阵过程中使用malloc函数分配了内存，在矩阵完成运算之后，一定要将矩阵的内存用free释放掉，以防止出现大量内存碎片，使程序崩溃。

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

需要注意的是释放的顺序必须是后申请的内存先释放，否则释放不彻底。
对于矩阵的元素设定需要设计一个函数，该函数将待设置矩阵和一维数组作为输入，将一维数组中的元素一次拷贝到矩阵的二维数组中。之前试过使用指向二维数组的指针，觉得不是很好，还是用一维数组吧。

/*
set datas to the matrix
*/
void set_mat_data(Matrix_t* mat,const float *data)
{
	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];
		}
	}
}

矩阵的加减实现方法区别不大，只需对应下标的元素相互加减即可，以下是矩阵相加的情形：

/*
mat=mat1+mat2
*/
Matrix_t add_mat(const Matrix_t* mat1, const Matrix_t* mat2)
{
	
	if (mat1->row != mat2->row)
	{
		printf("error, in add_mat: mat1->row != mat2->row\n");
		exit(1);
	}
	if (mat1->column != mat2->column)
	{
		printf("error, in add_mat: mat1->column != mat2->column\n");
		exit(1);
	}
	Matrix_t mat;
	uint16_t i, j;
	mat = create_mat(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;
}

矩阵转置时将矩阵的行和列对换；

/*
transpose the matrix, mat=mat'
*/
Matrix_t transpose_mat(const Matrix_t* mat)//mat'
{
	Matrix_t mat_T;
	mat_T = create_mat(mat->row, mat->column);
	uint16_t i, j;
	for (i = 0; i < mat->row; i++)
	{
		for (j = 0; j < mat->column; j++)
		{
			mat_T.data[i][j] = mat->data[j][i];
		}
	}
	return mat_T;
}

两个矩阵相乘，根据线性代数中的方法，第一个矩阵的第一行乘以第二个矩阵的每一列，得到结果矩阵的第一行的值。。。。。实现起来并不困难。

*
mat=mat1*mat2
*/
Matrix_t mult_mat(const Matrix_t *mat1, const Matrix_t* mat2)
{
	Matrix_t mat;
	if (mat1->column != mat2->row)
	{
		printf("error,In mult_mat: mat1->column != mat2->row\n");
		exit(1);
	}	
	else
	{
		mat = create_mat(mat1->row, mat2->column);
		clear_mat(&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;
}

举证行列式的值的计算是比较复杂的，在线性代数中，一般采用的方法是：将矩阵通过行变换（列主元消去法）变换为上三角矩阵，需要特别注意的是，矩阵两行互换时，行列式的值变号，编程时需要考虑。

/*
get matrix's derterminent value
*/
float det_mat(Matrix_t *m)
{
	uint16_t i, j, n, max_row;
	char swap_f;
	float max, k;
	float det=1;
	swap_f = 0;
	if (m->column != m->row)
	{
		printf("error:In det_mat (m->column != m->row)\n");
		exit(1);
	}
	Matrix_t mat = copy_mat(m);
	for (i = 0; i < mat.row-1; i++)
	{
		max = abs(mat.data[i][i]);
		max_row = i;
		for (j = i + 1; j < mat.row; j++)
		{
			if (max < abs(mat.data[j][i]))
			{
				max = abs(mat.data[j][i]);
				max_row = j;
			}
		}
		if (i != max_row)
		{
			swap_row_mat(&mat, i, max_row);
			swap_f++;
		}
		for (j = i + 1; j < mat.row; j++)
		{
			k = -mat.data[j][i]/mat.data[i][i];
			for (n= 0; n < mat.column; n++)
			{
				mat.data[j][n] = mat.data[i][n] * k + mat.data[j][n];
			}
		}
	}
	if (swap_f%2==1)swap_f = -1;
	else swap_f = 1;
	det = 1;
	for (i = 0; i < mat.column; i++)
		det *= mat.data[i][i];
	det *= swap_f;
	//show_mat("in det_mat: mat_upper",&mat);
	free_mat(&mat);
	return det;
}

矩阵的逆矩阵，这是难点所在，其实编程上容易实现的还是Gauss-Jordan列主元消去法，即：将 A|I 通过行变换为 I| $A^{-1}$

/*
get inverse matrix
use main column element of Gauss-Jordan algrithm: A|I  --->  I|A^(-1)
*/
Matrix_t inverse_mat(Matrix_t* m)
{
	if (det_mat(m) == 0)
	{
		printf("error: In inverse_mat(det_mat(mat) == 0)\n");
		exit(1);
	}

	Matrix_t mat = copy_mat(m);
	Matrix_t inv_mat = eye(m->row);
	
	int i, j, n,max_row;
	int swap_f = 0;
	float max,k;
	for (i = 0; i < mat.row - 1; i++)
	{
		max = abs(mat.data[i][i]);
		max_row = i;
		for (j = i + 1; j < mat.row; j++)
		{
			if (max < abs(mat.data[j][i]))
			{
				max = abs(mat.data[j][i]);
				max_row = j;
			}
		}
		if (i != max_row)
		{
			swap_row_mat(&mat, i, max_row);
			swap_row_mat(&inv_mat, i, max_row);
			swap_f++;
		}
		for (j = i + 1; j < mat.row; j++)
		{
			k = -mat.data[j][i] / mat.data[i][i];
			for (n = 0; n < mat.column; n++)
			{
				mat.data[j][n] = mat.data[i][n] * k + mat.data[j][n];
				inv_mat.data[j][n] = inv_mat.data[i][n] * k + inv_mat.data[j][n];
			}
		}
	}

	for (i = 0; i < mat.row; i++)
	{
		k = 1/mat.data[i][i];
		scale_row_mat(&mat,i, k);
		scale_row_mat(&inv_mat, i, k);
	}
	for (i = mat.row-1; i>0; i--)
	{
		for (j = i - 1; j >=0; j--)
		{
			k = -mat.data[j][i] / mat.data[i][i];
			for (n = 0; n < mat.column; n++)
			{
				mat.data[j][n] = mat.data[j][n] + k*mat.data[i][n];
				inv_mat.data[j][n] = inv_mat.data[j][n] + k*inv_mat.data[i][n];
			}
		}
		
	}
	//show_mat("in inverse_mat: mat", &mat);
	//show_mat("in inverse_mat: inv_mat", &inv_mat);
	free_mat(&mat);
	return inv_mat;
}

测试一下：

#include<iostream>
#include"myMatrix.h"

using namespace std;

int main()
{
	Matrix_t mat=create_mat(4,4);
	float data[16] = { 1, 2, 3, 4,
						2, 1, 4, 5,
						5, 4, 3, 5,
						7, 6, 5, 4 };
	set_mat_data(&mat, data);
	show_mat("mat", &mat);
	Matrix_t inv_mat=inverse_mat(&mat);
	show_mat("inv_mat", &inv_mat);

	system("pause");
	return 0;
}

这里写图片描述
以上是实现该矩阵运算库的核心代码了，感兴趣的可以下载本代码：
https://download.csdn.net/download/shuoyueqishilove/10433912

Alexon Xu

关注

13
点赞
踩
133

收藏

觉得还不错? 一键收藏
12
评论
基于C语言的矩阵运算库

最近本着锻炼自己编程能力的目的，捣鼓了一下矩阵运算，看到网上很多的矩阵运算库都是基于C++的，很少有基于C语言的，于是自己想要写一个纯C语言的矩阵运算库。大神看到了不要喷我，我只是个小白。个人感觉矩阵运算最重要的是数据结构的安排，我看到有些使用C语言编写的矩阵运算库里面的矩阵元素是使用以为数组存储的，这样的好处是确定了矩阵的行和列，直接将矩阵结构中的指针指向对应的数组即可；其形式如下：...
复制链接

扫一扫

专栏目录