C语⾔实现实数和复数矩阵及其各种运算(⼀)
⼀、前⾔
本连载⽂章主要内容是实现复数矩阵的各种运算,并利⽤matlab进⾏联写、联调,验证C语⾔编写的各个矩阵运算函数的正确性。其难点笔者认为在于矩阵的各种运算,C++中有Eigen库可⽤,以前在学slam和做课题时候在Ubuntu下⽤过⼀段时间的Eigen,功能强⼤,提供了各种典型的matrix计算接⼝,但是在C中,需要我们⾃⼰编写每个功能模块。谨以此连载⽂档记录⼀下算法开发的关键--复数矩阵运算的实现,纯粹为了学习交流。
⽹上关于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库⽂件的相关函数接⼝说明。
1.头⽂件
# 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类型相关的函数。
2.定义
这⾥我⽤⼀个结构体来完成对矩阵的定义和“模拟”。
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) 可以看到,两个不同类型的矩阵其结构体其实完全⼀样,只是内部指针变量类型不⼀样。
3.初始化内存分配
⾸先是写⼀个对⼀个复数进⾏初始化的函数,因为复数也是⼀个结构体,当定义⼀个复数变量时候,最好也对其进⾏初始化。此外,后续涉及到复数的累加累乘,因此也需要对其赋值和初始化:
/* 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有实部和虚部两个。
4.内存销毁
我们知道在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;
index复数}
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:

版权声明:本站内容均来自互联网,仅供演示用,请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198,我们将在24小时内删除。