第二回
2002/10/19

『シミュレーション概論』 講義メモ No.2 (10/21)

1.計算機の限界と数値計算

 計算機のメモリ(作業領域)と処理速度は有限であり,その制約のもとで計算を設計しなければならない。

<メモリ>
  普通のパソコンやワークステーションで作業に使用できるメモリは,現在100メガバイト程度である。このうち,OSや実行しようとしているプログラムそのもの,他の常駐プログラムなどの使用している分もあるから,計算作業で一時的に使用できるメモリはその一部分となる。

  代数方程式の解を求めるとか,1変数の関数の定積分を求めるような場合には,ほとんどメモリのことを気にする必要はないが,3次元領域での偏微分方程式の解を求めるような計算になるとメモリの制約が現れる。

  たとえば立方体状の連続領域を各次元方向を100等分したメッシュで計算する場合なら,記憶させておかなければならない要素数(変数の数)は1003=100万個である。1要素に対して実数変数1個(=4バイト:後述)が対応している場合なら,要するメモリは400万バイト,およそ4メガバイトである。これをもう少し精密に計算しようと思って各次元で4倍の400等分にした場合,要するメモリは64倍すなわち256メガバイトとなり,パンクしてしまう。

<計算速度>
  計算機の処理速度は主として使用されているCPUの性能で決まる。CPUの性能は,1ビットを書き換える速度で表され,例えば「ペンティアム4の1GHz」といえば,1秒間に0と1のビット変換を10億回実行できることを示す。数値計算では,足し算やかけ算など1回の基本演算にかかる時間の目安と考えてよい。

  パソコンを含めて最近のコンピュータなら,簡単な計算の場合,処理速度のことを気にするよりも,わかりやすいプログラム,間違いのないプログラムにすることの方に神経を使うべきである。しかしながら大がかりな数値計算の場合,プログラムの仕方によって数倍の差「半日で結果が出るものが1週間かかる」となると,卒論の締め切りに間に合うかどうか,という深刻なことにもなる。

  (問1)次の多項式の計算の場合,足し算とかけ算を何回必要とするか? 

      1 + 2*X + 3*X^2 + 4*X^3 + 5*X^4 + 6*X^5   (^はベキ乗,Fortranでは**)

    これを次のように書き換えたら何回になるか?

      1 + X* (2 + X* (3 + X* (4 + X* (5 + 6*X))))    (これならCでも書ける)


  さらに,計算の途中結果をディスプレイやハードディスクに書き出すような場合には,データの転送に多くの時間がかかることを知っておくべきであろう。

  大がかりな数値計算の場合,途中で予期しないエラーが起きると,それまでの計算がムダになる。このような場合,途中の結果をハードディスクに記録しておき,次はそこから再開するようにすればよい。しかしながらこれをあまり頻繁に行うと計算時間を大幅に遅らせてしまうことになる。

  (問2) [0,1)間の大量の実数値が次々に与えられるとき,(0,1)を100等分した区間に分けて,各区間の値が現れた頻度を整数配列 a(0:99)(要素数100)に記録しておき,あとでヒストグラムで描きたい。実数値 x(0≦x<1)が与えられたとき,これがどの区間に入るかを判定するにはどうしたらよいか? 

<並列計算>
  現在の最先端の大型計算機は CPU が多数組み込まれていて,これを並行して動作させることができるようになっている。アルゴリズムの中に,並行して同時に計算することができる部分を含んでいる場合には,並列処理できるようにプログラムすれば計算は格段に速くなる。

  次のような偏微分方程式(1次元拡散方程式)を解くことを考えてみよう。

     ∂a/∂t = ∂2a/∂x2

最も素朴な方法は,時間を dt,空間を dx の区間に分け,変数 a を配列とし

     d = dt/(dx**2)
     a(n) = a(n) + d * ( a(n-1) - 2*a(n) + a(n+1) )   (※注)

               ↑= { [a(n+1)-a(n)]/dx - [a(n)-a(n-1)]/dx }/dx

を n=0,1,2,...,N について実行することである。この場合,全ての位置(n)において,現在の両隣の値を見て次の瞬間の値を「いっせーのーで」決めるようにしなければならないから,この式をそのまま n=0,1,2,...,N と順番に(逐次的に)実行すると,要求されている計算とは異なることを遂行することになってしまう。つまり並列計算を必要とする例である。

  Fortran90ではこのことを意識して,並列処理できる部分を(バッファ配列を用意することなく)1行の文で書けるような文法が開発された。身の回りにあるパソコンやワークステーションは,たいていCPUを1個しかもたないものであり,実際に並列処理が実行されるわけではないが,将来おおがかりな計算をやらなければならなくなったときのために,並列計算の意味を理解しておこう。

