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

回复此楼

» 猜你喜欢

反铁磁体中的磁性切换：两种不同的机制已成功可视化已经有0人回复
求标准粉末衍射卡号 ICDD 01-076-1802 已经有0人回复
物理学I论文润色/翻译怎么收费? 已经有164人回复
新西兰Robinson研究所招收全奖PhD 已经有0人回复
石墨烯转移--二氧化硅衬底石墨烯已经有0人回复
笼目材料中量子自旋液体基态的证据已经有0人回复
数学教学论硕士可以读数学物理博士吗？已经有0人回复
德国亥姆霍兹Hereon中心汉堡分部招镁合金腐蚀裂变SCC课题方向2026公派博士生已经有4人回复
澳门大学应用物理及材料工程研究院潘晖教授课题组诚招博士后已经有11人回复

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

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

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

返回列表