如何用纯c语言优雅地实现一个矩阵运算库

最新推荐文章于 2025-11-18 06:15:50 发布

原创

最新推荐文章于 2025-11-18 06:15:50 发布 · 置顶 · 5.8k 阅读

150 ·

CC 4.0 BY-SA版权

文章标签：

#c语言矩阵运算库 #矩阵运算 #c语言矩阵求逆 #c语言矩阵库 #编写矩阵库

本文探讨了如何打造一个优雅易用的C语言矩阵运算库，涉及预定义数据类型、对象编程、内存管理、输入验证和编程细节。通过实例展示了矩阵创建、运算和分解等核心功能的实现。

文章目录

1.一个优雅好用的c语言库必须满足哪些条件
2.实现一个矩阵运算库的几点思考
3.完整c程序
参考资料

编程既是技术输出也是艺术创作。鉴赏高手写的程序，往往让人眼前一亮，他们思路、逻辑清晰，所呈现的代码简洁、优雅、高效，令人为之叹服。烂代码则像“屎山"一般让人作呕，软件难以维护最大的原因除了需求模糊的客观因素，最重要的主观因素还是代码写得烂。有生之年，愿能保持对编程的热情，不断提升编程能力，真正体会其乐趣，共勉！

1.一个优雅好用的c语言库必须满足哪些条件

这里给出的软件开发应遵循的一般原则，摘自Les Piegl和Wayne Tiller所著的一本非常经典的书The NURBS Book(Second Edition)。
(1)工具性(toolability):应该使用可用的工具来建立新的应用程序；
(2)可移植性(portability):应用程序应该容易被移植到不同的软件和硬件平台；
(3)可重用性(reusability):程序的编写应该便于重复使用代码段；
(4)可检验性(testability):程序应该简单一致，以便于测试和调试；
(5)可靠性(reliability):对程序运行过程中可能出现的各种错误应该进行合理、一致的处理，以使系统稳定、可靠；
(6)可扩展性(enhanceability):代码必须易于理解，以便可以容易地增加新的功能；
(7)可维护性(fixability):易于找出程序错误的位置；
(8)一致性(consistency):在整个库中，编程习惯应保持一致；
(9)可读性(communicability):程序应该易于阅读和理解;
(10)编程风格(style of programming):代码看起来像书上的数学公式那样以便于读者理解，同时遵循用户友好的编程风格;
(11)易用性(usability):应该使一些非专业的用户也能够方便地使用所开发的库来开发各种更高层次的应用程序;
(12)数值高效性(numerical efficiency):所有函数必须仔细推敲，保证其数值高效性;
(13)基于对象编程(object based programming):避免在函数间传递大量数据，并且使代码易于理解。

2.实现一个矩阵运算库的几点思考

(1)采用预定义的数据类型，避免直接使用编译器定义的数据类型

typedef  unsigned int ERROR_ID;
typedef int INDEX;
typedef short FLAG;
typedef int INTEGER;
typedef double REAL;
typedef char* STRING;
typedef void VOID;

使用预定义的数据类型，有利于程序移植，且可以提高可读性。例如，如果一个系统只支持单精度浮点数，那么只需修改数据类型REAL为float，如果需要更高精度，则只需要改为long double，达到一劳永逸的效果。定义INDEX与INTEGER数据类型是为了在编程时区分索引变量与普通整形变量，同样提高可读性。

(2)基于对象编程，定义矩阵对象

typedef  struct matrix
{
   
   
   INTEGER rows;
   INTEGER columns;
   REAL* p;
}MATRIX;

这里，用一级指针而非二级指针指向矩阵的数据内存地址，有诸多原因，详见博文:为什么我推荐使用一级指针创建二维数组？。
基于对象编程使函数接口更加稳定，不需要经常改动。例如，实现一个函数相乘的函数，你的函数接口可能是这样的：

void matrix_multiplication(double* A, int rowA, int colA, double* B, int rowB, int colB, double* C)

虽然输入参数较多，但一般情况下接口也无需改动。直到有一天，你需要计算维度很大的稀疏矩阵乘法，大量的零元素相乘消耗了大量时间，为了提升计算效率，你需要传入一个参数来表征矩阵是否为大型稀疏矩阵，对稀疏矩阵乘法做一些优化，这时，你的函数接口可能变成这样：

void matrix_multiplication(double* A, int rowA, int colA, int isSparseA, double* B, int rowB, int colB, int isSparseB, double* C)

函数接口的改变，必然迫使你修改所有的函数调用。使用基于对象编程很好地解决了这个问题，接口可以简洁地定义为：

ERROR_ID matrix_multiplication(_IN MATRIX* A, _IN MATRIX* B, _OUT MATRIX* C)

无论未来的业务需求需要增加什么功能，只需要对矩阵对象增加属性，修改函数的实现，接口可以永远不变！这对于矩阵运算库，太重要了，因为矩阵运算太常用了，为基础的函数库，接口一旦改变，调用者苦不堪言。

