C语言生成随机可逆方阵

1.前言

最近做并行计算作业的时候有一道题是让用并行的方式对一个矩阵求逆，这个实验的大致步骤是将一个写好的矩阵文件（一定格式）作为输入，使用一定的算法求出逆矩阵后再以文件的形式输出。因为在使用并行方式进行数值计算时，如果数据规模不够大，将很难体现并行方式的优越性，因此，拥有一个规模较大数据集就成为了这个实验成功的关键点之一。

关于数据集这件事儿，大神们纷纷使用matlab造出了随机二维数组文件，然而对于不会使用matlab的我（数学渣渣一枚）来说，想要得到一个大规模的矩阵文件真的好难啊！

2.关键问题与解决思路

有的同鞋可能会说：随机生成个二维数组还不简单么，几行代码的事嘛。但是请注意，这里我们需要的矩阵必须是随机且可逆的！矩阵可逆是要满足一定条件的，然而随机生成的矩阵可能压根就不是一个可逆矩阵，从而导致逆矩阵无法求出：

输入是这样的：

用它求逆后就变成这样了：

可以看到我的求逆程序报错了，显示not a number，这是因为这个随机生成的数组是不可逆的。因此我就想编写一个能生成任意大小方阵的程序。

首先我们要生成一个任意大小的可逆方阵，再把它按照一定的格式写入文件中。其实难点只有一个，就是如何构造可逆矩阵。我在网上找到了两种解决办法：

（1）生成随机数矩阵后，再根据行列式值或者秩等条件判断其是否可逆。

（2）先得到单位矩阵后在对其进行有限次的初等行变换。

运气不好的话，第一种方式的开销会很大，因此我们选择第二种方式。但是，在第二种方式中也存在着如下问题：

（1）如何进行初等行变换才能使最后的结果正确

（2）如何使初等变换后的矩阵看起来更匀称

下面我来解释这两点是啥意思。所谓正确就是在计算时没有出现越界，溢出等情况，因为当矩阵中的数字比较大的时候，进行初等行变换的结果可能会使某些数字过大从而超过计算机能表示的精度导致溢出；所谓匀称是指一个规模很大的矩阵中不会出现全是零或者某一大片扎堆是零另一部分非零。当然对于矩阵求逆来说，第一点是必须要满足的，第二点我们暂且作为功能性需求。

好了，下面我来介绍一下我的思路。

当我们得到单位矩阵之后，就要考虑对其进行初等行变换了，无非就是交换两行位置，将一行中所有的元素同乘以一个数，还有就是将一行中所有的元素同乘一个数后再加到另一行上，我的思路是，随机的选取一个主行（行数当然要在0~n-1之间，n是方阵的行列数），然后进行随机次数的行变换，但是行变换我仅采用“将一行中所有的元素乘以一个数之后再加到另一行上”这种方式，在每次初等行变换中执行乘法和加法之前要检测结果是否会溢出。

3.编码实施

3.1 程序整体框架

这本程序中，首先会提示用户输入方阵的行列数，然后就会为矩阵（二维数组）申请内存，这里我采用了动态数组的方式进行存储，要注意的是array本身是一个int**类型，即指针的指针，在array的第一维中存储的是一系列的int型指针，即int*类型的变量，而在其第二维中存储的才是int类型的变量。矩阵的内存申请完毕后，用户可以选择仅生成单位矩阵（输入0）还是非单位矩阵的可逆矩阵（输入1），由于生成非单位矩阵的可逆矩阵的逻辑包含生成单位矩阵的逻辑，因此我们简述生成非单位矩阵的可逆矩阵的过程：

（1）生成单位矩阵

（2）根据单位矩阵生成非单位矩阵的可逆矩阵

（3）将生成的矩阵写入文件

这几步在下面几个小节中详细介绍。

int main(int argc, char* argv[])
{
	int k = 0;
        int n = 0;
        int flag = 0;
        int** array;

        printf("please input the line and row number of the matrix:\n");
        scanf("%d",&n);

	//create matrix
        printf("start create a %d * %d matrix.\n",n,n);
        array = (int**)malloc(sizeof(int*)*n);

        for(k = 0; k < n; ++k)
        {
                array[k] = (int*)malloc(sizeof(int)*n);
        }

        printf("if you want get an identity matrix, please input 0, while if you want an invertible matrix, just input 1.\n");
        scanf("%d",&flag);

	if(flag == 0)
        {
                printf("your input word is %d, so you want an identity matrix:\n",flag);
                //get indentity matrix
                getidentitymatrix(n, array);

                //put indentity matrix into file
                putidentitymatrixintofile(n, array);

                //print identity matrix
                printidentitymatrix(n, array);
        }
	else if(flag == 1)
        {
                printf("your input word is %d, so you want an invertible matrix:\n",flag);
                //get indentity matrix
                getidentitymatrix(n, array);
                //get invertible matrix
                getinvertiblematrix(n, array);
                //put invertible matrix into file
                putinvertiblematrixintofile(n, array);
                //print invertible matrix
                printinvertiblematrix(n, array);
        }
	else
        {
                printf("error: you input a wrong number!\n");
                return -1;
        }

        //free matrix
        printf("free matrix.\n");
        freematrix(n, array);

        return 1;
}

注：代码中最后出现的freematrix函数是我专门为了释放数组内存写的，代码如下：

 

void freematrix(int n, int** array)
{
        int k = 0;
        for(k = 0; k < n; ++k)
        {
                free(array[k]);
        }
        free(array);
}

