fortran66のブログ

fortran について書きます。

メモ帳:週末小ネタ集

"you cannot do all of these at the same time."

パッと見、良さげな約束をしてくれますが、同時には成り立たないのは基本w

The problem with Common Lisp

Now I think that one of the main reasons for this is that while you can write scientific code in CL that will be (1) fast, (2) portable, and (3) convenient, you cannot do all of these at the same time.

Switching from Common Lisp to Julia
Switching from Common Lisp to Julia

COMMON LISP から Julia というところも・・・味わいが


ネスカフェ プレジデントCM 1983年 60秒 コーダー伯爵

1984CM NESCAFE プレジデント プリンス・デ・ラ・トレムイ
味わいは深く濃く


ネスカフェ プレジデントCM 1980年 30秒 ロドルフ殿下
ロドルフ殿下は全てに最高のものを求めるw


Retire Fortran?

A debate: Retire FORTRAN?: Yes (1984)

James R. McGraw
Lawrence Livermore National Laboratory, Livermore, California
http://physicstoday.scitation.org/doi/abs/10.1063/1.2916244

A debate: Retire FORTRAN?: No (1984)

David J. Kuck
University of Illinois, Urbana–Champaign
Michael Wolfe
Kuck & Associates, Inc., Champaign, Illinois
http://physicstoday.scitation.org/doi/abs/10.1063/1.2973891

Retire Fortran? A debate rekindled (1991)

Retire Fortran? A debate rekindled

SisalというLLNL謹製の聞いたことがないPASCAL系の関数型言語FORTRANの比較をしてその並列化での優位を説いています。
Sisal Language Tutorial Table of Contents

sourceforgeに廃墟が・・・
sourceforge.net

結果はKuckのOpenMPFortranが残り、LLNLの並列言語は忘れ去られました。

映画「ドリーム」

FORTRAN と IBM7090 が出てくる映画ということで見てきました。卓上計算機で解いたオイラー法に負けます。

三人の黒人女性がそれぞれ差別と偏見を乗り越えて NASA で活躍するという実話に基づく話です。教育テレビの道徳の時間に流れそうな全く興味のもてない筋書きですが我慢です。アメリカ映画らしく、途中で説教演説も聞かされますが我慢です。

話の細部はボロボロで少女マンガレベルですが、娯楽度とテンポはプリプリちぃちゃんの方が優ります。


映画『ドリーム』予告 IBM編



Using OpenMP 新刊

Using OpenMP の新刊が出るようです。
元々はベンダー毎にばらばらだった並列化 FORTRAN 用のコンパイラ指示行を統一する規格でしたが、20年の時を経て謎の進化を遂げています。

Using OpenMP -- The Next Step: Affinity, Accelerators, Tasking, and SIMD (Scientific and Engineering Computation)

Using OpenMP -- The Next Step: Affinity, Accelerators, Tasking, and SIMD (Scientific and Engineering Computation)

Using OpenMP: Portable Shared Memory Parallel Programming (Scientific and Engineering Computation)

Using OpenMP: Portable Shared Memory Parallel Programming (Scientific and Engineering Computation)

f:id:fortran66:20171015004714j:plain

1982 年の FORTRAN 誕生 25周年記念の集まりの記事

Annals of the History of Computing
Issue 1 • Jan.-March 1984

FORTRAN 誕生 25周年記念特集
IEEE Xplore: Annals of the History of Computing - ( Volume 6 Issue 1 )

このビデオが作られた会のようです。


The Beginnings of FORTRAN (Complete)

全部見たわけではありませんが、色々面白い逸話が紹介されています。

バッカスFORTRAN コンパイラは予定をズルズル遅れて、完成前は IBM 内部でも「できっこない」的な雰囲気だったのに、ひとたび出来上がって好評だと、手のひら返しとか。

FORTRAN コンパイラチームのメンバーの思い出で、本人はいつも遅刻ギリギリに来ていたが、他のメンバーは朝相当早くから揃っていて仕事がはかどっていた。不思議に思っていたら、ビルの向かいのアパートに住んでいる寝間着を着ないで寝ている若い女を見るためだった。その後、仕事場が引っ越してデパートの更衣室が見えるようになって、日がな窓の外を見る奴らばかりで、仕事がはかどらなくなった。最終的に、眺めの悪いところに引っ越すことになって、なんとかプロジェクトが完成したとか。

