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
      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

回复此楼

» 猜你喜欢

求助中医院临床中医中标的青年标书学习，谢谢已经有6人回复
浙大物院表面物理课题组长期招收普博及直博学生已经有259人回复
物理学I论文润色/翻译怎么收费? 已经有86人回复
【征稿】Springer Nature热电专题 Discover Materials / Discover Chemistry 已经有47人回复
25届硕士申请，材料计算方向已经有2人回复
课题组招1-2名硕士已经有1人回复
机时越来越贵之后，我把一部分 DFT 流程搬到了 AI 计算平台上已经有0人回复
斯坦福大学与Elsevier世界前2%学者新排名查询已经有1人回复
27申博推荐：南京航空航天大学国际前沿院光学方向招收博士生！已经有19人回复
华南师范大学+双一流学科物理学+广东省杰青团队2027年博士招生+机器学习驱动材料设计已经有1人回复
apl投稿状态已经有0人回复

1楼 2018-09-15 10:10:40

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

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

返回列表