我有一个有数千个双精度实数的网格。
它正在迭代,当它收敛到3位小数时,我需要它停止。
目标是让它尽可能快地运行,但每次(到3DP)都需要给出相同的结果。
在我做这样的事情的那一刻
REAL(KIND=DP) :: TOL = 0.001_DP
DO WHILE(.NOT. CONVERGED)
CONVERGED = .TRUE.
DO I = 1, NUM_POINTS
NEW POTENTIAL = !blah blah blah
IF (CONVERGED) THEN
IF (NEW_POTENTIAL < OLD_POTENTIAL - TOL .OR. NEW_POTENTIAL > OLD_POTENTIAL + TOL) THEN
CONVERGED = .FALSE.
END IF
END IF
OLD_POTENTIAL = NEW POTENTIAL
END DO
END DO
我在想,很多IF语句不会对性能有太大的影响。我想过在最后检查收敛性;找出平均值(对整个网格求和,除以num_points),然后检查它是否以与上面相同的方式收敛,但我不相信这总是准确的。
做这件事的最好方法是什么?
发布于 2010-12-17 12:20:48
如果我没理解错的话,你有某种时间步进,你在new_potential
中通过在old_potential
上的计算来创建值。然后让旧的等同于新的,继续前进。
您可以使用单个语句替换现有的收敛测试
converged = all(abs(new_potential - old_potential)<tol)
这样可能会更快。如果测试的速度是一个主要问题,你可以每隔一次测试一次(或者每三次或四次...)迭代
下面是一些评论:
1)如果您使用具有2个平面的潜在数组,而不是old_和new_potential,则可以通过在每次迭代结束时交换索引来将new_转换为old_。在您的代码中,会有大量的数据移动。
2)虽然在语义上使用While循环是正确的,但我总是使用迭代次数最多的do循环,以防永远不能满足收敛标准。
3)在您的声明REAL(KIND=DP) :: TOL = 0.001_DP
中,关于TOL数值的DP规范是多余的,REAL(KIND=DP) :: TOL = 0.001
就足够了。我也会把它设为一个参数,如果编译器知道它是不可变的,它可能会优化它的使用。
4)您实际上不需要在最外层的循环中执行CONVERGED = .TRUE.
,在第一次迭代之前设置它--这将为您节省一两纳秒。
最后,如果您的收敛标准是潜在数组中的每个元素都收敛到3dp,那么这就是您应该测试的。为你建议的平均值构造反例相对容易。然而,我担心的是你的系统永远不会在每个元素上收敛,你应该使用一些矩阵范数计算来确定收敛性。因此,在这个主题中不是一个教训的地方。
发布于 2010-12-17 18:02:19
收敛标准的计算方法是什么?除非它们比计算更糟糕,否则让IF语句尽快终止循环可能会更好,而不是猜测非常大量的迭代,以确保获得良好的解决方案。
Re High Performance Mark的建议#1,如果复制操作是运行时的重要部分,您也可以使用指针。
确定这些东西的唯一方法是测量运行时间……Fortran提供了测量CPU和时钟时间的内部函数。否则,你可能会修改你的代码的某一部分,使其更快,可能会使它更不容易理解,并可能引入错误,可能不会在运行时有太多改进……如果这一部分只占总运行时间的一小部分,那么再多的聪明也不会有多大的不同。
正如High Performance Mark所说,尽管当前的语义很优雅,但您可能希望防止无限循环。一种方法是:
PotentialLoop: do i=1, MaxIter
blah
Converged = test...
if (Converged) exit PotentialLoop
blah
end do PotentialLoop
if (.NOT. Converged) write (*, *) "error, did not converge"
发布于 2011-03-02 22:03:46
I = 1
DO
NEWPOT = !bla bla bla
IF (ABS(NEWPOT-OLDPOT).LT.TOL) EXIT
OLDPOT = NEWPOT
I = MOD(I,NUMPOINTS) + 1
END DO
也许更好
I = 1
DO
NEWPOT = !bla bla bla
IF (ABS(NEWPOT-OLDPOT).LT.TOL) EXIT
OLDPOT = NEWPOT
IF (I.EQ.NUMPOINTS) THEN
I = 1
ELSE
I = I + 1
END IF
END DO
https://stackoverflow.com/questions/4469847
复制