2013-01-08

Fortran 歴史的遺物

Fortran

Fortran は1950年代から過去との互換性を重視して、あまり機能を削らずに進化を続けたために、進化の過程で取り残された意味不明な命令がたくさんあります。メモ帳代わりに書いておきます。

DIM 関数 (DIMINISHING function)

DIM 関数とは、DIM(x, y) x - y >= 0.0 の時 DIM(x, y) = x - y; x - y < 0.0 の時 DIM(x, y) = 0.0 という謎な定義で、教科書の例題以外で使われているのも見たことがない関数です。自分で使ったこともありません。

そもそも DIM が何の略号なのかすら忘れ去られていて、IBM の FORTRAN IV のマニュアルにすら POSITIVE DIFFERENCE と書いているくらいです。最近、computer history museum のサイトの FORTRAN ARCHIVE COLLECTION にある IBM FORTRAN II マニュアルの初期の版(FORTRAN II for the IBM 704 Data Processing System : Reference Manual)に、DIM 関数の type of function を DIMINISHING としているのを見つけました。

同じ FORTRAN II のマニュアルでも後期の版では DIM 関数は positive difference と書いてあるので、なかなか見つけづらいものがあります。なお DIM 関数は FORTRAN I に存在しておらず、FORTRAN II で導入されてあっという間に名前の由来を忘れられてしまった関数のようです。。

CHAR ICHAR / ACHAR IACHAR EBCDIC と ASCII

文字の内部コードがらみ関数には機能がほとんど同じ二つの関数があります。CHAR(n) は、内部コード n の文字を返し、ACHAR は ASCII コード n の文字を返します。ICHAR、IACHARはそれぞれその逆関数になっています。

Fortran は作られたのが ASCII コードよりも古いので、元々文字コードはベンダー依存だったのですが、その後互換性問題から ASCII コード版の関数も用意されました。かつては IBM 社の EBCDIC が主要な文字コードだったものですが、今では ASCII コードが普遍的になったので、大抵の場合同じ挙動をするので、機能の重複する冗長な関数があるように見えています。

CSHIFT,　EOSHIFT 関数　　並列化（HPF) 機能としての配列シフト関数

配列シフト関数が並列化や HPF と関係して導入されたと言われると不思議な気がします。しかし、Fortran90の規格が80年代に形づくられたことを思い出して、70～80年代の並列計算機のアーキテクチャを考えるとこれが自然なものだと分かります。

例として以下に illiac のプロセッサのトポロジーを示します。プロセッサが二次元トーラス状のメッシュを形成していることが見て取れます。各プロセッサは二次元配列の１要素を受け持って計算を並列に行います。この場合には前後左右のプロセッサ間でメモリー内容を授受し合うと最も並列性良く動きます。80年代に一世を風靡したコネクション・マシンも似たようなアーキテクチャを取っていました。
f:id:fortran66:20130109014809p:plain
B. Svensson SIMD architectures http://www.hh.se/download/18.70cf2e49129168da0158000103477/1232956956122/SIMD+Architectures.pdf

2013-01-04

反復と射

Fortran

正月なので随想をば思いつくままに書きます。

Modern Fortran において留意すべき点として、FORTRAN77 まではソースコード上で区別されていなかった繰り返し構造を、概念的に大きく二つに分けて記述できる点にあると思います。その二つの区分の呼び名として適当な言葉が思い浮かばないので、とりあえず反復と射ということにしておきます。

反復と射

反復とは、普通の do loop 構造にあたるもので、loop 変数に従った順序付の繰り返し構造にあたります。

do i = 1, n
  a(i + 1) = a(i) + b(i)
end do

このように漸化式の形を取るものは、 i の順番に依存して実行結果が変わってくるので、i が 1 から始まって n まで順番に実行されることが意味を持ちます。

これに対して、ここで射と呼んでいるものは（写像のつもりなんですが）、本来対応を示しているもので、複数の操作をまとめて一気に実行できることを意味しているつもりです。

例として以下で indx が添え字の置換を表しているとして、置換にあたる操作は、

do i = 1, n
 　b(i) = a(indx(i))
