马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
n=3
real v(3,3),a(3,3)
u=1e-6
open(1,file='v.dat')
open(2,file='ev.dat')
write(2,*)'input stiffness information'
read(1,*)((a(i,j),j=1,n),i=1,n)
write(2,100)((a(i,j),j=1,n),i=1,n)
write(2,*)'input vector information'
do 10 i=1,n
v(i,i)=1
do 15 j=1,n
if(i.ne.j) v(i,j)=0
15 continue
10 continue
write(2,100)((v(i,j),j=1,n),i=1,n)
call eig(v,a,n,u)
write(2,*)'acqure vector information'
write(2,100)((v(i),j=1,n),i=1,n)
close(1)
close(2)
100 format(1x,3f6.2)
end
subroutine eig(v,a,n,u)
real v(n,n),a(n,n)
integer p,q
66 fm=0.0
do20 i=1,n
do20 j=1,n
if((i.ne.j).and.(abs(a(i,j)).gt.fm))then
fm=abs(a(i,j))
p=i
q=j
end if
20 continue
if(fm.le.u)then
stop
endif
x=-a(p,q)
y=(a(q,q)-a(p,p))/2.0
fai=x/sqrt(x*x+y*y)
if(y.lt.0)fai=-fai
sn=1.0+sqrt(1.0-fai*fai)
sn=fai/sqrt(2*sn)
cn=sqrt(1-sn*sn)
fm=a(p,p)
a(p,p)=fm*cn*cn+a(q,q)*sn*sn+a(p,q)*fai
a(q,q)=fm*sn*sn+a(q,q)*cn*cn-a(p,q)*fai
a(p,q)=0.0
a(q,p)=0.0
do30 i=1,n
fm=v(i,p)
v(i,p)=fm*cn+v(i,q)*sn
v(i,q)=-fm*sn+v(i,q)*cn
30 contiue
goto 66
end
主要提示我第二行,和57行,我很想不明白,为什么会错 ,没有语法错误的。这是个雅克比程序!先谢谢大家了 渴望得到大家的指点!! |