用 C 语言实现复数稀疏矩阵的基本运算
前言
在复习电力系统分析时,看到n年前清华大学PPT上比较有意思的思考题。这个思考题如下:
设计一种稀疏矩阵的存储方案,用 C 语言编程实现大电网导纳矩阵的排零存储和检索;针对采用的存储方案,设计算例,编写稀疏矩阵的加、减、乘、求逆等基本运算的程序,实现排零运算,分析稀疏技术对矩阵运算所需的存储空间、计算速度的改善情况。
这个问题比较有趣,大电网导纳矩阵本身就是一个复数稀疏对称矩阵。这里的对称并不是Hermitian Matrixs,而是对角元素相等。第二,就是C语言比较底层,并不像高级科学计算语言那样方便,实现这个功能比较有挑战性。之所以使用C语言,是因为自己积累了些数据结构于算法的部分知识,这些知识在设计C程序提供了一种有效的方法论。
稀疏矩阵的存储方式
1.三元组顺序表,存储稀疏矩阵的非零元和对应行和列的位置
(
i
,
j
)
(i,j)
(i,j)。反之,一个三元组
(
i
,
j
,
a
i
j
)
(i,j, a_{ij})
(i,j,aij)惟一确定了矩阵A的一个非零元。由此,稀疏矩阵可由表示非零元的三元组及其行列数惟一确定。
2.十字链表,在链表中,每个非零元可用一个含5个域的结点表示,其中
i
i
i、
j
j
j和
e
e
e这3个域分别表示该非零元所在的行,列和非零元的值,向右域right用以链接同一行中下一个非零元,向下域down用以链接同一列中下一个非零元。同一行的非零元通过right域链接成一个线性链表,同一列的非零元通过down域链接成一个线性链表,每个非零元既是某个行链表中的一个结点,又是某个列链表中的一个结点,整个矩阵构成了一个十字交叉的链表。
定义数据结构
1,复数的数据结构首选操作对象是复数,然而在CodeBlocks这个IDE当中,complex.h库当中,直接使用不了。直接复制一篇引用到complex.h库的C代码,直接给我报错了。也许是编译器的问题。因此那就定义Complex这个类吧,复数剧本功能的实现并不复杂。
typedef struct
{
float Realpart; // 实部
float Imagpart; // 虚部
}Complex;//自定义复数的基功能
Complex Create(float x,float y);
float GetReal(Complex C);
float GetImag(Complex C);
float Abs_complex(Complex C);
bool Equal(Complex c1, Complex c2);
Complex Add(Complex C1,Complex C2);
Complex Sub(Complex C1,Complex C2);
Complex Mul(Complex C1,Complex C2);
Complex Div(Complex C1,Complex C2);
void PrintC(Complex C);
2,稀疏矩阵的数据结构
/* 三元组稀疏矩阵元素类型 */
typedef Complex ElemType;
/* 三元组类型定义,主要用来存储非零元 */
typedef struct {
int i, j; // 该三元组非零元的行下标和列下标
ElemType e;
} Triple;
/* 三元组稀疏矩阵类型定义 */
typedef struct {
Triple data[MAXSIZE + 1]; // 非零元三元组表,data[0]未用
int mu, nu, tu; // 矩阵的行数、列数和非零元个数
} TSMatrix;
3,辅助定义的数据结构,这一部分用来求解稀疏矩阵的行列式和逆矩阵。定义的数据类似于MATLAB的full函数。
typedef struct {
int rows; // 矩阵的行数
int cols; // 矩阵的列数
Complex* data; // 矩阵元素的数组
} matrix_t;
稀疏矩阵乘法
两个复数稀疏矩阵相乘 ( Q = M × N ) (Q=M\times N) (Q=M×N)可以描述为:
Q初始化;
if (是非零矩阵) { // 逐行求积
for (arow = 1; arow <= M.mu; ++arow) { // 处理M的每一行
ctemp[] = 0;
// 累加器清零
计算中第arow行的积并存入ctemp[]中;
将ctemp[]中非零元压缩存储到Q.data;
} // for arow
} // if
对于稀疏矩阵的快速转置,很多的博客都有介绍。比如:
链接:稀疏矩阵的快速转置
求解复数稀疏矩阵的行列式
在求解之前,得先把三元组结构的复数稀疏矩阵转换成普通的矩阵。然后再matrix_t数据类型的矩阵定义函数求解行列式。
主要的方法有拉普拉斯展开定理(递归法)、LU分解法等等。这里只介绍递归法。拉普拉斯展开定理在线性代数也称作行列式按行或列展开。
一个n阶矩阵的行列式等于其按第
i
i
i行展开,对应元素与其代数余子式乘积的代数和。公式为:
D
=
∑
j
=
1
n
a
i
j
A
i
j
D=\sum_{j=1}^n{a_{ij}A_{ij}}
D=j=1∑naijAij
对于一阶矩阵来说,行列式就等于矩阵的元素,因此可以使用递归的方法求解矩阵的行列式。
核心代码
Complex determinant(matrix_t* m) {
int n = m->rows;
if(n == 1) {
return MAT(m, 0, 0, n);
}
Complex sum = Create(0.0, 0.0);
for(int k = 0; k < n; ++k) { // 对第一行进行展开,计算每个代数余子式的值,并累加到sum中
Complex factor = Mul(MAT(m, 0, k, n), Create(pow(-1, k), 0.0)); // 计算代数余子式的符号和系数
matrix_t* submatrix = matrix_create(n - 1, n - 1); // 创建一个子矩阵来存储余子式
for(int i = 1; i < n; ++i) { // 将原矩阵除去第一行和第k列的元素复制到子矩阵中
for(int j = 0, l = 0; j < n; ++j) {
if(j != k) {
MAT(submatrix, i - 1, l, submatrix->cols) = MAT(m, i, j, n);
l++;
}
}
}
Complex subdet = determinant(submatrix); // 递归地计算子矩阵的行列式值,即余子式的值
matrix_free(submatrix);
sum = Add(sum, Mul(factor, subdet));
}
return sum;
}
求解复数稀疏矩阵的逆矩阵
对于大电网导纳矩阵,它的逆矩阵对应着阻抗矩阵。求逆后,不保证矩阵的稀疏性。求解的方法也有很多。高斯-若当消元法等等,对于导纳矩阵,可以使用三角分解的方法求解逆矩阵。
首先,高斯-若当消元法对应于线性代数用增广矩阵,经过初等行变换和列变换,求解逆矩阵。
消元法的总结:
①消元的顺序是一列完成之后再进行下一列的,顺序不能乱。
②主元是每行中第一个不为
0
0
0的数,对于成功完成消元的方阵来说,主元在上三角阵的对角线上。
③消去
(
i
,
j
)
(i, j)
(i,j)位置元素的乘数记为
l
i
j
=
待消去的
(
i
,
j
)
位置的元素
(
j
,
j
)
位置的主元
l_{ij}=\frac{\text{待消去的}\left( i,j \right) \text{位置的元素}}{\left( j,j \right) \text{位置的主元}}
lij=(j,j)位置的主元待消去的(i,j)位置的元素
④消去
(
i
,
j
)
(i,j)
(i,j)位置的元素用公式:
第
i
行
−
l
i
j
×
第
j
行,
i
>
j
\text{第}i\text{行}-l_{ij}\times \text{第}j\text{行,}i>j
第i行−lij×第j行,i>j
⑤
A
x
j
=
e
j
Ax_{j}=e_{j}
Axj=ej中右边的
e
j
e_{j}
ej的元素也要参与消元,这样才能保证方程的解不变。
⑥消元的过程是自上而下的,求解的过程是自下而上的。
核心代码
int inverse_matrix(matrix_t* a, int n, matrix_t* b) {
int i, j, k; // 循环变量
int max_row; // 最大主元所在的行号
Complex temp; // 临时变量,用于交换和计算
Complex factor; // 消元因子
// 检查输入参数是否合法,即a和b是否都是n*n的方阵
if (a == NULL || b == NULL || a->rows != n || a->cols != n || b->rows != n || b->cols != n) {
return 0; // 参数不合法,返回0表示失败
}
// 将b初始化为单位矩阵,即对角线元素为1,其他元素为0
for (i = 0; i < n; i++) {
for (j = 0; j < n; j++) {
if (i == j) { // 对角线元素为1
MAT(b, i, j, n) = Create(1.0, 0.0);
} else { // 其他元素为0
MAT(b, i, j, n) = Create(0.0, 0.0);
}
}
}
// 对a进行高斯消元,同时对b进行同样的初等行变换,使得a变为上三角矩阵
for (k = 0; k < n - 1; k++) { // 遍历每一列
// 寻找第k列的最大主元所在的行号
max_row = k;
for (i = k + 1; i < n; i++) {
if (Abs_complex(Mul(MAT(a, i, k, n), MAT(a, i, k, n))) > Abs_complex(Mul(MAT(a, max_row, k, n), MAT(a, max_row, k, n)))) {
max_row = i;
}
}
// 如果最大主元为0,说明矩阵奇异,无法求逆
if (Abs_complex(Mul(MAT(a, max_row, k, n), MAT(a, max_row, k, n))) == 0.0) {
return 0; // 返回0表示失败
}
// 如果最大主元不在第k行,就交换第k行和第max_row行
if (max_row != k) {
for (j = 0; j < n; j++) {
// 交换a的两行元素
temp = MAT(a, k, j, n);
MAT(a, k, j, n) = MAT(a, max_row, j, n);
MAT(a, max_row, j, n) = temp;
// 交换b的两行元素
temp = MAT(b, k, j, n);
MAT(b, k, j, n) = MAT(b, max_row, j, n);
MAT(b, max_row, j, n) = temp;
}
}
// 用第k行的主元消去下面各行的第k列元素
for (i = k + 1; i < n; i++) {
// 计算消元因子
factor = Div(MAT(a, i, k, n), MAT(a, k ,k ,n));
// 更新a的各行元素
for (j = k; j < n; j++) {
MAT(a ,i ,j ,n) = Sub(MAT(a ,i ,j ,n), Mul(factor ,MAT(a ,k ,j ,n)));
}
// 更新b的各行元素
for (j = 0; j < n; j++) {
MAT(b ,i ,j ,n) = Sub(MAT(b ,i ,j ,n), Mul(factor ,MAT(b ,k ,j ,n)));
}
}
}
// 对a进行回代,同时对b进行同样的初等行变换,使得a变为单位矩阵
for (k = n - 1; k > 0; k--) { // 遍历每一列
// 如果对角线元素为0,说明矩阵奇异,无法求逆
if (Abs_complex(Mul(MAT(a ,k ,k ,n), MAT(a ,k ,k ,n))) == 0.0) {
return 0; // 返回0表示失败
}
// 用第k行的主元消去上面各行的第k列元素
for (i = 0; i < k; i++) {
// 计算消元因子
factor = Div(MAT(a ,i ,k ,n), MAT(a ,k ,k ,n));
// 更新a的各行元素
for (j = k; j >= 0; j--) {
MAT(a ,i ,j ,n) = Sub(MAT(a ,i ,j ,n), Mul(factor ,MAT(a ,k ,j ,n)));
}
// 更新b的各行元素
for (j = 0; j < n; j++) {
MAT(b ,i ,j ,n) = Sub(MAT(b ,i ,j ,n), Mul(factor ,MAT(b ,k ,j ,n)));
}
}
}
// 对角化a,同时对b进行同样的初等行变换,使得b变为逆矩阵
for (i = 0; i < n; i++) { // 遍历每一行
// 如果对角线元素为0,说明矩阵奇异,无法求逆
if (Abs_complex(Mul(MAT(a,i,i,n),MAT(a,i,i,n))) == 0.0) {
return 0; // 返回0表示失败
}
// 计算对角线元素的倒数
temp = Div(Create(1.0, 0.0), MAT(a,i,i,n));
// 更新a的对角线元素为1
MAT(a,i,i,n) = Create(1.0, 0.0);
// 更新b的各行元素
for (j = 0; j < n; j++) {
MAT(b,i,j,n) = Mul(temp, MAT(b,i,j,n));
}
}
return 1;
}
导纳矩阵的三角分解法(LDU)
矩阵的LDU分解来源于LU分解,LDU分解定理如下:
设
A
∈
R
n
×
n
A\in \mathbb{R}^{n\times n}
A∈Rn×n的所有顺序主子矩阵都非奇异,则
A
A
A存在单位下三角矩阵L和单位上三角矩阵
U
U
U,以及非奇异对角矩阵
D
D
D,使得
A
=
L
D
U
A=LDU
A=LDU,其中
L
L
L、
D
D
D、
U
U
U都是唯一的,反之,
A
A
A存在
L
D
U
LDU
LDU分解,则
A
A
A的所有顺序主子矩阵都非奇异。
各因子矩阵元素表达如下:
{
d
i
i
=
a
i
i
−
∑
k
=
1
i
−
1
l
i
k
u
k
i
d
k
k
u
i
j
=
(
a
i
j
−
∑
k
=
1
i
−
1
l
i
k
u
k
i
d
k
k
)
d
i
i
l
i
j
=
(
a
i
j
−
∑
k
=
1
j
−
1
l
i
k
u
k
j
d
k
k
)
d
j
j
\left\{ \begin{array}{l} d_{ii}=a_{ii}-\sum_{k=1}^{i-1}{l_{ik}u_{ki}d_{kk}}\ \ \ \ \ \ \ \ \ \\\ u_{ij}=\dfrac{\left( a_{ij}-\sum_{k=1}^{i-1}{l_{ik}u_{ki}d_{kk}} \right)}{d_{ii}}\ \ \\\ l_{ij}=\dfrac{\left( a_{ij}-\sum_{k=1}^{j-1}{l_{ik}u_{kj}d_{kk}} \right)}{d_{jj}}\\ \end{array} \right.
⎩
⎨
⎧dii=aii−∑k=1i−1likukidkk uij=dii(aij−∑k=1i−1likukidkk) lij=djj(aij−∑k=1j−1likukjdkk)
若A是对称矩阵,则有
A
T
=
A
=
L
D
U
=
U
T
D
T
L
T
A^{T}=A=LDU=U^{T}D^{T}L{T}
AT=A=LDU=UTDTLT
当
A
A
A的各阶主子式均不为零时,根据分解的唯一性,应有
L
T
=
U
L^{T}=U
LT=U或
U
T
=
L
U^{T}=L
UT=L,各因子的元素u和l的关系
u
i
j
=
l
j
i
u_{ij}=l_{ji}
uij=lji
最后分段求解即可。
核心代码
int inverse_SymmetricMatrix(matrix_t* a, int n, matrix_t* b) {
matrix_t* U = matrix_create(n, n);
matrix_t* D = matrix_create(n, n);
int i, j, k; // 循环变量
Complex sum; // 累加和
Complex temp; // 临时变量
if (a == NULL || b == NULL || a->rows != n || a->cols != n || b->rows != n || b->cols != n) {
return 0; // 参数不合法,返回0表示失败
}
// 对a进行分解,使得a = U^T * D * U
for (i = 0; i < n; i++) { // 遍历每一行
// 先计算D(i,i)
sum = Create(0.0f, 0.0f);
for (k = 0; k < i; k++) {
sum = Add(sum, Mul(Mul(MAT(U, k, i, n), MAT(D, k, k, n)), MAT(U, k, i, n))); // 累加U(k,i) * D(k,k) * U(k,i)
}
MAT(D, i, i, n) = Sub(MAT(a, i, i, n), sum); // 计算D(i,i) = a(i,i) - sum
// 再计算U(i,j)
for (j = i + 1; j < n; j++) {
sum = Create(0.0f, 0.0f);
for (k = 0; k < i; k++) {
sum = Add(sum, Mul(Mul(MAT(U, k, i, n), MAT(D, k, k, n)), MAT(U, k ,j ,n))); // 累加U(k,i) * D(k,k) * U(k,j)
}
MAT(U ,i ,j ,n) = Div(Sub(MAT(a ,i ,j ,n), sum), MAT(D ,i ,i ,n)); // 计算U(i,j) = (a(i,j) - sum) / D(i,i)
}
}
matrix_t* x = matrix_create(n ,n); // 用于存储U^T * x = I 的解
matrix_t* w = matrix_create(n ,n); // 用于存储w * D = x 的解
matrix_t* z = matrix_create(n ,n); // 用于存储U * z = w 的解
// 求解U^T * x = I
for (j = 0; j < n; j++) {
for (i = 0; i < n; i++) {
if (i < j) {
MAT(x ,i ,j ,n) = Create(0.0f ,0.0f);
} else if (i == j) {
MAT(x ,i ,j ,n) = Create(1.0f ,0.0f);
} else {
sum = Create(0.0f ,0.0f);
for (k = j; k < i; k++) {
sum = Add(sum, Mul(MAT(U ,k ,i ,n), MAT(x ,k ,j ,n)));
}
MAT(x ,i ,j ,n) = Sub(Create(0.0f ,0.0f), sum); // 计算x(i,j) = -sum
}
}
}
// 求解w * D = x
for (j = 0; j < n; j++) {
for (i = 0; i < n; i++) {
if (i < j) {
MAT(w ,i ,j ,n) = Create(0.0f, 0.0f); // 如果i < j,那么w(i,j) = 0
} else {
MAT(w, i, j, n) = Div(MAT(x, i, j, n), MAT(D, i, i, n)); // 如果i >= j,那么w(i,j) = x(i,j) / D(i,i)
}
}
}
// 求解U * z = w
for (j = 0; j < n; j++) {
for (i = n - 1; i >= 0; i--) {
sum = Create(0.0f, 0.0f);
for (k = i + 1; k < n; k++) {
sum = Add(sum, Mul(MAT(U, i, k, n), MAT(z, k, j, n)));
}
MAT(z, i, j, n) = Sub(MAT(w, i, j, n), sum); // 计算z(i,j) = w(i,j) - sum
}
}
// 将z复制到b,作为逆矩阵的结果
for (i = 0; i < n; i++) {
for (j = 0; j < n; j++) {
MAT(b, i, j, n) = MAT(z, i, j, n);
}
}
matrix_free(U);
matrix_free(D);
matrix_free(x);
matrix_free(w);
matrix_free(z);
return 1;
}
算例结果
分别用C语言计算的结果于MATLAB计算的结果比较
[1]:https://blog.csdn.net/weixin_43172803/article/details/103414500
[2]:https://zhuanlan.zhihu.com/p/123314715?utm_id=0