【用 C 语言实现复数稀疏矩阵的基本运算】

前言

在复习电力系统分析时，看到n年前清华大学PPT上比较有意思的思考题。这个思考题如下：

设计一种稀疏矩阵的存储方案，用 C 语言编程实现大电网导纳矩阵的排零存储和检索；针对采用的存储方案，设计算例，编写稀疏矩阵的加、减、乘、求逆等基本运算的程序，实现排零运算，分析稀疏技术对矩阵运算所需的存储空间、计算速度的改善情况。

这个问题比较有趣，大电网导纳矩阵本身就是一个复数稀疏对称矩阵。这里的对称并不是Hermitian Matrixs,而是对角元素相等。第二，就是C语言比较底层，并不像高级科学计算语言那样方便，实现这个功能比较有挑战性。之所以使用C语言，是因为自己积累了些数据结构于算法的部分知识，这些知识在设计C程序提供了一种有效的方法论。

稀疏矩阵的存储方式

1.三元组顺序表，存储稀疏矩阵的非零元和对应行和列的位置$(i,j)$。反之，一个三元组$(i,j,a_{ij})$惟一确定了矩阵A的一个非零元。由此,稀疏矩阵可由表示非零元的三元组及其行列数惟一确定。
2.十字链表，在链表中,每个非零元可用一个含5个域的结点表示,其中$i$、$j$和$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\times 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$行展开，对应元素与其代数余子式乘积的代数和。公式为：
$$D=\sum _{j=1}^n{a}_{ij}{A}_{ij}$$
对于一阶矩阵来说，行列式就等于矩阵的元素，因此可以使用递归的方法求解矩阵的行列式。
核心代码

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$的数，对于成功完成消元的方阵来说，主元在上三角阵的对角线上。
③消去$（i,j）$位置元素的乘数记为
$$l_{ij}=\frac{\text{待消去的}\left(i,j\right)\text{位置的元素}}{\left(j,j\right)\text{位置的主元}}$$
④消去$(i,j)$位置的元素用公式：
$$\text{第}i\text{行}-l_{ij}\times \text{第}j\text{行,}i>j$$
⑤$A{x}_j=e_j$中右边的$e_j$的元素也要参与消元，这样才能保证方程的解不变。
⑥消元的过程是自上而下的，求解的过程是自下而上的。

核心代码

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\in {\mathbb{R}}^{n\times n}$的所有顺序主子矩阵都非奇异，则$A$存在单位下三角矩阵L和单位上三角矩阵$U$,以及非奇异对角矩阵$D$，使得$A=LDU$,其中$L$、$D$、$U$都是唯一的，反之，$A$存在$LDU$分解，则$A$的所有顺序主子矩阵都非奇异。
各因子矩阵元素表达如下：
$$\{\begin{array}{l}{d}_{ii}=a_{ii}-{\sum }_{k=1}^{i-1}{l}_{ik}{u}_{ki}{d}_{kk}\\ u_{ij}=\frac{\left(a_{ij}-{\sum }_{k=1}^{i-1}{l}_{ik}{u}_{ki}{d}_{kk}\right)}{d_{ii}}\\ l_{ij}=\frac{\left(a_{ij}-{\sum }_{k=1}^{j-1}{l}_{ik}{u}_{kj}{d}_{kk}\right)}{d_{jj}}\end{array}$$

若A是对称矩阵，则有
$$A^T=A=LDU=U^T{D}^{T}LT$$
当$A$的各阶主子式均不为零时，根据分解的唯一性，应有$L^T=U$或$U^T=L$,各因子的元素u和l的关系$u_{ij}=l_{ji}$
最后分段求解即可。
核心代码

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