请假各位,如何用 乔列斯基分解获取矩阵的下三角矩阵啊。我查了何光渝老师的《Visual Fortran常用数值算法集》中有关对称方程组的乔列斯基(Cholesky)分解法的程序,但是不明白怎么得到矩阵分解后的下三角矩阵,求指点!
程序说明:
CHODCM(A,N,D,T)
N 整型变量,输入参数,矩阵A的阶数;
A NxN个元素的二维实型数组,输入、输出参数,输入时存放A的上三角部分,严格下三角部分可置为任意数,输出时,严格下撒娇部分,存放分解L的严格下三角 ,而上三角部分不变;
D N个元素的一维实型数组,输出参数,存放分解中对角均值D的对角元素;
T N个元数的一维实型数组,工作单元。
[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
|