プログラミングと数値計算

4倍精度演算(fortran90)

コメントする

fortran90において4倍精度で計算を行います。

メリット
情報落ちの回避
桁落ちの回避

デメリット
計算速度が遅い(倍精度の40~50倍!)
4倍精度に対応しているライブラリが少ない

計算式を工夫して、どうしても4倍精度でないと解けない場合のみ4倍精度を使うことを考えましょう。

  1. 4倍精度の存在理由
  2. 4倍精度演算の組み込み関数
  3. 4倍精度陽的ルンゲクッタ法(刻み幅制御)
  4. 4倍精度陰的ルンゲクッタ法(刻み幅制御)
  5. 4倍精度の一般行列の対角化(複素非エルミート)
  6. 4倍精度のLU分解による連立方程式

4倍精度の存在理由

4倍精度は倍精度演算を正しく行うための補助的存在です。
初めから最後まで4倍精度で行うことはほとんどありません。
4倍精度が効果を発揮する桁落ちの例を考えましょう。

計算
\(
1-(1+\frac{10^{-10}}{3})
\)
を考えます。
答えは\(-3.333333\cdots\times 10^{-11}\)です。
数値計算で倍精度と4倍精度で試してみますと、

program main
  implicit none
 
  write(6,*)1d0-(1d0+(1d-10)/3d0)
  write(6,*)1q0-(1q0+(1q-10)/3q0)
 
  stop
end program main

であり、実行結果は

$ ifort main.f90
$ ./a.out
 -3.333333609134570E-011
 -3.333333333333333333333342109654793E-0011
$

となります。計算過程で工夫もできないような場合、4倍精度を利用して

program main
  implicit none
  double precision::a
  real*16::b
 
  b=1q0-(1q0+(1q-10)/3q0)
  a=dble(b)
  write(6,*)a
 
  stop
end program main

として途中だけ4倍精度にして計算する、これが4倍精度の存在理由です。
上記プログラムを実行すると

-3.333333333333333E-011

となります。

想定

4倍精度演算は決まった書き方が無く、コンパイラに依存します(補足1)。
ここでは、
intel®fortran コンパイラ(ifort (IFORT) 16.0.2 20160204)
と
gnu fortran コンパイラ(GNU Fortran (Ubuntu 4.8.4-2ubuntu1~14.04.1) 4.8.4
)
の二つで動作することを考えます。

4倍精度変数の宣言方法

4倍精度実数、4倍精度複素数は

real*16::a
complex*32::c

と変数宣言をしましょう。
値を代入する際は、

a=1q0
c=cmplx(1q0,2q0,kind=16)

として”q”を用います。
※”c=cmplx(1q0,2q0)”としてはダメです。
複素数への代入はcmplx(実部、虚部、精度指定)でないといけません。
ifortだけであれば、qcmplxで自動的に精度が4倍精度指定されますが、gfortranではkindで指定しないと代入されません。実際、

program main
  implicit none
  complex*32::z

  z=cmplx(1q0/3q0,1q0)
  write(6,*)z  
 
  z=cmplx(1q0/3q0,1q0,kind=16)
  write(6,*)z

  stop
end program main

を実行すると、gfortranでもifortでも

$ ifort main.f90
$ ./a.out
 (0.333333343267440795898437500000000,1.00000000000000000000000000000000)
 (0.333333333333333333333333333333333,1.00000000000000000000000000000000)
$
$ gfortran main.f90
$ ./a.out
 ( 0.333333343267440795898437500000000000      ,  1.00000000000000000000000000000000000      )
 ( 0.333333333333333333333333333333333317      ,  1.00000000000000000000000000000000000      )
$

という結果を得て、正しく値が代入されません。

[adsense1]

4倍精度演算のための組み込み関数

基本的には総称名を用いるのが良いと思います[1,2]。

a=sin(4q0)
c=exp(cmplx(3q0,1q0,kind=16))

倍精度から4倍精度に変換する際は、

a=sin(4d0)*1q0
c=cmplx(exp(3q0),kind=16)

とすると良いでしょう。あとは他の精度の時と同じです。

