24小时热门版块排行榜

返回列表

chaoszero253

新虫 (初入文坛)

应助: 0 (幼儿园)
金币: 5
帖子: 5
在线: 49分钟
虫号: 6900181
注册: 2017-07-04
性别: GG
专业: 光学

[求助] 关于Jacobi法矩阵对角化

想用fortran语言通过Jacobi方法对角化一个对称矩阵，但是求出来的结果和matlab内置的对角化程序员结果不一样。请问我的对角化程序是哪里出问题了。

subroutine SOLVE(A,N,tezheng,tol)
   !jacobi法对角化
   implicit real*8(a-z)
   integer::N
   real*8::A(N,N),tezheng(N)
   real*8::A1(N,N),R(N,N),RT(N,N),U(N,N),thi(N,N)
   integer::i,j,k,m,c,b,z,t
   integer::p,q,x,y
   real*8 sum
   do i=1,N
      do j=1,N
            A1(i,j)=A(i,j)
      end do
   end do
   k=0
   p=0
   m=1
   do while(m==1)
      m=0
      do i=1,N-1
            do j=i+1,N
               R=0
               do q=1,N
                  R(q,q)=1
               end do
               if(abs(A1(i,j))>tol) then
                  m=1
!       判断非对角矩阵元是否大于给出的阈值，如果是则建立旋转矩阵R来归零矩阵元
                  if(A1(i,i)==A1(j,j))then
                        R(i,i)=1/1.41421356
                        R(j,j)=1/1.41421356
                        R(i,j)=1/1.41421356
                        R(j,i)=-1/1.41421356
                  else
                        thi(i,j)=0.5*atan(2*A1(i,j)/(A1(i,i)-A1(j,j)))
                        R(i,i)=cos(thi(i,j))
                        R(j,j)=cos(thi(i,j))
                        R(i,j)=-sin(thi(i,j))
                        R(j,i)=sin(thi(i,j))
                  end if
                  do x=1,N
                        do y=1,N
                           RT(x,y)=R(y,x)
                        end do
                  end do
                  U=matmul(RT,A1)
                  A1=matmul(U,R)
               end if
            end do
      end do
      print *,A1(9,9)
      p=p+1
      print *,p
      if (p==100) exit
   end do
   do i=1,N
      tezheng(i)=A1(i,i)
   end do
end subroutine SOLVE

回复此楼

» 猜你喜欢

大豆异黄酮分离已经有0人回复
湖南大学材料学院急招2026年博士生，临时增加一名博士联培指标已经有10人回复
物理学I论文润色/翻译怎么收费? 已经有295人回复
天津理工大学晶体材料全国重点实验室刘红军教授课题组招收博士生1-2名已经有1人回复
中国科学院物理研究所谌志国研究员团队招收2027年博士研究生已经有5人回复
2026年中德博士后交流项目 - 新型量子和磁性材料：材料制备表征和中子散射研究已经有12人回复
26申博推荐：南京航空航天大学国际前沿院光学方向招收博士生！已经有1人回复
如何从铁电相到顺电相。已经有1人回复

1楼 2018-07-16 13:41:46

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

相关版块跳转我要订阅楼主 chaoszero253 的主题更新

返回列表