Programming Practice

一個 n´n 的帶狀矩陣 (banded matrix) 包含下角和上三角帶狀區,帶寬 (bandwidth) r 的下三角帶狀區即是緊鄰正對角線 (diagonal) 之下的 r 條斜線的元素,帶寬為 s 的上三角帶狀區即是緊鄰正對角線 (diagonal) 之上的 s 條斜線的元素。一個 n´n 的帶狀矩陣只有下三角帶狀區、上三角帶狀區、和正對角線可以有非零的元素,其他的元素皆為零。例如,下列是一個下三角帶寬為 2,上三角帶寬為 3 的 8´8 矩陣:

帶狀 LU-分割 (banded LU-decomposition) 問題就是將一個 n´n 的帶狀矩陣 A 分割成兩個 n´n 的三角帶狀矩陣 LU。其中 L 是一個下三角帶狀矩陣 (lower triangular banded matrix) 且其對角線 (diagonal) 的元素為 1,U 是一個上三角帶狀矩陣 (upper triangular banded matrix) 且其對角線的元素都不為 0。 若 A 的下三角帶寬和上三角帶寬分別為 rs,則 L 的下三角帶寬為  rU 的上三角帶寬為 sL 的上三角帶寬和 U 的下三角帶寬可視為 0。

如果 A(k) 代表將 A 的前面 k 行和 k 列移除後的 (n-k)´(n-k) 帶狀矩陣,即

A(0) 開始到 A(n-1),帶狀  LU-分割的演算法,即是每個 A(k) 執行下列三個步驟:

  1. 計算矩陣 U 的第 k 列元素:uk,j = ak,j, for k£j£min(n-1,k+s),

  2. 計算矩陣 L 的第 k 行元素:li,k = ai,k / ak,k, for k£i£min(n-1,k+r),

  3. 計算矩陣 A(k+1) 的元素:ai,j = ai,j - li,k ´ uk,j, for k+1£i£min(n-1,k+r), and max(k+1,i-r)£j£min(n-1,min(i+s,k+s)),

寫一個 C 語言的程式以執行以下的步驟:

  1. 宣告 AA1、L、和 U 為 100´100 的二維陣列

  2. 輸入正整數 n 的值 (小於或等於 100);

  3. 使用副程式 rand() 隨機產生矩陣 A 元素的浮點數值 (假設 A 元素的值在 0 和 1 之間,並取小數點以下六位);同時,複製矩陣 A 的值到矩陣 A1;

  4. 輸出矩陣 A (至小數點以下四位);

  5. 計算 LU-分割的矩陣 L 和 矩陣 U 的值,以使得 A=LU

  6. 輸出矩陣 LU  (至小數點以下四位);

  7. 驗算 L U 的矩陣乘積是否等於 A1 (假設誤差值小於 0.0001)。

將帶狀 LU-分割寫成副程式 void lu_decompose_banded(int n, int r, int s, float A[100][100], float L[100][100], float U[100][100])。以下為程式執行範例:

Enter matrix size n: 8

Enter the lower bandwidth of matrix A, r:
3

Enter the upper bandwidth of matrix A, s:
4

Matrix A:
  0.0178   0.0305   0.0011   0.0239   0.0236
  0.0240   0.0194   0.0098   0.0282   0.0314   0.0106
  0.0260   0.0183   0.0189   0.0193   0.0270   0.0135   0.0312
  0.0109   0.0310   0.0011   0.0306   0.0140   0.0272   0.0040   0.0260
           0.0151   0.0304   0.0034   0.0132   0.0288   0.0084   0.0131
                    0.0178   0.0213   0.0153   0.0285   0.0250   0.0105
                             0.0279   0.0126   0.0229   0.0222   0.0143
                                      0.0053   0.0085   0.0276   0.0281

Matrix L:
  1.0000
  1.3483   1.0000
  1.4622   1.2098   1.0000
  0.6153  -0.5637   0.7216   1.0000
          -0.6955   5.0306   2.5716   1.0000
                    2.4692   2.2448   0.6228   1.0000
                             1.3006   0.1924   0.6681   1.0000
                                      0.1430  -1.0779   3.2024   1.0000

Matrix U:
  0.0178   0.0305   0.0011   0.0239   0.0236
          -0.0217   0.0084  -0.0040  -0.0004   0.0106
                    0.0072  -0.0108  -0.0069   0.0006   0.0312
                             0.0214   0.0042   0.0327  -0.0185   0.0260
                                      0.0368  -0.0511  -0.1011  -0.0537
                                              -0.0147   0.0524  -0.0144
                                                        0.0308   0.0004
                                                                 0.0189

The LU-decomposition program is correct.