第6回 配信講義　計算科学技術特論A （2023）

線形代数演算ライブラリBLAS とLAPACKの基礎と実践 (I) BLAS, LAPACK入門編 中田 真秀 理化学研究所 開拓研究本部 柚木研 2023-05-25 13:00 (JST) 計算科学技術特論A 13:00

BLAS, LAPACK入門編 講義内容 1. 線形代数は重要 2. あらゆる分野で標準化が重要 3. コンピュータでの数値計算 4. コンピュータでの線形代数 5. BLAS, LAPACKは線形代数の事実上の標準的ライブラリ 6. BLASのプログラム例: dgemm (行列-行列積) 7. LAPACKのプログラム例: dsyev (行列の対角化) 8. はじめての行列-行列積高速化

線形代数は重要 • いまからでも遅くないから線形代数ちゃんと勉強または復習しと こう • 線形代数 連立一次方程式を解く 行列-行列積 • 学問の基本的に全ては線形代数で書ける • なんでも(関数さえ)ベクトルに見えてくる→免許皆伝 • 抽象度が高い→自分の分野でイメージ作りするとよい • 中田の場合は量子力学; 非常にイメージしやすかった。 • 人間は何かを線形代数の言葉でしか考えられないのかもしれない a11 a12 a21 a22 A= ⋮ ⋮ am1 am2 ⋯ ⋯ ⋱ ⋯ a1n a2n ⋮ amn x1 x = x2 x3

ChatGPT4先生に聞いてみた

我々はコンピュータで 線形代数ばかりやってきた • 機械学習、AI • ChatGPT, Stable Di usion • 行列-行列積を大量にやっている • PlayStation5などでの3Dゲーム • 宇宙は巨大な量子コンピュータ ff https://prtimes.jp/main/html/rd/p/000000033.000035388.html

• https://prtimes.jp/main/html/rd/p/000000033.000035388.html

線形代数は紀元前までさかのぼれる • 九章算術 方程より • 中国、紀元前1世紀から紀 元後2世紀ころ、著者不詳 りゅう き • 263年頃魏の時代に劉徽に よって整理と注釈が加えら れた。 • https://ctext.org/nine-chapters/fang-cheng/zh • 人類史上初めての連立一次 • 方程式をGaussの消去法で 解いたと思う。 今でも1000年以上前の文 が何となく読めるのは凄い https://zh.wikisource.org/wiki/Page:Sibu̲Congkan0392-劉徽-九章算術-3-3.djvu/8 より

• https://zh.wikisource.org/wiki/Page:Sibu_Congkan0392-劉徽-九章算術-3-3.djvu/8

九章算術 方程より 問題 • 有上禾三秉，中禾二秉，下禾一秉，實三十九斗；上禾二秉，中禾三 秉，下禾一秉，實三十四斗；上禾一秉，中禾二秉，下禾三秉，實二 十六斗。問上、中、下禾實一秉各幾何？答曰:上禾一秉，九斗、四分 斗之一，中禾一秉，四斗、四分斗之一，下禾一秉，二斗、四分斗之 三。 • 方程術曰，置上禾三秉，中禾二秉，下禾一秉，實三十九斗，於右 方。中、左禾列如右方。以右行上禾遍乘中行而以直除。又乘其次， 亦以直除。然以中行中禾不盡者遍乘左行而以直除。左方下禾不盡 者，上為法，下為實。實即下禾之實。求中禾，以法乘中行下實，而 除下禾之實。餘如中禾秉數而一，即中禾之實。求上禾亦以法乘右行 下實，而除下禾、中禾之實。餘如上禾秉數而一，即上禾之實。實皆 如法，各得一斗。

九章算術 方程より 現代語訳 Powered by Google翻訳 • 問: 3束の上質のキビ、2束の中質のキビ、1 束の低質のキビが39個 のバケツに入っている。2束の上質のキビ、3束の中質のキビ、1束 の低質のキビが34個のバケツに入っている。1束の上質のキビ、2 束の中質のキビ、3束の低質のキビが26個のバケツに入っている。 上質、中質、低質のキビ1束はそれぞれバケツいくつになるか。 • 答: 上質 9 ¼ , 中質 4 ¼ , 低質 2 ¾ 個づつ • 上質のキビ3束、中質のキビ3束、低質のキビ1束を39バケツを右行 に置く。中行、左行も右のように並べる。右の上質を中行にかけ、 右行で引く。また左行にもかけて右行から引く。次に、中行の中質 のキビの余りを左行にかけて、中行で引く。左の低質に余りがある のでそして、割れば求まる(実を法で割る)。以下略

九章算術 方程より 現代語訳 Powered by Google翻訳 •問 3x + 2y + z = 39⋯(右) 2x + 3y + z = 34⋯(中) x + 2y + 3z = 26⋯(左) • (右)はそのまま、(中)は(中)を3倍したものから(右)を2倍したものを引き、 (左)を3倍して(左)から(右)を引く。 3(2x + 3y + z = 34) 2(3x + 2y + z = 39) 3(x + 2y + 3z = 39) 3x + 2y + z = 39 5y + z = 24⋯(中) 4y + 8z = 39⋯(左) • それから(左)を5倍する 3x + 2y + z = 39⋯(右) 5y + z = 24⋯(中) 20y + 40z = 195⋯(左) 36c = 99 あとは略