end do

のように、漸化式の反復の場合と同じような do loop で記述できますが、i の順番に依らずに同じ結果が得られます。したがって i が露わに記述される理由はないことになります。また代入の実行は順番を入れ替えても、並行して行っても構わないはずです。

Fortran90での代入

Fortran90 では、このような添え字に依存しない代入の繰り返しを記述すのに、

b = a(indx)

という記法が導入されました。

条件節のある場合

条件節による分岐がある場合、FORTRAN77では以下の様に記述されます。

do i = 1, n
  if (a(i) > 0) then
    b(i) =  sqrt( a(i) )
  else
    b(i) = -sqrt( a(i) )
  end if
end do

この場合も loop 変数 i には露わに依存性はありませんし、各々の i に対する代入は、順不同・同時・並列に実行してよいはずです。

where 構文

この場合には Fortran90 では where 構文を用いて記述します。

where (a > 0.0)
  b =  sqrt(a)
else where
  b = -sqrt(a)
end where

これによって、添え字への依存性がなく順不同・同時・並列に実行してよいことが明示されます。

forall 構文

似たような状況で loop 変数に露わに依存するものの、順不同に実行できるくりかえし構造もあります。

do i = 1, n
  ia(i) = i * (i - 1) / 2
end do

このような場合 Fortran90 では forall 構文を使います。

forall (i = 1:n)
  ia(i) = i * (i - 1) / 2
end forall

これは、代入が i に関して順不同・同時・並列に実行してよいことを明示的に表現しています。

mask による分岐

条件節による分岐がある場合を、次に示します。

do i = 1, n
  if ( abs( ia(i) ) > 1.0e3  ) ia(i) = ia(i) / i  
end do

forall 構文では mask が使えます。

forall (i = 1:n,  abs(ia(i)) > 1.0e3)
  ia(i) = ia(i) / i
end forall

forall 構文には mask 条件を満たさない場合の else 節にあたるものはありません。

do concurrent 構文

forall は代入操作しかできないため、Fortran2008 では類似の構文として do concurrent 構文が導入されました。

do i = 1, n
  x = i * 0.01
  a(i) = x * log(x) / (x - 1.0)
end do

は、以下の様に書き直されます。

do concurrent (i = 1:n)
  x = i * 0.01
  a(i) = x * log(x) / (x - 1.0)
end do

これによって、do loop 内の構造が順不同・同時・並列に実行できることが明示されます。

forall の注意点

ちなみに、forall 構文での複数行処理は、各行を完全に実行してから次の行に進みます。つまり

ia(1:n + 1) = -1
forall (i = 1:n)
  ia(i) = i
  ib(i) = i * ia(i + 1)
end forall

は

ia(1:n + 1) = -1
forall (i = 1:n) ia(i) = i
forall (i = 1:n) ib(i) = i * ia(i + 1)

と等価になります。したがって ib(i) = i * (i + 1) となり -i とはなりません。

まとめ

このように、FORTRAN77 の同じ do loop 構文の記述が、異なる複数の概念を表しています。これらを明示的に区別することで、プログラムの意味をよりはっきりできます。そしてこれはコンパイラによる最適化や並列化の助けとなります。

Fortran コンパイラは賢いので簡単な do loop などの並列性は簡単に見抜きますし、あえて明示的に区別しても対して利益はないかもしれませんが、将来的には必ず有用になると（正月気分で）想像します。

そもそも順番に実行すべき反復動作と、対応関係に基づいて複数の操作をまとめて実行することは別の抽象であるので、そのことを意識的・自覚的に認識して、別の構文で明示的に書き表すのが好ましいように思われます。

Modern Fortran には、そのための構文がいくつも用意されています。

2013-01-01

浮動小数フォーマット

Fortran

かつてはハードウェア毎に浮動小数フォーマットが異なっていて、汎用性のある数値計算プログラムを作るのが困難でした。Fortran90ではそのような問題を解決するために、浮動小数に関する情報を得るための命令がたくさん用意されました。もっともFortran90が普及する頃には、浮動小数もIEEE754にほぼ統一されてしまい、それらの命令が活用されることはあまりありませんでした。