IBMコンサルタントしていたフォン・ノイマンは、FORTRAN 開発提案に反対しなかったが無関心だった。当時すでに病気が進行していたからかもしれないとか。

当時の人々が、コンピュータ言語というものに懐疑的であった具体的なエピソードなども出ています。ベテランはアセンブラ固執しがちだったが、新参が続々 FORTRAN を使って普及とか。

f:id:fortran66:20171010222618p:plain


McCracken が FORTRAN 本を出版する頃の思い出もあります。

Guide to Fortran Programming

Guide to Fortran Programming

FORTRAN 〓@77E1 入門 (1969年)

FORTRAN 〓@77E1 入門 (1969年)

メモ帳:Modern Fortran の構造について

Fortran の構造

Fortran の分類

Fortran は、約 60 年間に渡って変化・発展してきたプログラム言語で、言語内容からみて FORTRAN I から FORTRAN77 までの 古典 FORTRAN (Classical FORTRAN) と、Fortran90 以降の現代 Fortran (Modern Fortran) に大きく分けられる。

古典 FORTRAN から受け継がれてきたもの

Fortran 言語には、最初期の FORTRAN I/II の頃から現在に至るまで、受け継がれてきた根幹部分がある。その中で

  1. 代入
  2. 配列
  3. SUBROUTINE / FUNCTION の別

が主要と思われる。

型とは二進列を整数や実数などの代数系で演算するためのパラメータ、代入は不可逆な過程で履歴を伴う作用、配列は記憶領域上で等質・連続な抽象データ構造、SUBROUTINE および FUNCTION は処理の抽象化の仕組みで戻り値・副作用の有無の別による。

現代 Fortran で加わったもの

Fortran の現代化改修には、二つの主要部分がある。

  1. 構造化プログラミング [ i) 抽象データ構造 (オブジェクト指向を含む) ii) 制御構造 ]
  2. 並列処理 [ i) 全配列操作 ii) coarray ]

構造化プログラミングは、N. Wirth 流のミニマリズムに則っており、オブジェクト指向関連の用語も Wirth に従っており、継承は拡張 (extention)、メソッドは派生型束縛手続き[?](type bound procedure) と呼ばれている。また例外構文・ジェネリックプログラミングも Wirth に従って、冗長機能として排されている。これは同じく Wirth に則った Google の GO 言語でも同様である。

現代 Fortran には、これらの他に

3. 浮動小数点数関係 [ i) 一般のフォーマット用 ii) IEEE754 用 ]
4. C 言語との相互運用 [ i) 古典 FORTRAN 水準 ii) 現代 Fortran 水準]

等の機能も付け加えられているが、付随的な部分とみなしてよいと思う。

構造化プログラミング

利用者定義派生型、モジュール、派生型束縛手続き[?]
ALGOL, PASCAL, MODULE-2, Oberon-2 (block, type, module, OO)
SUBROUTINE / FUNCTION, recursion
array

抽象データ構造

制御構造

(つづく)


並列化

対象ハードウェアの抽象化・アレイプロセッサ・パイプライン処理・メモリー共有型・分散型
implicit な(自動)並列化、explicit な並列化、コメント指示行、並列用構文
array

全配列操作

coarray

(つづく)


Modern Fortran Explained (Numerical Mathematics and Scientific Computation)

Modern Fortran Explained (Numerical Mathematics and Scientific Computation)

Classical Fortran: Programming for Engineering and Scientific Applications, Second Edition

Classical Fortran: Programming for Engineering and Scientific Applications, Second Edition

メモ帳: N. Wirth の Modula-2 の function

N. Wirth の Programming in Modula-2 1988年版はフォントが小さいので、function の章が 2 ページに収まって、Springer の preview で 全部読めるw

Programming in Modula-2 | Niklaus Wirth | Springer

  • 関数の副作用について

f:id:fortran66:20171005232214p:plain

  • まとめ

f:id:fortran66:20171005232326p:plain
3.副作用について
4.命名法について

Programming in Modula-2

Programming in Modula-2

function の side-effect について :追記あり

Fortran では、function の副作用は嫌われていて、文法的には許されていますが、道徳的に禁止になっています。Fortran95 以降では pure 指定子が出来たので、常に pure というのが行儀のよさだと思います。副作用のある場合は、サブルーチンを使います。

この点については、FORTRAN II で SUBROUTINE と FUNCTION が導入された時のマニュアルにも引数の値の変更はしないように書かれてあります。(厳密に言えば、昔の FORTRAN では変数は全て静的に確保されているので、静的変数に前回の履歴が残るという形で意図せずとも副作用がありますがw)