線形代数チートシート • • a11 a12 a21 a22 A= 行列 ⋮ ⋮ am1 am2 ベクトルと基底 b = (b1, b2, b3, ⋯, bn) x1 x = x2 x3 ⋯ ⋯ ⋱ ⋯ a1n a2n ⋮ amn ∀x = ∑ i ciei 1 e1 = 0 (0) • 行列、ベクトルの線形結合 a, b ∈ R, x, y ∈ V → ax + by ∈ V n • 内積 a·b = ∑ aibi i • 行列と行列の掛け算 (AB)ij = ∑ k Aik Bkj

線形代数チートシート • 連立一次方程式を解く • 行列の特異値問題 • 行列の固有値問題 −1 Ax = b ⟶ x = A b A = UΣV Ax = λx • エルミート行列とユニタリ行列 Hij = • 行列の固有値と固有ベクトル λ1 • 対角化 Hx = λx U †HU = 0 λ2 T −1 H* U ji ij 0 λ3 ⋱ λ∞ = U* ji

あらゆる分野で標準化が重要 13:12

あらゆる分野で標準化が重要 • International Organization for Standardization (ISO) • ISOは、1国1加盟を基本として、157カ国の国家標準化機関が加盟するネットワーク • 非常口ロゴマーク、カードサイズ、ねじ、食品安全、組織マネジメント... • Request For Comments (RFC) • Internet Engineering Task Forceが作る、インターネットのプロトコルの標準化。 個人参加が可能 • Système International d'unités (SI) 単位系 • IEEE Standards Association • 160ヶ国以上の共働志向のリーダーとともに、国際市場を創造・拡大し 健康と公共の安 全の保護を支援する合意形成に基づく組織 • 通信の規格 IEEE 802.11, • 数値計算の標準化: IEEE 754 • binary32, binary64, binary128などの精度 • de fact standard:事実上の標準。

あらゆる分野で標準化が重要 • 標準化がないとどうなるか • コンセントは各国で形状も電圧もばらばら • 自動車は右を走るべきか左を走るべきか？ • 日本特有:電気の交流は50Hz or 60Hz • 規格策定は各メーカー壮絶な戦い:ルール作る側はいつも強い • SIがあるのに、ヤードポンド法を使うアメリカ (ネジが合わない) • IEEE 754がない時 • 同じ数値計算プログラムの結果がコンピュータ、コンパイラのメーカーな どによって変わる事多々。 • 線形代数演算には特に標準化はされてない • しかし、往々にして規格というものは時代に追いつかない。

線形代数演算の標準化も重要 • 線形代数演算をコンピュータ上で行うなんらかの標準化はあっ たほうが良い • 事実上の標準ライブラリしかない • 商業ベースで開発は難しいか。 • LAPACKは大学で外部資金で主に開発。NAGなど • 似たライブラリが乱立すると極めてユーザーにとって不便 • BLASとLAPACKがあってよかった

コンピュータでの数値計算 13:16

コンピュータでの実数演算はどうするか • コンピュータは有限の整数しか扱えないため、特別な表記 (フォーマット)を使う。 • 浮動小数点数表記 2進数で32桁、64桁などのビット列を実数と みなす。 • 浮動小数点数は、符号、仮数部(fraction)、指数部(exponent)、 anは0 or 1から成る fraction n exponent ⏞ 1 ± 1+ an × 2m ∑ (2) n=1 k • 浮動小数点数を10進数で表す例。4ビット2進数−1.101 × 25 を10 進数になおす。 −1.101 × 25 = − (1 + 1 × 0.5 + 0 × 0.25 + 1 × 0.125) × 32 = 52

