Fortran Coder

标题: 对称正定矩阵分解，计算结果错误，实在找不到原因了 [打印本页]
作者: likm1110    时间: 2015-3-26 15:30
标题: 对称正定矩阵分解，计算结果错误，实在找不到原因了
书上给出的算法过程：

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

! 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.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'

Screenshot from 2015-03-26 15:27:02.png (263.33 KB, 下载次数: 315)

书上给出的算法过程

书上给出的算法过程
作者: pasuka    时间: 2015-3-26 18:14
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
pasuka 发表于 2015-3-26 18:14
matlab给出的结果，lz答案哪里错了？！
15.279446979521216                   0                   0     ...

但是老师让我把那本书所有程序都打一遍我也没办法。还得边打便给他改错我也是醉了。。。。。另外就是我手算验证过了，A=L*LT（LT是转置） 按照机器算的，和原矩阵对不上。难道是书上的算法给错了或者我写错了？
作者: pasuka    时间: 2015-3-27 12:09
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
估计你浮点数采用的是fast形式，这个可以设置的，要设置成严格类型，而且是IEEE扩展类型，结果精确才会超过6位
作者: likm1110    时间: 2015-3-28 10:51
pasuka 发表于 2015-3-27 12:09
单精度和双精度的问题，给你改了代码，自己慢慢琢磨吧

多谢！我好好看看！
作者: chiangtp    时间: 2015-3-29 00:52
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
chiangtp 发表于 2015-3-29 00:52
1. Single precision (~7 precision digits) 的計算環境下, 反算驗證, 兩個答案都滿足

多谢！！！
作者: mangix2010    时间: 2015-3-31 13:30
回复好专业，大神好～～
作者: chiangtp    时间: 2015-4-1 14:19
請心直口快的pasuka參考/評價: 數值方法與程式(林聰悟-林佳慧)
頁面擷取自-數值方法與程式(林聰悟-林佳慧)-1.rar (1.56 MB, 下载次数: 9)
頁面擷取自-數值方法與程式(林聰悟-林佳慧)-2.rar (1.15 MB, 下载次数: 9)




欢迎光临 Fortran Coder (http://bbs.fcode.cn/) Powered by Discuz! X3.2