纯c语言优雅地实现矩阵运算库的方法
短信预约 -IT技能 免费直播动态提醒
编程既是技术输出也是艺术创作。鉴赏高手写的程序，往往让人眼前一亮，他们思路、逻辑清晰，所呈现的代码简洁、优雅、高效，令人为之叹服。烂代码则像“屎山"一般让人作呕，软件难以维护最大的原因除了需求模糊的客观因素，最重要的主观因素还是代码写得烂。有生之年，愿能保持对编程的热情，不断提升编程能力，真正体会其乐趣，共勉！
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
#ifndef __COMMON_H__
#define __COMMON_H__
#include <math.h>
#include <stdio.h>
#include <malloc.h>
#include <stdlib.h>
#include <time.h>
#include <memory.h>
#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 //输入参数错误
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;
VOID init_stack(_IN_OUT STACKS* S);
VOID free_stack(_IN STACKS* S);
#endifmatrix.h
#ifndef __MATRIX_H__
#define __MATRIX_H__
#include "common.h"
void print_matrix(MATRIX* a, STRING string);
MATRIX* creat_matrix(_IN INTEGER rows, _IN INTEGER columns, _OUT ERROR_ID* errorID, _OUT STACKS* S);
MATRIX* creat_multiple_matrices(_IN INTEGER rows, _IN INTEGER columns, _IN INTEGER count, _OUT ERROR_ID* errorID, _OUT STACKS* S);
MATRIX* creat_zero_matrix(_IN INTEGER rows, _IN INTEGER columns, _OUT ERROR_ID* errorID, _OUT STACKS* S);
MATRIX* creat_eye_matrix(_IN INTEGER n, _OUT ERROR_ID* errorID, _OUT STACKS* S);
ERROR_ID matrix_add(_IN MATRIX* A, _IN MATRIX* B, _OUT MATRIX *C);
ERROR_ID matrix_subtraction(_IN MATRIX* A, _IN MATRIX* B, _OUT MATRIX* C);
ERROR_ID matrix_multiplication(_IN MATRIX* A, _IN MATRIX* B, _OUT MATRIX* C);
ERROR_ID matrix_inverse(_IN MATRIX* A, _OUT MATRIX* invA);
ERROR_ID matrix_transpose(_IN MATRIX* A, _OUT MATRIX* transposeA);
ERROR_ID matrix_trace(_IN MATRIX* A, _OUT REAL* trace);
ERROR_ID lup_decomposition(_IN MATRIX* A, _OUT MATRIX* L, _OUT MATRIX* U, _OUT MATRIX* P);
ERROR_ID solve_matrix_equation_by_lup_decomposition(_IN MATRIX* A, _IN_OUT MATRIX* B);
#endifcommon.c
#include "common.h"
VOID init_stack(_IN_OUT STACKS* S)
{
if (S == NULL)
{
return;
}
memset(S, 0, sizeof(STACKS));
}
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
#include "matrix.h"
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");
}
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;
}
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;
}
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;
}
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;
}
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;
}
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;
}
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;
}
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;
}
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;
}
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;
}
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;
}
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语言 矩阵运算内容请搜索编程网以前的文章或继续浏览下面的相关文章希望大家以后多多支持编程网！
免责声明：
① 本站未注明“稿件来源”的信息均来自网络整理。其文字、图片和音视频稿件的所属权归原作者所有。本站收集整理出于非商业性的教育和科研之目的，并不意味着本站赞同其观点或证实其内容的真实性。仅作为临时的测试数据，供内部测试之用。本站并未授权任何人以任何方式主动获取本站任何信息。
② 本站未注明“稿件来源”的临时测试数据将在测试完成后最终做删除处理。有问题或投稿请发送至: 邮箱/279061341@qq.com QQ/279061341