因为前面是使用动态数组的形式进行内存分配的，再加上是一个二维数组，所以在回收内存的时候还是略有复杂的。

 

3.2 生成单位矩阵

生成单位矩阵的代码逻辑很简单，就是行数和列数相等的位置上置1，其余的置0就可以了。

void getidentitymatrix(int n, int** array)
{
        int r = 0;
        int c = 0;

        for(r = 0; r < n; ++r)
        {
                for(c = 0; c < n; ++c)
                {
                        if(r == c)
                                array[r][c] = 1;
                        else
                                array[r][c] = 0;
                }
        }

}

3.3 生成可逆矩阵

生成可逆矩阵的步骤相对较复杂，要注意的有两点：

（1）初等行变换的次数随机（但不可无限制）

（2）初等行变换时的加法运算不可以溢出（因此要判断计算结果的合法性）

由于程序中需要用到随机数，我们首先使用srand函数设置随机数的种子，然后再使用rand函数生成我要进行初等行变换的次数（这里我将其设定在1000之内，如果你想了解关于更多随机数的东西，可以参看这篇文章：）；为了实现初等行变换中的加法，我需要一个中介数组存储矩阵中某一行中的元素，因此我采用动态的方法构建了一个名额日temparray的一维数组。接下来重头戏开始了，下面要进行一些循环，各位看官且听我慢慢道来。

最外圈的循环所控制的是初等行变换的次数，这个主循环里套着两个小循环：

（1）第一个循环开始之前有一句：

mainrownum = (int)(rand()%（n-1）);

这句话是随机的在矩阵0~n-1行中选择一个主行作初等行变换，第一个循环是将主行中的元素进行一定的处理

array[mainrownum][k])*((int)(rand()%5 - 10)

后存入中介数组temparray中（处理方式是我随便写的，也可以是别的运算方式）。

而这一句：

((uint16_max - (array[mainrownum][k])*((int)(rand()%5 - 10))) < 0) || ((uint16_max*(-1)) - (array[mainrownum][k])*((int)(rand()%5 - 10)) > temparray[k])

则是为了判断采用上述处理方式对矩阵元素进行运算后，元素数值是否会溢出。（如果你想了解更多，请看这里：）

第二个循环的作用是遍历矩阵中所有行，然后令处理过后主行和其他所有的行依次进行相加，这里也进行了溢出判断。

最后不要忘了释放内存哟~

void getinvertiblematrix(int n, int** array)
{
        int i = 0;
        int j = 0;
        int k = 0;
	int mainrownum = 0;

        int* temparray = null;

	srand((int)time(null));
	int transformtime = (int)(rand()%1000);
	printf("we will do %d times tansformation.\n",transformtime);

        temparray = (int*)malloc(sizeof(int)*n);

        for(i = 0; i < transformtime; ++i)
        {
		mainrownum = (int)(rand()%（n-1）);
                for(k = 0; k < n; ++k)
			if(((uint16_max - (array[mainrownum][k])*((int)(rand()%5 - 10))) < 0) || ((uint16_max*(-1)) - (array[mainrownum][k])*((int)(rand()%5 - 10)) > temparray[k]))
				temparray[k] = (array[mainrownum][k]);
			else
                        	temparray[k] = (array[mainrownum][k])*((int)(rand()%5 - 10));

                for(j = 0; j < n; ++j)
                        if(mainrownum != j)
                                for(k = 0; k < n; ++k)
				{
					if(((uint16_max - array[j][k]) < temparray[k]) || ((uint16_max*(-1)) - array[j][k] > temparray[k]))
						array[j][k] = array[j][k]/4;
					else
						array[j][k] = array[j][k] + temparray[k];
				}
        }

        free(temparray);
}

 

3.4 将可逆矩阵写入文件中

可逆矩阵算好之后，最后一步就是要将其存入文件中了，这里主要涉及到的是文件操作，其实本质还是矩阵输出，只不过这回是将元素输出到文件中罢了。

int putinvertiblematrixintofile(int n, int** array)
{
        file* fp = null;
        int i = 0;
        int j = 0;

        if((fp = fopen("input","w"))==null)
        {
                printf("error: writing file error!\n");
                return -1;
        }

        for(i = 0; i < n; ++i)
        {
                for(j = 0; j < n; ++j)
                {
                        if(j != (n-1))
                                fprintf(fp,"%d\t", array[i][j]);
                        else
                                fprintf(fp,"%d", array[i][j]);
                }
                fputs("\n",fp);
        }

        fclose(fp);

        return 1;

}

注意在输出时每行最后一个元素后面不应该有空格或者tab之类的分割，因此要做判断。对于每行中非末尾元素来说，最好使用\t，如果单纯使用空格的话，在数字很大而且既有负数又有正数的时候看起来会很乱。

4.实验结果

我是在64位的ubuntu 14.04上做的，因此使用gcc编译

gcc -o invertible invertiblematrix.c

我们以5*5的矩阵为例试一下

选择1，让其生成非单位矩阵的可逆矩阵

可以看出，做了780次行变换，只是矩阵有一列有很多0，看起来不太完美。

用vim看一下我们得到的文件（名为input）：

初步成功~

5. 小结

作为一个数学渣渣，这篇博文是我的一点愚见，算法并不完美，有很多瑕疵，还请大神们多多指正，如果有更好的方法也可以一同讨论并分享给更多的人。