矩陣除了基本的加法和乘法以外，還有很多種運算或分解。 其中，有不少運算都會用到 row 的基本操作(colunm 可自行類推)， 例如兩個 row 交換，某個 row 乘上一個常數加到另外一個 row 等等， 因此我們先來示範這兩種操作。 以下範例，會亂數產生一個矩陣，並進行上述的基本操作：
#include <stdio.h> #include <time.h> #include <stdlib.h> #define SIZE 5 void printArray(double array[][SIZE], int rowNum, int colNum){ int i, j; for(i=0;i<rowNum;i++){ for(j=0;j<colNum;j++){ printf("%.4f ", array[i][j]); } printf("\n"); } printf("\n"); } void rowExchange(double array[][SIZE], int size, int idx1, int idx2){ int i; double temp; for(i=0;i<size;i++){ temp = array[idx1][i]; array[idx1][i] = array[idx2][i]; array[idx2][i] = temp; } } void addRowToAnother(double array[][SIZE], int size, double c, int idx1, int idx2){ int i; for(i=0;i<size;i++){ array[idx2][i] += c * array[idx1][i]; } } int main() { double a[SIZE][SIZE]; int i, j, k; srand(time(NULL)); for(i=0;i<SIZE;i++){ for(j=0;j<SIZE;j++){ a[i][j] = rand()%10; } } printArray(a, SIZE, SIZE); printf("Swap row 1,4:\n"); rowExchange(a, SIZE, 1, 4); printArray(a, SIZE, SIZE); printf("2 * row 3 + row 0:\n"); addRowToAnother(a, SIZE, 2, 3, 0); printArray(a, SIZE, SIZE); return 0; }接著，我們來示範如何以高斯消去法(Gaussian Elimination)將矩陣進行三角化。 這個方法，我們其實在國中就學過，那時稱之為"加減消去法"； 高中時則告訴你"高斯消去法"這個名稱，以及引入矩陣的概念。 而也許你記得為了人工計算方便，我們會以好計算的數字所在的 row 當基準，來消去其他的 row； 在電腦裡面我們一樣可以如此操作，但還有另外一個理由，就是為了數值計算的穩定度。 以下範例，會示範如何將矩陣進行三角化，並展示中間每一步的運算過程：
#include <stdio.h> #include <time.h> #include <stdlib.h> #include <math.h> #define SIZE 5 void printArray(double array[][SIZE], int rowNum, int colNum){ int i, j; for(i=0;i<rowNum;i++){ for(j=0;j<colNum;j++){ printf("%7.4f ", array[i][j]); } printf("\n"); } printf("\n"); } void rowExchange(double array[][SIZE], int size, int idx1, int idx2){ int i; double temp; for(i=0;i<size;i++){ temp = array[idx1][i]; array[idx1][i] = array[idx2][i]; array[idx2][i] = temp; } } void addRowToAnother(double array[][SIZE], int size, double c, int idx1, int idx2){ int i; for(i=0;i<size;i++){ array[idx2][i] += c * array[idx1][i]; } } int main() { double a[SIZE][SIZE], maxVal; int i, j, maxRowIdx; srand(time(NULL)); for(i=0;i<SIZE;i++){ for(j=0;j<SIZE;j++){ a[i][j] = rand()%10; } } printf("Original matrix:\n"); printArray(a, SIZE, SIZE); for(j=0;j<SIZE-1;j++){ /* 在目前要處理的 column j 當中，找尋絕對值最大的 row i */ maxVal = fabs(a[j][j]); maxRowIdx = j; /* 先設對角線上是最大，再去比其他的 */ for(i=j+1;i<SIZE;i++){ if(fabs(a[i][j])>maxVal){ maxVal = fabs(a[i][j]); maxRowIdx = i; } } if(maxVal==0){ continue; } /* 如果最大值不在對角線上，則進行 row 的交換 */ if(maxRowIdx!=j){ printf("Swap row %d and %d:\n", j, maxRowIdx); rowExchange(a, SIZE, j, maxRowIdx); printArray(a, SIZE, SIZE); } /* 消去法 */ for(i=j+1;i<SIZE;i++){ printf("%.4f * row %d + row %d:\n", -a[i][j]/a[j][j], j, i); addRowToAnother(a, SIZE, -a[i][j]/a[j][j], j, i); printArray(a, SIZE, SIZE); } } return 0; }上述範例中的三角化，是把矩陣變成上三角矩陣(upper triangular matrix)， 當然如果你改變一下迴圈的方向， 很容易就可以得到下三角矩陣(lower triangular matrix)，甚至是單位矩陣(identity matrix)。上面提到可以透過列運算，把矩陣消去為單位矩陣，也許會讓你想起高中時的反矩陣求法： 先把原本矩陣跟單位矩陣，橫向並列為增廣矩陣； 而把左半邊原本矩陣透過列運算變成單位矩陣時，右半邊的矩陣就是反矩陣。 以下範例將示範基於此方法的反矩陣計算：
#include <stdio.h> #include <time.h> #include <stdlib.h> #include <math.h> #define SIZE 5 void printArray(double array[][SIZE], int rowNum, int colNum){ int i, j; for(i=0;i<rowNum;i++){ for(j=0;j<colNum;j++){ printf("%7.4f ", array[i][j]); } printf("\n"); } printf("\n"); } void rowExchange(double array[][SIZE], int size, int idx1, int idx2){ int i; double temp; for(i=0;i<size;i++){ temp = array[idx1][i]; array[idx1][i] = array[idx2][i]; array[idx2][i] = temp; } } void addRowToAnother(double array[][SIZE], int size, double c, int idx1, int idx2){ int i; for(i=0;i<size;i++){ array[idx2][i] += c * array[idx1][i]; } } void inv(double a[][SIZE], double out[][SIZE], double n){ int i, j, maxRowIdx; double maxVal, c; /* Downward */ for(j=0;j<SIZE-1;j++){ /* 在目前要處理的 column j 當中，找尋絕對值最大的 row i */ maxVal = fabs(a[j][j]); maxRowIdx = j; /* 先設對角線上是最大，再去比其他的 */ for(i=j+1;i<SIZE;i++){ if(fabs(a[i][j])>maxVal){ maxVal = fabs(a[i][j]); maxRowIdx = i; } } if(maxVal==0){ continue; } /* 如果最大值不在對角線上，則進行 row 的交換 */ if(maxRowIdx!=j){ rowExchange(a, SIZE, j, maxRowIdx); rowExchange(out, SIZE, j, maxRowIdx); } /* 消去法 */ for(i=j+1;i<SIZE;i++){ c = -a[i][j]/a[j][j]; addRowToAnother(a, SIZE, c, j, i); addRowToAnother(out, SIZE, c, j, i); } } /*Upward */ for(j=SIZE-1;j>0;j--){ /* 消去法 */ for(i=j-1;i>=0;i--){ c = -a[i][j]/a[j][j]; addRowToAnother(a, SIZE, c, j, i); addRowToAnother(out, SIZE, c, j, i); } } /* Diagonal */ for(i=0;i<SIZE;i++){ for(j=0;j<SIZE;j++){ out[i][j] /= a[i][i]; } } } int main() { double a[SIZE][SIZE], aInv[SIZE][SIZE]; int i, j; srand(time(NULL)); for(i=0;i<SIZE;i++){ for(j=0;j<SIZE;j++){ a[i][j] = rand()%10; aInv[i][j] = (i==j)?1:0; } } printf("Original matrix:\n"); printArray(a, SIZE, SIZE); printArray(aInv, SIZE, SIZE); inv(a, aInv, SIZE); printf("Inversed matrix:\n"); printArray(a, SIZE, SIZE); printArray(aInv, SIZE, SIZE); return 0; }關於行列式(determinant)的計算，也許你在線性代數等課會學到可以先進行三角化再來求解， 不過我們還是先來看看高中教的降階展開，要怎樣寫成程式。 在降階的方法中，要求 n-by-n 矩陣的行列式，需要計算 n 個 (n-1)-by-(n-1) 矩陣的行列式， 因此需要使用遞迴。範例如下：
#include <stdio.h> #include <time.h> #include <stdlib.h> #include <math.h> #define SIZE 5 void printArray(double array[][SIZE], int rowNum, int colNum){ int i, j; for(i=0;i<rowNum;i++){ for(j=0;j<colNum;j++){ printf("%7.4f ", array[i][j]); } printf("\n"); } printf("\n"); } double det(double array[][SIZE], int n){ double temp[SIZE][SIZE], result = 0; int sign = 1, i, ii, jj, iCounter; if(n==2){ return array[0][0] * array[1][1] - array[0][1] * array[1][0]; } for(i=0;i<n;i++){ /* copy */ iCounter = 0; for(ii=0;ii<n;ii++){ if(ii==i){ continue; } for(jj=1;jj<n;jj++){ temp[iCounter][jj-1] = array[ii][jj]; } iCounter++; } result += sign * array[i][0] * det(temp, n-1); sign = -sign; } return result; } int main() { double a[SIZE][SIZE]; int i, j; srand(time(NULL)); for(i=0;i<SIZE;i++){ for(j=0;j<SIZE;j++){ a[i][j] = rand()%10; } } printf("Original matrix:\n"); printArray(a, SIZE, SIZE); printf("Determinant: %.4f\n", det(a, SIZE)); return 0; }如果你在自己電腦上執行這個範例，可能會發現當矩陣大小增加到某個值後，程式會跑到天荒地老還沒有結果出現。 這是因為，若矩陣大小為 n-by-n，則降階法的執行時間，大約會與 nn 成比例(正式的說法，是"時間複雜度為O(nn)")。 因此，求行列式時，還是要透過三角化，會比較有效率。