C语言实现矩阵求秩和化约化阶梯形

阿表

(｡･ω･｡)

矩阵的秩是反映矩阵固有特性的一个重要概念，在线性代数中有非常重要的地位. 在本文中，笔者将依据自己对矩阵秩定义对理解，利用C语言实现矩阵求秩，并实现矩阵化约化阶梯形.

一、知识储备

• 一个矩阵 \boldsymbol{A}=\left[ \begin{array}{cccc} a_{11} &a_{12} &a_{13} & \cdots &a_{1n}\\ a_{21} & a_{22} & a_{23} & \cdots & a_{2n}\\ a_{31} & a_{32} & a_{33} & \cdots & a_{3n}\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ a_{m1} & a_{m2} & a_{m3} & \cdots & a_{mn}\\ \end{array} \right] 中不为零的子式的最大阶数称为这个矩阵的秩；

• 对矩阵进行初等行变换（包括互换、倍乘和倍加）不改变矩阵的秩，初等列变换也不改变矩阵的秩；下面的实现过程中，我们仅采用初等行变换；

• 求矩阵的秩通常采用这样的方法 ：将矩阵通过初等行（列）变换化为阶梯形，阶梯形矩阵的秩为其不全为零的行的个数，得到阶梯形矩阵的秩和原矩阵等秩，从而求出矩阵的秩.

二、基本思路

• 获取矩阵的行规模 r 和列规模 c ，同时获取矩阵；

• 将矩阵按行依次化为阶梯形，然后转化为约化阶梯形，最后求秩；笔者的思路如下：

（1）按矩阵的行开始循环，在每一行循环开始时， 找到该行第一个不为 0 的元素 \boldsymbol{a_{ij}} ，并利用 \boldsymbol{\rm{Gauss}} 消元法 将该元素所在列 j 上、该元素所在行 i 之后的元素都约为 0 ，这样在完成所有消元之后，我们可以获得一个 伪阶梯矩阵 ，该伪阶梯矩阵一定能通过行互换变为一个阶梯矩阵；

（2）把 伪阶梯矩阵 调整为阶梯矩阵，并把阶梯头化为 1 ，再把每个阶梯头上方均消为 0 ，这样便获得了约化阶梯形矩阵；

（3）统计含有不为 0 元素的行数，这个数就是矩阵的秩 r(\boldsymbol{A}) .

* 对于上文中定义的伪阶梯矩阵，举一个 3 阶方阵 \boldsymbol{B}=\left[ \begin{array}{cccc} a_{11}&a_{12}&a_{13}\\ a_{21}&a_{22}&a_{23}\\ a_{31}&a_{32}&a_{33} \end{array} \right] 来体会：

（1）如果元素 a_{11} 便非零，则将矩阵约化为 \left[ \begin{array}{cccc} a_{11}&a_{12}&a_{13}\\ 0&a_{22}'&a_{23}'\\ 0&a_{32}'&a_{33}' \end{array} \right] ；如果此时的元素 a_{22} 也非零，则矩阵约化为 \left[ \begin{array}{cccc} a_{11}&a_{12}&a_{13}\\ 0&a_{22}'&a_{23}'\\ 0&0&a_{33}'' \end{array} \right] ，矩阵成为了阶梯形；

（2）如果元素 a_{11} 为零而 a_{12} 非零，则找到矩阵第一行第一个非零元素，将矩阵约化为 \left[ \begin{array}{cccc} 0&a_{12}&a_{13}\\ a_{21}&0&a_{23'}\\ a_{31}&0&a_{33'} \end{array} \right] ；如果第二行第一列元素 a_{21} 非零，则矩阵约化为 \left[ \begin{array}{cccc} 0&a_{12}&a_{13}\\ a_{21}&0&a_{23}'\\ 0&0&a_{33}'' \end{array} \right] ，此时只要调换第一、二行，就可得到 \left[ \begin{array}{cccc} a_{21}&0&a_{23}'\\ 0&a_{12}&a_{13}\\0&0&a_{33} ''\end{array} \right] ，也能化为一个阶梯矩阵；

（3）其余状况类同.

三、实现方法

