錯誤編程FORTRAN中的Cholesky分解90

Im我的論文在波浪能設備上掙扎。由於我是FORTRAN 90的新手，我想提高自己的編程技能。因此，我剛從錯誤編程FORTRAN中的Cholesky分解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代碼主頁上沒有fortran的解決方案，所以任何暗示都是值得讚賞的。

非常感謝您提前。

來源

2014-03-05 user3385342

對於i和k循環的每次迭代，您都希望將sum1和sum2設置爲零。

來源

2014-03-05 22:15:27 francescalus

究竟應該在哪裏聲明sum1和sum2值？ – user3385342

在'do k = ...'和'do j = ...'之間調零。 – francescalus

非常感謝！我希望ress – user3385342

我終於找到了如何解決問題的更高階，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回答這麼快！

Regards，lbh

來源

2014-03-06 18:36:16 user3385342

錯誤編程FORTRAN中的Cholesky分解90

回答

相關問題