(3)除了特别编写的内存处理函数(使用栈链表保存、释放动态分配的内存地址)，不允许任何函数直接分配和释放内存

不恰当的分配、使用与释放内存很可能导致内存泄漏、系统崩溃等致命的错误。如果一个函数需动态申请多个内存，那么可能会写出这样啰嗦的程序:

double* x = NULL, * y = NULL, * z = NULL;

x = (double*)malloc(n1 * sizeof(double));
if (x == NULL) return -1;

y = (double*)malloc(n2 * sizeof(double));
if (y == NULL)
{
   
   
	free(x);
	x = NULL;
	return -1;
}

z = (double*)malloc(n3 * sizeof(double));
if (z == NULL)
{
   
   
	free(x);
	x = NULL;
	free(y);
	y = NULL;
	return -1;
}

为了优雅地实现动态内存分配与释放，Les Piegl大神分3步来处理内存申请与释放：
a)在进入一个新的程序时，一个内存堆栈被初始化为空；
b)当需要申请内存时，调用特定的函数来分配所需的内存，并将指向内存的指针存入堆栈中的正确位置；
c)在离开程序时，遍历内存堆栈，释放其中的指针所指向的内存。
程序结构大致如下：

STACKS S;
MATRIX* m = NULL;
INTEGER rows = 3, columns = 4;
ERROR_ID errorID = _ERROR_NO_ERROR;

init_stack(&S);

m = creat_matrix(rows, columns, &errorID, &S);
if (m == NULL) goto EXIT;
//do something
// ...

EXIT:
free_stack(&S);
return errorID;

(4)防御性编程，对输入参数做有效性检查，并返回错误号

例如输入的矩阵行数、列数应该是正整数，指针参数必须非空等等。

(5)注意编程细节的打磨

a)操作符(逗号，等号等)两边必须空一格；
b)逻辑功能相同的程序间不加空行，逻辑功能独立的程序间加空行；
c)条件判断关键字(for if while等)后必须加一空格，起到强调作用，也更清晰；
d)函数内部定义局部变量后，必须空一行后再编写函数主体。

3.完整c程序

本矩阵运算库只包含了矩阵的基本运算，包括创建任意二维/三维矩阵、创建零矩阵及单位矩阵、矩阵加法、矩阵减法、矩阵乘法、矩阵求逆、矩阵转置、矩阵的迹、矩阵LUP分解、解矩阵方程AX=B。
common.h

/*******************************************************************************
*     File Name :                        common.h
*     Library/Module Name :              MatrixComputation
*     Author :                           Marc Pony(marc_pony@163.com)
*     Create Date :                      2021/6/28
*     Abstract Description :             矩阵运算库公用头文件
*******************************************************************************/

#ifndef  __COMMON_H__
#define  __COMMON_H__

/*******************************************************************************
* (1)Debug Switch Section
*******************************************************************************/


/*******************************************************************************
* (2)Include File Section
*******************************************************************************/
#include <math.h>
#include <stdio.h>
#include <malloc.h>
#include <stdlib.h>
#include <time.h>
#include <memory.h>


/*******************************************************************************
* (3)Macro Define Section
*******************************************************************************/
#define _IN
#define _OUT
#define _IN_OUT
#define MAX(x,y) (x) > (y) ? (x) : (y)
#define MIN(x,y) (x) < (y) ? (x) : (y)

#define _CRT_SECURE_NO_WARNINGS
#define PI 3.14159265358979323846
#define POSITIVE_INFINITY 999999999
#define NEGATIVE_INFINITY -999999999
#define _ERROR_NO_ERROR		                                        0X00000000   //无错误
#define _ERROR_FAILED_TO_ALLOCATE_HEAP_MEMORY	                    0X00000001   //分配堆内存失败
#define _ERROR_SVD_EXCEED_MAX_ITERATIONS							0X00000002   //svd超过最大迭代次数
#define _ERROR_MATRIX_ROWS_OR_COLUMNS_NOT_EQUAL	                    0X00000003   //矩阵行数或列数不相等
#define _ERROR_MATRIX_MULTIPLICATION								0X00000004   //矩阵乘法错误(第一个矩阵的列数不等于第二个矩阵行数)
#define _ERROR_MATRIX_MUST_BE_SQUARE								0X00000005   //矩阵必须为方阵
#define _ERROR_MATRIX_NORM_TYPE_INVALID								0X00000006   //矩阵模类型无效
#define _ERROR_MATRIX_EQUATION_HAS_NO_SOLUTIONS						0X00000007   //矩阵方程无解
#define _ERROR_MATRIX_EQUATION_HAS_INFINITY_MANNY_SOLUTIONS	        0X00000008   //矩阵方程有无穷多解
#define _ERROR_QR_DECOMPOSITION_FAILED								0X00000009   //QR分解失败
#define _ERROR_CHOLESKY_DECOMPOSITION_FAILED						0X0000000a   //cholesky分解失败
#define _ERROR_IMPROVED_CHOLESKY_DECOMPOSITION_FAILED				0X0000000b   //improved cholesky分解失败
#define _ERROR_LU_DECOMPOSITION_FAILED								0X0000000c   //LU分解失败
#define _ERROR_CREATE_MATRIX_FAILED									0X0000000d   //创建矩阵失败
#define _ERROR_MATRIX_TRANSPOSE_FAILED								0X0000000e   //矩阵转置失败
#define _ERROR_CREATE_VECTOR_FAILED									0X0000000f   //创建向量失败
#define _ERROR_VECTOR_DIMENSION_NOT_EQUAL 							0X00000010   //向量维数不相同
#define _ERROR_VECTOR_NORM_TYPE_INVALID								0X00000011   //向量模类型无效
#define _ERROR_VECTOR_CROSS_FAILED									0X00000012   //向量叉乘失败
#define _ERROR_INPUT_PARAMETERS_ERROR								0X00010000   //输入参数错误


