fortran66のブログ

fortran について書きます。

EISPACK についての原作者講義

J. H. Wilkinson 理論、C. Moler FORTRAN 実装

Nick Higham のブログで、EISPACK についての原作者等による講義のビデオが紹介されていました。発掘の経緯についても書かれています。 The Argonne Tapesnickhigham.wordpress.com

3月に公開されて、現在まで 100 view 程度で草。

しかし、こういうリソースがあるのは白人様の文明の底力で恐れ入ります。ただで視聴乞食させていただきお靴を舐めたい所存であります。

f:id:fortran66:20190702021329p:plain
あいまいみー4巻

短いし、FORTRAN の話なので、とりあえず Moler の分だけ見てみました。

Moler の講義では、実対称行列を Hausholder 法で三重対角化するサブルーチン tred2 と、それをさらに QL 法で完全に対角化する tql2 を基本例としてソースコードを OHP で写しながら解説しています。 tql の途中で終わっているのが残念です。tred + tql は、かつて一世を風靡した numerical recipes にも、後継改良 77 版がそのまま収録されています。

tred で扱うのは実対称行列なので下三角に変換行列を貯めてゆくことでメモリーを節約していますが、Moler も巧妙だと称賛しています。元々の ALGOL 版からそうなっていたようです。

www.youtube.com Youtube 自動字幕では ALGOL が alcohol になっていて草。

今も開発の続く LAPACK は、線形方程式ライブラリの LINPACK と固有値問題ライブラリの EISPACK を合体させ FORTRAN77 に書き換えたものです。EISPACK は元々はドイツ人(C. Reinsch ?)が ALGOL でコツコツ集めていたのを、イギリス人の Wilkinson が加わって集大成し、それをさらにアメリカ人が FORTRAN66 に移植して成ったものです。

ALGOL はコンピュータ言語としてのほか、アルゴリズムを共有するための表記法としての目的を持っていました。

Wilkinson は誤差解析で有名ですが、固有値問題の本も書いています。その Wilkinson による EISPACK ライブラリの理論的背景講義の後、FORTRAN 移植を行った Moler の FORTRAN 実装に関する講義がなされています。Moler は MATLAB の開発者として有名ですが、その他にも G. E. Forsythe, M. A. Malcolm と数値計算レシピ本の嚆矢のような本も出しています。

Rounding Errors in Algebraic Processes

Rounding Errors in Algebraic Processes

基本的演算における丸め誤差解析 (1974年)

基本的演算における丸め誤差解析 (1974年)

The Algebraic Eigenvalue Problem (MONOGRAPHS ON NUMERICAL ANALYSIS)

The Algebraic Eigenvalue Problem (MONOGRAPHS ON NUMERICAL ANALYSIS)

Computer Solution of Linear Algebraic Systems (Automatic Computation)

Computer Solution of Linear Algebraic Systems (Automatic Computation)

Computer Methods for Mathematical Computations

Computer Methods for Mathematical Computations

Numerical Recipes in Fortran 77: The Art of Scientific Computing

Numerical Recipes in Fortran 77: The Art of Scientific Computing

講義では Machine Epsilon や二次方程式の桁落ち予防変形についても触れています。 また Moler は、元の ALGOL プログラムで、同一名の変数をプログラム中で全く異なる目的で何度も使っていることにちょっとだけ苦言を呈しています。

講義冒頭で触れられている J. H. Wilkinson and C. Reinsch, Handbook for Automatic Computation の Moler 本人による紹介記事が MATLAB のサイトで読めます。 blogs.mathworks.com

Handbook for Automatic Computation: Volume II: Linear Algebra (Grundlehren der mathematischen Wissenschaften)

Handbook for Automatic Computation: Volume II: Linear Algebra (Grundlehren der mathematischen Wissenschaften)

  • 作者: John H. Wilkinson,C. Reinsch,Friedrich L. Bauer,Alston S. Householder
  • 出版社/メーカー: Springer
  • 発売日: 1971/01/01
  • メディア: ハードカバー
  • この商品を含むブログを見る

tred & tql

昔、これらのルーチンにはお世話になりました。あの頃こういうリソースがあれば良かったろうにと思います。

Hausholder 法での三重対角化 tred のアルゴリズムは有限回のステップで終わりますが、tql の三重対角行列の対角化は反復法なので適当な閾値になるま繰り返します。5次方程式以上の解の公式がない以上、5次以上の行列の固有値問題にも有限ステップの解法が無いであろうからもっともな気がします。