コンピュータでの数の取扱いbinary64 (倍精度) • “754-2019 IEEE Standard for Floating-Point Arithmetic” • binary64 フォーマットは 10進16桁の有効桁がある (よく倍精度とよぶ) fraction n exponent 1 ± 1+ an × 2−1022∼1023 ∑ (2) n=1 52 • • binary32,128 などもある (単精度、四倍精度とよく呼ばれる)。 この規格に則って演算する場合がほとんど(最近の例外:PlayStation2) • 規格が無いときはコンピュータ毎に違う結果が得られたりした • b oat16等、半精度もでてきたが、IEEE 754規格外!! 混乱の予想 • FLOPS(フロップス）という単語が頻出 • 1秒間に1回浮動小数点数が計算できること=Floating point operation per second • 速さ:Core i7 (Broadwell, 10 cores, 3.5GHz): 560GFlops, NVIDIA fl TESLA H100: 34TFlops〜PFlops, 京コンピュータ10PFLOPS, HOKUSAI (3PFlops)、富岳：400PFlops (2022/5まで世界一) 世界一: Frontier@US 1.1ExaFlops

コンピュータでの線形代数 13:19

自作は避けたほうがいい • 線形代数演算をコンピュータでやるとき、プログラムを自作する場合が あるかもしれないが、自作は避けたほうがいい • 無駄な計算をしてしまうことがある。 • クラメールの公式を使うと宇宙時間待っても答えが出ない。 • 誤差が大きくなることがある。 • 逆行列は数値誤差を多く含む 数値アルゴリズムにおける精度と安定性 Accuracy and Stability of Numerical Algorithms, N. Higham 2002 • 行列の条件数が大きいと誤差も大きくなる • 速いコンピュータのはずがパフォーマンスが出ない • カタログに出てる理論性能値と比べて、1/100の速度しか出ない？ 計算科学のためのHPC技術1 まずは既存のライブラリを探すこと 安直だが

自作して失敗する例:連立一次方程式を クラメールの公式で解く • 以下のような連立一次方程式を解く。 a11x1 + a12 x2 + ⋯ + a1n xn = b1, a21x1 + a22 x2 + ⋯ + a2n xn = b2, ⋮ an1x1 + an2 x2 + ⋯ + ann xn = bn a11 a12 a21 a22 A := ⋮ ⋮ an1 an2 ⋯ ⋯ ⋱ ⋯ a1n a2n , ⋮ ann x1 x2 x := , ⋮ xn b1 b2 b := ⋮ bn Ax = b • 解は a1,1 ⋯ a1,i−1 b1 a1,i+1 ⋯ a1,n a2,1 ⋯ a2,i−1 b2 a2,i+1 ⋯ a2,n Ai := ⋮ ⋮ ⋮ ⋮ ⋮ an,1 ⋯ an,i−1 bn an,i+1 ⋯ an,n • として、行列Bに対する行列式det(B)を用いて、以下のようにかける。 det(Ai) xi = det(A)

自作して失敗する例:連立一次方程式を クラメールの公式で解く • 理論的にはクラメールの公式を使ってプログラムを作ると、 コンピュータで解ける。 • しかしながら、行列のサイズ n が大きくなるとすぐ解けな くなる。 • n=10で 36288000 回の演算 • n = 100 で 9.3 x 10157回の演算 • あっというまに宇宙時間より長い時間が必要になる。

自作して失敗する例:行列-行列積 ナイーブに作ると30-40倍くらい遅くなることがある̲ • 2000x2000の行列積AB＝Cを教科書どおりに計算するプログ ラムを作って • 0.4GFlops (AB)ij = ∑ k Aik Bkj •最適化されたBLASを用いた場合 $ octave ... 途中略... octave:1> n=4000; A=rand(n); B=rand(n); octave:2> tic(); C=A*B; t=toc(); octave:3> GFLOPS=2*n^3/t*1e-9 GFLOPS =84.1GFlops マシン：MacBookPro (2017) i5-7360U 2.3-3.6GHz 2 cores 理論性能：73-115GFlops

分野の違いと意識の違い (偏見あり) • 数学者の意識 : 原理的に可能, 解の存在のみ興味ある場合が多い • 情報系の数学より : アルゴリズムが多項式程度なものを考えたがる • 自然科学系研究者 : ともかく答えが求まる方法なら何でも良い • とりあえず求まればよい。問題が出るまで放置。よく指数関数的なアル ゴリズムを意識せずにゴリ押しする。 • HPC or 数値解析系研究者 : 1 clockでも速い方法 1bitでも転送量が少な い方法、1桁でも精度の良い方法などから選択 • ハード依存高め。さまざまな現実的な制限を考慮し、なるべく良い結果 を出す。

分野の違いと意識の違い (偏見あり) 連立一次方程式をどう解くか • 数学者の意識 : 行列式が0でないなら解ける。終わり • 情報系の数学より : LU分解も逆行列求める方法もオーダーは同じ O(n^3)なので違いはない、クラメールはO(n!)なんで論外。 • 自然科学系研究者 : とりあえず逆行列求めとこう。LU分解って何? 解 きたい問題は別だし他に考えるべきことが多いし後まわしにしよう。 • 最近はnumpyになげときゃいいや (実は正解) • HPC or 数値計算系研究者 : LU分解経由一択。基本であるdgemm だからキャッシュヒットも高いし、逆行列を使うより誤差もない少 よってLU分解経由以外あり得ない。クラメールは誤差も大きい。

BLAS, LAPACKは事実上の 標準の線形代数ライブラリ 13:28

線形代数演算ならBLAS, LAPACKを使おう • 行列、ベクトルなどの線形代数演算のライブラリはある? • BLASやLAPACKを用いましょう。 • コンピュータの仕組みを軽く知っておくと、手軽に高速な ライブラリの利用方法がわかる • 今回はBLAS, LAPACKをとりあえず使ってみる、という ことをやります

BLAS, LAPACK:世界最高の線形代数演算パッケージ • コンピュータで線形代数演算するならBLAS+LAPACKを使いましょう。 • 品質、信頼性がとても高く、無料で入手できる • http://www.netlib.org/lapack/ • 標準で高速版がインストールされていることが多い • http://www.openblas.net/ • 密行列のみ対応(疎行列は他のライブラリを利用する)。 • コンピュータでの行列線形代数演算の基礎中の基礎! BLAS, LAPACKは人類の至宝

BLAS, LAPACKを利用したソフトウェア • 著名な計算プログラムパッケージは大抵BLAS, LAPACKを利用して いる. • 機械学習系はほぼcuBLASを使ってる • 物理、化学ではGaussian, Gamess, ADF, VASP • 線形計画問題のCPLEX, NUOPT, GLPKなど.. • 高級言語からも利用可能 • Ruby, Python (numpy), Perl, Java, C, Mathematica, Maple, Matlab, R, octave, SciLab

BLAS : 基本的な線形代数サブプログラム群 •BLASはBasic Linear Algebra Subprogramsの略 •基礎的な線形代数の「サブ」プログラム Level 1 〜 3まである •Level 1 : ベクトル-ベクトルの内積など x ⋅ y •Level 2 : 行列-ベクトル積など y ← αAx + βy •Level 3 : 行列-行列積など C ← αAB + βC •Fortran90でさまざまなルーチンの仕様を提供している。 •参照実装の形で提供されている (Reference BLAS) •BLASのルーチンを「ブロック」にしてより高度なことをする。 •この実装を「お手本」とする •とても美しいコード! •高速版もある。

Level 1 BLAS • Level 1:ベクトル-ベクトル演算(+そのほか)のルーチン • ベクトルの加算 DAXPY y ← αAx + βy • 内積計算: DDOT < x, y > = ∑ i • 2-ノルム計算 ||x|| = N ∑ i xi yi | xi |2 • など15種類あり, さらに単精度, 倍精度, 複素単精度, 複素数倍精 度についての4通りの組み合わせがある.

Level 2 BLAS • Level 2:行列-ベクトル演算のルーチン • 行列-ベクトル積: DGEMV y ← αAx + βy • 上三角行列とベクトルの積:DTRMV x ← Ax • 上三角行列の連立一次方程式を解く:DTRSV x ← A −1x • 列ベクトルと行ベクトルの積: DGER A ← αxy t + A • など25種類あり, 同じように単精度、倍精度、複素数の4通 りの組み合わせがある。

Level 3 BLAS • Level 3 BLASは行列-行列演算のルーチン群 • 行列-行列積: DGEMM C ← αAB + βC • 対称行列-行列積: DSYMM C ← αAB + βC • 上(下)三角行列と行列の積: DTRMM B ← αAB • 対称行列の階数nの更新: DSYRK C ← αAA T + βC • 上三角行列の連立一次方程式を解く: DTRSM • など9種類ある。 B ← αA −1B

BLAS Quick Reference https://www.netlib.org/lapack/lug/node145.html

• https://www.netlib.org/lapack/lug/node145.html

BLAS Quick Reference https://www.netlib.org/lapack/lug/node145.html

• https://www.netlib.org/lapack/lug/node145.html

LAPACKとは? •LAPACK(Linear Algebra PACKage) : 線形代数パッケージ •BLASをビルディングブロックとして使いつつ、より高度な問題である連立一次方 程式、 最小二乗法固有値問題、固有値問題、特異値問題を解くことができる. •下請けルーチン群も提供する: 行列の分解(LU分解, コレスキー分解, QR分解, 特異 値分解, Schur分解, 一般化Schur分解)、条件数の推定ルーチン, 逆行列計算など。 •品質保証も非常に精密かつ系統的で、信頼がおける。 •パソコンからスーパーコンピュータまで様々なCPU、OS上で動く。 •Fortran 90で書かれ、3.11.0は1900以上のルーチンからなっている。 •webサイトはなんと約2億ヒット以上である •githubで開発が続いている https://github.com/Reference-LAPACK http://www.netlib.org/lapack

• http://www.netlib.org/lapack

LAPACK公式ドキュメント • http://www.netlib.org/lapack/lug/ : ユーザーガイ ド • http://www.netlib.org/lapack/faq.html : FAQ • http://www.netlib.org/lapack/lawns/index.html LAWN: LAPACK Working Notes : 実装の詳細、アル ゴリズム、パフォーマンスの比較など。 • 不断の努力の積み重ね

• http://www.netlib.org/lapack/lug/
• http://www.netlib.org/lapack/faq.html
• http://www.netlib.org/lapack/lawns/index.html

線形代数+コンピュータで最重要タスクたち •連立一次方程式問題 : Ax=b •最小二乗法 min ||b-Ax|| •固有値問題 Ax=λx •特異値問題 M = UΣV* 規模、精度、行列のタイプ、解き方に多様な応用がある

LAPACKのルーチンの種類 • Driver routines : 先程あげた固有値、連立一次方程式を解く • Simple driver: • Expert driver: Simple driverに比べて、条件数推定、解の改善、エラー バウンド、行列の平衡化などを行う • Computational routines • 上記タスクなどのために行うLU分解や三角行列のリダクションを行うが BLASよりは高級な処理を行う • Auxiliary routines • blockアルゴリズムのサブタスク、行列ノルム、スケーリングなどBLASの拡 張またはBLASに入れたほうがいいルーチンなど低レベル処理

LAPACKで連立一次方程式を解く simple driverたち http://www.netlib.org/lapack/lapackqref.ps

• http://www.netlib.org/lapack/lapackqref.ps

LAPACKで最小二乗法を解く simple driverたち http://www.netlib.org/lapack/lapackqref.ps

• http://www.netlib.org/lapack/lapackqref.ps

LAPACKで一般化固有値問題、一般化特異値問題を解く simple & divide and conquer driverたち http://www.netlib.org/lapack/lapackqref.ps

• http://www.netlib.org/lapack/lapackqref.ps

LAPACKで標準固有値問題、特異値問題を解く simple & dvide and conqure driverたち http://www.netlib.org/lapack/lapackqref.ps

• http://www.netlib.org/lapack/lapackqref.ps

様々な解法が存在していて、 様々なルーチンが存在する • たくさんLAPACKのルーチンを提示したが、これにそれ ぞれExpert driverや、RRR (relatively robust representation) 版などが存在する。 • simple/divide and conquer/RRR/Expertからどう やって選べばよいか? • これは問題に応じて個々人が選ぶ必要が出てくる。

LAPACKのルーチン構造 • 例えば実対称行列の固有値を求めるのにはdsyevを使ったが下 請けには34のルーチン群がある。 • dorgtr, dorgql, dorg2l, dorgqr, dlarfb • dlarf, dgemm, dcopy, dtrmv, dgemv, dger • dsyr2k, dlatrd, dsytd2, daxpy, dsymv, dlarfg, • dsyr2, dscal, dsteqr, dsterf, dlaev2, dlartg, • dlaset, swap, dlascl, dlasr, dlasrt, dlae2 dsyevルーチン相関図

LAPACKのルーチン構造 • 実特異値分解はもっと複雑。dgesvdだけでも 3503行あるが、殆どが総計46の下請けルー チンをコールしている。 dgesvdルーチン相関図

BLASプログラム例:dgemm 行列-行列積 13:40

BLAS、LAPACKのインストール • BLAS、LAPACKの開発環境を整える on Ubuntu 22.04 $ sudo apt-get update ; sudo apt-get upgrade $ sudo apt-get install gfortran g++ libblas-dev liblapack-dev liblapacke-dev ... g++ is already the newest version (4:11.2.0-1ubuntu1). gfortran is already the newest version (4:11.2.0-1ubuntu1). libblas-dev is already the newest version (3.10.0-2ubuntu1). liblapack-dev is already the newest version (3.10.0-2ubuntu1). liblapacke-dev is already the newest version (3.10.0-2ubuntu1). 0 upgraded, 0 newly installed, 0 to remove and 6 not upgraded. • こんな感じであればok

行列-行列の積 DGEMMを使ってみる • 行列-行列積DGEMMを使ってみる。ここでは 1 8 3 A = 2 10 8 9 −5 −1 9 8 3 B = 3 11 2.3 −8 6 1 α = 3,β = − 2 C ← αAB + βC • を計算するプログラムを書いてみる. • 答えは以下のようになる 21 336 70.8 −64 514 95 210 31 47.5 3 3 1.2 C= 8 4 8 6 1 −2

行列-行列の積:DGEMMの詳細 C ← αAB + βC • 今回はCBLASから、BLASを呼んでみる。 • BLASはFORTRANで書かれているが、Cから呼び出すことも可能。プロトタイプ宣言は以下のようになる void cblas̲dgemm(CBLAS̲LAYOUT layout, CBLAS̲TRANSPOSE TransA, CBLAS̲TRANSPOSE TransB, const int M, const int N, const int K, const double alpha, const double *A, const int lda, const double *B, const int ldb, const double beta, double *C, const int ldc); • layout : Column major or Row major • “transa”, “transb”, “transc” で行列を転置するか否かを指定。 • m, n, k は行列の次元。左図参照 • alpha, betaは行列積に対する掛けるスカラー • A, B, CはRow major形式で格納した行列へのポインタ • lda, ldb, ldc は 行列A, B, Cへのleading dimensionたち • 配列はゼロから始まるか1からはじまるか ここらへんがわかりにくいところ

FortranとCの違いの吸収 配列は0か1から始める? 多次元配列を手で展開 • Fortranでは配列は1からスタートするのが普通が, C, C++では, 0か らのみスタートする. 例えば • ループの書き方が一般的には1からNまで(FORTRAN)か, 0からn 未満か(C,C++). • ベクトルのxi要素へのアクセスはFORTRANではX(I)だが, Cでは x[i-1]となる. • 多次元配列はC/C++にはない。一次元配列に手で展開する。 • 行列のAij要素へのアクセスはFORTRANではA(I,J)だが, Cでは column majorとしてA[i-1+ (j-1)*lda]とするとよい

行列をColumn majorでメモリに格納する • 行列は2次元だが、コンピュータのメモリは1次元的である。次のような行 列を A= 1 2 3 (4 5 6) • 考えるとき、どのようにメモリに格納するか? column major式では、 • アドレスの小さい順から 1, 4, 2, 5, 3, 6 • のように格納する。 • FORTRANや、Matlab, octaveはcolumn majorである。

行列をRow majorでメモリに格納する 同じように、行列A 1 2 3 A= (4 5 6) • について、row majorでは、行方向に • 1, 2, 3, 4, 5, 6 • のように格納する。C/C++では普通row majorで格納。意 識的にcolumn majorにしたほうが良い • BLAS, LAPACKを呼び出すときは、行列格納方法に注意!

行列-行列の積 DGEMMを使ってみる • 行列-行列積DGEMMを使ってみる。ここでは 1 8 3 A = 2 10 8 9 −5 −1 9 8 3 B = 3 11 2.3 −8 6 1 α = 3,β = − 2 C ← αAB + βC • を計算するプログラムを書いてみる. • 答えは以下のようになる 21 336 70.8 −64 514 95 210 31 47.5 3 3 1.2 C= 8 4 8 6 1 −2

行列-行列の積のリスト
#include <stdio.h>
#include <cblas.h>
//Matlab/Octave format
void printmat(int N, int M, double *A, int LDA) {
double mtmp;
printf("[ ");
for (int i = 0; i < N; i++) {
printf("[ ");
for (int j = 0; j < M; j++) {
mtmp = A[i + j * LDA];
printf("%5.2e", mtmp);
if (j < M - 1) printf(", ");
} if (i < N - 1) printf("]; ");
else printf("] ");
} printf("]");
}
int main()
{
int n = 3; double alpha, beta;
double *A = new double[n*n];
double *B = new double[n*n];
double *C = new double[n*n];
A[0+0*n]=1; A[0+1*n]= 8; A[0+2*n]= 3;
A[1+0*n]=2; A[1+1*n]=10; A[1+2*n]= 8;
A[2+0*n]=9; A[2+1*n]=-5; A[2+2*n]=-1;
B[0+0*n]= 9; B[0+1*n]= 8; B[0+2*n]=3;
B[1+0*n]= 3; B[1+1*n]=11; B[1+2*n]=2.3;
B[2+0*n]=-8; B[2+1*n]= 6; B[2+2*n]=1;

C[0+0*n]=3; C[0+1*n]=3; C[0+2*n]=1.2;
C[1+0*n]=8; C[1+1*n]=4; C[1+2*n]=8;
C[2+0*n]=6; C[2+1*n]=1; C[2+2*n]=-2;
printf("# dgemm demo...\n");
printf("A
=");printmat(n,n,A,n);printf("\n");
printf("B
=");printmat(n,n,B,n);printf("\n");
printf(“C
=");printmat(n,n,C,n);printf("\n");
alpha = 3.0; beta = -2.0;
cblas_dgemm(CblasColMajor,CblasNoTrans,Cbl
asNoTrans, n, n, n, alpha, A, n, B, n,
beta, C, n);
printf("alpha = %5.3e\n", alpha);
printf("beta = %5.3e\n", beta);
printf("ans="); printmat(n,n,C,n);
printf("\n");
printf("#check by Matlab/Octave by:\n");
printf("alpha * A * B + beta * C =\n");
delete[]C; delete[]B; delete[]A;
}

行列-行列の積のコンパイルと実行 • 先ほどのリストを''dgemm_demo.cpp''などと保存する。 $ g++ dgemm_demo.cpp -o dgemm_demo -lblas -lapack • でコンパイルができる. 何もメッセージが出ないなら, コンパイルは成功である。実行 は以下のようになっていればよい。OctaveやMatlabにこの結果をそのままコピー& ペースとすれば答えをチェックできるようにしてある alpha = 3.000e+00 beta = -2.000e+00 ans=[ [ 2.10e+01, 3.36e+02, 7.08e+01]; [ -6.40e+01, 5.14e+02, 9.50e+01]; [ 2.10e+02, 3.10e+01, 4.75e+01] ] #check by Matlab/Octave by: alpha * A * B + beta * C

Block行列が便利になる例 • 行列の積を考える C = AB, Cij = ∑ k Aik Bkj • 行列積の定義は、要素ごとに積をとって和を取るだが、区 分行列にわけても、そのまま「数」のように積をとってよ い A A ⋯ A B B ⋯ B 11 12 1q ⋯ A2q C11 C12 C21 C22 AB = ⋮ ⋮ Cp1 Cp2 ⋯ ⋯ ⋱ ⋯ C1r C2r ⋮ Cpr ⋮ ⋮ Ap1 Ap2 ⋱ ⋯ 12 B21 B22 , B= ⋮ ⋮ ⋮ Bq1 Bq2 Apq A21 A22 A= 11 Cij = q ∑ k=1 ⋯ ⋱ ⋯ Aik Bkj 1r B2r ⋮ Bqr

Leading dimensionについて • 行列をさらに小さい行列に分けて考えることがある。 • これらを区分行列、小行列、ブロック行列などとよぶ。 • たとえば 以下のように、A, B, C, Dという行列を考えて 2 −1 5 A = −1 4 1 , (8 1 −2) −3 6 B= 1 3 , ( ) 4 1 C = (−4 2 6), D = (9 1) • それらを組み合わせてより大きな行列を作ることができ る。

Leading dimensionについて • 行列Aの区分行列A'にアクセスするにはど うしたらいい？ • leading dimensionを使うとよい。 • numpyのスライス、Fortran90の部分配 列に相当する。しかしCやC++にはない... • 関数をまたぐと、スライスでは対応で きなくなり、ldを渡す必要アリ • Aʼ[1,1] のアドレスから、Aʼ[P,Q]は Aʼ[1,1]+P*m+Qではなくて、 Aʼ[1,1]+P*LDA+Qとアクセスする。

たかが(?)行列積なのに 超めんどくさい • numpyならA*Bですむのになんで？ • 高速なライブラリをちゃんと設計するとこうなる • コンピュータへの分かりやすさ、人間への分かりやすさの違い

LAPACKプログラム例: 対称行列の対角化 (dsyev) 14:00

LAPACK実習: 行列の固有ベクトル、固有値を求める:DSYEV • 3x3 の実対称行列の固有ベクトル、固有値を求めよう。これらは三つあり、 1 2 3 A= 2 5 4 3 4 6 Avi = λivi (i = 1,2,3) という関係式が成り立つ。 • 固有値 λ1, λ2, λ3 は • で、 固有ベクトルは、 −0.40973,1.57715,10.83258 v1 = (−0.914357,0.216411,0.342225) v2 = (0.040122, − 0.792606,0.608413) v3 = (0.402916,0.570037,0.716042) となる

行列の固有ベクトル、固有値を求めるDSYEV •C/C++からLAPACKをコールするにはLAPACKEを用いる。プロトタイプは以下の通り •lapack_int LAPACKE_dsyev( int matrix_layout, char jobz, char uplo, lapack_int n, double* a, lapack_int lda, double* w ); •matrix̲layout : LAPACK̲COL̲MAJOR: column majorかrow majorか •jobz:固有値、固有ベクトルが必要か、固有値だけでよいか指定。 •uplo:行列の上三角、下三角を使うか。 •A, lda:行列Aとそのleading dimension •w:固有値を返す配列(昇順) (注意) 以下はFORTRANを直接呼ぶ場合は必要だったが、LAPACKEではなくなった • work, lwork:ワーク領域への配列またはポインタ、とそのサイズ • info: =0 正常終了。<0: INFO=-iではi番目の引数が不適当。>0: INFO=i 収束せず。

対称行列の対角化dsyevの例
#include <iostream>
#include <stdio.h>
#include <lapacke.h>
//Matlab/Octave format
void printmat(int N, int M,
double *A, int LDA) {
double mtmp;
printf("[ ");
for (int i = 0; i < N; i++) {
printf("[ ");
for (int j = 0; j < M; j++) {
mtmp = A[i + j * LDA];
printf("%5.2e", mtmp);
if (j < M - 1) printf(", ");
} if (i < N - 1) printf("]; ");
else printf("] ");
} printf("]");
}

int main()
{
int n = 3;
double *A = new double[n*n];
double *w = new double[n];
//setting A matrix
A[0+0*n]=1;A[0+1*n]=2;A[0+2*n]=3;
A[1+0*n]=2;A[1+1*n]=5;A[1+2*n]=4;
A[2+0*n]=3;A[2+1*n]=4;A[2+2*n]=6;
printf("A ="); printmat(n, n, A, n);
printf("\n");
LAPACKE_dsyev(LAPACK_COL_MAJOR,
'V', 'U', n, A, n, w);
//print out some results.
printf("#eigenvalues \n"); printf("w =");
printmat(n, 1, w, 1); printf("\n");
printf("#eigenvecs \n"); printf("U =");
printmat(n, n, A, n); printf("\n");
printf("#Check Matlab/Octave by:\n");
printf("eig(A)\n");
printf("U'*A*U\n");
delete[]w;
delete[]A;
}

対称行列の対角化dsyevの例 • 先ほどのリストを''eigenvalue̲demo.cpp''などと保存する。 •g++ dsyev_demo.cpp -o dsyev_demo -llapacke -lblas -llapack でコンパイルができる。何もメッセージが出ないなら, コンパイルは成功である。 • 実行は以下のようになっていればよい。同様にOctaveやMatlabにこの結果をそのままコピー&ペーストす れば答えをチェックできるようにしてある。 $ ./dsyev_demo A =[ [ 1.00e+00, 2.00e+00, 3.00e+00]; [ 2.00e+00, 5.00e+00, 4.00e+00]; [ 3.00e+00, 4.00e+00, 6.00e+00] ] #eigenvalues w =[ [ -4.10e-01]; [ 1.58e+00]; [ 1.08e+01] ] #eigenvecs U =[ [ -9.14e-01, 2.16e-01, 3.42e-01]; [ 4.01e-02, -7.93e-01, 6.08e-01]; [ 4.03e-01, 5.70e-01, 7.16e-01] ] #Check Matlab/Octave by: eig(A) U'*A*U

はじめての行列-行列積高速化 14:12

行列-行列積には高速化が必須 ナイーブに作ると30-40倍くらい遅くなることがある̲ • 2000x2000の行列積AB＝Cを教科書どおりに計算するプログ ラムを作って • 0.4GFlops (AB)ij = ∑ k Aik Bkj •最適化されたBLASを用いた場合 $ octave ... 途中略... octave:1> n=4000; A=rand(n); B=rand(n); octave:2> tic(); C=A*B; t=toc(); octave:3> GFLOPS=2*n^3/t*1e-9 GFLOPS =84.1GFlops マシン：MacBookPro (2017) i5-7360U 2.3-3.6GHz 2 cores 理論性能：73-115GFlops

歴史:高速なBLAS, LAPACKの実装について • Reference BLAS/LAPACKはある意味仕様書そのままなので、非常に低速である。 メモリの階層構造などは非常に意識して書かれているが、CPUに最適化は、各々がや ることになる。 • 2000年までは、だいたいSUN/IBM/DEC/Intel/Fujitsu/Hitachi/NECなどCPU、 OS,システムごとにベンダーの実装があり、高価で販売されていた、またはマシンを買 うとついてきたりした。 • ATLAS:R. Clint Whaley氏による, オートチューニング機構で高速化したBLAS。そ れまでの2001年より多くのコンピュータのBLAS環境を劇的に改善した, パイオニア 的存在。ハンドチューニングしたBLASよりは数%から数10%低速程度(オープンソー ス) • GotoBLAS/GotoBLAS2 : 後藤和茂氏が、アセンブラで様々なCPUに対応した BLASのソースコードを公開。マシンの性能の限界近くまで性能を追求(非オープンソー ス、最終バージョンはオープンソースとなった)。 2000年ころから高速なBLAS/LAPACKの環境が良くなってきた

現状: 高速なBLAS, LAPACKの実装について • OpenBLAS: Zhang Xianyi氏がGotoBLAS2の開発を引き継いだ。開発はアクティブ でSandyBridge以降のプロセッサ, OSにも対応している。ARM各種、AMD、Power, ICT Loongson-3A, 3Bにも対応。オープンソース https://www.openblas.net/ • Intel Math Kernel Library: Intelが開発している加速されたBLASおよびLAPACK。 x86系では最速と思われる。FFTなどもはいっている。https://software.intel.com/ en-us/mkl • GPU向けBLAS, LAPACK •b oat16はOpenBLAS, cuBLASが対応 !! 標準化されてない !! • MAGMAプロジェクトはCUDA, OpenCLなどGPUやアクセラレータ向けBLAS, LAPACKを開発している。 fl • https://icl.cs.utk.edu/magma/

• https://software.intel.com/en-us/mkl
• https://icl.cs.utk.edu/magma/

現在のコンピュータ

Top500 速いコンピュータのリスト • 1位はFrontier@US 1.2 Exa Flops = 一秒で10^18 回演算で きる 史上初Exa越え;2022/5 • 2位は富岳@JP (2021/11まで1 位) 440 PFlops • 3位はLUMI@FI • 4位はLeonardo@IT • 5位はSummit＠US • 神威・太湖之光@CNは7位 (2017/11まで1位) 入れ替わり激しい

binary64(倍精度)の計算スピードの比較 縦軸はFlops (一秒に何回計算できるか）

CPU黒歴史 • 遅いCPUはすべて市場から撤退させられる • i860, Itaniumnなど • 高速なCPUは多少使いにくくても人間が頑張る • NVIDIAのGPUなど • 理不尽に市場から撤退させられることもある • 政治的(V30、PEZYなど)、商業的、技術的問題

ボトルネック概念図 wikipediaより

プログラム全体の処理スピードを決める箇所 =ボトルネック

主な２つの コンピュータのボトルネック • 計算能力が高い/低い CPUの速度が速い/遅い • 計算バンド幅 • メモリの転送能力が高い/低い • メモリバンド幅 (フォン・ノイマン型ボトルネック) • 他にも • HDD/SSDの読み書き速度、SATA転送速度、レイテンシ(遅れ)=1bit欲しいと命令してからバスにの るまでの時間などなどなど…あらゆるところに存在する

2大コンピュータのボトルネック と高速化戦略 • やりたいことの特性を見極める • 計算バンド幅 < メモリバンド幅 • 計算能力が低く、メモリ転送が速い • 計算結果をなるべくメモリに残す • メモリバンド幅 < 計算バンド幅 • 計算能力が高く、メモリ転送は遅い • 少量のデータをなるべく使い回す

結局行列-行列積しか 高速化できない... • GEMM（行列-行列積）の高速化のみが主要な焦点。 • 行列-行列積はデータが使い回し＋演算量が多いので、高速化できる。 • 近い将来、メモリの大幅な高速化は困難。メモリは充放電器なので遅い。さらに光速 度も遅い。 • メモリ速度がボトルネックになるならメモリが高速なマシンを買うしかない。 • クロック速度の向上は期待できず、一方で演算器の微細化には可能性 • 計算機はまだ高速化する。さらにgemmに有利 • 全てのプログラムはGEMMを使用するように書き直すべし。 • メモリの転送は最小限に抑えること。むしろ計算の無駄を許容

なぜ行列積か • A, B, Cは n x n の行列として、積を計算する •C←A*B+ C 3回の浮動小数点演算 2n • • A * B : n3回積 ,n2(n-1)回和 • (A * B ) + C : n2回和 2回のデータの読み書き 4n • • A, B, C：3n2読み • ( A * B ) + C : n2書き • n→ でデータ読み書きより演算がボトルネックとなる。 2 / (2n3) * 8 -> 16/n 4n •

基本はblock化 SIMD+自明なマルチコア https://malithjayaweera.com/2020/07/blocked-matrix-multiplication/ B. Kargstrom et al., ACM Trans. Math. Soft., 24(3):268‒302, 1998. https://dl.acm.org/doi/pdf/10.1145/292395.292412 K. Goto et al., ACM Trans. Math. Soft. 34(3):1‒25, 2008 fl https://www.cs.utexas.edu/~ ame/pubs/GotoTOMS̲revision.pdf

• https://dl.acm.org/doi/pdf/10.1145/292395.292412
• https://www.cs.utexas.edu/~flame/pubs/GotoTOMS_revision.pdf
• https://malithjayaweera.com/2020/07/blocked-matrix-multiplication/

ChatGPT4先生に聞いてみた • ブロック化が基本 • 続けて、ブロック化はキャッシュ利用 の効率化を図るためのものであり、 キャッシュミスを減らすことが期待で きます。ただし、これは一例であり、実 際の最適化はハードウェア、特にCPU アーキテクチャの詳細に依存します。ま た、実際のパフォーマンスを向上させる ためには、並列化（マルチスレッディン グ）、SIMD命令の使用、適切なメモリ アライメントなど、さらなる最適化手 法を検討する必要があります。 と、ちゃんと教えてくれた ※こういうプログラムネタはすぐ検証でき るからChatGPTでも信頼できる(!?)

Block化は最も重要だが... • single core, single Blockのみ化をやってもあまり高速化しないように見える • 初心者にとって鬼門 • SIMD化もメモリアクセスを高速にするために重要。 • Kargstromの時代にはSIMD, マルチコアも無かった。 • パフォーマンスチューニングは、CPUやメモリのアーキテクチャを理解する必 要がある。 • ボトルネックの構造が分かりにくい。 • CPUのカウンタを使うのが面倒、専用ツールが必要か？

高速化という獣道への入り口 • 40行で結構高速なgemmを書くチュートリアル • https://en.algorithmica.org/hpc/algorithms/matmul/ • 100行でOpenBLASなみのコードを書くチュートリアル • https://cs.stanford.edu/people/shadjis/blas.html • CUDAのgemmを書くチュートリアル • https://siboehm.com/articles/22/CUDA-MMM • OpenCLでgemmを書くチュートリアル • https://cnugteren.github.io/tutorial/pages/page1.html 14:24

• https://en.algorithmica.org/hpc/algorithms/matmul/
• https://cs.stanford.edu/people/shadjis/blas.html
• https://siboehm.com/articles/22/CUDA-MMM
• https://cnugteren.github.io/tutorial/pages/page1.html

参考図書 • BLAS, LAPACKチュートリアル パート1 (基礎編) BLAS, LAPACKチュートリアル パート2 (GPU編) • http://nakatamaho.riken.jp/blas̲lapack̲tutorial̲part1.pdf • http://nakatamaho.riken.jp/blas̲lapack̲tutorial̲part2.pdf • LAPACK/BLAS入門 幸谷智紀 (https://www.morikita.co.jp/books/book/ 2226) • Matrix Computations, Gene H. Golub and Charles F. Van Loan • http://web.mit.edu/ehliu/Public/sclark/Golub%20G.H.,%20Van%20Loan%20C.F.-%20Matrix%20Computations.pdf • Accuracy and Stability of Numerical Algorithms, Nicholas J. Higham • LAPACK Working NOTE : 実装上の詳細ノートたちがある • http://www.netlib.org/lapack/lawns/index.html • 数値計算の常識, 伊理 正夫, 藤野 和建

• http://nakatamaho.riken.jp/blas_lapack_tutorial_part2.pdf
• http://www.netlib.org/lapack/lawns/index.html