编写主函数 ，实现用户输入矩阵的行规模、列规模和矩阵本身：

int main()
    float matrix[20][20];
    int i,j,r,c;
    printf("请输入矩阵的行规模:");
    scanf("%d",&r);
    printf("请输入矩阵的列规模:");
    scanf("%d",&c);
    printf("请输入矩阵:\n");
    for(i = 0;i < r;i ++)
        for(j = 0;j < c;j ++)
            scanf("%f",&matrix[i][j]);
    printf("矩阵的秩r = %d.\n",rank(matrix,r,c));
    printf("约化阶梯形矩阵为:\n");
    show_standard_echelon(matrix,r,c);
    return 0;
}

将矩阵化简为约化阶梯形、输出阶梯矩阵和求矩阵的秩，我们分别用函数

float standard_echelon(float matrix[20][20],int r,int c,int x,int y);
void show_standard_echelon(float matrix[20][20],int r,int c);
int rank(float matrix[20][20],int r,int c);

来实现.

2. 用 \boldsymbol{\rm{Gauss}} 消元法取得伪阶梯矩阵

说明：接下来的操作均是直接作用在二维数组本身上的. 应当认识到，我们向函数中传入的是数组的第一个元素的地址，因此在函数中改变数组的时候，原先的数组同样会发生变化，一次我们需要在函数中先备份一个二维数组的替身，最后在利用替身把二维数组恢复为传入时的模样. 以上就是向函数传入值和传入指针的区别，需要特别注意.

另外我们约定，函数中生成约化阶梯矩阵后，返回的是申请返回的行列位置上的元素（例如在调用函数时输入了 standard_echelon(matrix,r,c,2,3)，则函数返回其约化阶梯矩阵的第 3 行第 4 列的值而不是整个矩阵）. 之后的矩阵运算函数，均采用这种方式返回，这样可以规避用指针传递二维数组，便于调用时的数据统一；

接下来的第一步，我们用消元法把矩阵转化为伪阶梯形，先来看代码；

float times;
int i,j,r,c;
for(i = 0;i < r - 1;i ++)
    for(k = i + 1;k < r;k ++)
        j = 0;
        while(matrix[i][j] == 0)
            j ++;
        if(matrix[i][j] != 0)
            times = matrix[k][j] / matrix[i][j];
            for(j = 0;j < c;j ++)
                matrix[k][j] -= matrix[i][j] * times;
    }

按照我们之前的基本思路，在每一行的循环开始时，我们都去找该行的第一个非零元素，因此，可以通过

j = 0;
while(matrix[i][j] == 0)
    j ++;

实现非零元素的查找；

找到该行的非零元素之后，需要将该元素下方的全部元素消为0，这时候需要对元素所在行按照不同的倍率（记为 times）实行倍乘，然后减到相应的行中，实现消元. 对该元素所在行之后的行开始循环，这个倍率 times 就等于该元素所在列上其他元素 matrix[k][j] 与该元素 matrix[i][j] 的比值，则：

if(matrix[i][j] != 0)
    times = matrix[k][j] / matrix[i][j];
    for(j = 0;j < c;j ++)
        matrix[k][j] -= matrix[i][j] * times;
}

经过这一步，矩阵就被化为了伪阶梯矩阵. 再通过下面的过程：

for(i = 0;i < r;i ++)
    j = 0;
    while(matrix[i][j] == 0)
        j ++;
    if(matrix[i][j] != 0)
        times = matrix[i][j];
        for(j = 0;j < c;j ++)
            matrix[i][j] /= times;
}

我们就把全部的阶梯头化为了 1.

3. 调整为阶梯形： 要通过行变换将矩阵化为阶梯形，我们考虑到：阶梯矩阵中随着行数增加，每一行在 第一个非零元素出现之前的 0 零元素个数 一定是的增加（除非已经出现了元素全为 0 的行），因此我们只要数出 第一个非零元素出现之前的 0 零元素个数 ，把个数与行数 i 相对应起来，储存在数组 total 中：

int total[20] = {0};
for(i = 0;i < r;i ++)
    for(j = 0;j < c;j ++)
        if(matrix[i][j] == 0)
            total[i] ++;
            break;

