algorithm - 在 FORTRAN 90 中编程 Cholesky 分解时出错

Question

我正在为我关于波浪能设备的论文而苦苦挣扎。由于我是 FORTRAN 90 的新手，我想提高我的编程技能。因此，我只是从

http://rosettacode.org/wiki/Cholesky_decomposition

并尝试实现主页中解释的内容。基本上它是关于对 3x3 矩阵 A 的 Cholesky 分解进行编程。我知道已经有一些包可以为 Fortran 进行分解，但我想亲 body 验学习如何编程的努力。

编译没有错误，但是结果不匹配。尽管有元素 L(3,3)，但我基本上找出了所有元素。在附件中，您可以找到我在 Fortran 90 中从头开始创建的代码:

Program Cholesky_decomp

implicit none
!size of the matrix
INTEGER, PARAMETER :: m=3 !rows
INTEGER, PARAMETER :: n=3 !cols
REAL, DIMENSION(m,n) :: A, L

REAL :: sum1, sum2
INTEGER i,j,k

! Assign values to the matrix
A(1,:)=(/ 25,  15,  -5 /)   
A(2,:)=(/ 15,  18,   0 /)  
A(3,:)=(/ -5,   0,  11 /)

! Initialize values
L(1,1)=sqrt(A(1,1))
L(2,1)=A(2,1)/L(1,1)
L(2,2)=sqrt(A(2,2)-L(2,1)*L(2,1))
L(3,1)=A(3,1)/L(1,1)

sum1=0
sum2=0
do i=1,n
    do k=1,i
        do j=1,k-1
        if (i==k) then
            sum1=sum1+(L(k,j)*L(k,j))
            L(k,k)=sqrt(A(k,k)-sum1)    
        elseif (i > k) then
            sum2=sum2+(L(i,j)*L(k,j))
            L(i,k)=(1/L(k,k))*(A(i,k)-sum2)
        else
            L(i,k)=0
        end if
        end do
    end do
end do

!write output
do i=1,m
    print "(3(1X,F6.1))",L(i,:)
end do

End program Cholesky_decomp

你能告诉我代码中的错误是什么吗？当它应该是 L(3,3)=3 时，我得到了 L(3,3)=0。我完全迷路了，仅作记录:在 Rosetta 代码主页上没有 fortran 的解决方案，因此不胜感激。

非常感谢您。

Answer 1

对于 i 和 k 循环，您希望将 sum1 和 sum2 设置为零。