[线性代数] 对称正定矩阵分解，计算结果错误，实在找不到原因了

!---------------------------------------------------------------
! Description :  Module of Cholesky Decomposition Symmetric Positive Definite Matrix
!---------------------------------------------------------------
Module cholesky
Contains
  !-----------------------------------------------------
  ! Purpose : Subroutibe for Cholesky Decomposition
  ! Input parameters : 1. A     Symmetric Positive Definite Matrix
  !                    2. N     The Dimension
  ! Output parameters : L    Output the matrix : A=L*LT
  ! Post Script : Cholesky Decomposition is only for Symmetric Positive Definite Matrix
  !---------------------------------------------------
  Subroutine solve(A,L,N)
    Implicit None
    Integer :: N
    Real(8) :: A(N,N), L(N,N)
    Real(8) :: s
    Integer :: i, j, k
    L = 0.
    L(1,1) = dsqrt(A(1,1))
    L(2:,1) = A(2:,1)/L(1,1)
    Do j = 2,N
      s = 0.
      Do k = 1,j-1
        s = s + L(j,k)**2
      Enddo
      L(j,j) = dsqrt(A(j,j)-s)
      ! Warning : the range of i
      Do i = j+1,N
        s = 0.
        Do k = 1,j-1
          s = s+L(i,k)*L(j,k)
        Enddo
        L(i,j) = (a(i,j)-s)/L(j,j)
      Enddo
    Enddo
  End Subroutine solve 
End Module cholesky
!===============================================================
! Purpose : To Solve Symmetric Positive Definite Matrix(Cholesky Decomposition)
! Data files : 1. fin.dat    input file 
!              2. fout.txt   output file
!===============================================================
Program main
Use cholesky
  Implicit None
  Integer, Parameter :: N = 4
  Integer :: i, j
  Real(8) :: A(N,N), L(N,N)
  Open(unit=11,file='fin.txt')
  Open(unit=12,file='fout.txt')
  Read(11,*)
  Read(11,*) ((A(i,j),j=1,N),i=1,N)
  Write(*,*) ((A(i,j),j=1,N),i=1,N)
  Call solve(A,L,N)
  Write(12,21)
  21 Format (T5,'Cholesky Decomposition',/)
  Do i = 1,N
    Write(12,'(4(F10.6))') L(i,:)
  Enddo
End Program main