total 的下标与行数相同，然后利用类似 冒泡排序法 的基本思路将矩阵按行、以每行对应的 total 值由小到大重新排列. 由于编码需要，引入一个新的量 l 来表示行和列，代码如下；

int l;
float temp;
for(l = r - 1;l > 0;l --)
    for(i = 0;i < l;i ++)
        if(total[l] < total[i])
            for(j = 0;j < c;j ++)
                temp = matrix[l][j];
                matrix[l][j] = matrix[i][j];
                matrix[i][j] = temp;
            }

这样矩阵就被化成了阶梯形.

4. 化为约化阶梯形 ：基本思路不变，在每一行中找到第一个 1（阶梯头），然后把阶梯头上方的元素全部消为 0，就完成了约化，代码如下：

int l;
for(i = 0;i < r;i ++)
    j = 0;
    while(matrix[i][j] == 0)
        j ++;
    if(matrix[i][j] != 0)
        for(k = 0;k < i;k ++)
            times = matrix[k][j] / matrix[i][j];
            for(l = 0;l < c;l ++)
                matrix[k][l] -= times * matrix[i][l];
}

到此为止，本函数已经把矩阵化为了约化阶梯形，结果即为：

result = matrix[x][y];

最后考虑到传递指针可能带来的问题，在操作之前对二维数组进行备份，并在函数最后恢复备份：

float standard_echelon(float matrix[20][20],int r,int c,int x,int y)
    int i,j,k,l,total[20] = {0};
    float times,temp,result = 0,original_matrix[20][20];
    for(i = 0;i < r;i ++)
        for(j = 0;j < c;j ++)
            original_matrix[i][j] = matrix[i][j];
    for(i = 0;i < r - 1;i ++)
        for(k = i + 1;k < r;k ++)
            j = 0;
            while(matrix[i][j] == 0)
                j ++;
            if(matrix[i][j] != 0)
                times = matrix[k][j] / matrix[i][j];
                for(j = 0;j < c;j ++)
                    matrix[k][j] -= matrix[i][j] * times;
    for(i = 0;i < r;i ++)
        j = 0;
        while(matrix[i][j] == 0)
            j ++;
        if(matrix[i][j] != 0)
            times = matrix[i][j];
            for(j = 0;j < c;j ++)
                matrix[i][j] /= times;
    for(i = 0;i < r;i ++)
        for(j = 0;j < c;j ++)
            if(matrix[i][j] == 0)
                total[i] ++;
                break;
    for(l = r - 1;l > 0;l --)
        for(i = 0;i < l;i ++)
            if(total[l] < total[i])
                for(j = 0;j < c;j ++)
                    temp = matrix[l][j];
                    matrix[l][j] = matrix[i][j];
                    matrix[i][j] = temp;
    for(i = 0;i < r;i ++)
        j = 0;
        while(matrix[i][j] == 0)
            j ++;
        if(matrix[i][j] != 0)
            for(k = 0;k < i;k ++)
                times = matrix[k][j] / matrix[i][j];
                for(l = 0;l < c;l ++)
                    matrix[k][l] -= times * matrix[i][l];
    result = matrix[x][y];
    for(i = 0;i < r;i ++)
        for(j = 0;j < c;j ++)
            matrix[i][j] = original_matrix[i][j];
    return result;
}

5. 最后用 rank 函数计算矩阵的秩：

由于 standard_echelon 返回的是约化阶梯形矩阵中用户指定行列位置的元素，在计算秩时还需要先将矩阵还原出来：

float echelon_matrix[20][20];
int i,j;
for(i = 0;i < r;i ++)
    for(j = 0;j < c;j ++)
        echelon_matrix[i][j] = standard_echelon(matrix,r,c,i,j);

接下来，用一个量 none_zero 判断每一行元素是否含有非零元素，初始化其为 0，如果找到一个非零元素，则跳出列循环，并让秩的值 result 增加 1：

int none_zero = 0,result = 0;
for(i = 0;i < r;i ++)
    for(j = 0;j < c;j ++)
        if(echelon_matrix[i][j] != 0)
            none_zero = 1;
            break;
    if(none_zero == 1)
        result ++;
    none_zero = 0;
}

