24小时热门版块排行榜    

北京石油化工学院2026年研究生招生接收调剂公告
查看: 821  |  回复: 1

petersonic

铁虫 (初入文坛)

[求助] Fortan 反幂法求解矩阵按模最小特征值问题

本人用fortran反幂法求解最小特征值,其中求解线性方程组的部分由于是带状阵,因此采用了压缩空间的存储方式C(m,n)
但是子程序单独调试都能通过,放在一起就有问题了,哪位大神帮忙看看

program homework

implicit none
!定义
Integer , parameter :: p = SELECTED_REAL_KIND(8)  
real(p) ::        a(501),db=0.16,dc=-0.064,mc(5,501)
real(p)    :: Rs
integer :: i,n=501,m=5

!a初始化
do i=1,n
a(i)=(1.64-0.024*i)*sin(0.2*i)-0.64*exp(0.1/i)
end do
!给C赋初值5*501
mc(1,=dc
mc(2,=db
do i=1,n
mc(3,i)=a(i)
end do
mc(4,=db
mc(5,=dc
mc(1,1)=0.0
mc(1,2)=0.0
mc(2,1)=0.0
mc(5,n)=0.0
mc(4,n)=0.0
mc(5,n-1)=0.0


!求Rs(即按模最小的特征值),反幂法
call InPowoe(mc,m,n,Rs)
write(*,*)
write(*,*) "Rs===============",Rs
end program homework


!*****************************************************
!子程序反幂法求特征值(Inverse power of eigenvalue)         *
!*****************************************************
subroutine InPowoe(c,m,n,r)
implicit none
Integer , parameter :: p = SELECTED_REAL_KIND(8)
integer :: m,n
real(p) :: r,c(m,n)
real(p) :: u(n),y(n),t,max,flag
real(p) :: sum
integer :: i
sum=0.0
max=1.0
flag=1.0
u(1)=1.0
do i=2,n
u(i)=1.0
end do

!计算
do while(flag>=0.00001)
  do i=1,n
    sum=sum+u(i)*u(i)
  end do
   t=sqrt(sum)
   sum=0.0

  do i=1,n
    y(i)=u(i)/t
  end do

  !A*u=y,求解带状线性方程组
   call LU(c,u,y,m,n)
  !特征值
  do i=1,n
    sum=sum+y(i)*u(i)
  end do
    r=sum
        sum=0.0
  !结束判断
  flag=abs(r-max)/abs(r)
  max=r
end do
r=1.0/r
end subroutine         InPowoe


subroutine LU(a,x,b,m,n)
implicit none
Integer , parameter :: p = SELECTED_REAL_KIND(8)
integer :: i,j,k,r,s,t,m,n
real(p) :: a(m,n),x(n),b(n),c(m,n)
real(p) :: sum
sum=0.0
r=2
s=2

do i=1,m
  c(i,=a(i,
end do
!分解
do k=1,n
  do j=k,min(k+s,n)
    do t=max(1,k-r,j-s),k-1
          sum=sum+c(k-t+s+1,t)*c(t-j+s+1,j)
        end do
        c(k-j+s+1,j)=c(k-j+s+1,j)-sum
        sum=0.0
  end do
  
  if(k<n) then
  do i=k+1,min(k+r,n)
    do t=max(1,i-r,k-s),k-1
          sum=sum+c(i-t+s+1,t)*c(t-k+s+1,k)
        end do
        c(i-k+s+1,k)=(c(i-k+s+1,k)-sum)/c(s+1,k)
        sum=0.0
  end do
  end if
end do

!求解
do i=2,n
  do t=max(1,i-r),i-1
     sum=sum+c(i-t+s+1,t)*b(t)  
  end do
  b(i)=b(i)-sum
  sum=0.0
end do
x(n)=b(n)/c(s+1,n)
do i=n-1,1,-1
   do t=i+1,min(i+s,n)
      sum=sum+c(i-t+s+1,t)*x(t)
   end do
   x(i)=(b(i)-sum)/c(s+1,i)
   sum=0.0
end do
end subroutine LU
回复此楼
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

petersonic

铁虫 (初入文坛)

笑脸地方为  
2楼2014-11-09 12:42:50
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 petersonic 的主题更新
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 一志愿安徽大学计算机科学与技术学硕,331分求调剂 +5 蒋昌鹏qtj 2026-04-01 5/250 2026-04-02 08:10 by fxue1114
[考研] 298求调剂 +4 什么是胖头鱼 2026-03-30 6/300 2026-04-01 22:06 by 客尔美德
[考研] 290求调剂 +5 dfffsar 2026-03-29 5/250 2026-04-01 19:45 by 6781022
[考研] 0710生物学考研调剂 +3 李多米lee. 2026-03-27 4/200 2026-04-01 16:21 by zzchen2000
[考研] 材料0856 英一数二 323 求调剂 +9 袁sy 2026-04-01 9/450 2026-04-01 14:30 by wangjy2002
[考研] 材料与化工调剂一志愿大连海事085600,349 +9 吃的不少 2026-03-30 9/450 2026-04-01 11:24 by wangjy2002
[考研] 求调剂:一志愿:南京大学 专业:0705 总分320 ,本科985,四六级已过 +3 lfy760306 2026-03-31 3/150 2026-04-01 01:57 by Creta
[考研] 070300化学专业279调剂 +10 哈哈哈^_^ 2026-03-31 10/500 2026-03-31 23:13 by liu823948201
[考研] 299求调剂 +8 嗯嗯嗯嗯2 2026-03-27 8/400 2026-03-31 18:23 by lizhi8172
[考研] 318求调剂 +10 陈晨79 2026-03-30 10/500 2026-03-31 17:37 by 544594351
[考研] 考研调剂求助 +7 13287130938 2026-03-31 7/350 2026-03-31 16:39 by 690616278
[考研] 354求调剂 +3 lxb598 2026-03-31 4/200 2026-03-31 13:42 by sophie2180
[考研] 272求调剂,接受跨专业调剂! +3 闲鱼卢 2026-03-31 3/150 2026-03-31 13:00 by 替代品000
[考研] 一志愿浙江大学工科动力工程370,数一121,专业课135,现在能去哪里 +3 080700调剂 2026-03-30 4/200 2026-03-31 12:00 by KLMY666
[考研] 286求调剂 +5 丢掉懒惰 2026-03-27 8/400 2026-03-31 11:27 by Delta2012
[考研] 一志愿大连理工大学,机械工程学硕,341 +3 西瓜田的守望者 2026-03-30 3/150 2026-03-31 11:08 by asdfzly
[考研] 福建理工大学材料学院先进合金团队招收考研调剂学生 +3 大华金商都 2026-03-30 4/200 2026-03-31 01:04 by 方英俊602
[考研] 环境科学与工程334分求调剂 +6 王一一依依 2026-03-30 8/400 2026-03-30 11:52 by yjolah
[考研] 材料求调剂一志愿哈工大324 +7 闫旭东 2026-03-28 9/450 2026-03-28 08:51 by Xu de nuo
[考研] 315分求调剂 +7 26考研上岸版26 2026-03-26 7/350 2026-03-28 04:05 by fmesaito
信息提示
请填处理意见