1. 程式人生 > >C語言生成隨機可逆方陣

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. 小結

作為一個數學渣渣,這篇博文是我的一點愚見,演算法並不完美,有很多瑕疵,還請大神們多多指正,如果有更好的方法也可以一同討論並分享給更多的人。 最後,我已經把這個程式的原始碼傳到了我的資源中:http://download.csdn.net/detail/liu1075538266/9529058。感謝您閱讀我的部落格~