主题：重复，迭代，循环？
tianhy2010
[专家分：60] 发布于 2010-10-17 11:50:00
变量定义：
integer:: m=4,n=3
complex t0(2*m,2*m),t1(2*m,2*m),tt(2*m,2*m),t2(2*m,2*m),tt_(2*m,2*m),t3(2*m,2*m),t4(2*m,2*m) ,t(2*m,2*m),t5(2*m,2*m)
complex energy1(2*m,2*m),identity(2*m,2*m) ！w用energy1表示
complex h0(2*m,2*m),h2(2*m,2*m),h1(2*m,2*m) !都是已知矩阵，且h0表示H00,h1表示(H01+)，h2表示H01，tt表示ti,tt_表示ti_。
这样循环看下好不好：
t0=energy1-h0
t0=.i.t0
t0=t0.x.h2 !求t0,初始值
tt=t0 ！第一个t0值
t1=energy1-h0
t1=.i.t1
t1=t1.x.h1 !求t0_,初始值
t2=t1 ！将第一个t1赋给t2保存
do i=2,u ！循环，计算第i个t0值
t4=-t1.x.t0
t5=-t0.x.t1
t4=identity+t4+t5
t4=.i.t4
t0=t4.x.(t0**2) !ti
t2=t2.x.t0
tt=tt+t2
t1=t4.x.(t1**2) !t1_
! if(all(t1<eps*fi) .and. all(t0<eps*fi)) exit
t2=t2.x.t1
end do
t0=energy1-h0
t0=.i.t0
t0=t0.x.h2 !求t0,初始值
t2=t0
t1=energy1-h0
t1=.i.t1
t1=t1.x.h1 !求t0_,初始值
tt_=t1
do i=2,u
t3=-t1.x.t0
t5=-t0.x.t1
t3=identity+t5+t3
t3=.i.t3
t1=t3.x.(t1**2) !ti_
t2=t2.x.t1
tt_=tt_+t2
t0=t3.x.(t0**2) !ti
! if(all(t1<eps*fi) .and. all(t0<eps*fi)) exit
t2=t2.x.t0
end do
最后更新于：2010-10-17 12:54:00
回复列表 （共21个回复）
11 楼
tianhy2010 [专家分：60] 发布于 2010-10-18 10:21:00
try
12 楼
tianhy2010 [专家分：60] 发布于 2010-10-18 11:47:00
integer,parameter :: Max=1000
type cpx
complex,dimension(2*m,2*m) :: a
complex,dimension(2*m,2*m) :: b
end type cpx
type(cpx),dimension(Max) :: t
[em2][em2][em2][em2][em2][em2][em13][em13][em13][em13][em13][em13]
t是cpx类型，包含两种数组a,b分别存储ti和ti~.
t(1).a存储t0,t(1).b存储t0~,依次类推。
t0=energy1-h0
t0=.i.t0
t(1).a=t0.x.h2 !求t0,初始值
tt=t(1).a
t1=energy1-h0
t1=.i.t1
t(1).b=t1.x.h1 !求t0_,初始值
do i=2,u
t4=identity-t(i-1).a.x.t(i-1).b-t(i-1).b.x.t(i-1).a
t4=.i.t4
t(i).a=t4.x.(t(i-1).a.x.t(i-1).a) !ti
t(i).b=t4.x.(t(i-1).b.x.t(i-1).b) !ti_
t2=t2.x.t(i).b
if(i==2) then
tt=tt+t(i-1).b.x.t(i).a
else
tt=tt+t2.x.t(i).a
end if
! if(all(t(i).b<eps*hi) .and. all(t(i).a<eps*hi)) exit
end do
t0=energy1-h0
t0=.i.t0
t(1).a=t0.x.h2 !求t0,初始值
t2=t(1).a
t1=energy1-h0
t1=.i.t1
t(1).b=t1.x.h1 !求t0_,初始值
tt_=t(1).b
do i=2,u
t3=identity-t(i-1).a.x.t(i-1).b-t(i-1).b.x.t(i-1).a
t3=.i.t3
t(i).b=t4.x.(t(i-1).b.x.t(i-1).b) !ti_
t(i).a=t3.x.(t(i-1).a.x.t(i-1).a) !ti
t2=t2.x.t(i-1).a
if(i==2) then
tt_=tt_+t(1).a.x.t(2).b
else
tt_=tt_+t2.x.t(i).b
end if
! if(all(t1<eps*fi) .and. all(t0<eps*fi)) exit
end do
这样用type类型应该可以了吧？但是出现个错误：
f90: Fatal: There has been an internal compiler error (C0000005).
Error executing df.exe.
tianhynorand.obj - 1 error(s), 0 warning(s)
[em14][em14][em14][em14][em15][em15][em15][em15][em15][em15][em15][em15][em15]
13 楼
tianhy2010 [专家分：60] 发布于 2010-10-18 11:48:00
[quote]try[/quote]
tianhynorand.obj - 1 error(s), 0 warning(s)
[em14][em14][em14][em14][em15][em15][em15][em15][em15][em15][em15][em15][em15]
14 楼
tianhy2010 [专家分：60] 发布于 2010-10-18 11:48:00
[quote]
循环看上去很乱，费很大劲看不清楚，早建议你声明一个包含复数矩阵的type，再用type定义一个复数矩阵数组。这样去运算看起来就清晰多了。
integer,parameter :: Max=1000,i
type cpx
complex,dimension(8,8) :: a
complex,dimension(8,8) :: b
end type cpx
type(cpx),dimension(Max) :: t
引用的时候
t(:).a 用于存放一个序列的复数矩阵。
t(:).b 用于存放这个序列的复数矩阵的共轭。
再写程序，循环就是针对下标了，应该就清晰多了，试试吧。[/quote]
integer,parameter :: Max=1000
type cpx
complex,dimension(2*m,2*m) :: a
complex,dimension(2*m,2*m) :: b
end type cpx
type(cpx),dimension(Max) :: t
[em2][em2][em2][em2][em2][em2][em13][em13][em13][em13][em13][em13]
t是cpx类型，包含两种数组a,b分别存储ti和ti~.
t(1).a存储t0,t(1).b存储t0~,依次类推。
t0=energy1-h0
t0=.i.t0
t(1).a=t0.x.h2 !求t0,初始值
tt=t(1).a
t1=energy1-h0
t1=.i.t1
t(1).b=t1.x.h1 !求t0_,初始值
do i=2,u
t4=identity-t(i-1).a.x.t(i-1).b-t(i-1).b.x.t(i-1).a
t4=.i.t4
t(i).a=t4.x.(t(i-1).a.x.t(i-1).a) !ti
t(i).b=t4.x.(t(i-1).b.x.t(i-1).b) !ti_
t2=t2.x.t(i).b
if(i==2) then
tt=tt+t(i-1).b.x.t(i).a
else
tt=tt+t2.x.t(i).a
end if
! if(all(t(i).b<eps*hi) .and. all(t(i).a<eps*hi)) exit
end do
t0=energy1-h0
t0=.i.t0
t(1).a=t0.x.h2 !求t0,初始值
t2=t(1).a
t1=energy1-h0
t1=.i.t1
t(1).b=t1.x.h1 !求t0_,初始值
tt_=t(1).b
do i=2,u
t3=identity-t(i-1).a.x.t(i-1).b-t(i-1).b.x.t(i-1).a
t3=.i.t3
t(i).b=t4.x.(t(i-1).b.x.t(i-1).b) !ti_
t(i).a=t3.x.(t(i-1).a.x.t(i-1).a) !ti
t2=t2.x.t(i-1).a
if(i==2) then
tt_=tt_+t(1).a.x.t(2).b
else
tt_=tt_+t2.x.t(i).b
end if
! if(all(t1<eps*fi) .and. all(t0<eps*fi)) exit
end do
这样用type类型应该可以了吧？但是出现个错误：
f90: Fatal: There has been an internal compiler error (C0000005).
Error executing df.exe.
tianhynorand.obj - 1 error(s), 0 warning(s)
[em14][em14][em14][em14][em15][em15][em15][em15][em15][em15][em15][em15][em15]
15 楼
jstzhurj [专家分：4680] 发布于 2010-10-18 12:29:00
代码全部贴出。
16 楼
tianhy2010 [专家分：60] 发布于 2010-10-18 13:17:00
[em8][em8][em7][em7][em6][em6][em11][em11][em11][em12][em12][em12]
program main
use IMSL
implicit none
parameter pi=3.1415926
real e,eps,s !s表示小量，is是能量附加的小虚单位
complex:: fi=(0.0,1.0),conductance,hi=(1.0,0.0)
integer:: i,n,m,j,u
parameter(u=3,n=2**u,m=6)
integer,parameter :: Max=1000
type cpx
complex,dimension(2*m,2*m) :: a
complex,dimension(2*m,2*m) :: b
end type cpx
type(cpx),dimension(Max) :: t
complex h(2*m*n,2*m*n),h0(2*m,2*m),h2(2*m,2*m),h1(2*m,2*m) !哈密顿矩阵相关
complex t0(2*m,2*m),t1(2*m,2*m),tt(2*m,2*m),t2(2*m,2*m),tt_(2*m,2*m),t3(2*m,2*m),t4(2*m,2*m)!转移矩阵相关
complex energy1(2*m,2*m),identity(2*m,2*m)
s=0.0002
e=1.4
energy1=0.0
identity=0.0
do i=1,2*m
energy1(i,i)=e+s*fi
identity(i,i)=1.0*hi
end do
call hmax(h,h0,h1,h2,n,m)
!求转移矩阵tt,tt_
t0=energy1-h0
t0=.i.t0
t(1).a=t0.x.h2 !求t0,初始值
tt=t(1).a
t1=energy1-h0
t1=.i.t1
t(1).b=t1.x.h1 !求t0_,初始值
t2=t(1).b
do i=2,u
t4=identity-t(i-1).a.x.t(i-1).b-t(i-1).b.x.t(i-1).a
t4=.i.t4
t(i).a=t4.x.(t(i-1).a.x.t(i-1).a) !ti
t(i).b=t4.x.(t(i-1).b.x.t(i-1).b) !ti_
t2=t2.x.t(i).b
if(i==2) then
tt=tt+t(i-1).b.x.t(i).a
else
tt=tt+t2.x.t(i).a
end if
! if(all(t(i).b<eps*hi) .and. all(t(i).a<eps*hi)) exit
end do
t0=energy1-h0
t0=.i.t0
t(1).a=t0.x.h2 !求t0,初始值
t2=t(1).a
t1=energy1-h0
t1=.i.t1
t(1).b=t1.x.h1 !求t0_,初始值
tt_=t(1).b
do i=2,u
t3=identity-t(i-1).a.x.t(i-1).b-t(i-1).b.x.t(i-1).a
t3=.i.t3
t(i).b=t4.x.(t(i-1).b.x.t(i-1).b) !ti_
t(i).a=t3.x.(t(i-1).a.x.t(i-1).a) !ti
t2=t2.x.t(i-1).a
if(i==2) then
tt_=tt_+t(1).a.x.t(2).b
else
tt_=tt_+t2.x.t(i).b
end if
! if(all(t1<eps*fi) .and. all(t0<eps*fi)) exit
end do
write(*,*) tt_
stop
contains
subroutine hmax(h,h0,h1,h2,n,m)
real v
integer n,m,i,j
complex:: hi=(1.0,0.0)
complex::h(2*m*n,2*m*n)
complex h00(m,m),h01(m,m),h10(m,m)
complex h0(2*m,2*m)
complex h1(2*m,2*m),h2(2*m,2*m)
v=-3.0*hi
h00=0.0
h01=0.0
h=0.0
h0=0.0
h1=0.0
h2=0.0
h10=0.0
!y方向有m条横向链，x有n条带，所以一个大单胞总碳原子个数为2mn
do i=1,m
if(mod(i,2)==0) then
h00(i,i-1)=v
else
h00(i,i+1)=v
end if
end do
!右上角h01,左边列对右边列作用
do i=1,m
if(i==1) then
h01(i,m)=v
else
h01(i,i-1)=v
end if
end do
!左下角h10,同一条带右边列对左边列作用矩阵
do i=1,m
if(i==m) then
h10(i,1)=v
else
h10(i,i+1)=v
end if
end do
h0(1:m,1:m)=h00
h0(1:m,m+1:2*m)=h01
h0(m+1:2*m,1:m)=h10
h0(m+1:2*m,m+1:2*m)=h00
h1(m+1:2*m,1:m)=h10 !与右边带相互作用矩阵
h2(1:m,m+1:2*m)=h01 !与左边带相互作用矩阵
do i=1,2*m*(n-1)+1,2*m
h(i:i+2*m-1,i:i+2*m-1)=h0
if(i==1) then
h(1:2*m,2*m+1:4*m)=h1
! h(1:2*m,(n-1)*2*m+1:2*m*n)=h2
else if(i==(n-1)*2*m+1) then
h(i:2*m*n,i-2*m:i-1)=h2
! h(i:2*m*n,1:2*m)=h1
else
h(i:i+2*m-1,i+2*m:i+4*m-1)=h1
h(i:i+2*m-1,i-2*m:i-1)=h2
end if
end do
end subroutine
end program main
17 楼
jstzhurj [专家分：4680] 发布于 2010-10-18 16:32:00
.x.两边的复数矩阵加括号。
18 楼
tianhy2010 [专家分：60] 发布于 2010-10-18 16:40:00
该怎么感谢您呢[em4][em4][em4][em3][em3][em11][em11][em11][em2][em2][em2][em80][em80][em80][em75][em75][em76][em76][em76][em77][em77][em77][em74][em74][em74]
19 楼
doctorlive [专家分：800] 发布于 2010-10-18 17:43:00
如果已经解决了你的问题，那就结帖吧[em2]
20 楼
jstzhurj [专家分：4680] 发布于 2010-10-18 17:47:00
[quote]该怎么感谢您呢[em4][em4][em4][em3][em3][em11][em11][em11][em2][em2][em2][em80][em80][em80][em75][em75][em76][em76][em76][em77][em77][em77][em74][em74][em74][/quote]
用了type定义的，是不是清晰了？
我来回复
您所在位置：