对称正定矩阵分解，计算结果错误，实在找不到原因了
书上给出的算法过程：http://bbs.fcode.cn/data/attachment/forum/201503/26/152657loxkouoozihsthtp.png
代码：
!---------------------------------------------------------------
! Description :Module of Cholesky Decomposition Symmetric Positive Definite Matrix
!---------------------------------------------------------------
Module cholesky
Contains
!-----------------------------------------------------
! Purpose : Subroutibe for Cholesky Decomposition
! Input parameters : 1. A Symmetric Positive Definite Matrix
! 2. N The Dimension
! Output parameters : L Output the matrix : A=L*LT
! Post Script : Cholesky Decomposition is only for Symmetric Positive Definite Matrix
!---------------------------------------------------
Subroutine solve(A,L,N)
Implicit None
Integer :: N
Real(8) :: A(N,N), L(N,N)
Real(8) :: s
Integer :: i, j, k
L = 0.
L(1,1) = dsqrt(A(1,1))
L(2:,1) = A(2:,1)/L(1,1)
Do j = 2,N
s = 0.
Do k = 1,j-1
s = s + L(j,k)**2
Enddo
L(j,j) = dsqrt(A(j,j)-s)
! Warning : the range of i
Do i = j+1,N
s = 0.
Do k = 1,j-1
s = s+L(i,k)*L(j,k)
Enddo
L(i,j) = (a(i,j)-s)/L(j,j)
Enddo
Enddo
End Subroutine solve
End Module cholesky
!===============================================================
! Purpose : To Solve Symmetric Positive Definite Matrix(Cholesky Decomposition)
! Data files : 1. fin.dat input file
! 2. fout.txt output file
!===============================================================
Program main
Use cholesky
Implicit None
Integer, Parameter :: N = 4
Integer :: i, j
Real(8) :: A(N,N), L(N,N)
Open(unit=11,file='fin.txt')
Open(unit=12,file='fout.txt')
Read(11,*)
Read(11,*) ((A(i,j),j=1,N),i=1,N)
Write(*,*) ((A(i,j),j=1,N),i=1,N)
Call solve(A,L,N)
Write(12,21)
21 Format (T5,'Cholesky Decomposition',/)
Do i = 1,N
Write(12,'(4(F10.6))') L(i,:)
Enddo
End Program main
测试数据：
Input datas of Cholesky Decomposition
233.4615113.8423256.0623145.0697
113.8423 78.6033127.4298 95.3089
256.0623127.4298281.4721164.8676
145.0697 95.3089164.8676181.2339
问题在这里。结果与正确结果不同，但是我是在找不到错误了。。。。：
fout.txt 输出的计算结果在‘’的三个数上是错误的
计算结果（错误的）
15.2794470.000000 0.000000 0.000000
7.4506824.805272 0.000000 0.000000
16.7586100.534147‘0.579450’ 0.000000
9.4944345.112904’5.217081‘‘6.142468’
应为(正确的)
15.2794470.000000 0.000000 0.000000
7.4506824.805272 0.000000 0.000000
16.7586100.534147‘0.579464’ 0.000000
9.4944345.112904’5.216939‘‘6.142588’
——The program is from 'Fortran95/2003科学计算与工程-宋叶志，page 33-36'
matlab给出的结果，lz答案哪里错了？！
15.279446979521216 0 0 0
7.450681962022638 4.805271927871565 0 0
16.758610461700346 0.534147362099914 0.579449729096019 0
9.494433940864122 5.112903632921073 5.217080776328907 6.142467643894371
估计是双精度和单精度的浮点数误差
btw，宋某某那本书翻过几页，真是一本烂书，估计是混职称用的 pasuka 发表于 2015-3-26 18:14
matlab给出的结果，lz答案哪里错了？！
15.279446979521216 0 0 ...
但是老师让我把那本书所有程序都打一遍我也没办法。还得边打便给他改错我也是醉了。。。。。另外就是我手算验证过了，A=L*LT（LT是转置） 按照机器算的，和原矩阵对不上。难道是书上的算法给错了或者我写错了？
likm1110 发表于 2015-3-27 09:40
但是老师让我把那本书所有程序都打一遍我也没办法。还得边打便给他改错我也是醉了。。。。。另外就是我手 ...
单精度和双精度的问题，给你改了代码，自己慢慢琢磨吧
估计你浮点数采用的是fast形式，这个可以设置的，要设置成严格类型，而且是IEEE扩展类型，结果精确才会超过6位 pasuka 发表于 2015-3-27 12:09
单精度和双精度的问题，给你改了代码，自己慢慢琢磨吧
:-handshake多谢！我好好看看！ 1. Single precision (~7 precision digits) 的計算環境下, 反算驗證, 兩個答案都滿足
2. 矩陣A須以 double precision (15-16 precision digits) 表達, 方得見兩個相異的答案
3. 系統(矩陣A)有點 ill-conditioned, Conditionnumber (based on L2 norm) C~7000
chiangtp 发表于 2015-3-29 00:52
1. Single precision (~7 precision digits) 的計算環境下, 反算驗證, 兩個答案都滿足
{:3_59:}多谢！！！ 回复好专业，大神好～～ 請心直口快的pasuka參考/評價: 數值方法與程式(林聰悟-林佳慧)
页:
[1]