对称正定矩阵分解，计算结果错误，实在找不到原因了

likm1110 · 发表于 2015-3-26 15:30:15

书上给出的算法过程：

代码：

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

01!---------------------------------------------------------------
02! Description :  Module of Cholesky Decomposition Symmetric Positive Definite Matrix
03!---------------------------------------------------------------
04Module cholesky
05Contains
06  !-----------------------------------------------------
07  ! Purpose : Subroutibe for Cholesky Decomposition
08  ! Input parameters : 1. A     Symmetric Positive Definite Matrix
09  !                    2. N     The Dimension
10  ! Output parameters : L    Output the matrix : A=L*LT
11  ! Post Script : Cholesky Decomposition is only for Symmetric Positive Definite Matrix
12  !---------------------------------------------------
13  Subroutine solve(A,L,N)
14    Implicit None
15    Integer :: N
16    Real(8) :: A(N,N), L(N,N)
17    Real(8) :: s
18    Integer :: i, j, k
19    L = 0.
20    L(1,1) = dsqrt(A(1,1))
21    L(2:,1) = A(2:,1)/L(1,1)
22    Do j = 2,N
23      s = 0.
24      Do k = 1,j-1
25        s = s + L(j,k)**2
26      Enddo
27      L(j,j) = dsqrt(A(j,j)-s)
28      ! Warning : the range of i
29      Do i = j+1,N
30        s = 0.
31        Do k = 1,j-1
32          s = s+L(i,k)*L(j,k)
33        Enddo
34        L(i,j) = (a(i,j)-s)/L(j,j)
35      Enddo
36    Enddo
37  End Subroutine solve 
38End Module cholesky
39!===============================================================
40! Purpose : To Solve Symmetric Positive Definite Matrix(Cholesky Decomposition)
41! Data files : 1. fin.dat    input file 
42!              2. fout.txt   output file
43!===============================================================
44Program main
45Use cholesky
46  Implicit None
47  Integer, Parameter :: N = 4
48  Integer :: i, j
49  Real(8) :: A(N,N), L(N,N)
50  Open(unit=11,file='fin.txt')
51  Open(unit=12,file='fout.txt')
52  Read(11,*)
53  Read(11,*) ((A(i,j),j=1,N),i=1,N)
54  Write(*,*) ((A(i,j),j=1,N),i=1,N)
55  Call solve(A,L,N)
56  Write(12,21)
57  21 Format (T5,'Cholesky Decomposition',/)
58  Do i = 1,N
59    Write(12,'(4(F10.6))') L(i,:)
60  Enddo
61End Program main

测试数据：
Input datas of Cholesky Decomposition
233.4615  113.8423  256.0623  145.0697
113.8423 78.6033  127.4298 95.3089
256.0623  127.4298  281.4721  164.8676
145.0697 95.3089  164.8676  181.2339

问题在这里。结果与正确结果不同，但是我是在找不到错误了。。。。：
fout.txt 输出的计算结果在‘’的三个数上是错误的
   计算结果（错误的）
   15.279447  0.000000 0.000000 0.000000
   7.450682  4.805272 0.000000 0.000000
   16.758610  0.534147  ‘0.579450’ 0.000000
   9.494434  5.112904  ’5.217081‘  ‘6.142468’
   应为(正确的)
   15.279447  0.000000 0.000000 0.000000
   7.450682  4.805272 0.000000 0.000000
   16.758610  0.534147  ‘0.579464’ 0.000000
   9.494434  5.112904  ’5.216939‘  ‘6.142588’

                                                                                                            ——The program is from 'Fortran95/2003科学计算与工程-宋叶志，page 33-36'

pasuka · 发表于 2015-3-26 18:14:42

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，宋某某那本书翻过几页，真是一本烂书，估计是混职称用的

likm1110 · 发表于 2015-3-27 09:40:56

pasuka 发表于 2015-3-26 18:14
matlab给出的结果，lz答案哪里错了？！
15.279446979521216 0 0 ...

但是老师让我把那本书所有程序都打一遍我也没办法。还得边打便给他改错我也是醉了。。。。。另外就是我手算验证过了，A=L*LT（LT是转置）按照机器算的，和原矩阵对不上。难道是书上的算法给错了或者我写错了？

pasuka · 发表于 2015-3-27 12:09:07

likm1110 发表于 2015-3-27 09:40
但是老师让我把那本书所有程序都打一遍我也没办法。还得边打便给他改错我也是醉了。。。。。另外就是我手 ...

cholesky.zip (1.11 KB, 下载次数: 8)

单精度和双精度的问题，给你改了代码，自己慢慢琢磨吧

cholesky_test.f90 (1.27 KB, 下载次数: 7)

kerb · 发表于 2015-3-27 17:59:57

估计你浮点数采用的是fast形式，这个可以设置的，要设置成严格类型，而且是IEEE扩展类型，结果精确才会超过6位

likm1110 · 发表于 2015-3-28 10:51:35

pasuka 发表于 2015-3-27 12:09
单精度和双精度的问题，给你改了代码，自己慢慢琢磨吧

多谢！我好好看看！

chiangtp · 发表于 2015-3-29 00:52:23

1. Single precision (~7 precision digits) 的計算環境下, 反算驗證, 兩個答案都滿足

a.f90 (1.12 KB, 下载次数: 6)

2. 矩陣A須以 double precision (15-16 precision digits) 表達, 方得見兩個相異的答案

aa.f90 (2.61 KB, 下载次数: 6)

3. 系統(矩陣A)有點 ill-conditioned, Condition number (based on L2 norm) C~7000

ill-conditioned.pdf (120.86 KB, 下载次数: 8)

likm1110 · 发表于 2015-3-30 12:23:03

chiangtp 发表于 2015-3-29 00:52
1. Single precision (~7 precision digits) 的計算環境下, 反算驗證, 兩個答案都滿足

多谢！！！

mangix2010 · 发表于 2015-3-31 13:30:42

回复好专业，大神好～～

chiangtp · 发表于 2015-4-1 14:19:37

請心直口快的pasuka參考/評價: 數值方法與程式(林聰悟-林佳慧)

頁面擷取自-數值方法與程式(林聰悟-林佳慧)-1.rar (1.56 MB, 下载次数: 9)

頁面擷取自-數值方法與程式(林聰悟-林佳慧)-2.rar (1.15 MB, 下载次数: 9)

[线性代数] 对称正定矩阵分解，计算结果错误，实在找不到原因了

评分

评分

规矩勋章

水王勋章

元老勋章

热心勋章

贡献勋章

管理勋章

帅哥勋章

星光勋章