大晦日になって何でこんなこと書いてるんだと思うんだが、現状があまりに酷いので記事を書くことにした。
筆者の仕事というのは専らFEMによる数値解析で、あまり開発とかやっていなかったりする。どっちかというと計算用のモデルをくみ上げて計算機に放り込んで出てきた結果を纏めるといった程度であって、他の記事で書かれておられるような方々より高度な仕事をしているという訳ではない。
まぁ、実務であれ趣味であれ世の中広いからいろいろ見聞する必要だけはあるんだろうが、井の中の蛙にとってそれは難しい。で、その井の中の蛙からしてみても驚愕というか、唖然と言うか名状し難いものに出くわすことが稀によく起こる。
我が国に有限要素法が導入されて半世紀以上が立つ。初期はプログラムの自作から始まり、現在はソフトウエアを使いこなすことが求められるように変わっていった。
ここで取り上げるのは現在もしつこく続いているらしい、そのFEMプログラムの作成についてである。
平たく言うと50年も経てばプログラミングに求められるものは変わり果てていると言っても過言ではない。初期のころ、つまり自分が垣間見ることすらなかった70年代のコンピューターの能力は今よりずっと低く、その中で効率的な手法を考えざるを得なかった。その最たるものが自己書き換えプログラムというもので、プログラム実行中にプログラムの実行コードを書き換えて別の動作をさせるのである。アスキー版「ハッカー英和辞典」にはその例として「必ず抜ける無限ループ」というのがある。
この例を見てお前は何を言ってるんだという人は正しい。無限ループから抜けるわけないじゃないか。その通り。無限ループというのはあるところをずっと周回するものだから抜けるわけがなく、またどこかで実行位置がループの先頭に必ず戻る。それがずっと繰り返すから無限ループというのである。
このプログラムの作者はアセンブラどころか素の機械語で直接プログラミングをしていた。この時点で今なら正気を疑う人がいるだろうが、当時は当たり前だった。それどころかそうしないと「効率がいい」プログラムを書くことが出来なかった。ここで言う効率と言うのはプログラムサイズと速度の両方を言う。こんなことをするには当然自分がプログラムをしている機械について熟知していないと不可能である。最適化を煮詰めようとする場合などは今でもそうである。
当時のコンピューターの能力はせいぜい80年代の16ビット機程度と思われるので、メモリも多くなかったから、なるべくコンパクトに書く必要があった。速度も当然遅いので最高効率を叩き出すべくあらゆる手段を講じる必要があった。それを達成するための手段の一つが自己書き換えコードである。
で、この必ず抜ける無限ループはどうやったかというと、ループ回数に達するとメモリ上のNOPをジャンプ命令に書き換えるのである。今なら絶対にこんなことは不可能である。多分というか必ず一般保護例外を出してプログラムは終了する。
で、このプログラムをどうやって解読するんだという話になる。コメントなんてものは存在しない。というか入れられない。機械語だけで書かれたコードにコメントを入れることが出来ないからである。解読する方も書いたものと同程度の知識が必要である。そうしないとこのようなものを解読することは不可能だろう。
これが全てでないにしても70年代のプログラム開発というのは、そして少し遅れた80年代のパソコンでの開発というのはこのようなものだった。流石にパソコンで自己書き換えコードは実装されなかったかもしれないが、断言はできない。何しろ「真のプログラマ」はそのようなことを平然と行うからである。それは置くとしても、当時としては実行効率を考えるならアセンブラを使った方が遥かに早かったからである。
註1)アスキー版「ハッカー英和辞典」では「高級言語を軟弱とみなし、スクリーンエディタを軟弱者の道具と見做して自分のコードをデバッガで編集する手合いのことである。
この状況が変わったのが、パソコン上では90年代半ばである。MS-Windows95が登場した辺りから、コンピューターの能力は一気に増大した。まぁ、GUIなんてものを入れたらメモリをバカ食いするし、それ以前にグラフィック自体が資源を食うのでコンピューターの能力を上げないと使い物にならないという現実が横たわっているのだが。
と、ここまで前置きをして、旧Fortranでの開発話に移る。今ならエディタを使用してプログラムを書くのだろうが、当時は手で書いていた。ここでも「お前は何を言っとるんじゃ」と思うだろうが殆ど事実である。何しろ当時のFORTRAN入門書でこう書かれていたのだから事実というしかない。まずコーディングシートと呼ばれる専用の紙を用意する。ここにプログラムを書くのである。書き終わったら機械にかけると紙に書いた内容を読み取ってパンチカードに穿孔する。で、このカードを壱行ごとに機械から吐き出すのである。出力は当然この穴を開けた状態のパンチカードである。それをさらに機械に読み取らせることで、コンピューター本体のメモリにパンチカードの内容が「転送される」。そしてようやく本来の処理を始めるという訳である。
一見するとめんどくさい方法だが、こんなことをしていた理由など単純で、パンチカードの山が当時のコンピューターの外部記憶(今でいうハードディスク等である)に相当したからである。一行につき一枚のカードが使われているが書き誤りなどがあれば、そのカードを廃棄して差し替えることはできる。当時もハードディスクは存在したが、容量が低く常時ユーザーのデーターを置くなんてことは出来なかった(と思われる)。従ってこれらの高価な機材には必要な時にデーターを書き込み、使い終わったらデーターを消去するようになっていたはずである。と言うかそのようしないとデーターでHDDの容量が直ぐ一杯になるから、直ちに消せと言われまくった人もいるはずである。
これが遅くとも80年代の終わりまでの大学や企業におけるプログラムの開発状況である。もちろんFORTRANも例外ではない。
個人向けのコンピューターが出てきたのは70年代末期である。といってもCPUから始まり各素子を基盤にはんだ付けする必要があったので、電子工作の知識と実践が出来ないといけない代物だった。その後パーソナルコンピューターが販売されるに至る。今使っているパソコンの直接の起源がこれである。当時と今で異なるのは外部記憶などはなくあってもテープレコーダーを使うしか無い(ハードディスク？それって何？)事だったり、フロッピーディスクが出始めたりした時期だった。ディスプレイもブラウン管式のものしかなくテレビを使ったりしていた時代である。プログラムに使えるものといったBASICだけだった時代でもある。CやPascalが出始めたのは80年代半ばであり、実はFortranコンパイラも登場していた。嘘と思われるだろうが今のIntel FortranはMS-Fortranが起源である。
90年代になってPC用のHDDが出たが40Mbyteといったものが多かった。これでも当時としては大容量である。何事もだが当時のことを今の目では見てはいけないのである。
勿論このような外付けのソフトを動かすにはOSが必要である。8ビット機用には(私は触ったことがないが)、CP/Mというのがあった。16ビット機用には同じくCP/MやMS-DOSがあった。これらの言語製品はこの上で実行するようになっていて単独では動かなかった。今では不思議でも何でもないが、8ビット機時代のBasicが直接立ち上がって処理待ちをする時代から、明示的にOSを先に起動しそれからこれらのプログラムを実行する時代へと変わっていったのである。
プログラミングの方はどうだというと、BASICの時代には近代的プログラムの萌芽すらなかった。そりゃそうだろう。制御文はif文の他while文があったが後者が使われていたとは言い難い。どちらかというと割り当て型goto文が多用されていた。これには理由があって、外部機器を使うときの割り込みの為に使われていた。勿論計算型goto文も使われている。サブルーチンには引数なんてものはない。従って変数はこれ全て大域変数でプログラムの管理どころではない状態になる。構造体なんてのは勿論ない。
当時のgoto文の使われ方は単純なジャンプの他に割り当て型(イベントドリブン)や計算型goto文というのがあり、多重分岐にはそれらしか使えなかった。モノによってはgotoがgosubになっていたりするが同じことである。
これが当時の8ビット機における開発状況である。速度優先ならアセンブラがあったが、これも恐らく大同小異である。
それがOS上でのソフト開発が一変すると状況が変わってくる。開発には高級言語が不可欠であり、そのためにC言語が普及しだした。当時のPC上だとCの開発の主力を担っていたのはMicrosoft C Compilerである。OSベンダーの出す開発ソフト故である。これと互角の戦いを繰り広げていたのがTurbo PascalやTurbo Cを出していたBorlandだった。前者は高い最適化とベンダー提供の製品が強みだったが高価だった。対するBorlandは低価格と統合開発環境が売りだった。最も値段についてはMicrosoftだけの専売特許ではなく実は他も同じようなものである。
C言語には当時のプログラミングに必要なものが一式揃っていた。つまりこれ一つあれば現代的プログラムが組めるわけである。C言語はATTのベル研のハッカーたちがプログラマ向けに作ったのだから当然と言える。ポインタ関連が面倒くさい点を除けばだが。
こうしてC言語はプログラム言語の主流にのし上がった。これらのプログラミング言語が齎したのは高速な処理を行いうるプログラムの作成だけではない。データー構造・可読性といったソフトウエア工学の賜物である。
だがFortranはこの流れから完全に取り残された。旧規格最後のFORTRAN77の最初のコンパイラが出たのは80年代になってからである。この後の大改訂になるFortran90、最初のmodern fortranのコンパイラが出たのは93年である。ここに至ってC言語とタメを張れないことないレベルになった。
しかし猫も杓子もmodern fortranに移行とはならなかった。今問題とされているのはプログラムの可読性であり、開発効率の向上である。最早旧規格の風土病などに付き合っている暇は誰にもないのである。それもこの風土病を直そうともせずただ漫然としているだけである。要するにソフトウエア工学的見地からすれば彼らは「何もしなかった」のである。これが正しい歴史的解釈である。大学教授だのといった人種が彼らが習ったのはFORTRAN66である。80年代前半まで学生だった連中ならばこう断言しても間違いなかろう。存在しない言語処理系の教育など不可能である。
こんな連中が教鞭を取る立場になったらどうなるか。答えは明白である。要するに下手糞プログラマー(と言うのも烏滸がましい連中)が出来るだけである。更に言うと、今どきの大学は家が建つほどの授業料を取っているにもかかわらず初心者未満の連中しかできないのかという話になる。
更に問題なのは言語の機能すら熟知していない連中が量産されることである、これには根拠がある。いや、この根拠こそが問題なのだが、丸善の「Fortran90/95による有限要素法プログラミング」(藤井元岐阜大教授以下2名)のP79を見ればわかる。この部分には看過できない問題がある。要するに車輪の最発明である。
この部分の説明をすると、ファイル名に連番を入れることを考える。まぁ、hoge.001.dat、hoge.002.datという具合に最大でhoge.999.datという風にファイル名をつけるわけである。このようなファイル名をプログラム中で生成するにはどうするかである。Fortran77以降であればFORMAT文で(a,".",i3.3,".",a)を指定し、それぞれwrite文で書き込み先にファイル機番ではなく文字列を与え、引数にファイル名(この場合hoge)、番号、ファイル名を与えればよい。この部分はFORTRAN77ではこのようになる。
character(100) fname
100 format(a,".",i3.3,".",a)
write(fname,100) 'hoge',number,'dat'
書けばいいのは、たったのこれだけである。あと必要なのは桁数のチェックだけだろう。つまり書籍中にある10行以上がこれで消せることになる。それを何を考えたかだらだらと数値から文字への変換作業なるものを桁数分実行している。この方法だと桁数が多くなったときに書き直しが生じるしプログラムが無駄に長くなる。上記の方法なら3桁から4桁に増やす場合i3.3をi4.4に変更するだけである。そういえば内部ファイル機能を使う必要があるが、こんなことも知らないのが丸解りである。
modern fortranではこの部分はこうなる。
character(100):: fname
write(fname,’(a,".",i3.3,".",a)') 'hoge',number,'dat'
だがもっと重大な点がある。このような機能はANSI FORTRAN66には実はない。FORTRAN77で追加されたのである。つまり藤井以下本書の査読を行った全員がこの機能を知らなかったと言われても反論は出来ない。現実にこんな書き方しているのだからそうなる。何見てたんだというレベルではない。彼ら教授やら研究員やらはこの機能について知らなければならないはずである。使わないから知らないという言い訳も通用しない。研究員の方は兎に角、教授というのは教育者である。知らなくてどうするんだ？
この素晴らしいプログラム原文をここに書いてもいいが、いい加減指が腐るのでやめた。
更に危惧するのはどうやら日本計算工学会なるものが後援しているらしいことである。この本のカバー末尾に同会出版書籍の一覧が載せられている。このプログラムを見たはずの一人、生出佳もここの会員である。なんでも調布の方のメカニカルデザインの研究者とあるが、ここの研究員ですらこれでは絶望しかないのが解っているのだろうか。豊田中研の連中諸共解かってはおるまいな。それとも藤井が駄々でも捏ねまくったのだろうか。共著者と査読者全員の顔に泥どころか糞尿を塗りたくる行為にしかならないのだが。
これには後日談がある。勤め先でも同じようなものを見たのである。こちらは書式設定にi0を用いて数値を文字列化し、行頭に0を追加するという方法をとっているが、面倒は同じである。最大桁数が固定ならば上の方法をとったほうが余程スマートである。それは兎に角、これを書いた研究員が岐阜大出身とは思えないから、話を際限なく大きくするとこのような知っておいても損しない機能についてはどこも教えていないということになる。
ここまでくると、今まで何を教えていたんだと言いたいが、教育の目的として置くべきはセミプロ並みのプログラムを書けるようにするべきである。金をふんだくって出来上がるのが3流未満は許されないのである。この調子で(退役したが)京や富岳を使わせる気かといいたくなる。効率の悪い書き方をするのが目に見えているので、コンパイラは一杯仕事するしかなく、最新のアーキテクチャーに関係ないコードを生成するのだけは分かる。
この書籍のプログラムに横溢しているのは結果としての開発効率の軽視である。後の方に出てくる要素剛性行列の計算部分を見ればわかると思うが、何かの修行かといいたくなる内容である。同じ問題が日刊工業新聞社の「有限要素法のつくり方！」にもある。この本の内容も大概で京大の副学長様が前書きを書いているが騙されてはいけない。これらの本は一言で言って有害図書である。この書籍の内容はその偉い人の顔に泥を塗りたくる作業である。それともこいつら全員特殊な趣味でもあるんですかね？ああ、自分の顔じゃないからいいのか。そうなのか。
問題その2は初心者向けに書いた気になっているころだろう。それすら怪しいのだが。第一初心者を初心者のままにしてどうする。初心者は初心者から脱出できない内容である。この調子で教えていたらそりゃ世界で通用するソフトウエアの作成など只の妄想である。だって設計すら教えないんだもの。できるわけがない。
と思ったら、下には下がいた。いて欲しくなかったんですが。東大の中島研吾教授の熱伝導解析プログラムである。まず冒頭に燦然と輝くimplicit real*8が存在する。不適切な命名ばかりが並ぶ変数表を見ると暗然たる気持ちになる。何で節点総数が二つあるのか意味不明であるし、その名前もよりによってNとNPである。一文字の大域変数など悪夢の始まりそのものである。特にNのような変数はループで使われやすい。こんなもんを出した時点で、またデーターモジュール名が全く意味をなさないpfem_utilである時点でどうしようもない。
つまり変数名の付け方が不適切である。例えば最大反復回数にはITERが使われているが、MXITERにでもしたらよかっただろう。ただのITERにはそんな意味はない。このほかにITER_ACTUALとかという変数があってこれが只の反復回数である。ならばそれぞれ単純にMAX_ITER・ITERにすればよい。他ではこうしているのに、なんでこんな命名をしているのか理解できない。この世界では長いものに巻かれたり寄らば大樹の陰は大正義だったりする。この長いものや大樹は世間一般では標準というが、それに従っておけば問題ないという意味である。他にも節点セットを配列化しているが、構造体にするべきである。MPP版ではそうしているのに何故やらない。これも名前の配列、セット当たりの節点の個数をそれぞれ配列にしているが、なら肝心の節点リストはどうなっているかというと、全部一つの配列に押し込んでいる。こんな手間暇かけるくらいなら構造体にして番号・名前・節点数・節点リストをひとまとめにするのがよい。というか。そうしろやと怒鳴りつけたくなる。よくこんなものをYouTubeに出したと感心するが(中島以下関係者全員にとっては)どこも違和感がないのだろう。
また、読み込んだ変数をプログラム中で破壊している。これもコネクティビティであれば、読み込んだ節点番号と実際のデーターの位置が異なることは日常茶飯事なので内部インデックスに置き換える必要があるから仕方ないともいえるが、このプログラムでは反復計算における許容値を破壊している。混乱の元である。この変数はRESIUという名前だがこれも適切ではない。共益勾配法の許容値なんだから、CG_TOLとした方がよかっただろう。RESIUの元になったresiualという単語はどちらかというと残りの意味であって自分の知る限り許容値の意味ではない。それはtoleranceである。構造解析ではこのようになるが熱伝導解析では異なるので主か？それとも英語も残念なんですか？
慣用的に使われる単語を慣用的に使わないのに何か意味があるのだろうか。
因みにresidualには「説明のつかない」「除去できない」の意味もある。普通は「残りの」「残余の」という意味で使われている。いずれにしても許容値という意味で使うのは間違っているとしか言いようがない。
で、なんでこんなことしなきゃならんのだ
構造解析における要素剛性行列の計算に相当する要素熱容量行列の計算を行っているJACOBという副プログラムだが、パラメーターが31個も並んでいる。藤井もだが、抑々こんなにだらだらと並べたらサブルーチンを呼び出すのに時間がかかるだろうが。この連中は今の関数呼び出しのシーケンスも知らないと見える。
まず関数を呼び出すにあたってこれらのパラメーターはスタックに積まれる。Intelベースのアーキテクチャーだと最初の4個がレジスタ行になる。勿論レジスタ内のデーターを上書きするので、これらもスタックに退避する必要が生じるが、こうすることで、実行時の負担を低くしているのである。中島の場合、31個のうち24個が座標データーなので一つの配列にまとめることが出来る。これで8個にまで減った。藤井の場合引数は10個あるが8個がコネクティビティだから、これも配列にすれば全部スタックではなく、レジスタに格納されることになり、ただのジャンプ命令同然となる。更に最初の2個はヤング率とポアソン比だから、これも初めからDマトリックスにしておけばパラメータは僅か2個にまで減少するし、例題中にあるように要素ごとにDマトリックスを計算する手間を省ける。要するに設計の不備である。
中島の場合最初の一つがヤコビアンの行列式、次の3個が∂N/∂X・∂N/∂Y・∂N/∂Zでさらに∂N/∂r1・∂N/∂r2・∂N/∂r3が続く。これらを纏めると最終的に引数は4個となり全部レジスタに格納できる。
どうしてこんな書き方をするのかというと、こちらが早いと考えているんだろう。実際70年代のコンピューターであれば早い方法だった。コンピューターにはスタックなんてものがないことがザラである、嘘と思うだろうがIBMのマシンにはそんなものが無かったことが知られる。次にFortranの副プログラムの引数は変数のアドレスを渡している。つまり値をコピーするのではなく、変数の在処を副プログラムに教えているわけだ。だから呼び出した先で引数の値をうっかり変更しようものなら、意味不明なエラーが生じることになる。それは置くとして、では、これらの副プログラムを一度しか呼ばない場合はどうなるかである。そのような場合、引数のアドレスを呼び出し先のアドレスに置き換えればいい。そうすれば処理はジャンプ分を挟むにしてもシームレスなものとなるし、実行速度も上がることが期待できる。
そして藤井諸共X1からX8のようにだらだらと変数を個別に宣言しているのも同じ理由である。スタックなぞ無いから変数はこれ全て大域領域に置かれることになる。通常コンパイラは出てきた順に変数をメモリ上に配置するだろうから、それを期待した最適化も行いうる。そんなところだろう。であればこれらの変数を用いた内積の計算もいちいち手で展開しているのも説明がつく。ループを使えばいいのにと思うし、今のFortranならdot_productがある。これを使えば、ループを使わずに済むし、コンパイラによってはもっといいコードを生成する可能性だってある。
しかし、ループを使えばその分処理時間がかかるから手で展開したほうが早い、そういうことなのだろう。何しろ昔はコンピューターそのものの処理が遅かったのである。処理速度を稼ぐにはループを使わなければいいということになる。
しかし、今は違う。コンピューターの処理速度は大幅に上がった。処理速度は機械によるが桁違いなはずである。それよりもソフトウエアの大規模化や複雑化に伴て可読性やコード開発の効率化が重視されるようになった。それらを無視しているのは計算力学だけではないのか。
そして、なぜ一体であるべき変数をばらすと駄目なのかを説明すると、これらはスタックに積まれることになる。そうであれば、これらは全部スタックポインタ相対でアクセスされる。期待が持てるとしたら変数によっては連続した領域に配置されるだろうから、それを基にした最適化を行ってくれるという希望である。だが、そういう風にはならない。配列であればどこにあっても連続した領域に配置されるという仮定を置くことが出来るから、ベクトル化といった最適化を容易に行いうる。だが個別の変数とした場合、コンパイラがその仮定はしないとしても文句は言えない。しなかったらどうなるかというと昔懐かしい書き方をしたコードが吐き出されることになる。嘘だと思うんなら一度アセンブラ出力をしてみることを薦める。本当にそうなるのだ。
更に今のCPUのレジスタは64ビットとは限らない。モノによってはこの8倍の512ビットなんてのもある。倍精度実数なら8個のデーターを格納できるサイズである。この二人の書いたプログラムで使われている要素は6面体一次要素といい、節点数は8個である。つまり各データーが8個ある。SIMD計算をフルに使えた場合、ヤコビアンの計算では、格納6回、並列乗算9回、水平加算9回、ストア9回で済むが、昔懐かしい方法では一つの要素の計算で格納16回、乗算8回、加算7回、ストア1回である。格納については工夫すれば48回に収まるがレジスタも48個必要である。因みにIntel CPUのレジスタは16か32個の筈だから全然不足である。さらに工夫してもストア32回だが、これだけでレジスタを使い切ることになる。
その中島が教えているのがHPCである。なんでもHigh Performance Computingの略らしいが、このプログラムを見る限りそんなものはどこにも見えない。
ついでに下らない部分があったので指摘しておく。このプログラムでは全体剛性の計算に先にも書いたように共益勾配法が使われており、全体剛性行列を作る際に0でない部分だけを記憶するようにしている。この中で各行の処理を行う際に、必要な部分だけ記憶するという作業がある。この位置を記憶するのに際して、何を思ったか大きい順に並べ直し、格納時には小さい順に並べ直す(要するに逆順に代入)しているのである。こんなことをするくらいなら初めから専用のルーチンを作成し、データーに各行の情報を渡せば一発である。どこがハイパフォーマンス何だろうか。
このほかにも56回のおんなじサブルーチンの呼び出しがあったりする。これも渡される変数名を見ると要素内部の節点番号を配列ではなくばらしたものを渡している。同じものを渡す場合は処理しないようにしているが、そんなものはサブルーチン内で直ちに帰るようにすればそんな面倒なことはしなくて済む。引数も配列にしておけば二重ループを書くだけで済む。突っ込みどころ満載である。
で、中島研吾のプロフィールを見たら日本計算工学会会員とあった。藤井の本といい、こいつの本といい、ここは下手糞なプログラムを曝す趣味でもあるんだろうか。
要するに昔々学生の頃習った方法に固執しているだけの輩が教授になると、こうなるということである。そこには他言語で取り入れられている流儀やいい点を取るということはしない。その結果古臭い方法が目立つものとなり、この結果がFortran撲滅論に繋がるわけである。
Fortran撲滅論を潰すいい方法は一つしかない。それは現在使われている主要言語と遜色のない書法を使うことであり、プログラムやデーターの設計といったものに気を配ることである。この部分が今@本的に欠如しているのが問題なのだが。それにもかかわらず最新のソフトウエア作成論に従ったものを書いて、今のFortranが「古くて新しい」ものであることを証明するしかないだろう。
