純c語言優(yōu)雅地實現(xiàn)矩陣運算庫的方法
更新時間：2021年08月17日 11:25:17 作者：Marc Pony
本文主要介紹了純c語言優(yōu)雅地實現(xiàn)矩陣運算庫，文中通過示例代碼介紹的非常詳細，具有一定的參考價值，感興趣的小伙伴們可以參考一下
編程既是技術輸出也是藝術創(chuàng)作。鑒賞高手寫的程序，往往讓人眼前一亮，他們思路、邏輯清晰，所呈現(xiàn)的代碼簡潔、優(yōu)雅、高效，令人為之嘆服。爛代碼則像“屎山"一般讓人作嘔，軟件難以維護最大的原因除了需求模糊的客觀因素，最重要的主觀因素還是代碼寫得爛。有生之年，愿能保持對編程的熱情，不斷提升編程能力，真正體會其樂趣，共勉！
1.一個優(yōu)雅好用的c語言庫必須滿足哪些條件
這里給出的軟件開發(fā)應遵循的一般原則，摘自Les Piegl和Wayne Tiller所著的一本非常經(jīng)典的書The NURBS Book(Second Edition)。
(1)工具性(toolability):應該使用可用的工具來建立新的應用程序；
(2)可移植性(portability):應用程序應該容易被移植到不同的軟件和硬件平臺；
(3)可重用性(reusability):程序的編寫應該便于重復使用代碼段；
(4)可檢驗性(testability):程序應該簡單一致，以便于測試和調(diào)試；
(5)可靠性(reliability):對程序運行過程中可能出現(xiàn)的各種錯誤應該進行合理、一致的處理，以使系統(tǒng)穩(wěn)定、可靠；
(6)可擴展性(enhanceability):代碼必須易于理解，以便可以容易地增加新的功能；
(7)可維護性(fixability):易于找出程序錯誤的位置；
(8)一致性(consistency):在整個庫中，編程習慣應保持一致；
(9)可讀性(communicability):程序應該易于閱讀和理解;
(10)編程風格(style of programming):代碼看起來像書上的數(shù)學公式那樣以便于讀者理解，同時遵循用戶友好的編程風格;
(11)易用性(usability):應該使一些非專業(yè)的用戶也能夠方便地使用所開發(fā)的庫來開發(fā)各種更高層次的應用程序;
(12)數(shù)值高效性(numerical efficiency):所有函數(shù)必須仔細推敲，保證其數(shù)值高效性;
(13)基于對象編程(object based programming):避免在函數(shù)間傳遞大量數(shù)據(jù)，并且使代碼易于理解。
2.實現(xiàn)一個矩陣運算庫的幾點思考
(1)采用預定義的數(shù)據(jù)類型，避免直接使用編譯器定義的數(shù)據(jù)類型
typedef unsigned int ERROR_ID; typedef int INDEX; typedef short FLAG; typedef int INTEGER; typedef double REAL; typedef char* STRING; typedef void VOID;
使用預定義的數(shù)據(jù)類型，有利于程序移植，且可以提高可讀性。例如，如果一個系統(tǒng)只支持單精度浮點數(shù)，那么只需修改數(shù)據(jù)類型REAL為float，達到一勞永逸的效果。定義INDEX與INTEGER數(shù)據(jù)類型是為了在編程時區(qū)分索引變量與普通整形變量，同樣提高可讀性。
(2)基于對象編程，定義矩陣對象
typedef struct matrix
{
INTEGER rows;
INTEGER columns;
REAL* p;
}MATRIX;
這里，用一級指針而非二級指針指向矩陣的數(shù)據(jù)內(nèi)存地址，有諸多原因，詳見博文:為什么我推薦使用一級指針創(chuàng)建二維數(shù)組？。
(3)除了特別編寫的內(nèi)存處理函數(shù)(使用棧鏈表保存、釋放動態(tài)分配的內(nèi)存地址)，不允許任何函數(shù)直接分配和釋放內(nèi)存
不恰當?shù)姆峙?、使用與釋放內(nèi)存很可能導致內(nèi)存泄漏、系統(tǒng)崩潰等致命的錯誤。如果一個函數(shù)需動態(tài)申請多個內(nèi)存，那么可能會寫出這樣啰嗦的程序:
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;
}
為了優(yōu)雅地實現(xiàn)動態(tài)內(nèi)存分配與釋放，Les Piegl大神分3步來處理內(nèi)存申請與釋放：
a)在進入一個新的程序時，一個內(nèi)存堆棧被初始化為空；
b)當需要申請內(nèi)存時，調(diào)用特定的函數(shù)來分配所需的內(nèi)存，并將指向內(nèi)存的指針存入堆棧中的正確位置；
c)在離開程序時，遍歷內(nèi)存堆棧，釋放其中的指針所指向的內(nèi)存。
程序結(jié)構(gòu)大致如下：
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)防御性編程，對輸入?yún)?shù)做有效性檢查，并返回錯誤號
例如輸入的矩陣行數(shù)、列數(shù)應該是正整數(shù)，指針參數(shù)必須非空等等。
(5)注意編程細節(jié)的打磨
a)操作符(逗號，等號等)兩邊必須空一格；
b)邏輯功能相同的程序間不加空行，邏輯功能獨立的程序間加空行；
c)條件判斷關鍵字(for if while等)后必須加一空格，起到強調(diào)作用，也更清晰；
d)函數(shù)內(nèi)部定義局部變量后，必須空一行后再編寫函數(shù)主體。
3.完整c程序
本矩陣運算庫只包含了矩陣的基本運算，包括創(chuàng)建任意二維/三維矩陣、創(chuàng)建零矩陣及單位矩陣、矩陣加法、矩陣減法、矩陣乘法、矩陣求逆、矩陣轉(zhuǎn)置、矩陣的跡、矩陣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 //分配堆內(nèi)存失敗
#define _ERROR_SVD_EXCEED_MAX_ITERATIONS 0X00000002 //svd超過最大迭代次數(shù)
#define _ERROR_MATRIX_ROWS_OR_COLUMNS_NOT_EQUAL 0X00000003 //矩陣行數(shù)或列數(shù)不相等
#define _ERROR_MATRIX_MULTIPLICATION 0X00000004 //矩陣乘法錯誤(第一個矩陣的列數(shù)不等于第二個矩陣行數(shù))
#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 //創(chuàng)建矩陣失敗
#define _ERROR_MATRIX_TRANSPOSE_FAILED 0X0000000e //矩陣轉(zhuǎn)置失敗
#define _ERROR_CREATE_VECTOR_FAILED 0X0000000f //創(chuàng)建向量失敗
#define _ERROR_VECTOR_DIMENSION_NOT_EQUAL 0X00000010 //向量維數(shù)不相同
#define _ERROR_VECTOR_NORM_TYPE_INVALID 0X00000011 //向量模類型無效
#define _ERROR_VECTOR_CROSS_FAILED 0X00000012 //向量叉乘失敗
#define _ERROR_INPUT_PARAMETERS_ERROR 0X00010000 //輸入?yún)?shù)錯誤
/*******************************************************************************
* (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: 創(chuàng)建矩陣 Input: 矩陣行數(shù)rows，列數(shù)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: 創(chuàng)建多個矩陣 Input: 矩陣行數(shù)rows，列數(shù)columns，個數(shù)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: 創(chuàng)建零矩陣 Input: 矩陣行數(shù)rows，列數(shù)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: 創(chuàng)建單位矩陣 Input: 矩陣行數(shù)rows，列數(shù)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: 矩陣轉(zhuǎn)置 Input: 矩陣A Output: 矩陣A的轉(zhuǎn)置 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: 創(chuàng)建矩陣
Input: 矩陣行數(shù)rows，列數(shù)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非空才能執(zhí)行指針操作
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: 創(chuàng)建多個矩陣
Input: 矩陣行數(shù)rows，列數(shù)columns，個數(shù)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: 創(chuàng)建零矩陣
Input: 矩陣行數(shù)rows，列數(shù)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: 創(chuàng)建單位矩陣
Input: 矩陣行數(shù)rows，列數(shù)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: 矩陣轉(zhuǎn)置
Input: 矩陣A
Output: 矩陣A的轉(zhuǎn)置
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語言優(yōu)雅地實現(xiàn)矩陣運算庫的方法的文章就介紹到這了,更多相關c語言 矩陣運算內(nèi)容請搜索腳本之家以前的文章或繼續(xù)瀏覽下面的相關文章希望大家以后多多支持腳本之家！
相關文章
這篇文章主要介紹了C和C++混合編程問題,需要的朋友可以參考下2015-10-10
C++超詳細講解隱藏私有屬性和方法的兩種實現(xiàn)方式
為了避免因為將類庫中的私有成員開放給類的使用方而導致的軟件邏輯外泄，因此需要將對外代碼中的私有成員隱藏起來，下面我們來了解一下隱藏私有屬性和方法的兩種實現(xiàn)方式2022-05-05
C++實現(xiàn)LeetCode(111.二叉樹的最小深度)
這篇文章主要介紹了C++實現(xiàn)LeetCode(111.二叉樹的最小深度),本篇文章通過簡要的案例,講解了該項技術的了解與使用,以下就是詳細內(nèi)容,需要的朋友可以參考下2021-07-07
今天小編就為大家分享一篇關于C語言數(shù)組和指針的問題一道非常值得深思的筆試題，小編覺得內(nèi)容挺不錯的，現(xiàn)在分享給大家，具有很好的參考價值，需要的朋友一起跟隨小編來看看吧2018-12-12
這篇文章主要介紹了c++ 動態(tài)轉(zhuǎn)換的相關資料,幫助大家更好的理解和使用c++編程，感興趣的朋友可以了解下2020-09-09
這篇文章主要介紹了C++之Qt5雙緩沖機制案例教程,本篇文章通過簡要的案例,講解了該項技術的了解與使用,以下就是詳細內(nèi)容,需要的朋友可以參考下2021-07-07
這篇文章主要介紹了C++11 并發(fā)指南之Lock 詳解，文中通過示例代碼介紹的非常詳細，對大家的學習或者工作具有一定的參考學習價值，需要的朋友們下面隨著小編來一起學習學習吧2020-02-02
這篇文章主要介紹了C語言中的初階基礎，介紹了其相關概念，具有一定參考價值。需要的朋友可以了解下，希望能夠給你帶來幫助2021-11-11
這篇文章主要介紹了C++實現(xiàn)合并排序的方法,實例分析了合并排序的原理與相關實現(xiàn)技巧,需要的朋友可以參考下2015-07-07
這篇文章主要介紹了VC++角色游戲中的人物初始化模塊,對大家學習VC++有一定的幫助,需要的朋友可以參考下2014-08-08