f:id:fortran66:20171004040611p:plain
FORTRAN II for the IBM 704 Data Processing System : Reference Manual | 102653989 | Computer History Museum
(p. 28)
fortran66.hatenablog.com

また、Fortran95 の impure 指定子に類する ABNORMAL FUNCTION 指定子が、1960年ごろの UNIVAC 1107 FORTRAN IV にあったそうです。
fortran66.hatenablog.com

Fortran で function の副作用が嫌われている理由は、最適化のためだと思います。副作用がなければ、関数は sin(x) + sin(x) が 2*sin(x) とまとめられるように、f(x) + f(x) を 2 * f(x) と出来ます。また FORTRAN では、式の項の評価順序が決まっておらず、コンパイラが自由に並べ替えて最適化できるので、副作用のある関数が項として複数現れる式では、結果が不定になってしまいます。

なお Fortran の組み込み関数や組み込みサブルーチンでもこの原則は守られていて、他言語では関数になっている乱数や時間もサブルーチンとして与えられて、引数で値を得るようになっています。

Wirth の場合

一方 N. Wirth の著述をみると PROGRAMMING IN MODULA-2 では、13章(p.52)に副作用のある関数の例が出てきて、引数の変更を含めて、特に禁止するようなことを書いていません。

Programming in Modula-2 (Monographs in Computer Science)

Programming in Modula-2 (Monographs in Computer Science)

しかし、のちの Programming in Oberon: Steps Beyond Pascal and Modula では 6.6.1 Side-effect の節で、副作用を好ましくもなく、良いプログラミング・スタイルでもないとしています。しかし、例として疑似乱数を挙げて、限定的使用も認めています。しかしその場合でも、許されるのは静的な変数を通じての依存性で、引数への副作用は戒めています。

Programming in Oberon: Steps Beyond Pascal and Modula (ACM Press)

Programming in Oberon: Steps Beyond Pascal and Modula (ACM Press)

f:id:fortran66:20171004032527p:plain

メモ帳:Wirth は上記 MODULA-2 本で、一般の関数名は名詞、論理型を返す関数名は形容詞、サブルーチン名は動詞にすることを勧めています。(上記 MODULA-2 本は、ほぼ同内容のものを Oberon 版にも書き直しており(上記本とは別のもの)、多分元々は PASCAL 版もしくは ALGOL 版で、内容は共通しているものと推測します。)

追記:
PASCAL report 初版では、グローバル変数も引数も代入するのはまかりならぬ(PURE であれ)となっていますが、改訂版では、その記述がごっそり削られています。 

f:id:fortran66:20171004213925p:plain
N. Wirth, The Programming Language Pascal (Revised Report) Acta INformatica, 1 (1971), 35-63.
http://oberoncore.ru/_media/library/wirth_the_programming_language_pascal.pdf

f:id:fortran66:20171004213641p:plain
N. Wirth, The Programming Language Pascal (Revised Report) (1972)
http://www.eah-jena.de/~kleine/history/languages/Wirth-PascalRevisedReport.pdf

Kernighan の場合

元々 C 言語は サブルーチンと関数の区別がなく、かつ関数の返り値をエラーコードとして利用して、引数の方の副作用で結果を返すことを常用しているので、関数の純粋性の概念は薄いと思います。

RATFOR を使った Software Tools でも、しょっぱなから getc(c) 関数を定義して、返り値のほかに引数にも同じ値を入れて返しています。Kernighan は、これは文法に違反していないと本文でわざわざ断っていますが、FORTRAN 的には行儀がよくない作法になります。

なぜこのようにしたかといえば、DO WHILE () 構造のためです。DO WHILE 構造は、Fortran90 から落選しそうになるくらい Fortran では微妙な構文です。Fortran では代入をはじめ実行文は値を持たないので、DO WHILE 構文の条件式には代入文や実行文を書けせん。したがって、DO WHILE ( (c=getc()) /= EOF) のような、1文字読んで変数に代入してその値が EOF かどうか調べる、といった書き方が出来ません。このようなスタイルは C 言語系のプログラムではよく見ますが、Fortran では使えません。DO WHILE 構文の使い勝手があまり良くないのです。

Software Tools の巻頭一番で、副作用付きの関数 getc(c) を導入して DO WHILE (getc(c) /= EOF) と getc 関数の返り値で EOF のチェックをするとともに、関数引数の副作用で変数 c に返り値と同じ値を代入するしかなくなっています。

