最近在企業實習，花了三個星期的時間肝了一個室內定位演算法（雛形，還未優化），實現matlab和C的聯寫、聯調，其難點筆者認爲在於矩陣的各種運算，C++中有Eigen庫可用，以前在學slam和做課題時候在Ubuntu下用過一段時間的Eigen，功能強大，提供了各種典型的matrix計算介面，但是在C中，需要我們自己編寫每個功能模組。演算法核心程式碼matlab只有大概100行（這裏容我再次感嘆matlab真是爲matrix計算而生，逆天強大），但是C我寫了大概有1500多行。謹以此連載文件記錄一下演算法開發的關鍵–複數矩陣運算的實現。
網上關於C語言實現矩陣運算的demo較多，但是大體都不夠完整和準確，很多function直接git下來各種問題，且用起來沒有統一的標準，另外，只在實數域（Double）內討論，但是我們知道在演算法領域，很多矩陣運算都涉及到複數（Complex）運算，複數域內矩陣運算其實與實數域還是有很大區別的（PS：這也是我在寫C程式時候遇到的坑，之後本文會一一介紹）。
因此，在本連載文章中，我將完全使用C語言完成對矩陣的定義、記憶體初始化、存取索引、記憶體銷燬的基本功能，更重要的是編寫求取代數餘子式、行列式、逆矩陣、協方差矩陣、特徵值、特徵向量等功能函數。所有函數demo我都單獨拿出來仔細講解，並結合matlab，程式設計對比驗證demo正確性，大家可以放心閱讀和使用。
1. 知識儲備：熟練C語言，精通指針、陣列等C語言「靈魂」知識；
2. IDE工具：Visual Studio2019Community; MATLAB2019a;
3. C相關標頭檔案：complex庫
注意到社羣內有小夥伴用二維陣列來定義一個矩陣，但是矩陣並不總是二維的，在matlab裏面，矩陣可以是N維的，若是試圖用更高維度的陣列去完成定義，在後續的存取矩陣/陣列元胞時候顯得非常冗餘。雖說C是結構化程式設計語言，但一定的「魯棒性」也是需要的。另外我是用C++和matlab居多，所以總是或多或少在程式設計時候有一定的C++的思想。後者物件導向，具有類別範本和模板函數。
PS：在之後貼出來的code中，我都會分別給出Double和Complex兩份demo，以及呼叫到的complex庫檔案的相關函數介面說明。
# include <complex.h> // Complex Numbers
# include <stdlib.h> // malloc/free, et.al
說明：
關於複數標頭檔案，C11開始提供了大量函數介面，這裏我只對用到的單獨拿出來講一下。有興趣的可以點開原始檔看看都有什麼內容。首先是定義，原始檔見下：
#ifndef _C_COMPLEX_T
#define _C_COMPLEX_T
typedef struct _C_double_complex
{
double _Val[2];
} _C_double_complex; // This is what I need
typedef struct _C_float_complex
{
float _Val[2];
} _C_float_complex;
typedef struct _C_ldouble_complex
{
long double _Val[2];
} _C_ldouble_complex;
#endif
typedef _C_double_complex _Dcomplex; // This is what I need
typedef _C_float_complex _Fcomplex;
typedef _C_ldouble_complex _Lcomplex;
簡而言之，complex.h標頭檔案通過定義一個結構體來表示一個複數。具體說來，結構體內部成員是一個一維陣列，儲存兩個數據，由該陣列來儲存複數的實部（real）和虛部（imaginary），根據實/虛部浮點精度不同定義了三種陣列，分別是double、float、long double，對應定義了三種結構體型別，再對應每個function都定義了三個函數。謹記，這三種常規型別指的是複數的實部和虛部的數據型別。演算法中，我用到的是double型別。所以以下涉及到複數運算的我都只調用了complex.h標頭檔案中跟double型別相關的函數。
這裏我用一個結構體來完成對矩陣的定義和「模擬」。
typedef _Dcomplex ComplexType; // Complex Variables
typedef double DoubleType; // Double Variables
/* Complex Cell */
typedef struct // Matrix Structor
{
int row, column; // Size
ComplexType* arrayComplex; // Pointer to Cell
}Matrix;
/* Double Cell */
typedef struct
{
int row, column;
DoubleType* arrayDouble;
}Matrix2Double;
/* bool */
typedef enum
{
False = 0, True = 1
}Bool;
說明：
(1) 雖然complex.h標頭檔案內已經對複數型別_C_double_complex取了別名_Dcomplex，爲了區別double和complex，我對二者再次取了一個別名，遵守Google命名規則，隨後的定義矩陣的兩個結構體型別和bool變數的列舉型別同理；
(2) 重點來了，兩種不同元素型別的矩陣，我定義了兩個結構體，一個是元胞是實數型別的，一個是複數型別的。結構體內部兩個int型別的member代表的是矩陣的行數（row）和列數（column），指針member用於指向儲存元素/元胞（cell）的一段記憶體，因此矩陣元胞可以用指針來進行存取和遍歷，此時內部數據結構指針是1D的，但卻可以表示2D及以上的矩陣，從這裏可以看到指針代替陣列的優越性；
(3) 最後，筆者還定義了一個列舉型別，後續的用於安全機制 機製的函數有用到；
(4) 可以看到，兩個不同類型的矩陣其結構體其實完全一樣，只是內部指針變數型別不一樣。
首先是寫一個對一個複數進行初始化的函數，因爲複數也是一個結構體，當定義一個複數變數時候，最好也對其進行初始化。此外，後續涉及到複數的累加累乘，因此也需要對其賦值和初始化：
/* Initiation of Complex Number */
void InitComplex(ComplexType* Complex) // Transfer of Pointer
{
Complex->_Val[0] = 0.0;
Complex->_Val[1] = 0.0;
}
然後是對矩陣進行初始化：
/* Initiation of Complex Matrix */
void InitComplexMatrix(Matrix* matrix, int row, int column) // Transmiss of Pointer
{
int size = row * column * sizeof(ComplexType);
if (size <= 0)
{
puts("ERROE: An invalid matrix!\n");
return;
}
matrix->arrayComplex = (ComplexType*)malloc(size); // initiate pointer
if(matrix->arrayComplex)
{
matrix->row = row; // initiate row and size
matrix->column = column;
for (int index = 0; index < row * column; index++) // initiate cell
{
InitComplex(matrix->arrayComplex + index); // call InitComplex() function
}
}
}
/* Initiation of Double Matrix */
void InitDoubleMatrix(Matrix2Double* matrix, int row, int column)
{
int size = row * column * sizeof(DoubleType);
if (size <= 0)
{
puts("ERROE: An invalid matrix!\n");
return;
}
matrix->arrayDouble = (DoubleType*)malloc(size);
if (matrix->arrayDouble)
{
matrix->row = row;
matrix->column = column;
for (int row_i = 0; row_i < row; ++row_i)
{
for (int column_j = 0; column_j < column; ++column_j)
{
matrix->arrayDouble[row_i * matrix->column + column_j] = 0.0;
}
}
}
}
說明：
(1) 矩陣初始化涉及到確定矩陣row/column的大小、初始化元胞，也就是對結構體變數進行初始化賦值，其核心是對指針動態申請一段記憶體，採用的是malloc()函數，包含在標頭檔案stdlib.h中；
(2) 用指針存取/遍歷（iteration）矩陣元胞，其實也就是指針指向的每個單元內儲存的數，在initiate cell環節，這裏其實有三種表示方式，主要是指針存取元素、矩陣存取元胞具有不同的方式。後續存取元素/元胞時候我都是隨機選一種寫的（sorry，我當時寫程式碼時候這些都不是重點，所以想到哪種寫哪種），但是爲了避免讀者閱讀混亂，這裏，筆者都給出來：
// Solution_1
for (int index = 0; index < row * column; index++) // initiate cell
{
InitComplex(matrix->arrayComplex + index);
}
// Solution_2
for (int index = 0; index < row * column; index++)
{
InitComplex(&(matrix->arrayComplex[index]));
}
// Solution_3
for (int row_i = 0; row_i < row; ++row_i)
{
for (int column_j = 0; column_j < column; ++column_j)
{
InitComplex(matrix->arrayComplex + row_i * matrix->column + column_j);
}
}
(3) 另外，爲了安全起見，我還加入了兩個if語句，可以進一步減少最後compile時出現的各種堆疊溢位等記憶體bug和exception，比如free掉已經不存在的記憶體或者malloc失敗的記憶體，也就是C語言使用指針出現的野指針和堆錯誤；
(4) double型別的矩陣類同complex型別，甚至在初始化值時候還要簡單一點，因爲每個double元胞只有一個元素，而complex有實部和虛部兩個。
我們知道在Linux內，malloc()和free()系統函數總是成對出現，有記憶體申請必有記憶體的釋放，同樣在VS下，利用一個結構體變數定義了一個矩陣之後，呼叫了malloc()函數申請了一段memory，故而需要釋放：
/* Validity of Complex Matrix */
Bool IsNullComplexMatrix(const Matrix* matrix)
{
int size = matrix->row * matrix->column;
if (size <= 0 || matrix->arrayComplex == NULL)
{
return True;
}
return False;
}
/* Free Memory of Complex Matrix */
void DestroyComplexMatrix(Matrix* matrix)
{
if (!IsNullComplexMatrix(matrix)) // Nested Call of IsNullComplexMatrix()
{
free(matrix->arrayComplex); // Pair: malloc--free
matrix->arrayComplex = NULL;
}
matrix->row = matrix->column = 0;
}
/* Validity of Double Matrix */
Bool IsNullDoubleMatrix(const Matrix2Double* matrix)
{
int size = matrix->row * matrix->column;
if (size <= 0 || matrix->arrayDouble == NULL)
{
return True;
}
return False;
}
/* Free Memory of Double Matrix */
void DestroyDoubleMatrix(Matrix2Double* matrix)
{
if (!IsNullDoubleMatrix(matrix))
{
free(matrix->arrayDouble);
matrix->arrayDouble = NULL;
}
matrix->row = matrix->column = 0;
}
說明：
(1) 首先，定義一個IsNullComplexMatrix()函數用於判斷矩陣是否初始化成功，主要判斷結構體內部的member是否初始化成功（體現在矩陣，就是行列數和記憶體）。屬於保障記憶體安全的操作。看到沒，這裏我就呼叫了前文說到的Bool列舉型別成員；
(2) 之後，定義一個DestroyComplexMatrix()函數釋放結構體內部的member，即巢狀呼叫free()釋放掉指針成員指向的記憶體，並置爲NULL（這步雖然不是必須的但是筆者建議各位讀者每次都加上），再將兩外兩個int型別的member變數置爲0；
(3) 關於double型別，同樣給出，不再贅述。
(4) 補充一點，在我的演算法框架中，出現了大量矩陣的定義，也就是多次記憶體的申請和釋放，甚至在幾個核心的演算法函數內部，出現了二十個以上的矩陣，要是每個矩陣呼叫一次初始化和銷燬函數，會嚴重影響程式設計美觀和效率。因此，我又寫了幾個初始化（initiate）和記憶體銷燬（destroy）函數，其目的是一次性只調用一個函數將我需要的所有的矩陣進行記憶體的申請和釋放。當然，這個功能只在我的具體演算法場景下用得到，所以我只是做一個補充拓展貼出來，另外，我只給出記憶體銷燬函數，後續我會依次寫demo驗證一下這個function：
/* Free Memory of Complex Matrice Array */
void DestroyComplexMatrixArray(Matrix matrixArray[], int num) // Array Transfer--->Pointer Transfer
{
if (num) // if no matrix
{
for (int i = 0; i < num; i++)
{
DestroyComplexMatrix(&matrixArray[i]); // Nested Call of DestroyComplexMatrix()
}
}
}
/* Free Memory of Double Matrice Array */
void DestroyDoubleMatrixArray(Matrix2Double* matrixArray, int num)
{
if (num) // if no cell
{
for (int i = 0; i < num; i++)
{
DestroyDoubleMatrix(&matrixArray[i]);
}
}
}
這裏我也定義了三個函數用於返回相應的值：
/* Return Matrix Row Size */
int MatrixRow(const Matrix* matrix)
{
return matrix->row;
}
/* Return Matrix Column Size */
int MatrixColumn(const Matrix* matrix)
{
return matrix->column;
}
/* Return Complex Matrix Size */
int MatrixSize(const Matrix* matrix)
{
return MatrixRow(matrix) * MatrixColumn(matrix) ; // Size Refers to Numbers of Cell of Matrix
}
說明：
(1) 由於這裏的矩陣筆者用的結構體來~~「模擬」~~ 實現，我們爲了讓矩陣看起來」更像「我們數學中那樣的「格式」，數學美學要用計算機實現，還是定義了相關函數來返回我們要的量，因此事實上，可以直接採用如下格式直接獲取：
matrix->row;
matrix->column;
matrix->row * matrix->column;
這是第一章節的內容，整體上來講，只是爲了完成矩陣的定義、初始化、記憶體銷燬等，但是非常重要，因爲後續的所有矩陣運算函數和介面都基於此，事實上，我遇到的很多問題都是出現在這一部分。爲後續方便讀者閱讀，或者幫助讀者程式設計時候排雷，筆者對個人程式設計習慣和這部分遇到的坑做個總結：