/*******************************************************************************
* (4)Struct(Data Types) Define Section
*******************************************************************************/
typedef  unsigned int ERROR_ID;
typedef int INDEX;
typedef short FLAG;
typedef int INTEGER;
typedef double REAL;
typedef char* STRING;
typedef void VOID;

typedef  struct matrix
{
   
   
	INTEGER rows;
	INTEGER columns;
	REAL* p;
}MATRIX;

typedef struct matrix_node
{
   
   
	MATRIX* ptr;
	struct matrix_node* next;
} MATRIX_NODE;

typedef struct matrix_element_node
{
   
   
	REAL* ptr;
	struct matrix_element_node* next;
} MATRIX_ELEMENT_NODE;

typedef struct stacks
{
   
   
	MATRIX_NODE* matrixNode;
	MATRIX_ELEMENT_NODE* matrixElementNode;

	// ...
	// 添加其他对象的指针
} STACKS;

/*******************************************************************************
* (5)Prototype Declare Section
*******************************************************************************/
/**********************************************************************************************
Function: init_stack
Description: 初始化栈
Input: 无
Output: 无
Input_Output: 栈指针
Return: 无
Author: Marc Pony(marc_pony@163.com)
***********************************************************************************************/
VOID init_stack(_IN_OUT STACKS* S);


/**********************************************************************************************
Function: free_stack
Description: 释放栈
Input: 栈指针
Output: 无
Input_Output: 无
Return: 无
Author: Marc Pony(marc_pony@163.com)
***********************************************************************************************/
VOID free_stack(_IN STACKS* S);

#endif

matrix.h

/*******************************************************************************
*     File Name :                        matrix.h
*     Library/Module Name :              MatrixComputation
*     Author :                           Marc Pony(marc_pony@163.com)
*     Create Date :                      2021/6/28
*     Abstract Description :             矩阵运算库头文件
*******************************************************************************/

#ifndef  __MATRIX_H__
#define  __MATRIX_H__

/*******************************************************************************
* (1)Debug Switch Section
*******************************************************************************/


/*******************************************************************************
* (2)Include File Section
*******************************************************************************/
#include "common.h"


/*******************************************************************************
* (3)Macro Define Section
*******************************************************************************/


/*******************************************************************************
* (4)Struct(Data Types) Define Section
*******************************************************************************/


/*******************************************************************************
* (5)Prototype Declare Section
*******************************************************************************/

最低0.47元/天解锁文章

31 条评论

Marc Pony 2023.08.16
哪一个

保持理智+相信未来 2023.08.16
对于行列不相等的情况不能求吗?

2004v2004 2023.05.12
矩阵求逆，更改矩阵里一个值就会变化很大，你这个是没有问题的。

2004v2004 2023.05.09
REAL tmp_a[16 * 16] = { 1,2,3,6,5,5,8,7,2...... }; 空间有限省略下我这么赋值可以，为什么用程序 REAL tmp_matrix_d[16 * 16] = {0}; int ik = 0; for (i = 0; i < 16; i++) // 16 { for (j = 0; j < 16; j++) // 16 { sum = 0; for (k = 0; k < 4; k++) // 4 { sum += a[i][k] * Transpose_matrix[k][j]; } tmp_matrix_d[ik] = sum; ik++; } } 打印输出 tmp_matrix_d和tmp_a结果是一样的，调用 init_stack(&S); A = creat_matrix(16, 16, &errorID, &S); // A->p = tmp_a;// 直接赋值就可以， A->p = tmp_matrix_d; // 使用这个得到的就不正确这俩差到那里了？？？？？？ print_matrix(A, "A"); invA = creat_matrix(16, 16, &errorID, &S); errorID = matrix_inverse(A, invA); // 矩阵LUP分解 printf("\n"); printf("矩阵LUP分解 :"); print_matrix(invA, "invA"); 矩阵LUP分解得到的值是不一样的， tmp_a; 得到的是正确的，而tmp_matrix_d 就是错的，是为什么呢