技術標籤：陣列對稱

在科學和工程計算中，經常需要求解形如

而科學和工程的實際計算中，經常遇到係數矩陣

為正定陣，則有如下三角陣

使

將矩陣關係式

據此可逐行求出矩陣

基於矩陣分解式

可順序計算出

而由

可逆序求得

由於矩陣分解時公式含有開方運算，所以該演算法稱為平方根法，又叫Cholesky分解法。

根據上述公式，編寫程式即可對方程進行求解：

subroutine cholesky_full(n, a, y)
    implicit none
    
    integer, intent(in) :: n
    real, intent(inout) :: a(n, n), y(n)
    
    integer :: i, j, k
    real :: temp
    
    ! 分解矩陣，生成下三角陣L
    ! 工程問題中的很多矩陣非常龐大，所以，計算過程中的資料應該直接存放在原始陣列a中，
    ! 而不是新建立一個數組
    do i = 1, n
        ! 公式中，j的取值範圍為1到j-1，此處換成1到j，可以將分解式統一起來，省去一次判斷。
        ! 因為j=i時，j迴圈雖然會執行錯誤的操作、生成錯誤的a(i,j)結果，
        ! 但a(i,j)馬上就會被最外層的i迴圈生成的正確資料替換
        do j = 1, i
            temp = a(i, j)
            do k = 1, j-1
                temp = temp - a(i, k) * a(j, k)
            end do
            a(i, j) = temp / a(j, j)
            ! a(j, i) = 0.  ! 對角線上方d的元素賦0，可有可無
        end do
        a(i, i) = sqrt(temp)
    end do
    
    ! 根據Ly=b求解出y
    do i = 1, n
        temp = y(i)
        do j = 1, i-1
            temp = temp - a(i, j) * y(j)
        end do
        y(i) = temp / a(i, i)
    end do
    
    ! 求解出x
    do i = n, 1, -1
        temp = y(i) / a(i, i)
        y(i) = temp
        ! 公式中k的範圍為i+1到n，此處為1到i-1，因為下方a(i,k)的下標和公式中交換了順序
        do k = 1, i-1
            y(k) = y(k) - temp * a(i, k)
        end do
    end do
end subroutine
 

以上程式碼的Cholesky分解部分與前文公式基本上一致，很好理解，但引入了一個臨時變數temp，用於儲存資料。而如果我們將j、k兩層迴圈交換一下位置，再稍微調整一下迴圈計數器的取值範圍，就可以不借助臨時變數直接完成分解操作。程式碼如下：

do i = 1, n
    do k = 1, i - 1
        a(i, k) = a(i, k) / a(k, k)
        do j = k + 1, i
            a(i, j) = a(i, j) - a(i, k) * a(j, k)
        end do
    end do
    a(i, i) = sqrt(a(i, i))
end do
 

參考資料：

王能超. 高等學校教材, 數值分析簡明教程, （第2版）[M]. 2003.
吳建平, 王正華, 李曉梅. 稀疏線性方程組的高效求解與平行計算[M]. 湖南科學技術出版社, 2004.