纯c语言优雅地实现矩阵运算库的方法

2024-10-24 14:03:17

(1)采用预定义的数据类型，避免直接使用编译器定义的数据类型
(2)基于对象编程，定义矩阵对象
(3)除了特别编写的内存处理函数(使用栈链表保存、释放动态分配的内存地址)，不允许任何函数直接分配和释放内存
(4)防御性编程，对输入参数做有效性检查，并返回错误号
(5)注意编程细节的打磨

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，达到一劳永逸的效果。定义INDEX与INTEGER数据类型是为了在编程时区分索引变量与普通整形变量，同样提高可读性。

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

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

这里，用一级指针而非二级指针指向矩阵的数据内存地址，有诸多原因，详见博文:为什么我推荐使用一级指针创建二维数组？。

(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
*******************************************************************************/

void print_matrix(MATRIX* a, STRING string);

/**********************************************************************************************
Function: creat_matrix
Description: 创建矩阵
Input: 矩阵行数rows，列数columns
Output: 错误号指针errorID，栈指针S
Input_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);

/**********************************************************************************************
Function: creat_multiple_matrices
Description: 创建多个矩阵
Input: 矩阵行数rows，列数columns，个数count
Output: 错误号指针errorID，栈指针S
Input_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);

/**********************************************************************************************
Function: creat_zero_matrix
Description: 创建零矩阵
Input: 矩阵行数rows，列数columns
Output: 错误号指针errorID，栈指针S
Input_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);

/**********************************************************************************************
Function: creat_eye_matrix
Description: 创建单位矩阵
Input: 矩阵行数rows，列数columns
Output: 错误号指针errorID，栈指针S
Input_Output: 无
Return: 矩阵指针
Author: Marc Pony(marc_pony@163.com)
***********************************************************************************************/
MATRIX* creat_eye_matrix(_IN INTEGER n, _OUT ERROR_ID* errorID, _OUT STACKS* S);

/**********************************************************************************************
Function: matrix_add
Description: 矩阵A + 矩阵B = 矩阵C
Input: 矩阵A,矩阵B
Output: 矩阵C
Input_Output: 无
Return: 错误号
Author: Marc Pony(marc_pony@163.com)
***********************************************************************************************/
ERROR_ID matrix_add(_IN MATRIX* A, _IN MATRIX* B, _OUT MATRIX *C);

/**********************************************************************************************
Function: matrix_subtraction
Description: 矩阵A - 矩阵B = 矩阵C
Input: 矩阵A,矩阵B
Output: 矩阵C
Input_Output: 无
Return: 错误号
Author: Marc Pony(marc_pony@163.com)
***********************************************************************************************/
ERROR_ID matrix_subtraction(_IN MATRIX* A, _IN MATRIX* B, _OUT MATRIX* C);

/**********************************************************************************************
Function: matrix_multiplication
Description: 矩阵乘法C = A * B
Input: 矩阵A,矩阵B
Output: 矩阵C
Input_Output: 无
Return: 错误号
Author: Marc Pony(marc_pony@163.com)
***********************************************************************************************/
ERROR_ID matrix_multiplication(_IN MATRIX* A, _IN MATRIX* B, _OUT MATRIX* C);

/**********************************************************************************************
Function: matrix_inverse
Description: 矩阵求逆
Input: 矩阵A
Output: 矩阵A的逆矩阵
Input_Output: 无
Return: 错误号
Author: Marc Pony(marc_pony@163.com)
***********************************************************************************************/
ERROR_ID matrix_inverse(_IN MATRIX* A, _OUT MATRIX* invA);

/**********************************************************************************************
Function: matrix_transpose
Description: 矩阵转置
Input: 矩阵A
Output: 矩阵A的转置
Input_Output: 无
Return: 错误号
Author: Marc Pony(marc_pony@163.com)
***********************************************************************************************/
ERROR_ID matrix_transpose(_IN MATRIX* A, _OUT MATRIX* transposeA);

/**********************************************************************************************
Function: matrix_trace
Description: 矩阵的迹
Input: 矩阵A
Output: 矩阵A的迹
Input_Output: 无
Return: 错误号
Author: Marc Pony(marc_pony@163.com)
***********************************************************************************************/
ERROR_ID matrix_trace(_IN MATRIX* A, _OUT REAL* trace);

/**********************************************************************************************
Function: lup_decomposition
Description: n行n列矩阵LUP分解PA = L * U
Input: n行n列矩阵A
Output: n行n列下三角矩阵L，n行n列上三角矩阵U，n行n列置换矩阵P
Input_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);

/**********************************************************************************************
Function: solve_matrix_equation_by_lup_decomposition
Description: LUP分解解矩阵方程AX=B,其中A为n行n列矩阵，B为n行m列矩阵，X为n行m列待求矩阵(写到矩阵B)
Input: n行n列矩阵A
Output: 无
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);

#endif

common.c

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

/*******************************************************************************
* (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
*******************************************************************************/

/*******************************************************************************
* (6)Global Variable Declare Section
*******************************************************************************/

/*******************************************************************************
* (7)File Static Variable Define Section
*******************************************************************************/

/*******************************************************************************
* (8)Function Define Section
*******************************************************************************/

/**********************************************************************************************
Function: init_stack
Description: 初始化栈
Input: 无
Output: 无
Input_Output: 栈指针
Return: 无
Author: Marc Pony(marc_pony@163.com)
***********************************************************************************************/
VOID init_stack(_IN_OUT STACKS* S)
{
 if (S == NULL)
 {
  return;
 }

 memset(S, 0, sizeof(STACKS));
}

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

 if (S == NULL)
 {
  return;
 }

 while (S->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 :                      2021/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("\n");
 for (i = 0; i < a->rows; i++)
 {
  for (j = 0; j < a->columns; j++)
  {
   printf("%f  ", a->p[i * a->columns + j]);
  }
  printf("\n");
 }

 printf("\n");
}

/**********************************************************************************************
Function: creat_matrix
Description: 创建矩阵
Input: 矩阵行数rows，列数columns
Output: 错误号指针errorID，栈指针S
Input_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;

 *errorID = _ERROR_NO_ERROR;
 if (rows <= 0 || columns <= 0 || errorID == NULL || S == NULL)
 {
  *errorID = _ERROR_INPUT_PARAMETERS_ERROR;
  return NULL;
 }

 matrix = (MATRIX*)malloc(sizeof(MATRIX));
 matrixNode = (MATRIX_NODE*)malloc(sizeof(MATRIX_NODE));
 matrixElementNode = (MATRIX_ELEMENT_NODE*)malloc(sizeof(MATRIX_ELEMENT_NODE));
 if (matrix == NULL || matrixNode == NULL || matrixElementNode == NULL)
 {
  free(matrix);
  matrix = NULL;
  free(matrixNode);
  matrixNode = NULL;
  free(matrixElementNode);
  matrixElementNode = NULL;

  *errorID = _ERROR_FAILED_TO_ALLOCATE_HEAP_MEMORY;
  return NULL;
 }

 matrix->rows = rows;
 matrix->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_matrices
Description: 创建多个矩阵
Input: 矩阵行数rows，列数columns，个数count
Output: 错误号指针errorID，栈指针S
Input_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;

 *errorID = _ERROR_NO_ERROR;
 if (rows <= 0 || columns <= 0 || count <= 0 || errorID == NULL || S == NULL)
 {
  *errorID = _ERROR_INPUT_PARAMETERS_ERROR;
  return NULL;
 }

 matrix = (MATRIX*)malloc(count * sizeof(MATRIX));
 matrixNode = (MATRIX_NODE*)malloc(sizeof(MATRIX_NODE));
 if (matrix == NULL || matrixNode == NULL)
 {
  free(matrix);
  matrix = NULL;
  free(matrixNode);
  matrixNode = NULL;

  *errorID = _ERROR_FAILED_TO_ALLOCATE_HEAP_MEMORY;
  return NULL;
 }

 for (i = 0; i < count; i++)
 {
  p = creat_matrix(rows, columns, errorID, S);
  if (p == NULL)
  {
   free(matrix);
   matrix = NULL;
   free(matrixNode);
   matrixNode = NULL;

   *errorID = _ERROR_FAILED_TO_ALLOCATE_HEAP_MEMORY;
   return NULL;
  }

  matrix[i] = *p;
 }

 matrixNode->ptr = matrix;
 matrixNode->next = S->matrixNode;
 S->matrixNode = matrixNode;

 return matrix;
}

/**********************************************************************************************
Function: creat_zero_matrix
Description: 创建零矩阵
Input: 矩阵行数rows，列数columns
Output: 错误号指针errorID，栈指针S
Input_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;

 *errorID = _ERROR_NO_ERROR;
 if (rows <= 0 || columns <= 0 || errorID == NULL || S == NULL)
 {
  *errorID = _ERROR_INPUT_PARAMETERS_ERROR;
  return NULL;
 }

 matrix = creat_matrix(rows, columns, errorID, S);
 if (*errorID == _ERROR_NO_ERROR)
 {
  memset(matrix->p, 0, rows * columns * sizeof(REAL));
 }

 return matrix;
}

/**********************************************************************************************
Function: creat_eye_matrix
Description: 创建单位矩阵
Input: 矩阵行数rows，列数columns
Output: 错误号指针errorID，栈指针S
Input_Output: 无
Return: 矩阵指针
Author: Marc Pony(marc_pony@163.com)
***********************************************************************************************/
MATRIX* creat_eye_matrix(_IN INTEGER n, _OUT ERROR_ID* errorID, _OUT STACKS* S)
{
 INDEX i;
 MATRIX* matrix = NULL;

 *errorID = _ERROR_NO_ERROR;
 if (n <= 0 || errorID == NULL || S == NULL)
 {
  *errorID = _ERROR_INPUT_PARAMETERS_ERROR;
  return NULL;
 }

 matrix = creat_matrix(n, n, errorID, S);
 if (*errorID == _ERROR_NO_ERROR)
 {
  memset(matrix->p, 0, n * n * sizeof(REAL));
  for (i = 0; i < n; i++)
  {
   matrix->p[i * n + i] = 1.0;
  }
 }

 return matrix;
}

/**********************************************************************************************
Function: matrix_add
Description: 矩阵A + 矩阵B = 矩阵C
Input: 矩阵A,矩阵B
Output: 矩阵C
Input_Output: 无
Return: 错误号
Author: Marc Pony(marc_pony@163.com)
***********************************************************************************************/
ERROR_ID matrix_add(_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 < rows; i++)
 {
  for (j = 0; j < columns; j++)
  {
   C->p[i * columns + j] = A->p[i * columns + j] + B->p[i * columns + j];
  }
 }
 return errorID;
}

