我在波浪能设备上的论文挣扎。由于我是Fortran 90的新手,因此我想提高自己的编程技能。因此,我刚刚从
中摘了一个例子http://rosettacode.org/wiki/cholesky_decomposition
试图实施主页中解释的内容。基本上,它即将对3x3矩阵A进行的cholesky分解。我知道已经有一些软件包可以使Fortran进行分解,但是我想体验自己在学习如何编程方面的努力。
编译没有错误,但结果不匹配。尽管有元素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 Code主页上,没有解决方案的解决方案,因此任何提示都将不胜感激。
非常感谢您。
您要为i
和k
循环的每次迭代设置sum1
和sum2
。
我终于找到了如何解决问题,4x4矩阵等。这是最终代码:
Program Cholesky_decomp
!*************************************************!
!LBH @ ULPGC 06/03/2014
!Compute the Cholesky decomposition for a matrix A
!after the attached
!http://rosettacode.org/wiki/Cholesky_decomposition
!note that the matrix A is complex since there might
!be values, where the sqrt has complex solutions.
!Here, only the real values are taken into account
!*************************************************!
implicit none
INTEGER, PARAMETER :: m=3 !rows
INTEGER, PARAMETER :: n=3 !cols
COMPLEX, DIMENSION(m,n) :: A
REAL, DIMENSION(m,n) :: 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 /)
!!!!!!!!!!!!another example!!!!!!!
!A(1,:) = (/ 18, 22, 54, 42 /)
!A(2,:) = (/ 22, 70, 86, 62 /)
!A(3,:) = (/ 54, 86, 174, 134 /)
!A(4,:) = (/ 42, 62, 134, 106 /)
! Initialize values
L(1,1)=real(sqrt(A(1,1)))
L(2,1)=A(2,1)/L(1,1)
L(2,2)=real(sqrt(A(2,2)-L(2,1)*L(2,1)))
L(3,1)=A(3,1)/L(1,1)
!for greater order than m,n=3 add initial row value
!for instance if m,n=4 then add the following line
!L(4,1)=A(4,1)/L(1,1)
do i=1,n
do k=1,i
sum1=0
sum2=0
do j=1,k-1
if (i==k) then
sum1=sum1+(L(k,j)*L(k,j))
L(k,k)=real(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
期待听到有关评论,更好的编程方法,更正和任何形式的反馈。感谢2 Francescalus这么快回答!
,lbh