問題によっては、最初の三重対角化に Lanczos 法を使って、非対角項が十分小さくなった適当なところで打ち切ると、小さい行列の対角化になって(tql は O(N3) のアルゴリズムなので)計算量を大幅に減らせます。 

f:id:fortran66:20190702024123p:plain

あいまいみー コミック 1-8巻 セット

あいまいみー コミック 1-8巻 セット

あいまいみー 9 (バンブーコミックス WINセレクション)

あいまいみー 9 (バンブーコミックス WINセレクション)

【寝言】最近の面白かった番組

NHK ビバ合唱

この土曜日の朝にやっていた、フィンランドのラヤトン合唱団の巻は結構良かったです。梅雨空に合う。

www4.nhk.or.jp

ラヤトン 無限の森へ フィンランド・アカペラの響き 音楽CDアルバム全15曲と絵本

ラヤトン 無限の森へ フィンランド・アカペラの響き 音楽CDアルバム全15曲と絵本

北欧の森の物語

北欧の森の物語

映画「君の名は。

テレビで放送していたので、ついながらみしてしまいました。これ以前の新海誠の作品は、どれも耐えられなくて5分と見てられなかったのですが、これは楽しめました。

夢の中で大事な出会いをしていながら、目が覚めた途端あっという間に忘れて薄れてしまい、後に焦燥感だけが残るという感じが良かったです。何を探しているのかも分からないまま探す感は不毛。

そういえばこの前見た夢の中で、高層ビルの見晴らしの良い上層階の、数フロアをぶち抜いたおしゃれ空間のスタバの様なカフェの中を、右手にカップラーメン、左手にお~いお茶を持って、空いてる席を探してひたすらさまよい歩きました。このままでは麺が伸びてしまうという、やるせない焦燥感をよく覚えています。

伊藤園 おーいお茶 緑茶 525ml×24本

伊藤園 おーいお茶 緑茶 525ml×24本

【寝言】最近のオモシロ記事

クリントン

民主党  哀れなり、クリントン家は蚊帳の外 本誌紹介 | ニューズウィーク日本版 オフィシャルサイト

ヒラリーが次期大統領選を完全に諦めたのは、先ごろモニカ・ルインスキー事件に関してビル・クリントンをかばってルインスキーを非難するという政治センスの無さを披露し、プロの支持者層からも見放された事情があるようです。パワハラ+セクハラ案件の典型例で、今の時代ではシビアさが上がっているのですが、貧困調査かなw

ビル・クリントンは、葉巻をモニカ・ルインスキーの局部に突っ込んだ後、口にくわえて 'It tastes good.' とのたまいました。当時の私は葉巻がしけって火がつきにくかろうと思ったのですが、よくよく葉巻の対称性を考えると、咥える方を突っ込んだ可能性もあると思うのであります。

コンサルw

www.foreignaffairsj.co.jp

土人国の生産単位どもは、相変わらず白人のご主人様をありがたがっているのかと思ったら、そういうことも無くみんな馬鹿にしているようです。春秋戦国の頃の鶏鳴狗盗の類か。

【メモ帳】flang Video

EuroLLVM Developers’ Meeting

視聴者の少なさw 内容の乏しさw  f18の話も少し。


2019 EuroLLVM Developers’ Meeting: S. Scalpone “Flang Update”


2018 EuroLLVM Developers’ Meeting: S. Scalpone “Flang -- Project Update”


2018 LLVM Developers’ Meeting: S. Scalpone “Flang Update”

【新刊】Fulton ヤング・タブロー

フルトン「ヤング・タブロー」

本屋をのぞいたら Fulton の Young Tableaux の訳本が新刊で出ていました。

表題は別に普通に「ヤング盤」でいいんじゃないかと思いましたが、訳者前書きに思う所が書いてありました。そのうちゆっくり立ち読みしたいです。

LMSST: 35 Young Tableaux (London Mathematical Society Student Texts)

LMSST: 35 Young Tableaux (London Mathematical Society Student Texts)

ヤング・タブロー: 表現論と幾何への応用

ヤング・タブロー: 表現論と幾何への応用

【メモ帳】U. Fano と G. Racah はいとこ同士

最近になって、今頃知ったことw (&昔のソース付け)

Giulio Racah, spectroscopy and group theoretical methods in physics

  •  U. Fano と G. Racah はいとこ同士

裕福なイタリア系ユダヤ人。Racah は E. Fermi と同期?

