[Fortran] 纯文本查看 复制代码
SUBROUTINE chodcm(a,n,d,t) !乔列斯基分解子程序
REAL a(n,n),d(n),t(n)
INTEGER i,j
do i=1,n
sum=a(i,i)
do j=1,i-1
t(j)=a(j,i)
do k=1,j-1
t(j)=t(j)-t(k)*a(j,k)
end do
if (d(j)==0.) then
if (t(j)/=0.) then
pause 'no cholesky decomposition'
else
a(i,j)=1.0
endif
else
a(i,j)=t(j)/d(j)
endif
sum=sum-t(j)*a(i,j)
end do
d(i)=sum
end do
END SUBROUTINE chodcm
PROGRAM D1R10 !主程序
PARAMETER(N=2)
DIMENSION A(N,N),B(N),D(N),T(N)
DIMENSION C(N,N),X(N)
DATA A/1,-0.4179,-0.4179,1/
!Save matrix a for later testing
DO L=1,N
DO K=1,N
C(K,L)=A(K,L)
END DO
END DO
!DO Cholesky decomposition
CALL CHODCM(C,N,D,T)
!Solve equations for each right-hand vector
DO K=1,N
X(K)=B(K)
END DO