(Fortran编程)数值分析之幂法

module POWER
!乘幂法求主特征值及特征向量
contains
subroutine SOLVE(A,N,tezheng,u,tol)
implicit real*8(a-z)
integer::N,i,k
real*8::A(N,N)
real*8::u(N),u0(N),v(N)

do i=1,N
u0(i)=1d0
end do

u=u0
m0=0
write(12,200)
200 format(/,T5,"k",T25,"u(k)规范化向量",T60,"max{v(k)}",//)
do k=1,500
v=matmul(A,u)
call MAX_ROU(v,N,m1)
! write(12,*)((A(i,j),j=1,N),i=1,N )
! write(12,220) u
!220 format("u=",3f12.4)
u=v/m1
! write(12,250) v
!250 format("v=",3f12.4,/)
if (dabs(m1-m0)m0=m1
write(12,300) k,u,m1
300 format(I5,T15,"(",F6.4,2F12.4,")",T55,F16.6)
end do
tezheng=m1
end subroutine SOLVE


subroutine MAX_ROU(r,N,ma)
implicit real*8(a-z)
integer::N,i,k
real*8::r(N)

ma=dabs(r(1))
k=1
do i=2,N
if(dabs(r(i))>ma) then
ma=dabs(r(i))
k=i !用k记录下标，但ma不取r(i)
end if
end do

ma=r(k)

end subroutine MAX_ROU
end module POWER


program mifa
use POWER
implicit real*8(a-z)
integer,parameter::N=3
real*8::A(N,N),u(N)
open(unit=11,file="juzhenA.txt")
open(unit=12,file="jieguo.txt")
read(11,*)
read(11,*)((A(i,j),j=1,N),i=1,N )
!write(12,*)((A(i,j),j=1,N),i=1,N )
call SOLVE(A,N,tezheng,u,1D-7)
write(12,100) tezheng,u
100 format(/,T3,"主特征值为：",/,f12.7,//,T3,"特征向量为：",/,T3,"(",f9.7,2(f12.7),")t")
end program mifa