这里给出的软件开发应遵循的一般原则，摘自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 : 2023/6/28* Abstract Description : 矩阵运算库公用头文件*******************************************************************************/#ifndef __COMMON_H__#define __COMMON_H__/******************************************************************************** (1)Debug Switch Section*******************************************************************************//******************************************************************************** (2)Include File Section*******************************************************************************/#include #include #include #include #include #include /******************************************************************************** (3)Macro Define Section*******************************************************************************/#define _IN#define _OUT#define _IN_OUT#define MAX(x,y) (x) > (y) ? (x) : (y)#define MIN(x,y) (x) matrixNode != NULL){matrixNode = S->matrixNode;S->matrixNode = matrixNode->next;free(matrixNode->ptr);matrixNode->ptr = NULL;free(matrixNode);matrixNode = NULL;}while (S->matrixElementNode != NULL){matrixElementNode = S->matrixElementNode;S->matrixElementNode = matrixElementNode->next;free(matrixElementNode->ptr);matrixElementNode->ptr = NULL;free(matrixElementNode);matrixElementNode = NULL;}// ...// 释放其他指针}matrix.c
/******************************************************************************** File Name : matrix.c* Library/Module Name : MatrixComputation* Author : Marc Pony(marc_pony@163.com)* Create Date : 2023/2/24* Abstract Description : 矩阵运算库源文件*******************************************************************************//******************************************************************************** (1)Debug Switch Section*******************************************************************************//******************************************************************************** (2)Include File Section*******************************************************************************/#include "matrix.h"/******************************************************************************** (3)Macro Define Section*******************************************************************************//******************************************************************************** (4)Struct(Data Types) Define Section*******************************************************************************//******************************************************************************** (5)Prototype Declare Section*******************************************************************************//******************************************************************************** (6)Global Variable Declare Section*******************************************************************************//******************************************************************************** (7)File Static Variable Define Section*******************************************************************************//******************************************************************************** (8)Function Define Section*******************************************************************************/VOID print_matrix(MATRIX* a, STRING string){INDEX i, j;printf("matrix %s:", string);printf(" ");for (i = 0; i rows; i++){for (j = 0; j columns; j++){printf("%f ", a->p[i * a->columns + j]);}printf(" ");}printf(" ");}/**********************************************************************************************Function: creat_matrixDescription: 创建矩阵Input: 矩阵行数rows，列数columnsOutput: 错误号指针errorID，栈指针SInput_Output: 无Return: 矩阵指针Author: Marc Pony(marc_pony@163.com)***********************************************************************************************/MATRIX* creat_matrix(_IN INTEGER rows, _IN INTEGER columns, _OUT ERROR_ID* errorID, _OUT STACKS* S){MATRIX* matrix = NULL;MATRIX_NODE* matrixNode = NULL;MATRIX_ELEMENT_NODE* matrixElementNode = NULL;if (errorID == NULL){return NULL;}*errorID = _ERROR_NO_ERROR;if (rows columns = columns;matrix->p = (REAL*)malloc(rows * columns * sizeof(REAL)); //确保matrix非空才能执行指针操作if (matrix->p == NULL){free(matrix->p);matrix->p = NULL;free(matrix);matrix = NULL;free(matrixNode);matrixNode = NULL;free(matrixElementNode);matrixElementNode = NULL;*errorID = _ERROR_FAILED_TO_ALLOCATE_HEAP_MEMORY;return NULL;}matrixNode->ptr = matrix;matrixNode->next = S->matrixNode;S->matrixNode = matrixNode;matrixElementNode->ptr = matrix->p;matrixElementNode->next = S->matrixElementNode;S->matrixElementNode = matrixElementNode;return matrix;}/**********************************************************************************************Function: creat_multiple_matricesDescription: 创建多个矩阵Input: 矩阵行数rows，列数columns，个数countOutput: 错误号指针errorID，栈指针SInput_Output: 无Return: 矩阵指针Author: Marc Pony(marc_pony@163.com)***********************************************************************************************/MATRIX* creat_multiple_matrices(_IN INTEGER rows, _IN INTEGER columns, _IN INTEGER count, _OUT ERROR_ID* errorID, _OUT STACKS* S){INDEX i;MATRIX* matrix = NULL, *p = NULL;MATRIX_NODE* matrixNode = NULL;if (errorID == NULL){return NULL;}*errorID = _ERROR_NO_ERROR;if (rows matrixNode = matrixNode;return matrix;}/**********************************************************************************************Function: creat_zero_matrixDescription: 创建零矩阵Input: 矩阵行数rows，列数columnsOutput: 错误号指针errorID，栈指针SInput_Output: 无Return: 矩阵指针Author: Marc Pony(marc_pony@163.com)***********************************************************************************************/MATRIX* creat_zero_matrix(_IN INTEGER rows, _IN INTEGER columns, _OUT ERROR_ID* errorID, _OUT STACKS* S){MATRIX* matrix = NULL;if (errorID == NULL){return NULL;}*errorID = _ERROR_NO_ERROR;if (rows rows != B->rows || A->rows != C->rows || B->rows != C->rows|| A->columns != B->columns || A->columns != C->columns || B->columns != C->columns){errorID = _ERROR_MATRIX_ROWS_OR_COLUMNS_NOT_EQUAL;return errorID;}rows = A->rows;columns = A->columns;for (i = 0; i p[i * columns + j] + B->p[i * columns + j];}}return errorID;}/**********************************************************************************************Function: matrix_subtractionDescription: 矩阵A - 矩阵B = 矩阵CInput: 矩阵A,矩阵BOutput: 矩阵CInput_Output: 无Return: 错误号Author: Marc Pony(marc_pony@163.com)***********************************************************************************************/ERROR_ID matrix_subtraction(_IN MATRIX* A, _IN MATRIX* B, _OUT MATRIX* C){INDEX i, j;INTEGER rows, columns;ERROR_ID errorID = _ERROR_NO_ERROR;if (A == NULL || B == NULL || C == NULL){errorID = _ERROR_INPUT_PARAMETERS_ERROR;return errorID;}if (A->rows != B->rows || A->rows != C->rows || B->rows != C->rows|| A->columns != B->columns || A->columns != C->columns || B->columns != C->columns){errorID = _ERROR_MATRIX_ROWS_OR_COLUMNS_NOT_EQUAL;return errorID;}rows = A->rows;columns = A->columns;for (i = 0; i p[i * columns + j] - B->p[i * columns + j];}}return errorID;}/**********************************************************************************************Function: matrix_multiplicationDescription: 矩阵乘法C = A * BInput: 矩阵A,矩阵BOutput: 矩阵CInput_Output: 无Return: 错误号Author: Marc Pony(marc_pony@163.com)***********************************************************************************************/ERROR_ID matrix_multiplication(_IN MATRIX* A, _IN MATRIX* B, _OUT MATRIX* C){INDEX i, j, k;REAL sum;ERROR_ID errorID = _ERROR_NO_ERROR;if (A == NULL || B == NULL || C == NULL){errorID = _ERROR_INPUT_PARAMETERS_ERROR;return errorID;}if (A->columns != B->rows || A->rows != C->rows || B->columns != C->columns){errorID = _ERROR_MATRIX_MULTIPLICATION;return errorID;}for (i = 0; i rows; i++){for (j = 0; j columns; j++){sum = 0.0;for (k = 0; k columns; k++){sum += A->p[i * A->columns + k] * B->p[k * B->columns + j];}C->p[i * B->columns + j] = sum;}}return errorID;}/**********************************************************************************************Function: matrix_inverseDescription: 矩阵求逆Input: 矩阵AOutput: 矩阵A的逆矩阵Input_Output: 无Return: 错误号Author: Marc Pony(marc_pony@163.com)***********************************************************************************************/ERROR_ID matrix_inverse(_IN MATRIX* A, _OUT MATRIX* invA){INDEX i;INTEGER n;MATRIX* ATemp = NULL;ERROR_ID errorID = _ERROR_NO_ERROR;STACKS S;if (A == NULL || invA == NULL){errorID = _ERROR_INPUT_PARAMETERS_ERROR;return errorID;}if (A->rows != A->columns){errorID = _ERROR_MATRIX_MUST_BE_SQUARE;return errorID;}init_stack(&S);n = A->rows;ATemp = creat_matrix(n, n, &errorID, &S);if (errorID != _ERROR_NO_ERROR) goto EXIT;memcpy(ATemp->p, A->p, n * n * sizeof(REAL));memset(invA->p, 0, n * n * sizeof(REAL));for (i = 0; i p[i * n + i] = 1.0;}errorID = solve_matrix_equation_by_lup_decomposition(ATemp, invA);EXIT:free_stack(&S);return errorID;}/**********************************************************************************************Function: matrix_transposeDescription: 矩阵转置Input: 矩阵AOutput: 矩阵A的转置Input_Output: 无Return: 错误号Author: Marc Pony(marc_pony@163.com)***********************************************************************************************/ERROR_ID matrix_transpose(_IN MATRIX* A, _OUT MATRIX* transposeA){INDEX i, j;ERROR_ID errorID = _ERROR_NO_ERROR;if (A == NULL || transposeA == NULL){errorID = _ERROR_INPUT_PARAMETERS_ERROR;return errorID;}if (A->rows != transposeA->columns || A->columns != transposeA->rows){errorID = _ERROR_MATRIX_TRANSPOSE_FAILED;return errorID;}for (i = 0; i rows; i++){for (j = 0; j columns; j++){transposeA->p[j * A->rows + i] = A->p[i * A->columns + j];}}return errorID;}/**********************************************************************************************Function: matrix_traceDescription: 矩阵的迹Input: 矩阵AOutput: 矩阵A的迹Input_Output: 无Return: 错误号Author: Marc Pony(marc_pony@163.com)***********************************************************************************************/ERROR_ID matrix_trace(_IN MATRIX* A, _OUT REAL *trace){INDEX i;ERROR_ID errorID = _ERROR_NO_ERROR;if (A == NULL || trace == NULL){errorID = _ERROR_INPUT_PARAMETERS_ERROR;return errorID;}if (A->rows != A->columns){errorID = _ERROR_MATRIX_MUST_BE_SQUARE;return errorID;}*trace = 0.0;for (i = 0; i rows; i++){*trace += A->p[i * A->columns + i];}return errorID;}/**********************************************************************************************Function: lup_decompositionDescription: n行n列矩阵LUP分解PA = L * UInput: n行n列矩阵AOutput: n行n列下三角矩阵L，n行n列上三角矩阵U，n行n列置换矩阵PInput_Output: 无Return: 错误号Author: Marc Pony(marc_pony@163.com)参考：https://zhuanlan.zhihu.com/p/84210687***********************************************************************************************/ERROR_ID lup_decomposition(_IN MATRIX* A, _OUT MATRIX* L, _OUT MATRIX* U, _OUT MATRIX* P){INDEX i, j, k, index, s, t;INTEGER n;REAL maxValue, temp;ERROR_ID errorID = _ERROR_NO_ERROR;if (A == NULL || L == NULL || U == NULL || P == NULL){errorID = _ERROR_INPUT_PARAMETERS_ERROR;return errorID;}if (A->rows != A->columns){errorID = _ERROR_MATRIX_MUST_BE_SQUARE;return errorID;}n = A->rows;memcpy(U->p, A->p, n * n * sizeof(REAL));memset(L->p, 0, n * n * sizeof(REAL));memset(P->p, 0, n * n * sizeof(REAL));for (i = 0; i p[i * n + i] = 1.0;P->p[i * n + i] = 1.0;}for (j = 0; j = j) that maximizes | U(i, j) |index = -1;maxValue = 0.0;for (i = j; i p[i * n + j]);if (temp > maxValue){maxValue = temp;index = i;}}if (index == -1){continue;}//Interchange rows of U : U(j, j : n) < ->U(i, j : n)for (k = j; k p[s];U->p[s] = U->p[t];U->p[t] = temp;}//Interchange rows of L : L(j, 1 : j - 1) < ->L(i, 1 : j - 1)for (k = 0; k p[s];L->p[s] = L->p[t];L->p[t] = temp;}//Interchange rows of P : P(j, 1 : n) < ->P(i, 1 : n)for (k = 0; k p[s];P->p[s] = P->p[t];P->p[t] = temp;}for (i = j + 1; i p[s] = U->p[s] / U->p[j * n + j];for (k = j; k p[i * n + k] -= L->p[s] * U->p[j * n + k];}}}return errorID;}/**********************************************************************************************Function: solve_matrix_equation_by_lup_decompositionDescription: LUP分解解矩阵方程AX=B,其中A为n行n列矩阵，B为n行m列矩阵，X为n行m列待求矩阵(写到矩阵B)Input: n行n列矩阵AOutput: 无Input_Output: n行m列矩阵B(即n行m列待求矩阵X)Return: 错误号Author: Marc Pony(marc_pony@163.com)***********************************************************************************************/ERROR_ID solve_matrix_equation_by_lup_decomposition(_IN MATRIX* A, _IN_OUT MATRIX* B){INDEX i, j, k, index, s, t;INTEGER n, m;REAL sum, maxValue, temp;MATRIX* L = NULL, * U = NULL, * y = NULL;ERROR_ID errorID = _ERROR_NO_ERROR;STACKS S;if (A == NULL || B == NULL){errorID = _ERROR_INPUT_PARAMETERS_ERROR;return errorID;}if (A->rows != A->columns){errorID = _ERROR_MATRIX_MUST_BE_SQUARE;return errorID;}init_stack(&S);n = A->rows;m = B->columns;L = creat_matrix(n, n, &errorID, &S);if (errorID != _ERROR_NO_ERROR) goto EXIT;U = creat_matrix(n, n, &errorID, &S);if (errorID != _ERROR_NO_ERROR) goto EXIT;y = creat_matrix(n, m, &errorID, &S);if (errorID != _ERROR_NO_ERROR) goto EXIT;memcpy(U->p, A->p, n * n * sizeof(REAL));memset(L->p, 0, n * n * sizeof(REAL));for (i = 0; i p[i * n + i] = 1.0;}for (j = 0; j = j) that maximizes | U(i, j) |index = -1;maxValue = 0.0;for (i = j; i p[i * n + j]);if (temp > maxValue){maxValue = temp;index = i;}}if (index == -1){continue;}//Interchange rows of U : U(j, j : n) < ->U(i, j : n)for (k = j; k p[s];U->p[s] = U->p[t];U->p[t] = temp;}//Interchange rows of L : L(j, 1 : j - 1) < ->L(i, 1 : j - 1)for (k = 0; k p[s];L->p[s] = L->p[t];L->p[t] = temp;}//Interchange rows of P : P(j, 1 : n) < ->P(i, 1 : n), C = P * B，等价于对B交换行for (k = 0; k p[s];B->p[s] = B->p[t];B->p[t] = temp;}for (i = j + 1; i p[s] = U->p[s] / U->p[j * n + j];for (k = j; k p[i * n + k] -= L->p[s] * U->p[j * n + k];}}}for (i = 0; i p[i * n + i]) p[i * m + j] - sum;}}//U * x = yfor (j = 0; j = 0; i--){sum = 0.0;for (k = i + 1; k p[i * n + k] * B->p[k * m + j];}B->p[i * m + j] = (y->p[i * m + j] - sum) / U->p[i * n + i];}}EXIT:free_stack(&S);return errorID;}test_matrix.c
#include "matrix.h"void main(){REAL a[3 * 3] = { 1,2,3,6,5,5,8,7,2 };REAL b[3 * 3] = {1,2,3,6,5,4,3,2,1};MATRIX *A = NULL, * B = NULL, * C = NULL, * D = NULL, * E = NULL, * Z = NULL, * invA = NULL, *m = NULL;ERROR_ID errorID = _ERROR_NO_ERROR;REAL trace;STACKS S;init_stack(&S);Z = creat_zero_matrix(3, 3, &errorID, &S);print_matrix(Z, "Z");E = creat_eye_matrix(3, &errorID, &S);print_matrix(E, "E");A = creat_matrix(3, 3, &errorID, &S);A->p = a;print_matrix(A, "A");B = creat_matrix(3, 3, &errorID, &S);B->p = b;print_matrix(B, "B");C = creat_matrix(3, 3, &errorID, &S);D = creat_matrix(3, 3, &errorID, &S);invA = creat_matrix(3, 3, &errorID, &S);errorID = matrix_add(A, B, C);errorID = matrix_subtraction(A, B, C);errorID = matrix_multiplication(A, B, C);errorID = matrix_transpose(A, D);print_matrix(D, "D");errorID = matrix_trace(A, &trace);errorID = matrix_inverse(A, invA); print_matrix(invA, "invA");m = creat_multiple_matrices(3, 3, 2, &errorID, &S);m[0].p = a;m[1].p = b;free_stack(&S);} 参考资料The NURBS Book(Second Edition). Les Piegl,Wayne Tiller