http://www.fortran.bcs.org/2007/jubileeprog.php
Bugs I Have Known and Loved - Ron Bell, AWE Aldermaston
http://www.fortran.bcs.org/2007/jubilee/bugs.pdf
常数も変数のようにメモリー上に置かれ、かつ最適化すると共有されてしまいます。ところが、引数で常数を書き換えるようなサブルーチンに値が渡されると、共有されていた値も変わってしまって、色々面倒を起こします。
昔の FORTRAN では、同様にパラメータ文で定義した常数もサブルーチン側で書き換わります。
この仕様は昔引っかかっりました。配列とスカラーを区別するためにコロンをつけておいたらスピード低下、原因は部分配列ではコピーが作られることで、コロンがつくとたとえ全配列渡しであってもコピーが作られてしまうという罠でした。最近のコンパイラでは起きないような気がします。
日本 NAG のサイトでやっている Fortran 検定(http://www.nag-j.co.jp/fortran/exam/index.html)をやってみました。文字列がらみの問題を軒並み間違ってしまいました。色々勘違いしている。ﾔｳﾞｧｲｗｗｗ
文字変数は読み込み時は左詰め、書き出し時は右詰めｗ
FORTRAN77 規格には倍精度複素数はないが、Fortran90 以降にはある。
いつもメールで通知が来る少し前に上がっているような気がします。今回は主に bug fix で、新しい機能などは入らなかったようです。
updateも１０まで来たので、そろそろヴァージョン番号が上がって Fortran2003 完全対応になる日も近いでしょうかね？
暇つぶしにネットにあったコネクションマシンのCM-Fortranのマニュアルを見ていたら、array constructor の節に A = [1:10:2] に類する記述がありました。よもやと思って Intel Fortran で試したところ通るｗｗｗ 前々から出来たらいいなと思っていた機能なので、今まで気づかなかったのかと焦りました。しかし調べてみるとF2008でも非標準の機能のようです。
ただDEC Fortran(Ver.5)の時代から非標準の拡張として備わっていたようです。この三連符はarray constructorの表記としてきわめて自然な気がするので、是非とも標準規格に入っていただきたいところです。
コネクションマシンのCM-FORTRANのF77からの拡張は、ほぼFortran90の配列操作の部分で、例題などを見てもそのままF90で通りそうな感じです。
http://www.sthmuseum.org/downloads/CM5/GettingStartedinCMFortran.pdf
http://people.csail.mit.edu/bradley/cm5docs/CMFortranProgrammingGuide.pdf
参考： Bryan Carpenter The Development of Data-Parallel Programming
http://www.hpjava.org/talks/beijing/hpf/introduction/introduction.html
program test implicit none integer :: a(6) = [1:12:2], n print *, a print * n = 10 print *, product([1:n]) stop end program test
前回のものに比べて、必要なワーク配列を減らしました。ここでは最悪の場合に必要になるサイズ、すなわち並べ替えるデータと同じ大きさの整数配列二つを用意しています。しかし、実験的に調べてみると。ランダムなデータの場合、並べ替える数が十分に増えると２０％以下しか実際は使われません。
（最悪の場合とは、ソートの途中で[2,1,4,3,6,5,....N,N-1]のような並びが現れる状況だと思います。）
クイックソートの分割が十分に進んだ段階で、挿入ソートなり別の方法を使うようにすれば、必要なワーク配列を2N/nにできると思われます。ここでNは並べ替えるデータの総数、nはソート方法を切り替えるデータ数の閾値。
再帰版に比べて１０倍くらい速いようです。90で書いていますが77にもすぐ直せる命令しか使っていません。
module m_quick implicit none contains subroutine qsort(x) real, intent(in out) :: x(:) integer, allocatable :: indx(:), iwk(:) integer :: n allocate( indx(size(x)) ) allocate( iwk (size(x)) ) n = size(x) call qsort2(n, x, indx, iwk) return end subroutine qsort subroutine qsort2(n, x, indx, iwk) integer, intent(in) :: n real, intent(in out) :: x(n) integer, intent(out) :: indx(n), iwk(n) integer :: i, j, kk0, kk1, k0, k1, ix0, ix1, iindx, nindx real :: tmp if (n < 2) return indx(1) = 1 indx(2) = n !size(x) nindx = 2 do if (nindx == 0) exit iindx = 0 do j = 1, nindx - 1, 2 kk0 = indx(j) kk1 = indx(j + 1) k0 = kk0 + 1 k1 = kk1 do if (k0 > k1) exit do ix0 = k0, k1 ! L if (x(ix0) > x(kk0)) exit end do do ix1 = k1, k0, -1 ! R if (x(ix1) < x(kk0)) exit end do if (ix0 >= ix1) exit tmp = x(ix0) !swap x(ix0) = x(ix1) x(ix1) = tmp k0 = ix0 + 1 k1 = ix1 - 1 end do tmp = x(kk0) ! CSHIFT(x(kk0:ix1), 1) do i = kk0, ix1 - 1 x(i) = x(i + 1) end do x(ix1) = tmp if (kk0 < ix1 - 1) then ! L iindx = iindx + 1 iwk(iindx) = kk0 iindx = iindx + 1 iwk(iindx) = ix1 - 1 end if if (kk1 > ix0) then ! R iindx = iindx + 1 iwk(iindx) = ix0 iindx = iindx + 1 iwk(iindx) = kk1 end if end do nindx = iindx indx(1:nindx) = iwk(1:nindx) end do return end subroutine qsort2 end module m_quick module m_kahan implicit none contains pure function sum_kahan(xx) result(s) real, intent(in) :: xx(:) real :: c, s, x, t integer :: i c = 0.0 s = 0.0 do i = 1, size(xx) x = xx(i) - c t = s + x c = (t - s) - x s = t end do return end function sum_kahan end module m_kahan program main use m_quick use m_kahan implicit none integer, parameter :: nmax = 10**7 real :: x(nmax), s0, s1 REAL(8) :: t0, t1 integer :: i call RANDOM_SEED() do i = 1, 1000 call RANDOM_NUMBER(x) s0 = sum_kahan(x) call CPU_TIME(t0) call qsort(x) call CPU_TIME(t1) s1 = sum_kahan(x) print *, t1 - t0, ANY( x(:size(x) - 1) > x(2:) ), s0 - s1, s0 end do stop end program main
NaNを入力とする時の関数の出力は、規格上の仕様が定まっていないらしくベンダー依存のようです。Intel Fotranの場合、ほとんどの関数はNaNが出力になるようですが、SIGN関数は符号ビットだけをみているようで、NaNから実数へよみがえることができるようです。
追記：ifort は負号と関係なし、gfortran は符号の逆。sign(a, b) の定義は、b の符号に従ってその符号をつけた a を返すというものなので、b が NaN でも a のなにがしかを返すのでしょう。しかし NaN は符号を持たないのでベンダー依存で適当な符号をつけているようです。
±無限大を引数とする場合、ATANは予期する答えを与えてくれるようです。
IEEE745規格にのっとるようになってから、実数ではゼロ割りで憤死！ということが少なくなったわけですが、整数ではまだゼロ割りで憤死します。（当然といえば当然ですが・・・）
program test implicit none real, parameter :: pi = 4.0 * atan(1.0) real :: x x = 1.0 / 0.0 print *, 'x =', x, 'atan(x) = ', atan(x) / pi x = 1.0 / -0.0 print *, 'x =', x, 'atan(x) = ', atan(x) / pi print * x = 0.0 / 0.0 print *, 'x =', x, 'sign(1.0, x)=', sign(1.0, x) print * x = 0 / 0 print *, 'x = 0/0 ', x stop end program test
本の出版はまだですが、例題プログラムのページが未完成ながらできているようです。
http://flibs.sourceforge.net/examples_modern_fortran.html