/**********************************************************************************************
Function: matrix_subtraction
Description: 矩阵A - 矩阵B = 矩阵C
Input: 矩阵A,矩阵B
Output: 矩阵C
Input_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 < rows; i++)
 {
  for (j = 0; j < columns; j++)
  {
   C->p[i * columns + j] = A->p[i * columns + j] - B->p[i * columns + j];
  }
 }
 return errorID;
}

/**********************************************************************************************
Function: matrix_multiplication
Description: 矩阵乘法C = A * B
Input: 矩阵A,矩阵B
Output: 矩阵C
Input_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 < A->rows; i++)
 {
  for (j = 0; j < B->columns; j++)
  {
   sum = 0.0;
   for (k = 0; k < A->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_inverse
Description: 矩阵求逆
Input: 矩阵A
Output: 矩阵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 < n; i++)
 {
  invA->p[i * n + i] = 1.0;
 }

 errorID = solve_matrix_equation_by_lup_decomposition(ATemp, invA);

EXIT:
 free_stack(&S);
 return errorID;
}

/**********************************************************************************************
Function: matrix_transpose
Description: 矩阵转置
Input: 矩阵A
Output: 矩阵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 < A->rows; i++)
 {
  for (j = 0; j < A->columns; j++)
  {
   transposeA->p[j * A->rows + i] = A->p[i * A->columns + j];
  }
 }
 return errorID;
}

/**********************************************************************************************
Function: matrix_trace
Description: 矩阵的迹
Input: 矩阵A
Output: 矩阵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 < A->rows; i++)
 {
  *trace += A->p[i * A->columns + i];
 }
 return errorID;
}

/**********************************************************************************************
Function: lup_decomposition
Description: n行n列矩阵LUP分解PA = L * U
Input: n行n列矩阵A
Output: n行n列下三角矩阵L，n行n列上三角矩阵U，n行n列置换矩阵P
Input_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 < n; i++)
 {
  L->p[i * n + i] = 1.0;
  P->p[i * n + i] = 1.0;
 }

 for (j = 0; j < n - 1; j++)
 {
  //Select i(>= j) that maximizes | U(i, j) |
  index = -1;
  maxValue = 0.0;
  for (i = j; i < n; i++)
  {
   temp = fabs(U->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 < n; k++)
  {
   s = j * n + k;
   t = index * n + k;
   temp = U->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 < j; k++)
  {
   s = j * n + k;
   t = index * n + k;
   temp = L->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 < n; k++)
  {
   s = j * n + k;
   t = index * n + k;
   temp = P->p[s];
   P->p[s] = P->p[t];
   P->p[t] = temp;
  }

  for (i = j + 1; i < n; i++)
  {
   s = i * n + j;
   L->p[s] = U->p[s] / U->p[j * n + j];
   for (k = j; k < n; k++)
   {
    U->p[i * n + k] -= L->p[s] * U->p[j * n + k];
   }
  }
 }
 return errorID;
}

/**********************************************************************************************
Function: solve_matrix_equation_by_lup_decomposition
Description: LUP分解解矩阵方程AX=B,其中A为n行n列矩阵，B为n行m列矩阵，X为n行m列待求矩阵(写到矩阵B)
Input: n行n列矩阵A
Output: 无
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 < n; i++)
 {
  L->p[i * n + i] = 1.0;
 }

 for (j = 0; j < n - 1; j++)
 {
  //Select i(>= j) that maximizes | U(i, j) |
  index = -1;
  maxValue = 0.0;
  for (i = j; i < n; i++)
  {
   temp = fabs(U->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 < n; k++)
  {
   s = j * n + k;
   t = index * n + k;
   temp = U->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 < j; k++)
  {
   s = j * n + k;
   t = index * n + k;
   temp = L->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 < m; k++)
  {
   s = j * m + k;
   t = index * m + k;
   temp = B->p[s];
   B->p[s] = B->p[t];
   B->p[t] = temp;
  }

  for (i = j + 1; i < n; i++)
  {
   s = i * n + j;
   L->p[s] = U->p[s] / U->p[j * n + j];
   for (k = j; k < n; k++)
   {
    U->p[i * n + k] -= L->p[s] * U->p[j * n + k];
   }
  }
 }

 for (i = 0; i < n; i++)
 {
  if (fabs(U->p[i * n + i]) < 1.0e-20)
  {
   errorID = _ERROR_MATRIX_EQUATION_HAS_NO_SOLUTIONS;
   goto EXIT;
  }
 }

 //L * y = C
 for (j = 0; j < m; j++)
 {
  for (i = 0; i < n; i++)
  {
   sum = 0.0;
   for (k = 0; k < i; k++)
   {
    sum = sum + L->p[i * n + k] * y->p[k * m + j];
   }
   y->p[i * m + j] = B->p[i * m + j] - sum;
  }
 }

 //U * x = y
 for (j = 0; j < m; j++)
 {
  for (i = n - 1; i >= 0; i--)
  {
   sum = 0.0;
   for (k = i + 1; k < n; k++)
   {
    sum += U->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

到此这篇关于纯c语言优雅地实现矩阵运算库的方法的文章就介绍到这了,更多相关c语言矩阵运算内容请搜索我们以前的文章或继续浏览下面的相关文章希望大家以后多多支持我们！

C语言实现矩阵运算案例详解

C语言实现矩阵运算给定一个n×n的方阵,本题要求计算该矩阵除副对角线.最后一列和最后一行以外的所有元素之和.副对角线为从矩阵的右上角至左下角的连线. 输入格式: 输入第一行给出正整数n(1<n≤10):随后n行,每行给出n个整数,其间以空格分隔. 输出格式: 在一行中给出该矩阵除副对角线.最后一列和最后一行以外的所有元素之和. 输入样例: 4 2 3 4 1 5 6 1 1 7 1 8 1 1 1 1 1 输出样例: 35 #include <stdio.h> #include <
纯c语言优雅地实现矩阵运算库的方法

目录 1.一个优雅好用的c语言库必须满足哪些条件 2.实现一个矩阵运算库的几点思考 (1)采用预定义的数据类型,避免直接使用编译器定义的数据类型 (2)基于对象编程,定义矩阵对象 (3)除了特别编写的内存处理函数(使用栈链表保存.释放动态分配的内存地址),不允许任何函数直接分配和释放内存 (4)防御性编程,对输入参数做有效性检查,并返回错误号 (5)注意编程细节的打磨 3.完整c程序参考资料编程既是技术输出也是艺术创作.鉴赏高手写的程序,往往让人眼前一亮,他们思路.逻辑清晰,所呈现的代码简洁
C语言实现用户态线程库案例

轮子年年有人造,我们也来凑热闹,参考协程实现,大概有以下几种方法: 1)利用setjmp,longjmp 2)利用ucontext接口函数 3)汇编 (线程无非就是多了个抢占功能,由定时器触发,而非自愿让出运行权限) 因为我写的时候还没看到其他帖子,如果看到了,铁定会用最直观的ucontext接口写的(注意,在macOSX中已经标注为废除,头文件得换做sys/ucontext.h),结果就是我用了汇编来写,但是尽量不用汇编来写整个switch_to调度函数(这样有个明显的坏处,那就是用gas/n
R语言学习RcppEigen进行矩阵运算

目录创建cpp文件代码示例其他矩阵操作命名基础用法定义矩阵对矩阵的一些基础操作1 基础操作2 矩阵基础运算1 矩阵基础运算2 求最小最大值.迹等点乘等特征值与特征向量形式转换矩阵初始化0 Map等操作求解Ax = b 前面我们介绍了一些基本的Rcpp的用法:让你的R代码更快--Rcpp入门,但用基础的Rcpp来进行矩阵运算还是非常麻烦,没有现成的函数来让我们使用. 这时我们就想到:是否可以调用别的库来解决矩阵运算的一些问题呢?这就需要我们的RcppEigen包,也就是C+
R语言技巧Rcpp与Eigen库之间的相互转换

当我们在使用Rcpp时,进行矩阵运算最简单的是使用Eigen库进行相关操作,可以很轻松地讲R中向量化与矩阵化的思想应用到C++代码上,从而对代码进行加速.可参考前面的博客:利用RcppEigen进行矩阵运算. 但有时,我们却必须使用Rcpp进行DataFrame,List等对象格式的处理.或者如果我们涉及到缺失值的处理,也需要使用Rcpp中的函数来做. 所以,如何在两种矩阵或向量格式,NumericVector/Matrix与VectorXd/MatrixXd之间相互转化就变得非常重要. 我们可
一文带你了解Go语言实现的并发神库conc

目录前言 worker池 Stream ForEach和map ForEach map 总结前言哈喽,大家好,我是asong:前几天逛github发现了一个有趣的并发库-conc,其目标是: 更难出现goroutine泄漏处理panic更友好并发代码可读性高从简介上看主要封装功能如下: 对waitGroup进行封装,避免了产生大量重复代码,并且也封装recover,安全性更高提供panics.Catcher封装recover逻辑,统一捕获panic,打印调用栈一些信息提供一个并发
易语言去除多余的支持库方法

易语言不得不说是一个很方便的编程工具,其强大的支持库.模块能快速的帮助我们编程,不过,一般完整版的易语言都带了很多的支持库,我们其实不需要,反倒让我们找我们需要的支持库的时候感到麻烦,怎么去除多余的支持库呢? 1.方法一: 比较推荐的一种方法,临时性的去除,想用的时候还可以加载出来. 打开易语言. 2.点击上方菜单中的工具,选择弹出菜单里面的支持库配置,如图: 3.会弹出一个窗口,里面会列出你电脑里面保存的各种支持库,可以看到有好多种: 4.把你不想用的支持库前面的对扣去掉,易语言就不会加载这个
Nginx纯配置实现日志实时上报的思路与方法

目录前言实现思路实现步骤 1. 编译 Nginx 2. 配置文件如下总结前言 Nginx 作为常用的负载均衡网关. 会产生大量的日志. 但是由于 Nginx 的配置文件是一种声明式的编程范式, 不方便描述流程控制, 因此不能通过简单的指令实现日志的上报. 通常 Nginx 的日志上报是需要写一个shell脚本或其他语言的脚本来定时解析 Nginx 的 log 文件, 然后进行上报. 利用 NJS 模块, 可以实现实时的日志上报. 但是由于 NJS 模块支持的指令的限制, 无法通过单一指
go语言通过odbc操作Access数据库的方法

本文实例讲述了go语言通过odbc操作Access数据库的方法.分享给大家供大家参考.具体如下: 这里需要用到go-odbc库,下载地址为:https://github.com/weigj/go-odbc 复制代码代码如下: package main; import ( "fmt" "database/sql" _"odbc/driver" ) func main(){ conn,err := sql.Open("odbc&q
Angular 如何使用第三方库的方法

Angular 的组件与模块看似好像与现有各种第三方类库(例如:lodash.moment 等)使用上有点格格不入,这很大的原因是 TypeScript 造成的假象.三足鼎立的前端其实都是雷同的,不管是哪种前端框架都可以使用到这些第三方类库. 以下我将从另一个视角解释 Angular 如何使用第三方类库的一种经验做法. 一.写在前面在开始之前,需要先了解一下 TypeScript 模块系统 --模块是指在其自身作用域里执行,而不是在全局作用域里:模块间是依靠 export 和 import 建
vue-cli3 DllPlugin 提取公用库的方法

vue 开发过程中,保存一次就会编译一次,如果能够减少编译的时间,哪怕是一丁点,也能节省不少时间.开发过程中个人编写的源文件才会频繁变动,而一些库文件我们一般是不会去改动的.如果能把这些库文件提取出来,就能减少打包体积,加快编译速度.本文主要讲述在 vue-cli3 中利用 DllPlugin 来进行预编译. 1.安装相关插件 yarn add webpack-cli@^3.2.3 add-asset-html-webpack-plugin@^3.1.3 clean-webpack-plugin