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