(※注) このまま鵜呑みにして使用しないこと! dt/(dx)^2 が大きすぎると,あっという間に発散する。このため,dx をあまり小さくできないという制約があり,「何でも詳しければよいというものでもない」例である。


2.計算機におけるデータ表現

(1) 文字データ
    アルファベット大小文字,数字,記号などを,0-255 の間に番号付け(コード化),1文字=1バイトで表現する。
    漢字体系は1文字2バイト(全角文字)

(2) 論理データ
    {True, False} の2値しかとらないが,普通,4バイトデータで表現される。

(3) 整数データ
    Fortranでは標準で(Cではlongint),32桁の2進数として4バイトで表現する。
    負の数には「232の補数表現」が用いられ,引き算を足し算で行う。

  したがって,表現できる整数の範囲は  -231≦I≦231-1 である。

  この範囲を超えた整数の入力は普通拒否されるが,計算の途中で越えても(オーバーフロー)エラーメッセージは出ず,単に上位にオーバーフローしたビットを切り捨てるだけで作業が進むことも多いので注意が必要である。

<整数の内部表現を求める>Fortranプログラム Windows用実行ファイル

INTEGER :: j, ia=1
DO WHILE(ia /= 0)
   PRINT *,'Input an integer (or 0 to stop):'
   READ *, ia
   PRINT '(5X, 4(8I1, 1X) )', (IBITS(ia, j, 1), j=31, 0, -1)
END DO
END

IBITS(i,j,k) は,整数 i の第 j ビット目(下位から数える)から k ビットを取り出す関数。

(4) 実数(浮動小数点数)
  0.0以外の実数を,f×βn (1≦f<β)の形で表して,{符号,指数 n,仮数部 f} を4バイト(Fortranの標準実数,Cのfloat)で表すが,システムによって表現方法が異なる。

<実数の内部表現を求める>Fortranプログラム Windows用実行ファイル

REAL :: a=1.0
CHARACTER :: ic*11
INTEGER :: ia, j
DO WHILE(a /= 0.0)
   PRINT *, 'Input a real (or 0.0 to stop):'
   READ *, a
   WRITE(ic,'(I11)') a;  READ(ic,'(I11)') ia  ! ※
   PRINT '(5X,I1,1X,8I1,1X,23I1)', (IBITS(ia,j,1), j=31,0,-1)
END DO
END

※ FortranのI編集記述子(I変換)を用いて,実数のビット表現と同じビット表現を持つ整数(のアスキー表現)を内部ファイル(=文字変数ic)に無理矢理に書きこみ,これを整数変数として読み出してから,ビット表現を出力する,というアクロバットな方法。

これで調べてみると,ここのパソコンやワークステーションでは,x = ±f×2n-127 として

    先頭の1ビットで符号(正なら0,負なら1)
    次の8ビットで 2n-127 の指数 n を
    残りの23ビットで仮数部 f の2進数小数部 を

