Fortran 2023
Intel Frotran 202.0.0 は、少しですが Fortran 2023 対応を始めたので、ちょっと試してみることにします。
度数による三角関数
従来の Fortran の三角関数は弧度によるものでした。古代メソポタミア文明より伝わる 60 進法の度数を用いる場合、円周率を使ってラジアンに変換する必要がありました。これは大した手間ではないですが、数値計算的には度数用のルーチンを作った方が計算量も少なく最後の一ビットまで正確に計算できることが知られていて、よく数値計算ライブラリに度数用の三角関数が用意されておりました。
今回、Fortran 組み込み関数として採用されることになりました。
以下では、角度 90 度のところで、度数と弧度とで計算した結果を比べたいと思います。
program test2024 implicit none real, parameter :: pi = 4.0 * atan(1.0) real :: deg, rad deg = 90.0 rad = deg / 180.0 * pi print *, sind(deg), sin(rad) print *, cosd(deg), cos(rad) print *, tand(deg), tan(rad) end program test2024
1.000000 1.000000
0.0000000E+00 -4.3711388E-08
Infinity -2.2877332E+07
計算結果を見ると、従来型の三角関数では本来 0 と ∞ になるべきものが、有効桁の分のずれが出ていることが分かります。
do concurrent の reduction 演算
reduction 演算とは配列の次元が減るような演算、たとえば sum や product のような演算を指す言葉です。
do concurrent 命令は、同時並列平行に実行しても結果が変わらないことを指す反復命令で、反復される命令がループ変数に対して互いに独立していることが要求されていました。しかし OpenMP の parallel loop で reduction ができることに鑑みて、Fortran の標準文法で、これに相当すものを実現することになりました。
反復される命令が独立にならないので do concurrent の後ろに reduce 指定子を置きます。OpenMP のような塩梅で、実行も OpenMP に還元されるようです。
収束の遅い$\zeta(2) = {\pi2\over 6}$を計算してみることにします。
program reduc use, intrinsic :: iso_fortran_env implicit none real(real64), parameter :: pi = 4 * atan(1.0_real64) integer, parameter :: n = 10**7 real(real64) :: s, x(n) forall(integer :: i = 1:n) x(i) = 1.0_real64 / i s = 0.0_real64 do concurrent(integer :: i = 1:n) reduce(+:s) s = s + x(i)*x(i) end do print *, s, pi**2 / 6, s - pi**2 / 6 end program reduc
コンパイル時のオプションとして OpenMP を指定する必要があります。またこの時、stack overflow が出るので Linker のオプションで stack を 10**8 位に大きくとることにしました。
1.64493396684830 1.64493406684823 -9.999992989229156E-008