これは、DO WHILE 構文をすてて、素直に DO..END DO の無限ループを EXIT 命令で抜ける形にすれば解決できますが、脱出が LOOP の中途から起こることになります。REPEAT..UNTIL 型の構文は滅多に使われることがないので、Fortran では採用されていませんので、個人的には更に一歩進んで DO WHILE 構文も捨てた方が良かったろうと思います。

Software Tools

Software Tools

ソフトウェア作法

ソフトウェア作法


  

行列積を Cannon の algorithm で

昔、コネクションマシンの本を読んだときに、各並列プロセッサは 1bit 演算で、浮動小数点数も 1bit 演算を繰り返して計算すると知って、これは駄目ですわと思ったのですが、グスタフソンの動的可変長浮動小数点数フォーマットのことを思うと、任意長に出来る 1bit 演算もありかなと思ってパラパラと斜め読みしなおしてみました。


日本語版には訳者がまとめた CM-2 の資料がありますが、その中に Fortran8x に基づく CM-FORTRAN のサンプルプログラムがあったので、今の Fortran で動かしてみました。CM-FORTRAN はおおむね Fortran-95 水準なので、ほぼそのまま動きます。(1bit 演算とは関係ありませんw)

コネクションマシン―65,536台のプロセッサから構成される超並列コンピュータ

コネクションマシン―65,536台のプロセッサから構成される超並列コンピュータ

サンプルプログラムは、サブルーチンで、Cannon のアルゴリズムによる行列積の並列計算です。元々二次元アレイプロセッサ用のものなので、あまり意味はありません。CSHIFT/EOSHIFT 関数が 並列計算用組み込み関数に分類されるのは、二次元トーラス上のアレイプロセッサに配列要素を割り当てる前提があったためです。


Cannon のアルゴリズムの説明動画を貼っておきます。

Cannon's algorithm for matrix multiplication

実行結果

乱数で生成した行列 A,B の積 C=A*B を、組み込み関数の MATMUL と Cannon 法の二通りで求めて、各行列要素の差の絶対値の最大値を求めています。また、各ステップごとの経過時間を秒単位で書き出しています。

Intel Fortran で自動並列化をかけましたが、たいして効果はありません、MATMUL は MKL のルーチンへの自動置換オプションを選んだので並列化が効いています。

行列要素ではなく、ブロック行列に適用しなければいけなかろうかと思います。

  80610.66     :start
 2.1000000E-02 :matmul
  3.666000     :Cannon
 0.0000000E+00 :maxval
max difference  5.6457520E-04
続行するには何かキーを押してください . . .

ソース

サブルーチンは、原文そのままを書こうと思ったのですが、コメント文と CSHIFT の分を直してあります。CSHIFT は引数の順番の定義が一部変わっていました。

    module m_mod
    contains
      SUBROUTINE CANNON(M1, M2, R, N)
      INTEGER N
      REAL M1(N,N), M2(N,N), R(N,N)
!C    Skew argument arrays.
      M1 = CSHIFT(M1, (/0:N-1/), 2)
      M2 = CSHIFT(M2, (/0:N-1/), 1)
!C    Multiply matrices
      R = 0.0
      DO I = 1,N
           R = R + M1*M2
           M1 = CSHIFT(M1,1,2)
           M2 = CSHIFT(M2,1,1)
      END DO
!C    Restore argument arrays.
      M1 = CSHIFT(M1, -(/0:N-1/), 2)
      M2 = CSHIFT(M2, -(/0:N-1/), 1)
      RETURN
      END 
    end module m_mod

    program CM
      use m_mod
      implicit none
      integer, parameter :: n = 1000 
      real :: a(n, n), b(n, n), c(n, n), d(n, n)
      call RANDOM_NUMBER(a)
      call RANDOM_NUMBER(b)
      call stamp('start')
      d = matmul(a, b)
      call stamp('matmul')
      call CANNON(a, b, c, n)
      call stamp('Cannon')
      call stamp('maxval')
      print *, 'max difference', maxval(abs(c - d))
    contains
      subroutine stamp(text)
        character(*), intent(in) :: text
        integer :: ic0 = 0, ic, irate
        call SYSTEM_CLOCK(ic, irate)
        print *, (ic - ic0) / real(irate), ':', text
        ic0 = ic
      end subroutine stamp  
    end program CM