如何用乔列斯基分解求下三角矩阵？

cdf5041353

请假各位，如何用 乔列斯基分解获取矩阵的下三角矩阵啊。我查了何光渝老师的《Visual Fortran常用数值算法集》中有关对称方程组的乔列斯基（Cholesky）分解法的程序，但是不明白怎么得到矩阵分解后的下三角矩阵，求指点！

程序说明：
CHODCM(A,N,D,T)
N     整型变量，输入参数，矩阵A的阶数；
A      NxN个元素的二维实型数组，输入、输出参数，输入时存放A的上三角部分，严格下三角部分可置为任意数，输出时，严格下撒娇部分，存放分解L的严格下三角 ，而上三角部分不变；
D   N个元素的一维实型数组，输出参数，存放分解中对角均值D的对角元素；
T    N个元数的一维实型数组，工作单元。

[Fortran] 纯文本查看 复制代码

查看代码

打印?

01	SUBROUTINE chodcm(a,n,d,t)                 !乔列斯基分解子程序

02	REAL a(n,n),d(n),t(n)

03	INTEGER i,j

04

do i=1,n

05	sum=a(i,i)

06	do j=1,i-1

07	t(j)=a(j,i)

08	do k=1,j-1

09	t(j)=t(j)-t(k)*a(j,k)

10

end do

11	if (d(j)==0.) then

12	if (t(j)/=0.) then

13	pause 'no cholesky decomposition'

14

else

15	a(i,j)=1.0

16

endif

17

else

18	a(i,j)=t(j)/d(j)

19

endif

20	sum=sum-t(j)*a(i,j)

21

end do

22

d(i)=sum

23

end do

24	END SUBROUTINE chodcm

25

26

27

28	PROGRAM D1R10              !主程序

29	PARAMETER(N=2)

30	DIMENSION A(N,N),B(N),D(N),T(N)

31	DIMENSION C(N,N),X(N)

32	DATA A/1,-0.4179,-0.4179,1/

33	!Save matrix a for later testing

34

DO L=1,N

35

DO K=1,N

36	C(K,L)=A(K,L)

37

END DO

38

END DO

39	!DO Cholesky decomposition

40	CALL CHODCM(C,N,D,T)

41	!Solve equations for each right-hand vector

42

DO K=1,N

43

X(K)=B(K)

44

END DO

kerb

你的程序分解的是A的复制品C，你输出A当然没有变化，你应该输出C看看

cdf5041353

但是我将上述程序运行完，输出结果a，仍然是/1,-0.4179,-0.4179,1/,与Matlab采用chol(A)的结果（下三角矩阵）：1,0，-4.179，0.985有所差别，所以不晓得问题出在哪

楚香饭

程序说明：
CHODCM(A,N,D,T)
N     整型变量，输入参数，矩阵A的阶数；
A      NxN个元素的二维实型数组，输入、输出参数，输入时存放A的上三角部分，严格下三角部分可置为任意数，输出时，严格下三角部分，存放分解L的严格下三角 ，而上三角部分不变；
D   N个元素的一维实型数组，输出参数，存放分解中对角均值D的对角元素；
T    N个元数的一维实型数组，工作单元。