LDL分解的Fortran代码

LDL分解的Fortran代码

LDL分解是将一个对称矩阵分解为一个对角矩阵 $D$ 和两个互为转置的三角阵 $L$, $L^T$ 的乘积 $A=LD{L}^T$。
具体做法可以看成将不断重复以下计算直到$C_i$的规格为$1\times 1$，对于$n\times n$ 的矩阵要重复n-1次

$$\left[\begin{array}{cc}a& v^T\\ v& C_i\end{array}\right]=\left[\begin{array}{cc}1& o\\ v/a& E\end{array}\right]\left[\begin{array}{cc}a& o\\ o& C_{i+1}\end{array}\right]\left[\begin{array}{cc}1& v^T/a\\ o& E\end{array}\right]$$
其中$C_{i+1}=C_i-v{v}^T/a$
Fortran代码如下

subroutine ldl_decompose(a,l,d,n,erro)
    implicit none
    integer,intent(in)::n
    real,intent(inout)::a(n,n)
    real,intent(out)::l(n,n),d(n,n)
    integer,intent(out)::erro
    integer::i,j,k
    do i=1,n
        do j=1,i-1
            if (a(i,j)/=a(j,i)) then
                erro=0
                stop
            end if
        end do
    end do
    do i=1,n
        if (a(i,i)+1.0==0.0) then
            erro=0
            stop
        end if
        do j=i+1,n
            a(i,j)=a(i,j)/a(i,i)
            a(j,i)=a(i,j)
        end do
        do j=i+1,n
            a(j,j)=a(j,j)-a(j,i)*a(j,i)*a(i,i)
            do k=i+1,j-1
                A(j,k)=A(j,k)-A(k,i)*A(j,i)*A(i,i)
                A(k,j)=A(j,k)
            end do
        end do
    end do
    do i=1,n
        d(i,i)=A(i,i)
        l(i,i)=1
        do j=1,i-1
            l(i,j)=A(i,j)
            l(j,i)=0
            d(i,j)=0
            d(j,i)=0
        end do
    end do
    erro=1 
end subroutine