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

回复此楼

» 猜你喜欢

【复旦大学】二维材料方向招收2026年博士研究生1名已经有0人回复
北京纳米能源与系统研究所王中林院士/曹南颖研究员课题组2026级硕/博/博后招生已经有10人回复
物理学I论文润色/翻译怎么收费? 已经有68人回复
荷兰Utrecht University超快太赫兹光谱王海教授课题招收2026 CSC博士生已经有23人回复
反铁磁体中的磁性切换：两种不同的机制已成功可视化已经有0人回复
求标准粉末衍射卡号 ICDD 01-076-1802 已经有0人回复
新西兰Robinson研究所招收全奖PhD 已经有0人回复
石墨烯转移--二氧化硅衬底石墨烯已经有0人回复

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

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

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

返回列表