主题：jstzhurj老师，我又将程序简化了，能否看下？
tianhy2010
[专家分：60] 发布于 2010-08-28 22:18:00
还是上次那代码。我想了下，由于相邻的两个特征值之间可能有好几个p，比如在eigenvalue(4)和eigenvalue(5)之间有8个p，分别是23、24、25、。。。、30，他们的conr值相同，我只要算一个就可以了，然后把这个值赋给其他的conr，比如得到了conr(23),可以把conr(23)直接赋给conr(24、25、26、...、30),这样就不用再重复计算了，我改了下程序，帮忙看下吧？呵呵
q=1
do i=1,4*m*n-1 ！每个i找满足的p,先对i循环
conduc=0.0
k=0
do p=q,450000
if((eigenvalue(i)<=Ef(p)).and.(eigenvalue(i+1)>Ef(p))) then ！如果p满足条件，就执行下面的操作，再看下满足这个条件的p有几个
q=p+1 !p+1看下还能否满足条件，如果能那就k+1,直到不能满足条件跳出if语句
k=k+1 !计算总的能满足i条件的p的个数
end if
do temp=1,i
do j=i+1,4*m*n
ea=matmul(traneigv(temp,:),hx)
ec=matmul(traneigv(temp,:),hy)
eb=matmul(traneigv(j,:),hx)
ed=matmul(traneigv(j,:),hy)
conduc=conduc+(dot_product(ea,eigenvector(:,j))*dot_product(ed,eigenvector(:,temp))
end do
end do
conr(p)=2*pi*aimag(conduc)/S
！将能满足i条件的所有conr(p)都赋值,如果k>1，那就说明不止一个能满足条件了
if(k>1) then
do ii=1,k
conr(p+ii)=conr(p)
end do
end if
q=p
exit ！退到计算下一个i的地方
end do
end do
最后更新于：2010-08-30 22:21:00
回复列表 （共58个回复）
31 楼
tianhy2010 [专家分：60] 发布于 2010-08-30 21:15:00
matlab运算得到了比较好的结果。
换成complex*16又试了，fortran不行，跟matlab得到的结果不同，我看了，还是特征向量的问题啊。数学上倒是说过，对于一个矩阵来说，特征向量有很多个,我曾经把公式带进去试过：Hr=Ar,H,A,r分别是矩阵、矩阵特征值和特征向量，在matlab下，二者严格相等，但是在fortran下，二者只在三位小数范围内相等，后面的就不相等了。今天把fortran得到的矩阵和matlab得到的矩阵拿出来看了下，基本是相同的矩阵，至少矩阵元素前10位数是相等的。真不明白为什么会得到这种结果出来。如果运行几次就会发现，每次得到的矩阵特征向量还都差不多，尽管跟matlab得到的结果差距较大。
[em69][em69][em60][em60][em30][em30][em32][em32][em41][em41][em51][em51][em11][em11][em12][em12][em13][em13][em14][em14][em9][em9][em80][em80][em72][em72][em72][em19][em19][em16][em16][em4][em4][em3][em3][em2][em2][em1][em1][em7][em7][em7]
32 楼
jstzhurj [专家分：4680] 发布于 2010-08-30 21:36:00
[quote]matlab运算得到了比较好的结果。
换成complex*16又试了，fortran不行，跟matlab得到的结果不同，我看了，还是特征向量的问题啊。数学上倒是说过，对于一个矩阵来说，特征向量有很多个,我曾经把公式带进去试过：Hr=Ar,H,A,r分别是矩阵、矩阵特征值和特征向量，在matlab下，二者严格相等，但是在fortran下，二者只在三位小数范围内相等，后面的就不相等了。今天把fortran得到的矩阵和matlab得到的矩阵拿出来看了下，基本是相同的矩阵，至少矩阵元素前10位数是相等的。真不明白为什么会得到这种结果出来。如果运行几次就会发现，每次得到的矩阵特征向量还都差不多，尽管跟matlab得到的结果差距较大。
[em69][em69][em60][em60][em30][em30][em32][em32][em41][em41][em51][em51][em11][em11][em12][em12][em13][em13][em14][em14][em9][em9][em80][em80][em72][em72][em72][em19][em19][em16][em16][em4][em4][em3][em3][em2][em2][em1][em1][em7][em7][em7][/quote]
Hr=Ar在Fortran下不严格相等，这是库函数的问题。要看这种误差是多少？我以为若是误差很小（小于百分之几）就可以暂且不管了，本人对这方面涉及很少，不知有没有更好的库函数。除此之外，程序是否还有其他更大的问题呢？
既然每次算的特征向量相等，loop=10这种无意义的循环也就没有必要了。[em18][em18][em67][em67][em77][em77]
33 楼
jstzhurj [专家分：4680] 发布于 2010-08-30 22:06:00
do j=1,4*m*n
do i=1,450000-1
if((Ef(i)<=eigenvalue(j)) .and. (eigenvalue(j)<Ef(i+1))) then
ds(i)=ds(i)+1
exit
end if
end do
end do
这段程序怎么这么让人费解？你明白循环嵌套执行的过程么？！
34 楼
tianhy2010 [专家分：60] 发布于 2010-08-30 22:07:00
能做到这样就很高兴了，可能程序还有更大的问题，谢谢大姐。自己再检查公式，检查程序看下哪里出问题了吧。我看到徐士良程序集上有给出的解实矩阵特征向量的程序，不知道有没有解虚数矩阵特征向量的程序？
[em32][em32][em35][em35][em36][em36][em37][em37][em38][em38][em30][em30][em72][em72][em62][em62][em68][em68][em80][em80][em78][em78][em77][em77][em4][em4][em4][em4][em4][em8][em8][em8][em8][em8]
35 楼
tianhy2010 [专家分：60] 发布于 2010-08-31 10:19:00
[quote]
do j=1,4*m*n
do i=1,450000-1
if((Ef(i)<=eigenvalue(j)) .and. (eigenvalue(j)<Ef(i+1))) then
ds(i)=ds(i)+1
exit
end if
end do
end do
这段程序怎么这么让人费解？你明白循环嵌套执行的过程么？！[/quote]
确实很费解。应该是这样吧：
!计算态密度
do i=1,450000-1
do j=1,4*m*n
if((Ef(i)<=avereig(j)) .and. (avereig(j)<Ef(i+1))) then
ds(i)=ds(i)+1
exit
end if
end do
end do
这地方我计算态密度，就是两个相邻的Ef之间有几个本征值。先从i=1开始检查，看本征值有没有落在Ef(1)和Ef(2)之间的，如果有avereig(j)的话，ds(i)+1，直到有某个avereig(j)不在这个区间，那么退出循环，就到i循环继续下一个循环。
36 楼
jstzhurj [专家分：4680] 发布于 2010-08-31 10:40:00
[quote][quote]
do j=1,4*m*n
do i=1,450000-1
if((Ef(i)<=eigenvalue(j)) .and. (eigenvalue(j)<Ef(i+1))) then
ds(i)=ds(i)+1
exit
end if
end do
end do
这段程序怎么这么让人费解？你明白循环嵌套执行的过程么？！[/quote]
确实很费解。应该是这样吧：
!计算态密度
do i=1,450000-1
do j=1,4*m*n
if((Ef(i)<=avereig(j)) .and. (avereig(j)<Ef(i+1))) then
ds(i)=ds(i)+1
exit
end if
end do
end do
这地方我计算态密度，就是两个相邻的Ef之间有几个本征值。先从i=1开始检查，看本征值有没有落在Ef(1)和Ef(2)之间的，如果有avereig(j)的话，ds(i)+1，直到有某个avereig(j)不在这个区间，那么退出循环，就到i循环继续下一个循环。[/quote]
exit?满足条件就退出？这样ds最多等于1，还能累加么？exit不是乱用的！！！
37 楼
tianhy2010 [专家分：60] 发布于 2010-08-31 10:55:00
加个else应该可以了吧？满足条件就累计，不满足就退出。还有个耗时问题。由于我这里Ef(i)和
avereig(j)都是单调递增的，我也参照以前给的做法，设一个q变量吧
q1=1
q2=1
do i=q1,450000-1
do j=q2,4*m*n
if((Ef(i)<=avereig(j)) .and. (avereig(j)<Ef(i+1))) then
ds(i)=ds(i)+1
else
q1=i+1 !不用让i从1开始循环了，从i+1开始循环，前面那些已经检查过了
q2=j ！如果j不满足这个条件，就让j的循环从这个数开始
exit
end if
end do
end do
38 楼
jstzhurj [专家分：4680] 发布于 2010-08-31 11:09:00
[quote]加个else应该可以了吧？满足条件就累计，不满足就退出。还有个耗时问题。由于我这里Ef(i)和
avereig(j)都是单调递增的，我也参照以前给的做法，设一个q变量吧
q1=1
q2=1
do i=q1,450000-1
do j=q2,4*m*n
if((Ef(i)<=avereig(j)) .and. (avereig(j)<Ef(i+1))) then
ds(i)=ds(i)+1
else
q1=i+1 !不用让i从1开始循环了，从i+1开始循环，前面那些已经检查过了
q2=j ！如果j不满足这个条件，就让j的循环从这个数开始
exit
end if
end do
end do
[/quote]
你编程是东一榔头西一拐！i变量还要你操心吗？退出之后它自动加1执行了，什么条件退出？！哎......[em1][em6][em7][em21][em21]
39 楼
tianhy2010 [专家分：60] 发布于 2010-08-31 11:25:00
早饭没吃，吃完回来再看[em13][em13][em13][em13][em13][em13][em72][em72][em72][em72][em21][em21][em21][em21][em21][em21][em21][em21][em21][em21][em11][em11][em11][em11][em11]
40 楼
tianhy2010 [专家分：60] 发布于 2010-08-31 12:09:00
当avereig(j)不在Ef(i)和Ef(i+1)之间的时候就退出。
比如Ef(12)=0.01,Ef(13)=0.21,
avereig(1)=0.05,avereig(2)=0.1,avereig(3)=0.15,avereig(4)=0.20,avereig(5)=0.25,当j=5的时候，avereig(5)就不在Ef(12)和Ef(13)之间，就退出，看下avereig(5)在不在Ef(13)和Ef(14)之间，在的话就统计ds。
q1=1
do i=q1,450000-1
do j=q2,4*m*n
if((Ef(i)<=avereig(j)) .and. (avereig(j)<Ef(i+1))) then
ds(i)=ds(i)+1
else
q1=i+1 !不用让i从1开始循环了，从i+1开始循环，前面那些已经检查过了
exit
end if
end do
end do
我来回复
您所在位置：