ここでは、メモ帳代わりに各種フォーマットをメモしておきます。
$\pm f\times\beta^m$
f は仮数部、β進表現で
$f={c_1\over\beta}+{c_2\over\beta^2}+\dots+{c_n\over\beta^t}$
[tex:0

Computer	$\beta$	t	L	U	$\epsilon_M=\beta^{1-t}$
IEEE754(single)	2	23	-126	127	$1.19\times10^{-7}$
IEEE754(double)	2	52	-1022	1023	$2.22\times10^{-16}$
IEEE754(quadruple)	2	112	-16382	16383	$1.93\times10^{-34}$

Computer	$\beta$	t	L	U
NEC SX(float0);IEEE	2	52	-1022	1023
NEC SX(float1);IBM	2	56	-64	63
NEC SX(float2);CRAY	2	48	-16384	16383

(from SENAC Vol.30-4,1997)

Computer	$\beta$	t	L	U	$\epsilon_M=\beta^{1-t}$
Univac 1108	2	27	-128	127	$1.49\times10^{-8}$
Honeywell 6000	2	27	-128	127	$1.49\times10^{-8}$
PDP-11	2	24	-128	127	$1.19\times10^{-7}$
Control Data 6600	2	48	-976	1070	$7.11\times10^{-15}$
Cray-1	2	48	-16384	8191	$7.11\times10^{-15}$
Illiac-IV	2	48	-16384	16383	$7.11\times10^{-15}$
Stun(Russian)	3	18	?	?	$7.74\times10^{-9}$
Burroughs B5500	8	13	-51	77	$1.46\times10^{-11}$
Hewlett Packard HP-45	10	10	-98	100	$1.00\times10^{-9}$
Texas Instruments SR-5x	10	12	-98	100	$1.00\times10^{-11}$
IBM360 and 370	16	6	-64	63	$9.54\times10^{-7}$
IBM360 and 370	16	14	-64	63	$2.22\times10^{-16}$
Telefunken TR440	16	$9{1\over2}$	-127	127	$5.84\times10^{-11}$
Maniac II	65536	$2{11\over16}$	-7	7	$7.25\times10^{-9}$

(from G.E.Forsythe, M.A.Malcolm and C.B.Moler, Computer methods for mathematical computations, 1977)

Intel Fortran (IEEE754) $x=s\times b^e\times\sum_{k=1}^pf_k\times b^{-k}$

Computer	$\beta$	$p$	$e_{\rm min}$	$e_{\rm max}$
real(4)	2	24	-125	128
real(8)	2	53	-1021	1024
real(16)	2	113	-16381	16384

(from Intel Fortran User Reference Guide)

ソースコード

    program Console5
      implicit none
      integer, parameter :: ks = 4, kd = 8, kq = 16
      real(ks) :: s
      real(kd) :: d
      real(kq) :: q
      print *, 'Radix        b', radix(s)
      print *, 'Digits       t', digits(s)
      print *, 'Min Exponent L', minexponent(s)
      print *, 'Max Exponent U', maxexponent(s)
      print *, 'Epsilon       ', epsilon(s)
      print *
!
      print *, 'Radix        b', radix(d)
      print *, 'Digits       t', digits(d)
      print *, 'Min Exponent L', minexponent(d)
      print *, 'Max Exponent U', maxexponent(d)
      print *, 'Epsilon       ', epsilon(d)
      print *
!
      print *, 'Radix        b', radix(q)
      print *, 'Digits       t', digits(q)
      print *, 'Min Exponent L', minexponent(q)
      print *, 'Max Exponent U', maxexponent(q)
      print *, 'Epsilon       ', epsilon(q)
      stop
    end program Console5

実行結果

L,Uの定義が本文と違っていたけど仕方ないね。脳内で変換してみてね。

Radix b 2
Digits t 24
Min Exponent L -125
Max Exponent U 128
Epsilon 1.1920929E-07
Radix b 2
Digits t 53
Min Exponent L -1021
Max Exponent U 1024
Epsilon 2.220446049250313E-016
Radix b 2
Digits t 113
Min Exponent L -16381
Max Exponent U 16384
Epsilon 1.925929944387235853055977942584927E-0034
続行するには何かキーを押してください . . .

2012-12-18

A FORTRAN Coloring Book, by Roger E. Kaufman

Fortran

「FORTRAN ぬりえ」という本がかつてあったようです。実際は塗り絵ではなく、Dr.Suess（スース先生）風の絵入り本のようです。wikipedia http://en.wikipedia.org/wiki/Dr._Seuss
出版年が70年代末なので、FORTRAN66で書かれていると思われます。

http://www.seas.gwu.edu/~kaufman1/FortranColoringBook/ColoringBkCover.html

中身は結構真面目なようで、ラプラス方程式を解いていたりするようです。参考ページ。

2012-12-10

Tutorial: Introduction to Intel Visual Fortran

Fortran

http://software.intel.com/en-us/articles/tutorial-introduction-to-intel-visual-fortran

2012-11-23

Springerの電子書籍が安い！

Fortran

今週末に米国感謝祭の黒字の金曜セール、日本の初売りみたいなのがありますが、その週明けの月曜はリアル店舗で買いそびれた人がネット店舗に押し掛けるサイバー月曜とされています。元々はドイツの癖に Springer も参戦するようです。

Steve Lionel のTweetによると、Ian Chivers and Jane Sleightholme 著の Introduction to Programming with Fortran With Coverage of Fortran 90, 95, 2003, 2008 and 77 2nd Edition が定価 $90 のところ PDF版だけですが $15になるようです。まぁどっちかというとf90時代から書き足していって、拡張しすぎの温泉旅館みたいになってるしょーもない本なんですがｗ買うか。

http://www.apress.com/springer/computer-science/9780857292322

http://www.apress.com/customer-support/

それにしても、日本の新嘗祭を進駐軍のアホどもをだまくらかして、勤労感謝の日にカモフラージュさせて守った先人たちは偉いですね！

2012-11-21

四半世紀前の並列プログラムの問題点が今と同じ

Fortran

A Comparison of 12 Parallel Fortran Dialects http://dl.acm.org/citation.cfm?id=624571.624759 を斜め読みしたのですが、四半世紀前に並列プログラミングの問題点として挙げられたことが、今とさして変わらないことが分かりました。
共有メモリー型におけるキャッシュ・コヒーレンスの問題、自動並列化やコンパイラ・ディレクティブの記法、種々のメモリーロックの方法。
分散メモリー型におけるメモリー同期まわりの簡単さと、双方向メッセージ・パッシングの面倒さ。コンパイラによる自動並列化の困難さ。アセンブラとコアダンプ・デバッグを思わせるプログラマへの負荷。
また、キャッシュの効果でコア数に対してスピードアップが傾き１以上で増えるスーパーリニア効果についても書いてありました。

fortran66のブログ

fortran について書きます。

Fortran 歴史的遺物

DIM 関数 (DIMINISHING function)

CHAR ICHAR / ACHAR IACHAR EBCDIC と ASCII

CSHIFT,　EOSHIFT 関数　　並列化（HPF) 機能としての配列シフト関数

反復と射

反復と射

Fortran90での代入

条件節のある場合

where 構文

forall 構文

mask による分岐

do concurrent 構文

forall の注意点

まとめ

浮動小数フォーマット

ソースコード

実行結果

A FORTRAN Coloring Book, by Roger E. Kaufman

Tutorial: Introduction to Intel Visual Fortran

Springerの電子書籍が安い！

四半世紀前の並列プログラムの問題点が今と同じ

DIM 関数 (DIMINISHING function)

CHAR ICHAR / ACHAR IACHAR EBCDIC と ASCII

CSHIFT, EOSHIFT 関数 並列化（HPF) 機能としての配列シフト関数

反復と射

Fortran90での代入

条件節のある場合

where 構文

forall 構文

mask による分岐

do concurrent 構文

forall の注意点

まとめ

実行結果

CSHIFT,　EOSHIFT 関数　　並列化（HPF) 機能としての配列シフト関数