如此一来，便实现了利用约化阶梯形计算矩阵的秩的过程.

把上述过程全部综合起来并适当简化，得到全过程的代码如下：

#include <stdio.h>
#include <math.h>
int rank(float matrix[20][20],int r,int c);
float standard_echelon(float matrix[20][20],int r,int c,int x,int y);
void show_standard_echelon(float matrix[20][20],int r,int c);
int main()
    float matrix[20][20];
    int i,j,r,c;
    printf("请输入矩阵的行规模:");
    scanf("%d",&r);
    printf("请输入矩阵的列规模:");
    scanf("%d",&c);
    printf("请输入矩阵:\n");
    for(i = 0;i < r;i ++)
        for(j = 0;j < c;j ++)
            scanf("%f",&matrix[i][j]);
    printf("矩阵的秩r = %d.\n",rank(matrix,r,c));
    printf("约化阶梯形矩阵为:\n");
    show_standard_echelon(matrix,r,c);
    return 0;
int rank(float matrix[20][20],int r,int c)
    float echelon_matrix[20][20];
    int i,j,none_zero = 0,result = 0;
    for(i = 0;i < r;i ++)
        for(j = 0;j < c;j ++)
            echelon_matrix[i][j] = standard_echelon(matrix,r,c,i,j);
    for(i = 0;i < r;i ++)
        for(j = 0;j < c;j ++)
            if(echelon_matrix[i][j] != 0)
                none_zero = 1;
                break;
        if(none_zero == 1)
            result ++;
        none_zero = 0;
    return result;
float standard_echelon(float matrix[20][20],int r,int c,int x,int y)
    int i,j,k,l,total[20] = {0};
    float times,temp,result = 0,original_matrix[20][20];
    for(i = 0;i < r;i ++)
        for(j = 0;j < c;j ++)
            original_matrix[i][j] = matrix[i][j];
    for(i = 0;i < r - 1;i ++)
        for(k = i + 1;k < r;k ++)
            j = 0;
            while(matrix[i][j] == 0)
                j ++;
            if(matrix[i][j] != 0)
                times = matrix[k][j] / matrix[i][j];
                for(j = 0;j < c;j ++)
                    matrix[k][j] -= matrix[i][j] * times;
    for(i = 0;i < r;i ++)
        j = 0;
        while(matrix[i][j] == 0)
            j ++;
        if(matrix[i][j] != 0)
            times = matrix[i][j];
            for(j = 0;j < c;j ++)
                matrix[i][j] /= times;
    for(i = 0;i < r;i ++)
        for(j = 0;j < c;j ++)
            if(matrix[i][j] == 0)
                total[i] ++;
                break;
    for(l = r - 1;l > 0;l --)
        for(i = 0;i < l;i ++)
            if(total[l] < total[i])
                for(j = 0;j < c;j ++)
                    temp = matrix[l][j];
                    matrix[l][j] = matrix[i][j];
                    matrix[i][j] = temp;
    for(i = 0;i < r;i ++)
        j = 0;
        while(matrix[i][j] == 0)
            j ++;
        if(matrix[i][j] != 0)
            for(k = 0;k < i;k ++)
                times = matrix[k][j] / matrix[i][j];
                for(l = 0;l < c;l ++)
                    matrix[k][l] -= times * matrix[i][l];
    result = matrix[x][y];
    for(i = 0;i < r;i ++)
        for(j = 0;j < c;j ++)
            matrix[i][j] = original_matrix[i][j];
    if(fabs(result) <= 0.0005)
        result = 0;
    return result;
void show_standard_echelon(float matrix[20][20],int r,int c)
    float echelon_matrix[20][20];
    int i,j;
    for(i = 0;i < r;i ++)
        for(j = 0;j < c;j ++)
            echelon_matrix[i][j] = standard_echelon(matrix,r,c,i,j);
    for(i = 0;i < r;i ++)
        for(j = 0;j < c;j ++)
            if(fabs(echelon_matrix[i][j]) < 0.0005)