Following a visit to the National Bureau of Standards in 1954, Racah developed a great interest in the application of computers to theoretical spectroscopic analysis. During two years he developed additional programs (see ref. 57) to diagonalize matrices, comparing the resulting diagonalizations with experimental data and to compute line strengths.

Weizmann 研究所に建設せられた WEIZAC Computer は 一語 40bit の二進コンピュータで、スピードは IBM704 の約半分。対角化アルゴリズムは Jacoby-von Neumann method らしいが、多分いまでいうヤコビ法。

G. Racah. The use of the Weizac in theoretical spectroscopy, Bulletin of the Research Council of Israel 8F 1 (1959).

ワイズマンのテーマ

ワイズマンのテーマ

 J. Backus がアセンブラやってらんないと FORTRAN を提案し IBM の公式になったのが 1954 年頃。当時 IBM の顧問をしていて、末期がんだった J. von Neumann は興味示さず。

FORTRAN Anecdotes, January-March 1984, pp. 59-64, vol. 6. CSDL | IEEE Computer Society

https://archive.org/details/fortran-anecdotes/page/n3archive.org

It is impossible to probe the beginnings of an event of the magnitude of FORTRAN without the name of John von Neumann. Cuthbert Hurd and John Backus have reasonably established that the FORTRAN project got a green light in about December 1953. At that time, von Neumann was a consultant to IBM and worked closely with Hurd, who had hired him. Daniel Leeson first suggested that von Neumann had attended a meeting at which FORTRAN was discussed. JAN and I asked Hurd whether or not von Neumann had any link to the FORTRAN project in its early stages and/or if he had even expressed his feelings about the need or importance for such a project. Hurd did remember an early meeting on the subject of language development which von Neumann attended, as had John Backus, Frank Beckman, and John Greenstadt, but he could remember no specifics on any aspect. On reflection, he thought that von Neumann had not been impressed.

fortran66.hatenablog.com

その数式、プログラムできますか?

その数式、プログラムできますか?

From Mathematics to Generic Programming (English Edition)

From Mathematics to Generic Programming (English Edition)

Fortran2018 の 16 進フォーマット

Hexadecimal Format

J.Reid の The New Features of Fortran 2018 中の Hexadecimal input/output を intel fortran 19.1 の beta 版で試してみます。出力は文中の possible output と少し違って、桁がずれて出力されました。

N2145 wg5-fortran.org

16 進数浮動小数点用フォーマットといえば、昔の IBM360 のものが有名ですが、なぜ今頃これが導入されるのかよく分かりません。

ソース・プログラム

とりあえず二進表記にして、4桁づつまとめて 16 進数にすれば対応がつかめます。

    program Console3
        implicit none

        print *, 1.375, '= 1.0 + 0.25 + 0.125 = 1+1/4+1/8 = 1.0110 = 0X1.6'
        print *, int(b'10110')/16.0 
        print '(a, ex0.1)', '   1.375 = ', 1.375  
        print *
!
        print *, -15.625, '= -(15.0 + 0.5 + 0.125 = 15+1/2+1/8) = -1111.1010 = -0XF.A' 
        print *, -int(b'11111010000000000')/16.0**3/2
        print '(a, ex14.4e3)', '  -15.625 = ', -15.625
        print *
!
        print *, 1048580.0, '= 16^4 + 4 = 1 0000 0000 0000 0000 0100 = 1000.0000 0000 0000 0010*2^17 = 8.0002P17'
        print *, b'100000000000000000100'
        print '(a, ex0.0)', '   1048580.0 = ', 1048580.0
        
    end program Console3

出力結果

Reid の与えた possible output は小数点位置を常に 1.xxx の形に位取りしていますが、intel fortran は 4 ビット分の 15 以下の数で位取りしているようです。

   1.375000     = 1.0 + 0.25 + 0.125 = 1+1/4+1/8 = 1.0110 = 0X1.6
   1.375000
   1.375 = 0XB.0P-3

  -15.62500     = -(15.0 + 0.5 + 0.125 = 15+1/2+1/8) = -1111.1010 = -0XF.A
  -15.62500
  -15.625 = -0XF.A000P+000

   1048580.
 = 16^4 + 4 = 1 0000 0000 0000 0000 0100 = 1000.0000 0000 0000 0010*2^17 = 8.0002P17
     1048580
   1048580.0 = 0X8.0002P+17

Modern Fortran Explained: Incorporating Fortran 2018 (Numerical Mathematics and Scientific Computation)

Modern Fortran Explained: Incorporating Fortran 2018 (Numerical Mathematics and Scientific Computation)