ykb.f90

用雅可比方法求矩阵的特征值和特征向量的FORTRAN程序

!用雅可比方法求矩阵的特征值和特征向量的FORTRAN程序
program YaKeBi 
parameter(n=5) 
integer p,q,loop 
double precision a(n,n),b(n,n),u(n,n),s,t2s,ss,cs,temp1,temp2 
open(100,file='in.txt',status='old') 
read(100,*) ((a(i,j),j=1,n),i=1,n) 
close(100) 
loop=0 
open(200,file='process.txt',status='unknown') 
do  i=1,n 
	do  j=1,n 
		u(i,j)=0.0d0 
	end do 
end do 
do  i=1,n 
	u(i,i)=1.0d0 
end do 
4 loop=loop+1 
if(loop.gt.100) then 
	write(*,*) 'error' 
	goto 90 
endif 
p=1 
q=1 
do i=1,n 
	do j=1,n 
		b(i,j)=a(i,j) 
	end do 
end do 
do  i=1,n 
	do  j=1,n 
		if(dabs(a(i,j)).gt.dabs(a(p,q))) then 
			p=i 
			q=j 
		endif 
	end do 
end do 
if(b(p,q).lt.1.0d-10) then 
goto 46 
endif 
t2s=2.0d0*a(p,q)/(a(p,p)-a(q,q)) 
s=0.5d0*atan(t2s) 
cs=cos(s) 
ss=sin(s) 
a(p,q)=0.0d0 
a(q,p)=0.0d0 
do   j=1,n 
	if((j.ne.p).and.(j.ne.q)) then 
		a(p,j)=b(p,j)*cs+b(q,j)*ss 
		a(j,p)=a(p,j) 
		a(q,j)=-b(p,j)*ss+b(q,j)*cs 
		a(j,q)=a(q,j) 
	endif 
	temp1=u(j,p) 
    temp2=u(j,q) 
	u(j,p)=temp1*cs+temp2*ss 
	u(j,q)=temp1*ss+temp2*cs 
end do 
a(p,p)=b(p,p)*(ss**2)+2.0d0*b(p,q)*cs*ss+b(q,q)*(cs**2) 
a(q,q)=b(p,p)*(ss**2)-2.0d0*b(p,q)*cs*ss+b(q,q)*(cs**2) 
write(200,*) 'Monitor Loop:',loop 
write(200,*) 'a(n,n):' 
do  i=1,n 
	write(200,45) (a(i,j),j=1,n) 
END DO 
write(200,*) 'u(n,n):' 
do  i=1,n 
	write(200,45) (u(i,j),j=1,n) 
END DO 
45 format(100d30.14) 
46 write(200,*) '--------------------------------------------------' 
goto 4 
open(300,file='out.txt',status='unknown') 
write(300,45) (a(i,i),i=1,n) 
write(300,*) '--------------------归一化之前--------------------' 
do i=1,n 
	write(300,45) (u(i,j),j=1,n) 
END DO
do   j=1,n 
	temp1=0.0d0 
	do   i=1,n 
		temp1=temp1+u(i,j)**2 
	end do 
	temp1=sqrt(temp1) 
	do i=1,n 
		u(i,j)=u(i,j)/temp1 
	end do 
end do 
write(300,*) '--------------------归一化之后--------------------' 
do  i=1,n 
	write(300,45) (u(i,j),j=1,n) 
end do 
90 close(200) 
close(300) 
end