Fortran Coder

查看: 6545|回复: 3
打印 上一主题 下一主题

[线性代数] 如何用乔列斯基分解求下三角矩阵?

[复制链接]

15

帖子

5

主题

0

精华

入门

F 币
100 元
贡献
50 点
跳转到指定楼层
楼主
发表于 2015-8-29 16:10:03 | 只看该作者 回帖奖励 |倒序浏览 |阅读模式
请假各位,如何用 乔列斯基分解获取矩阵的下三角矩阵啊。我查了何光渝老师的《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

分享到:  微信微信
收藏收藏 点赞点赞 点踩点踩

723

帖子

4

主题

0

精华

大师

农村外出务工人员

F 币
650 元
贡献
333 点

新人勋章爱心勋章水王勋章元老勋章热心勋章

沙发
发表于 2015-8-30 08:01:15 | 只看该作者
程序说明:
CHODCM(A,N,D,T)
N     整型变量,输入参数,矩阵A的阶数;
A      NxN个元素的二维实型数组,输入、输出参数,输入时存放A的上三角部分,严格下三角部分可置为任意数,输出时,严格下三角部分,存放分解L的严格下三角 ,而上三角部分不变;
D   N个元素的一维实型数组,输出参数,存放分解中对角均值D的对角元素;
T    N个元数的一维实型数组,工作单元。

15

帖子

5

主题

0

精华

入门

F 币
100 元
贡献
50 点
板凳
 楼主| 发表于 2015-8-31 11:20:34 | 只看该作者
但是我将上述程序运行完,输出结果a,仍然是/1,-0.4179,-0.4179,1/,与Matlab采用chol(A)的结果(下三角矩阵):1,0,-4.179,0.985有所差别,所以不晓得问题出在哪

59

帖子

2

主题

0

精华

大师

F 币
810 元
贡献
476 点
地板
发表于 2015-11-18 16:12:30 | 只看该作者
你的程序分解的是A的复制品C,你输出A当然没有变化,你应该输出C看看
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

捐赠本站|Archiver|关于我们 About Us|小黑屋|Fcode ( 京ICP备18005632-2号 )

GMT+8, 2024-11-22 12:01

Powered by Tencent X3.4

© 2013-2024 Tencent

快速回复 返回顶部 返回列表