[thread display] [new arrival display] [word search] [past log] [管理用]
| Subject | ： Re^2: SAD Update. V1.0.10.4.17a40 Bug in TrackParticles with RAD and NPARA >1 |
| Date | ： 2011/08/10(Wed) 11:00:14 |
| Contributor | ： Akio Morita |
> Dear Users,
>
> 1. Another strange behavior was found in the logic of EPS in BEND with RAD. The division number did not necessarily increase/decrease for smaller/larger values opt EPS.
以下、重箱の隅
この部分のコードで、分割数の上限(旧4/新1024)が設定されていますが、根拠は有るのでしょうか？
tbend()自身に再帰しているので、最終的には ndiv == 1を満足するところまで再帰分割されるので、
call stack depthを考えると一度に分割した方が効率が良いように思えます。
もし、上限を設定するとすれば、Integer型の表現域の上限が根拠になると思います
Fortran95的には変数から、 上限 RADIX(ndiv)*(RADIX(ndiv)**(DIGITS(ndiv)-1) - 1)+1を導出できます
# RADIX(ndiv)**DIGITS(ndiv) - 1が途中の計算でオーバーフローするのを避けるために面倒な表現に…
# コンパイル時に確定する定数なので、まともなコンパイラなら実行時の負荷は増えません
また、上限の制約を int()で変換された後に行いますが、Real型の値
(1 + (abs(phi0) * max(h0/100.d0,1.d0/rphidiv) / eps1))がInteger型の
表現域を越えている場合、どのような値に変換されるか保証はされていなかったと思います
したがって、max()はReal型でとる必要が有ります int(max(dble(ndivmax), ...))
さらに、例外状況を考えると、eps1で除算した項がNaNである可能性が排除出来ないので、
NaNを含む入力に max()が何を返すか保証がなければ、int()にNaNが渡る可能性を考慮する必要があります。
現状の実装で起こりうるケースとしては、EPSが 十分小さな正の値の時、
(... / eps1)が Integer型の表現域を越える大きな値 or +Infinityになり、
int()での変換結果が INT_MAXに張り付かない実装では、意図通りに
再帰分割が働かないことになります。
例えば FreeBSD/amd64環境の gfortran 4.5.4 20110721では、
* int(1.d0 / 1.d-300)は -2147483648になります
* int(NaN)は、-2147483648になります
同様に、BEND B = (EPS = 0/0);で EPSに NaNを設定可能です。
意図的に設定しなくても、何らかの計算結果から EPSパラメータをセットする場合に NaNになる可能性を排除出来ません