表現しているようである。(0.0 だけは,全てのビットが0) そうすると,この表現で近似的に表される実数は

    有効数字は6桁 (∵ 2-23=1.19×10-7
    絶対値の範囲は10-37〜1037

である。

 表現の基本は2進数であるから,10進数で短い有限小数であっても,2進数表現では無限小数になることを知っておく必要がある。例えば10進数の 「ちょうど 0.1」 でも,2進数では 0.0001100110011001.... と循環小数になるから有限2進数では近似的にしか表現できない。このため0.1 を10回足しても決して「ちょうど1.0」にはならない。

 もっと有効数字の多い数値あるいは絶対値の大きい(小さい)数を扱う必要がある場合には,「倍精度(8バイト長)」(Fortran では REAL(8),Cでは double)の実数が用意されている。各型の実数で扱える絶対値の範囲,バイト長,有効数字の桁数はシステムによって異なるが,Fortran90にはこれらの値を得る関数が用意されている。(変数 x に4バイトまたは8バイトの型宣言が行われているとき,それぞれ,RANGE(x),KIND(x),PRECISION(x) で与えられる。)

  (問) 倍精度実数で扱われる絶対値の範囲,有効数字の長さを調べて見よ。
      倍精度は「REAL(8):: x 」または「DOUBLE PRECISION x」で宣言する。
   
3.計算機の特性に起因する誤差についての注意

  以上で述べたように実数値は近似的にしか扱われない。これによる誤差を丸め誤差という。特に以下のような場合にはこのことが如実に現れ,所定の有効数字よりも遙かに粗い,予期せぬ誤差を生じるので注意が必要である。さらに詳しい説明,特にC言語のプログラムについては,早川先生の講義ノート(第2講)を参照のこと。

(1)絶対値が大幅に異なる二つの数の加減算
  有効数字が6桁のときに,例えば 1.23456 に 0.0000001 を加えても 1.2345601 にはならない。手計算の場合には有効数字の考えを無視したこういう計算をすることはないが,プログラムを組んで何回も足し算を繰り返すような場合(例えば数値積分など)に,うっかりするとやってしまう。

試しに,トリビアルな例として 1.0/n を n 回,順番に足しあわせていくプログラムを実行してみよう。

REAL :: s, x
INTEGER :: i, k, n
DO k=1, 8
   n=10**k
   x=1.0/n
   s=0.0
   DO i=1, n
      s= s + x
   END DO
   PRINT *, n, s
END DO
END
           n          s 
          10   1.000000    
         100  0.9999993    
        1000  0.9999907    
       10000   1.000054    
      100000   1.000990    
     1000000   1.009039    
    10000000   1.064767    
   100000000  0.2500000    

1/1000 くらいでも結果の有効数字は怪しくなっているであろう。

最適化  コンパイルの際にプログラムの流れを先読みし,このような誤差を最小限に押さえるよう改良してしまうオプションを「最適化」という。システムによっては最適化がデフォルトになっている場合がある。たとえばセンターのパソコンに入っているDVFはそうである。オプションを何も指定しないでコンパイルすると,上のような結果にはならず,この範囲ではほとんど誤差は見られない。上のプログラム名を「test.f90」として,コンパイルの際に「実数処理の最適化」を解除するため

   DF test.f90 -Op (オーピィ)

としてから実行してみること。プログラムどおりの処理をしてほしいときには,このように指定する。
  ワークステーションの場合は何もオプション指定しないと実数の最適化はされないようである。

   f90 test.f90

でコンパイルして実行すると,上のような結果が得られる。

  最適化によって誤差の回避や高速化などが行われるが,どのような改良がされるのかについては,よほど習熟しないと理解できないものである。原理的には,プログラムに書いたとおりに処理され,原理的には以上述べたような誤差を生じるものと覚悟しておくべきである。

(2)非常に近い二つの数の引き算(桁落ち
  例えば 0.123456 と 0.123455 の引き算を行うと,差は 0.000001 で,当然のことながら有効数字は6桁ではなくなる。筆算では分かっていても,プログラム上は気が付かないことが起こり得る。
  教訓的な例は2次方程式の解の計算である。

     x2 - 2bx + c = 0   (簡単のため,b, c >0 とする)

の解は
     x = b ± (b2-c)1/2

であるが,もし b≫c であれば,小さい方の解 x1 非常に小さくなり,相対的な誤差が大きくなる。すなわち,大きい方の解 x2 の有効数字との間にアンバランスを生じる。この場合,解と係数の関係を用いて

     x1 = c / x2

とすれば,小さい方の解の有効数字が失われることはない。
  同じ種類の誤差は,数値微分の場合にも生じる。連続関数 f(x) の点 x における導関数の値 f'(x) は

     [f(x+h/2) - f(x-h/2)] / h 

で近似的に与えられるが,数値計算においては h を小さくすればするほど近似が良くなるとは言えないことが理解できよう。
  ものは試し,f(x) = x^2 / 2 の x=1 における導関数の値(理論上は1である)を,h=1/n として計算してみよう。

REAL :: f, h, df
INTEGER :: k, n
f(x) = x**2 / 2
DO k=1, 8
   n=10**k
   h=1.0/n
   df=( f(1.0 + h/2.0) - f(1.0 - h/2.0) ) / h
   PRINT *, n, df
END DO
END

やはり「最適化」しないでコンパイルした結果は以下のとおりである。やはり n=1000 くらいから有効数字が大幅に失われていていることがわかる。
           n         df
          10  0.9999996    
         100  0.9999990    
        1000  0.9999871    
       10000  0.9995699    
      100000   1.001358    
     1000000  0.9536743    
    10000000  0.5960464    
   100000000  0.0000000E+00

※ DVFでデフォルトの最適化コンパイルでは,以上の範囲ではすべて結果は1.0であった。(Fortranを使える人は両方を試して見よ。)

  このように計算機を用いた数値計算では,なんでもかんでも精密に計算すればよいというものではないことを理解しておく必要がある。