ルンゲクッタ法の4倍精度刻み幅制御

4倍精度でルンゲクッタ法を行います。(刻み幅制御ルンゲクッタ法についてはルンゲクッタ法の説明と刻み幅制御をご覧ください。)
対象とする微分方程式は
\(
\displaystyle \frac{d^2y}{dx^2}=\left(-\frac{1}{4x^2}+\frac{i}{x}-1\right)y
\)
を初期条件
\(
\begin{align}
\left. y(x)\right|_{x=1}=\exp(i) \\
\left. \frac{dy(x)}{dx}\right|_{x=1}=\left(\frac{1}{2}+i\right)\exp(i)
\end{align}
\)
の条件下で解くものです。解析解は\(y(x)=\exp(ix)\sqrt{x},\frac{dy(x)}{dx}=\left(\frac{1}{2x}+i\right)\exp(ix)\sqrt{x} \)となります。
メインのコードは

program main
  use qRKmod
  implicit none
  integer::i,N,Nx,info
  real*16::dx,tol,xbound,x0,tx
  real*16,allocatable::x(:)
  complex*32,allocatable::y(:),ss(:)
  complex*32::rkfc
  external::rkfc
 
  N=2
  allocate(y(1:N))
 
  x0=1q0; dx=1000q0; Nx=2
  allocate(x(0:Nx-1))
  do i=0,Nx-1
     x(i)=x0+dble(i)*dx
  enddo

  y(1)=exp(cmplx(0q0,1q0,kind=16))
  y(2)=cmplx(0.5q0,1q0,kind=16)*y(1)

  call qrk_preparation("rkf45")

  tol=1q-18
  write(10,'(3e30.17e4)')x(0),real(y(1)),imag(y(1))
  do i=1,Nx-1
     info=0; tx=x(i-1); xbound=x(i)
     call qrkf45_e(rkfc,tx,y,xbound,info,tol)
     if(info.eq.-1)write(6,'(A,f10.5,A,f10.5)')"strange point between ",x(i-1)," and ",x(i)
 
     write(10,'(3e30.17e4)')x(i),real(y(1)),imag(y(1))
  enddo
  write(6,*)x(0),y(1)
 
  call qrk_deallocation("rkf45")
 
  stop
end program main
!----------------
function rkfc(N,x,y,s)
  implicit none
  integer,intent(in)::N,s
  real*16,intent(in)::x
  complex*32::y(1:N)
  complex*32::rkfc
 
  rkfc=0q0
  if(s.eq.1)then
     rkfc=y(2)
  elseif(s.eq.2)then
     rkfc=(-1q0/(4q0*x*x)+cmplx(0q0,1q0,kind=16)/x-1q0)*y(1)
  else
     write(6,*)"unknown s, program stop"
     stop
  endif
 
  return
end function rkfc

で、4倍精度用のルンゲクッタ法のモジュールはこちら。

1ステップ当りの精度を18桁に設定し、実行すると
実行すると、

$ ifort qrkmod.f90 qmain.f90
$ ./a.out
   1.00000000000000000000000000000000      
 (-12.4004402201339396388318226191018,29.1071998369283977905998988164457)
$
$ gfortran qrkmod.f90 qmain.f90
$ ./a.out
   1.00000000000000000000000000000000000       ( -12.4004402201339396388318226207936776      ,  29.1071998369283977905998988157304999      )
$

という結果を得ます。
ここに表示されているのは、解析解の\(x=1001\)での値で、
-12.400440220133925015780697384759199505100325758422733906… +
29.107199836928403740368952842025341239069669658649312423… i
が厳密値となります。
x=1001で数値計算解と比較すると15桁一致していることが分かります。

計算時間は倍精度計算の約40倍かかりますが、どうしても16桁近い精度が欲しい場合や桁落ちが激しいことが分かっている場合には、途中を4倍精度で演算して後で倍精度に戻す、と言う方法が良いでしょう。

[adsense2]

計算時間について

ルンゲクッタ法で、倍精度、4倍精度を比較します。
中身が複雑なのであまり良い比較方法ではありませんが、目安としてどの程度か分かるでしょう。
CPU時間を計測します。

