修改版 LDLᵀ 分解演算法實作與 C++ 程式碼詳解

編輯：關於C語言

前言

LDLᵀ 分解是求解對稱正定線性方程組的經典數值方法，屬於丘列斯基分解的變體。本文介紹一個修改版的 LDLᵀ 分解實作，包含完整的 C++ 原始碼、輸入輸出範例，以及關鍵步驟的詳細註解。程式碼展示了如何從檔案讀取矩陣與向量、執行分解、前向/回代求解，並輸出分解後的 L、D 矩陣與解向量。

核心演算法流程

1. 記憶體配置與資料讀入

程式啟動時提示使用者輸入輸入/輸出檔案路徑，開啟檔案串流並讀取矩陣維度 n。接著動態配置一維陣列儲存 n×n 矩陣 a、右側向量 b、對角矩陣對角元 g 以及解向量 x。採用列主序（Column-major）將二維矩陣壓扁為一維，存取公式為 a[i*n + j]。

2. 分解主迴圈（Step 1–10）

Step 1：檢查 a[0]（即 A₁₁）是否接近零，避免除以零。
Step 2–3：初始化第一列，計算 g[1] 並更新 a[n]、a[n+1]。
Step 4–10：對 i = 3 … n 迭代：
- Step 5：取得第 i 列第一個元素到 g[1]。
- Step 6：內層迴圈 j = 2 … i-1 計算 g[j]（扣除前面已計算項的貢獻）。
- Step 7–8：將 g 除以對應對角元，寫回 a 形成 L 矩陣的非對角元。
- Step 9：計算第 i 行對角元更新量 t。
- Step 10：更新對角元 a[i*n-n+i-1] 並檢查正定性。

3. 前向代入與回代（Step 11–13）

Step 11：前向代入求解 L y = b，結果暫存於 b。
Step 12：求解最後一個未知數 x[n]。
Step 13：回代求解剩餘未知數 x[i]。

輸入輸出範例

輸入檔案

10
1 1 1 1 1 1 1 1 1 1
1 2 2 2 2 2 2 2 2 2
1 2 3 3 3 3 3 3 3 3
1 2 3 4 4 4 4 4 4 4
1 2 3 4 5 5 5 5 5 5
1 2 3 4 5 6 6 6 6 6
1 2 3 4 5 6 7 7 7 7
1 2 3 4 5 6 7 8 8 8
1 2 3 4 5 6 7 8 9 9
1 2 3 4 5 6 7 8 9 10
8 2 -12 4 -45 6 -5 9 10 30

輸出結果

g(1)=1
g(2)=1
g(3)=1
g(4)=1
g(5)=1
g(6)=1
g(7)=1
g(8)=1
g(9)=1
g(10)=-6.27744e+066
y=(8,-6,-14,16,-49,51,-11,14,1,20);
diag(D)=(1,1,1,1,1,1,1,1,1,1);
L=
1 0 0 0 0 0 0 0 0 0
1 1 0 0 0 0 0 0 0 0
1 1 1 0 0 0 0 0 0 0
1 1 1 1 0 0 0 0 0 0
1 1 1 1 1 0 0 0 0 0
1 1 1 1 1 1 0 0 0 0
1 1 1 1 1 1 1 0 0 0
1 1 1 1 1 1 1 1 0 0
1 1 1 1 1 1 1 1 1 0
1 1 1 1 1 1 1 1 1 1
The roots are:(14,8,-30,65,-100,62,-25,13,-19,20);

程式碼關鍵片段

// Step 4–10 主分解迴圈
for(i=3; i<w; i++) {
    g[1] = a[i*n-n];                 // Step 5
    for(j=2; j<i; j++) {            // Step 6
        for(t=0, k=1; k<j; k++)
            t += g[k] * a[j*n-n+k-1];
        g[j] = a[i*n-n+j-1] - t;
    }
    a[i*n-n] = g[1] / a[0];          // Step 7
    for(j=2; j<i; j++)              // Step 8
        a[i*n-n+j-1] = g[j] / a[j*n-n+j-1];
    for(t=0, k=1; k<i; k++)         // Step 9
        t += g[k] * a[i*n-n+k-1];
    a[i*n-n+i-1] -= t;               // Step 10
    if(a[i*n-n+i-1] < 1e-8) { ... }
}

注意事項與改進方向

程式使用舊式 C++ 標頭（<iostream.h> 等），現代編譯器建議改用 <iostream>、<fstream> 並加上 using namespace std;。
動態記憶體以 new/delete 管理，建議改用 std::vector 或智慧指標避免記憶體洩漏。
閾值 1e-8 硬編碼，實際應用可視矩陣條件數調整。
輸出格式固定，可進一步封裝為函式庫供其他專案呼叫。