上記微分方程式を10回解いたとき、
倍精度では
0.298 [秒]
4倍精度では
11.268[秒]
でした。
約40倍遅いことが分かります。
計算が工夫できるならするべきです。むやみに4倍精度を使ってはなりません。

補足1)コンパイラに依存してしまうとは？

例えば、
\(1/3\)を計算してみましょう。

program main
  implicit none

  write(6,*)qext(1d0/3d0)
 
  stop
end program main

”qext”は倍精度→4倍精度の型変換をする組み込み関数です。
このコードをifortで行うと計算が実行され、
実行結果

$ ifort main.f90
$ ./a.out
  0.333333333333333314829616256247391  
$

を得ます。しかし、同じコードをgfortranでコンパイルしようとすると、

$ gfortran main.f90
main.f90:4.12:

  write(6,*)qext(1d0/3d0)
            1
Error: Function 'qext' at (1) has no IMPLICIT type
$

となり、コンパイルが出来ません。
この場合、型変換用の組み込み関数を使わないで、

program main
  implicit none

  write(6,*)(1d0/3d0)*1q0
 
  stop
end program main

とするのが良いでしょう。
異なる型の演算は高い方の型に合わせられることを利用します。
すると、どちらのコンパイラでもokで、

$ gfortran main.f90
$ ./a.out
  0.333333333333333314829616256247390993      
$ ifort main.f90
$ ./a.out
  0.333333333333333314829616256247391      
$

となります。
gfortranとifortで表示されている桁数が違うのは標準出力の違いだけのようです。
なぜなら、組み込み関数[3]
———————————————————–
precision:10進数で表現できる桁数
epsilon:機械精度(1にこの値を足しても1のままになる最大の値)
tiny:指数表現での最小値
huge:指数表現での最大値
———————————————————–
を使って調べると、

program main
  implicit none
  real*16::a
 
  write(6,*)precision(a)
  write(6,*)epsilon(a)
  write(6,*)tiny(a)
  write(6,*)huge(a)

  stop
end program main

実行結果

$ ifort main.f90
$ ./a.out
          33
  1.925929944387235853055977942584927E-0034
  3.362103143112093506262677817321753E-4932
  1.189731495357231765085759326628007E+4932
$
$ gfortran main.f90
$ ./a.out
          33
   1.92592994438723585305597794258492732E-0034
   3.36210314311209350626267781732175260E-4932
   1.18973149535723176508575932662800702E+4932
$

という結果が得られるからです。

4倍精度の陰的ルンゲクッタ法

4倍精度の陰的ルンゲクッタ法については
陰的Runge-Kutta法
に載せてありますので、そちらをご覧ください。

4倍精度の一般行列の対角化

複素非エルミート行列の対角化については
4倍精度の一般行列の対角化
に載せてありますので、そちらをご覧ください。

4倍精度のLU分解による連立方程式の解法

4倍精度のLU分解を用いた連立方程式の解法は
九州大学の渡部 善隆様が公開なさっているGECP(Gaussian Elimination with Complete Pivoting, 一般実行列に対する連立1次方程式の数値解を完全ピボット選択付き Gauss の消去法によって求める Fortran サブルーチン)
によって実現できます[4]。
http://yebisu.cc.kyushu-u.ac.jp/~watanabe/RESERCH/GECP/index.html。
再配布可能です。

参考文献

[1]9.58 CMPLX — Complex conversion function -The GNU Fortran Compiler(gfortran complex精度指定)
[2]第 3 章 FORTRAN 77 および VMS 組み込み関数 -Sun Studio 12: Fortran ライブラリ・リファレンス
[3]浅岡 香枝、平野 彰雄、”はじめてのFortran90”
[4]渡部 善隆, 一般実行列に対する連立1次方程式の数値解を完全ピボット選択付き Gauss の消去法http://yebisu.cc.kyushu-u.ac.jp/~watanabe/RESERCH/GECP/index.html

関連記事

このエントリーをはてなブックマークに追加

コメントを残す

メールアドレスが公開されることはありません。 ※ が付いている欄は必須項目です