--- title: 第7回 配信講義　計算科学技術特論A （2023） tags: author: [R-CCS　計算科学研究推進室](https://www.docswell.com/user/2300203199) site: [Docswell](https://www.docswell.com/) thumbnail: https://bcdn.docswell.com/page/5JGXGN1Q7L.jpg?width=480 description: 第7回 6月1日　線形代数演算ライブラリ BLAS と LAPACK の基礎と実践２ 二回目は、講演者は任意精度のBLAS, LAPACKを開発しており、それについても利 用方法、応用について紹介する。また、大きなライブラリを開発するときにどの ような点に注意しているかも解説する published: May 29, 23 canonical: https://www.docswell.com/s/2300203199/5VV61J-2023-05-29-130040 --- # Page. 1 ![Page Image](https://bcdn.docswell.com/page/5JGXGN1Q7L.jpg) 線形代数演算ライブラリBLASとLAPACKの 基礎と実践 (II) MPLAPACKについて (高精度BLAS, LAPACK) 中田 真秀 理化学研究所 開拓研究本部 柚木研究室 2023-06-01 計算科学技術特論A # Page. 2 ![Page Image](https://bcdn.docswell.com/page/47Q8XQDWEP.jpg) BLAS, LAPACK高精度計算編 講義内容 1. MPLAPACKの開発に影響を与えた考え方:富豪的プログラミング 2. 私は如何にして心配するのをやめて高精度LAPACKを作るようになったの か (半正定値計画法) 3. コンピュータにおける(高精度)浮動小数点数演算 4. MPLAPACKの誕生 (旧名MPACK) 5. 8年の停滞のあとまさかのブレイクスルー！MPLAPACK 2.0.1のリリース ＋ベンチマーク結果 6. LAPACK (MPLAPACK)のテストは何をしてるのか。 7. MPLAPACKを使ってみる 8. 巨大なライブラリの書き直しで何を学んだか。 # Page. 3 ![Page Image](https://bcdn.docswell.com/page/KE426NZ171.jpg) MPLAPACK開発に 影響を与えた考え方 # Page. 4 ![Page Image](https://bcdn.docswell.com/page/L71MV5R5JG.jpg) 富豪的プログラミング • http://www.pitecan.com/fugo.html • 「メモリや実行効率を気にしないでお気楽にプログラムを作る」 • 「効率を重視したプログラムは作るのが大変ですし、 ちゃんと動かすには デバッグも大変です。 富豪的プログラミングでは一番単純で短いアルゴリ ズムを使います。」 • https://web.archive.org/web/20030902021235/http:// www.pitecan.com/fugo.html • 初出2003年9月頃？ 増井俊之による • 知ったのは2006年頃か # Page. 5 ![Page Image](https://bcdn.docswell.com/page/G7W3L51W72.jpg) コロンブスの卵 • 自明なことをやってみたかった • ぶっちゃけLAPACKをC++に変えただけやん • ちがうちがうちがう 数値線形代数に新しい価値を加える # Page. 6 ![Page Image](https://bcdn.docswell.com/page/4JZVMNP273.jpg) 私は如何にして心配するのを やめて高精度LAPACKを作 るようになったのか # Page. 7 ![Page Image](https://bcdn.docswell.com/page/YE65N846JV.jpg) 半正定値計画法の紹介 primalとdualという等価な問題がある。 X, Y は半正定値(固有値が0以上)の行列。半正定値を保ったまま線 形関数の最小値を行う 最適化条件は、XY=0 (X, Y少なくともどちらかは条件数が発散) 特徴: primal は最適解の上界、dualは最適解の下界を出し、最適になると一致 Max-cut (NP-完全の問題)に対して0.878近似を与える(Williamson-Goeman) 内点法で効率的に解ける -> 藤澤らのSDPAが実装では最高速 # Page. 8 ![Page Image](https://bcdn.docswell.com/page/GE5Y59Q5E4.jpg) 量子化学から半正定値計画問題へ • 量子化学の基底状態のエネルギーなどは縮約密度行列を 変分すれば求まる。Eg = min Γ2∈N−rep. ̂ 2 HΓ • 中田ら J. Chem. Phys 2001(https://  aip.scitation.org/doi/10.1063/1.1360199) で、 これが半正定値計画問題に帰着できて、実際にSDPAとい うソルバで説いて、分子、原子でかなり良い結果を得た。 # Page. 9 ![Page Image](https://bcdn.docswell.com/page/972GN2PG7R.jpg) 色々悩んだ挙げ句(2000年代前半) •量子化学的には有効数字10桁くらいほしい。あと１­２桁く らい足りない •最適解付近で行列の条件数が発散するため、解を精度 良く求めるのは難しい。 多倍長精度計算が必要 # Page. 10 ![Page Image](https://bcdn.docswell.com/page/DJYVPY5GJM.jpg) 色々悩んだ挙げ句(2000年代前半) •量子化学的には有効数字10桁くらいほしい。あと１­２桁くらい足りない •最適解付近で行列の条件数が発散するため、解を精度良く求めるのは難しい。 •半正定値計画ソルバーSDPA(藤澤ら)を使っていた。 •BLAS, LAPACKを沢山使っていた。プログラムはC++で書かれており、かなり複雑 •初期バージョンはmesarchを使っておりBLAS, LAPACKを用いていなかった。 •精度は８桁で十分とみな思っていた。工場の在庫が1万個と1万1個に差はない。 •精度保証は2000年代前半にはあまり発達しておらず •LU分解など部品のいくつかは精度保証はできたが、全部やるのは途方もなく感じた。 •多倍長精度による計算は時間かかるが、それをやってみるかなと。 •すでにBNCpackという幸谷 智紀先生による多倍長精度の線形代数ライブラリがあった。 •CでGMP/MPFRを直に叩いていた。 •SDPAはC++で複雑。BNCpackを使うように書き直すのは複雑だしバグが混入しそうでイヤ # Page. 11 ![Page Image](https://bcdn.docswell.com/page/V7NM5PN8E8.jpg) 高精度半正定値計画ソルバ SDPA-GMPの開発を行った • 富豪的プログラミング: C++でdoubleの感覚で多倍長精度計 算をすべし。 • SDPAはC++で書かれていた。 • 確実に解があるといいたい。 • 細くても(遅くても)切れる(高い精度の)プロ用の道具がほしい。 # Page. 12 ![Page Image](https://bcdn.docswell.com/page/YJ9NY3RXJ3.jpg) 日本人的に陥りやすい罠(私見) • 富豪的プログラミングは中田に衝撃を与えた。 • 日本人は超絶技巧が好き • 発展的な、複雑な理論、計算手法開発。 • 日本は常にリソースが不足している • 1bitでも少なく、1clockでも速く、1pJでも消費電力が少なくしたい。 • パラダイムシフトや、コロンブス卵的研究は多数派ではない。 • BLAS, LAPACKのようなAPIを定められない • BLAS, LAPACKは科学に大きな影響を与えたプログラムの一つ。 • https://www.nature.com/articles/d41586-021-00075-2 • 線形代数は自分が今思っているよりはるかに重要なんじゃないか？と思い始めた。 • 人類の歴史に古くから存在する。 # Page. 13 ![Page Image](https://bcdn.docswell.com/page/GJ8K6MW9ED.jpg) MPLAPACK前史 • BLAS, LAPACKをなるべくそのまま移植することにした。 • SDPAに使われている関数を実装していった。 • コレスキー分解、対称行列の固有値を求めるまでok • dsyevまでたどり着くまで… 50個くらい移植 • https://github.com/nakatamaho/sdpa-gmp/tree/master/ • GNUやオープンソースを信奉してたのですぐリリース。 # Page. 14 ![Page Image](https://bcdn.docswell.com/page/LJLZ93NM7R.jpg) 結果 • J. Chem. Phys. 2008 • Hubbard modelでU/t-> でちゃんと値が出てきた。 • SDPの問題で不安定、最適解がわかってない問題へも応用 できた。 • https://ieeexplore.ieee.org/document/5612693 • 81回も引用していただき感謝 (2023-05-24;google scholar) # Page. 15 ![Page Image](https://bcdn.docswell.com/page/47MVZDX67W.jpg) SDPの収束の様子 Gpp250 内点がないill conditionな問題 100000 SDPA-GMP SDPA-QD SDPA-DD SDPA 1 maxG11 (Slater条件を満足) 100000 SDPA-GMP SDPA-QD SDPA-DD SDPA 1 1e-05 1e-05 1e-10 µ µ 1e-10 1e-15 1e-20 1e-15 1e-25 1e-20 1e-25 1e-30 1e-35 0 20 40 60 80 # of iterations 100 120 0 5 10 15 20 25 # of iterations 30 35 40 # Page. 16 ![Page Image](https://bcdn.docswell.com/page/P7R264NNJ9.jpg)   高精度半正定値計画法の応用例 • 三次元イジングモデル; 共形場理論を用いた結果が厳密解と等しい予想がある。 • https://arxiv.org/pdf/1406.4858.pdf : SDPA-GMPを用いている https://arxiv.org/pdf/1502.02033.pdf 独自ソルバ • # Page. 17 ![Page Image](https://bcdn.docswell.com/page/PJXDM2N8EX.jpg) SDPの精度保証 J. Chem. Phys. 2001の我々の結果、誤差が大きい場合があったことも検証してくれた # Page. 18 ![Page Image](https://bcdn.docswell.com/page/3JKXGMNLED.jpg) コンピュータにおける (高精度)浮動小数点数演算 # Page. 19 ![Page Image](https://bcdn.docswell.com/page/LE35NYV6J5.jpg) 浮動小数点演算と多倍長精度演算ライブラリ • binary64(倍精度) Binary128 • 精度が必要な問題に対しては、単に仮数部を増やして力づくで解く方法を採用。 • GMP: 任意精度ライブラリ、高速。 • MPFR: GMP+IEEE754ライクな丸めつき、つまりより誠実な演算をする。 • MPC: MPFRに基づく複素数演算ライブラリ。初等関数も実装されている。 • DD, QD: 倍精度変数を複数もつことでほぼ4, 8倍精度を実現。いわゆるIBM方式 ユーザのニーズに合わせて、高精度演算ライブラリはいくつも存在する # Page. 20 ![Page Image](https://bcdn.docswell.com/page/8EDDW58MEG.jpg)  binary64 • binary64(倍精度) 10進数16桁程度 • 64bit=8bytesで１つの実数を表現する • いわゆる「倍精度」 • もっともポピュラーな浮動小数点数。ほぼすべてのプロ セッサでハードウェアのサポートがあるため高速。 # Page. 21 ![Page Image](https://bcdn.docswell.com/page/V7PX9G8QE8.jpg) ̲Float64x • ̲Float64x : 80bit浮動小数点数。IEEE 754 2008で規格外になった • (拡張倍精度) 10進数19桁程度 • 80bit=8bytesで１つの実数を表現する; 128bitで格納するか? 80bitで格納するか? • いわゆる「拡張倍精度」exponentのrangeがbinary64より大きい。 • x87/68881でサポートされいる/いた浮動小数点数。もはやSSEやAVXが主流なのでほぼ使われ てない。後方互換のためだけに存在する。遅い。 • CではTS 18661-3で(今更)初めてサポートができた。C++のサポートはなし(GNU拡張はあり) • 32bitのx86ではx87を使っていてlong doubleが̲Float64であったこともあった。 # Page. 22 ![Page Image](https://bcdn.docswell.com/page/2JVLM6NP7Q.jpg) binary128 • binary128はIEEE 754 2008で定義された • 10進数33桁ほどの精度 • いわゆる四倍精度。 • ハードウェアサポートはIBM z13程度。他はすべてソフト ウェア実装。非常に遅くなる。 # Page. 23 ![Page Image](https://bcdn.docswell.com/page/5EGXGNKQJL.jpg) binary128の状況は混沌としている • binary128はIEEE 754-2008で規格化され、gcc4.3から̲̲ oat128という形で使えている(GNU拡張) • Cでは、TS 18661-3でbinary128の実装である、̲Float128が定義された。 • C++ではまだ何も定まっていない(!) が、GNU拡張で̲Float128が使える場合もある。 • libquadmathやglibcのサポートは必須。 • Gcc10.xでもaarch64のように̲̲ oat128/libquadmathが使えない環境はある。わざと潰してある。 • 少なくとも6種類存在。 • long doubleのみ存在そしてこれはbinary128。̲ oat128,̲Float128は使えない(CentOS7, 富岳AArch64) • ̲Float128はgcc/glibcでサポートされていて、long doubleと同じ(CentOS8 AArch64) • ̲Float128はgccで̲̲ oat128とlibquadmathとして存在している(MacOS x86̲64, mingw64) • ̲Float128はglibcでサポートされていて、binary128型は一つしか利用できない(Intel C/C++ 多分バグ) • ̲Float128はgccとglibcでサポートされていて、̲Float128, ̲̲ oat128両方使える(Ubuntu20.04, amd64) • ̲Float128はgccとglibcでサポートされていて、̲Float128のみ使える(Ubuntu22.04, amd64) https://qiita.com/mod̲poppo/items/8a61bdcc44d8afb5caed fl fl fl fl fl は大変参考にした。 # Page. 24 ![Page Image](https://bcdn.docswell.com/page/4JQ8XQNW7P.jpg) • • binary64をふたつ並べて10進32桁ほどの精度を実現 • Knuth-Dekkerの方法を使うと無誤差でbinary64の数同士の足し算、 掛け算ができる。 • IBM方式ともいわれ、IBMのマシンやコンパイラ、PowerMac、 PowerPCのgccのlgong doubleは倍々精度である/あった(obsolated) • Y. HidaらのQD libraryを用いた。   倍々精度 # Page. 25 ![Page Image](https://bcdn.docswell.com/page/K7426NG1J1.jpg) 8倍精度 • binary64を4つ並べて10進64桁ほどの精度を実現 • Knuth-Dekkerの方法を使うと無誤差でbinary64の数同 士の足し算、掛け算ができる。 • Y. HidaらのQD libraryを用いた。 # Page. 26 ![Page Image](https://bcdn.docswell.com/page/LJ1MV5D5EG.jpg) The GNU Multiple Precision Arithmetic Library • 整数の配列を浮動小数点とみなす。 • 速度に主眼を置いており、fractionは64bit単位でないと 増やせない、丸めモードがない、exponentのrangeが固 定という制約がある。 • sqrtはあるが、exp, logなどがない→MPLAPACKでatan # Page. 27 ![Page Image](https://bcdn.docswell.com/page/GJW3L5YWJ2.jpg) The GNU MPFR Library • 整数の配列を浮動小数点とみなす。 • GMPをベースに、IEEE 754ライクな丸めモードを搭載+exponentが 可変 (精度を意図的に低くできる)。真面目に計算するときに重要。 • 初等関数(sin, cos, tan, log, exp, arctan, sqrt...)などがある • 複素数ライブラリのMPCも存在する(MPLAPACKで利用) • GMPと比べると少し遅い。 # Page. 28 ![Page Image](https://bcdn.docswell.com/page/4EZVMNX2J3.jpg) MPLAPACKの誕生 (旧名MPACK) # Page. 29 ![Page Image](https://bcdn.docswell.com/page/Y765N8P6EV.jpg) 線形代数: 現行BLAS, LAPACKではAPI定義、規模、スピー ドという価値のみ。精度という新しい価値を加えたい。 • 50個BLAS, LAPACKの関数を実装していって、どうせならBLAS は全部、LAPACKもいくつかやってみてもいいじゃないか。でき るんじゃないかと思い始めた。 • BLAS, LAPACKは素直で綺麗なコード • FORTRAN77からCに書き直すのは良いことなんではないか？ • 線形代数は歴史を紐解くと長いし、重要度はますます増すばかり。 • 精度という価値はあんまり重要視されてない。線形代数に一つ、貢 献できるかも。 # Page. 30 ![Page Image](https://bcdn.docswell.com/page/G75Y59K574.jpg) MPLAPACK (旧MPACK)へ • 多倍長精度のBLAS, LAPACKはちゃんと作っておこう。 • C++の勉強を兼ねて使ってみた。 • 体力気力がどこまで続くかやってみようと思ってやってみた。BLASはすぐ終わるが、 LAPACKは見えない。完璧を目指すより何かを出すことを目標。 • 開発方針 • C++で書いて、doubleや oatと同じ感覚で使えるようにすること。ただし、C＋＋は 便利なCであること。 • FORTRAN77から脱却し、バグが入りにくいライブラリを。 • なるべく多くの多倍長精度ライブラリを取り込むこと。 • 標準化と、実装の最適化を分けること。 fl • SDPAとはリンクできること。実用例が重要。 # Page. 31 ![Page Image](https://bcdn.docswell.com/page/9J2GN2WGJR.jpg) プログラミングモデル • 関数の名前はPre xを変えた。“R” や “C”を用いた。あとは小文字。 • daxpy -> Raxpy, zgemm -> Cgemm • INTEGER, REAL, COMPLEXという仮想的な型を用いて実装した。 • C++のクラスを使ってREALなどを実現。 • typedefを用いて REAL -> mpf̲class, qd̲real, ̲̲ oat128など実際の型を指定した。 • テンプレートは用いず。 • Cgemm<mpf̲real, mpf̲complex>はかっこ悪い • C++も進化してきたしテンプレート版作成予定 fl fi • 初等関数sin, cos, log, exp…などはあれば使うが、なければdoubleで代用。 # Page. 32 ![Page Image](https://bcdn.docswell.com/page/DEYVPYLGEM.jpg) 実装方法 • BLAS, LAPACKをf2cを用いてCのコードに変換 • 初期の頃は手探り状態で、目で見ながら変換してた。 • Sedをヘビーに使って修正 • C++に変換。 • パーサに手を入れようとしたが殆ど読めず。yacc, lexなどを学んだことな くて、理解できず。 • ひたすら手で直す。 • ランダムな値を入れてBLAS, LAPACKのルーチンと比較。差が小さければok # Page. 33 ![Page Image](https://bcdn.docswell.com/page/VJNM5P4878.jpg) FORTRANとC++の違いの吸収方法 • 配列は0から始まる(C++) 1から始まる (FORTRAN) • 二次元配列がC++にはないので、leading dimensionを陽に用いて一次元配列に開いた。どの ルーチンにも必ずleading dimensionの記述はある。 • 一次元配列はA(I) -> a(i-1)と1引くことにした。 • 二次元配列はa(i,j) -> a[(i-1) + (j-1)*lda]と変換した。 • Do loop は 1から始まり、Nで終わることが多い。 • Do loopは+1/-1されるかが分からない場合がある。以下のイディオムが役に立った。 • for (i=p; inc>0 ? i<=n : i>=n ; i=i+inc) { • COMMON, FORMATもWRITEもBLAS, LAPACKには存在しなかった。 • これはいいこと。I/Oを使わないので楽 • ただしテストプログラムでFORMAT, WRITEは多用している。 • GOTOは残した。しかし手でcontinue, breakなどにしたこともあった。 • 初期の変換は目で行っていて「自然な」ループ (for (i=0; i<n; i++))にしていたこともある。 # Page. 34 ![Page Image](https://bcdn.docswell.com/page/YE9NY3QXE3.jpg) dgemmの変換を見てみよう # Page. 35 ![Page Image](https://bcdn.docswell.com/page/GE8K6MG97D.jpg) BLAS: dgemmの抜粋 # Page. 36 ![Page Image](https://bcdn.docswell.com/page/LELZ93GMJR.jpg) dgemmをf2cに通した場合 # Page. 37 ![Page Image](https://bcdn.docswell.com/page/4JMVZDQ6JW.jpg) 多倍長精度Rgemm # Page. 38 ![Page Image](https://bcdn.docswell.com/page/PJR2648NE9.jpg) Rgemm(旧版) # Page. 39 ![Page Image](https://bcdn.docswell.com/page/PEXDM2887X.jpg) 「正しい値がでているか」 が難しい。 • どうしたら正しい値が出てるか、のチェックは本質的に難しい。 • 「品質保証」Quality assurance (QA) と呼ばれる • まず、MPFR版とオリジナルのBLASに、ランダムな値を色々入れてみて、コード を隈なく通っているかどうかを確認しつつ、十分結果が近ければMPFR版は正し い、とした。 • GMP版、QD版....はMPFR版と比較してより厳しい閾値でチェックした。 • これは本質的にはQAではないが、BLAS, LAPACKは正しいことを根拠に、正 しいことにした。研究での利用には耐えているので、いいと思う。 • 時間もものすごくかかる。Rgemmだと、n, m, k 以外にlda, ldb, ldcもあるの で6重ループになる。さらにtransposeで4通り…と積み重なる。 # Page. 40 ![Page Image](https://bcdn.docswell.com/page/3EKXGMKL7D.jpg) Rgemmのチェックプログラム # Page. 41 ![Page Image](https://bcdn.docswell.com/page/L735NYZ6E5.jpg) MPACK 0.8.0 (2012-12-25) • 多倍長精度のBLAS, LAPACK • API提供を目的 • MBLASはすべて実装、MPALACPKは100ルーチン実装。 • 実対称、エルミート固有値問題、LU, コレスキー分解、逆行列、条件数推定対応 • Linux/BSD/Mac/Win対応 • GMP, MPFR, binary128, quad-double, double-double, double, long double 対応 • Rgemm (double-double) はCUDAでC2050対応 • Intel Xeon Phi対応 • C++で書かれていて、double/ fl • 2-条項BSDライセンス oatのようにプログラム可能。 # Page. 42 ![Page Image](https://bcdn.docswell.com/page/87DDW5RM7G.jpg) MBLAS 76ルーチン 古いリリース # Page. 43 ![Page Image](https://bcdn.docswell.com/page/VJPX9GWQ78.jpg)  MLAPACK 100ルーチン 古いリリース • # Page. 44 ![Page Image](https://bcdn.docswell.com/page/2EVLM68PJQ.jpg) 8年の停滞のあと まさかのブレイクスルー MPLAPACK 2.0.1のリリース! # Page. 45 ![Page Image](https://bcdn.docswell.com/page/57GXGN5QEL.jpg) • 燃え尽きた 2012を最後にリリースが できなくなった。 • 高精度な非対称行列の固有値問題、特異値分解のルーチンはコミュニティから多数要望があった が、あまりの複雑さに音を上げてしまった。 • 品質保証をどうしたらいいかわからない • 100個書いたが、これ以上書くモチベーションが無くなった。 • ランダムな値を入れて結果をLAPACKと比較して、という戦略だけでは、品質保証は難しい。 • Rsteqrなど、QR法で対角化する場合、かなり場合分けが多い。隈なくコードをチェック できているか？ • 何をしているのか理解するのに時間がかかる、または分からない場合。 • 自分の興味が移ってきた。 • 機械学習の勃興により、久しぶりに量子化学をやり始めたら面白かった。 • https://nakatamaho.riken.jp/pubchemqc.riken.jp/ NeurIPSでみんなに競っていただきま • 2021,2022と機械学習での予想コンテストKDD, した。 • https://ogb.stanford.edu/docs/lsc/pcqm4mv2/ # Page. 46 ![Page Image](https://bcdn.docswell.com/page/4EQ8XQZWJP.jpg) 2012に残された課題 • mpackだと他のプログラムパッケージと名前がかぶっていて、mplapackに直したかった。 • F2cで変換したコードをもとに、手でC++に変換する、には限界があった。 • Sedを使って直しても、可読性が低い。一日に変換できる量に限界があり、また、集中度によりバ グがどう混入してるかわからない。 • Cで書かれたパーサは書き直せなかった。出力されたコードの可読性は重要だとわかった。 • 実際の半正定値計画問題が解けるとかなり満足して、モチベーションが下がった。 • コード変換より、品質保証の手間が大きすぎた。 • 600個程度あるLAPACKのルーチンのうち100個程度は変換できた。隈なくコードをなめるよう なテストケースを書く必要がある。特にあるルーチンが理解できない場合、品質保証が難しい。 • LAPACKのTESTING以下を移植するのは良いとはわかっていたが、FORTRAN+f2cでの format文を手で直すのは困難。 • C++の変化が意外と激しくて、すぐエラーが出てくる。 • 開発環境の維持にコストが掛かりすぎる。 • 仮想マシンも維持が結構めんどくさい。リソースも多くとられる。 • OSやアーキテクチャはたくさんある。 • IEEE 754-2008でbinary128が入ったが、コンパイラ、OSの状況は未だ混沌としている。 # Page. 47 ![Page Image](https://bcdn.docswell.com/page/KJ426N31E1.jpg) 特異値分解、非対称行列の固有値問題ソルバ は特に非常に複雑 # Page. 48 ![Page Image](https://bcdn.docswell.com/page/LE1MV5157G.jpg) dgesvd(特異値分解)のcall graph dgeev(非対称行列の対角化)のcall graph dgesvdだけでも3000行以上あるのに dgesvdより長いし、依存する関数もより多い さらに30個以上の品質保証が必要。 f2cと 気力だけでは 品質保証が 不可能 # Page. 49 ![Page Image](https://bcdn.docswell.com/page/GEW3L58WE2.jpg) どうやってデバッグおよび 品質保証するか? # Page. 50 ![Page Image](https://bcdn.docswell.com/page/47ZVMN82E3.jpg) LAPACK/TESTINGのC++への移行 f2cでは事実上不可能 • BLASは代数的な演算だけなので、大雑把にBLASと MPBLASの値がだいたい合っていれば、実装がうまく行っ ているとは言える。 • LAPACKはどうすればいい? • LAPACKのTESTING/で良いテストを行っている。ので これに乗っかりたい。 • しかし、FORTRANのwrite, formatを手で直すのも事 実上不可能。 # Page. 51 ![Page Image](https://bcdn.docswell.com/page/YJ65N8K67V.jpg) f2c変換されたチェックルーチンの format文の例 • LAPACKのTESTINGで良いテストを。 • # Page. 52 ![Page Image](https://bcdn.docswell.com/page/GJ5Y59P5J4.jpg) 課題の克服 • 中田からは何もしてない。待っているだけであった。 • Fable: FORTRANからC++へ、可読性のあるコンバータの登場。 • Pythonで書かれているため、パーサへの手入れが簡単。 • 2012年頃はpython知らなかった。2015年ころから使い始めた。 • FEM (Fortran EMulator) の搭載 : write, formatがほぼそのまま使えた。 • Clang-format: C++の整形ツール。sedに渡しやすい。 • Docker: OSやアーキテクチャを手軽に試せた。 • Github: webベースの存在がありがたい。 • デパス: 鈍い頭痛、肩こりに20年位悩まされていたのが消えた。 • コロナ: 通勤しなくて良くなった。 # Page. 53 ![Page Image](https://bcdn.docswell.com/page/LE35NY35J5.jpg) fableの登場(2012) • FORTRANからC++へのコードのコンバーター • コードを生成するときの方針としては: • 厳密に同じ動きをさせる(f2c) • 人間が読めて、たいていそのままで動く。さらに直す場合もな るべく手間が最小限(fable) … dsyevまではそのままで動くと論 文にはあった。 http://cci.lbl.gov/fable/ # Page. 54 ![Page Image](https://bcdn.docswell.com/page/8EDDW54YEG.jpg) fableの登場(2012) • Fotran -> C++ への変換部分はすべてpythonで書かれているため、可読性が 高い。 • 気に入らない変換があれば、直接手に入れられる。ライブラリも多く、正規 表現も使える。 • FEM (fortran emulation library) • I/Oやintrinsic(内在関数)はコンバートではなくライブラリに。 • LAPACKのTESTINGでwrite/format文を処理することが可能になった • ←まさかのC++のコード! # Page. 55 ![Page Image](https://bcdn.docswell.com/page/V7PX9GM2E8.jpg) Dockerの登場 • 仮想技術を使ったLinuxコンテナを作れる。 • OSなどの違いを吸収するのは意外と面倒。 • 仮想マシンだと、管理が必須。リソースも多く消費する。 • Dockerだと軽量なので気軽に環境の構築、破棄が可能。 • 気が向いたときに環境構築だけに時間を費やして時間切れ、はもう嫌だ。 • 他の人と成功、失敗を共有しやすい。再現度が高くなる • pythonの環境は無茶苦茶なので、仮想マシンを割り当てるのが妥当。 • Binary128のサポートは過渡期にあり、OSにより随分違うため試せる環境が必 要。 • PowerPC64le (IBM Power9?)を普通のインテルマシンで試せる # Page. 56 ![Page Image](https://bcdn.docswell.com/page/2JVLM69X7Q.jpg) 8年の停滞とv1.0.0への道のり • FableでLAPACKのTESTINGをC++に移植できるかも!! • もしかしたら今回は非対称行列の固有値問題や、特異値分解 を実装できるかもしれないと思い、2021-03くらいから開 発を再開してみた。 • https://github.com/nakatamaho/mplapack/tree/ v1.0.0 開発の再開 # Page. 57 ![Page Image](https://bcdn.docswell.com/page/5EGXGNZRJL.jpg) 開発方針 • Fableに手を加えてなるべく自動変換を行う • BLAS, LAPACKともに自動変換で対処 • できないところはparserに手を加えたり、sedで対処 • 最後に手で直す # Page. 58 ![Page Image](https://bcdn.docswell.com/page/4JQ8XQLY7P.jpg) MPLAPACK特有の手入れ • dconjg, dsqrtなどのトークンを後にsedで変換できるように文字列挿入。 • 文字列をFEM特有からcharの配列に直した。 • カッコをつけたり外したり、は正規表現で対処した部分もある。 • 複素数の宣言 COMPLEX a = COMPLEX(0, 1);のように初期化対応。 • 二次元配列、多次元配列は、fableでは関数になってた。これらはすべて一次元配列に直し た。 • leading dimensionは決め打ち(!) • a(i,j) -> a[(i-1) + (j -1)* lda] にパーサの段階で変換。 • DO LOOPはforにする。 • Incrementが正、負が分かる場合は+1, -1 に。 • Incrementが変数の場合はソースを読んで手で直す。 • Incrementが正負とりうる場合は、イディオムを手で挿入(そんなに多くない)。 # Page. 59 ![Page Image](https://bcdn.docswell.com/page/K7426NDZJ1.jpg) MPBLASはほぼ自動 生成できるようになった。 • BLASは、C++、多倍長精度化も一箇所を除き自動で生 成できるようになった! • conj(a) * a は実数だがこれは自動では変換できなかっ たのでパッチを当てている。 • ヘッダも自動生成するようにした。 # Page. 60 ![Page Image](https://bcdn.docswell.com/page/LJ1MV5ZDEG.jpg) MPLAPACKでも、かなりの 部分が自動的に変換できた • さすがにMPLAPACKの全ては自動変換できなかった • それでもLAPACK 3.9.1 の約1000個ルーチンのうち、約350個は自動変換 できた(無修正でコンパイルできた) • その他も軽微な修正で出来た • 2021-04-04から2021-04-20までですべてのルーチンの手変換が終わっ た。 • mpackのころは、何ヶ月もかかった。しかもできたルーチンは挙動不審。 • 混合精度は今回はターゲットとせず(むしろしたほうがいいか)。 • XBLAS利用による解の精度向上もターゲットとせず。 # Page. 61 ![Page Image](https://bcdn.docswell.com/page/GJW3L598J2.jpg) MPLAPACKの自動変換の 品質保証 • 2012年までの品質保証はまず行ってみた。 • LAPACKとMPLAPACKに同時にランダムに生成した 行列を入れて、その値を比較するということを今回も 行った。 • ほぼ無修正でテストにパスした。 • これまでのルーチンは実装は変わったが、使えること がわかった。 # Page. 62 ![Page Image](https://bcdn.docswell.com/page/4EZVMN99J3.jpg) LAPACKのTESTINGの移植 • Fable + FEMで、format文はほぼそのまま使えることが判明 したため、移植を試みた。 • 2021-05-12にとりあえず、 • LIN/xlintst{R,C}̲mpfr : MPFR版の実、複素線形方程式 ルーチンのテスター • EIG/xlintst{R,C}̲mpfr : MPFR版の実、複素固有値問題 程式ルーチンのテスター のビルドに成功。 # Page. 63 ![Page Image](https://bcdn.docswell.com/page/Y765N8KDEV.jpg) LAPACKのTESTINGの移植 • ほぼ自動変換したプログラムのどこにバグが存在するか? • チェッカープログラムの移植不具合により適切に検出で きてない場合がほとんど。 • それ以外は • maxlocの仕様がわかってなかった • Do loopで増減が両方の場合が存在(イディオムで対処) # Page. 64 ![Page Image](https://bcdn.docswell.com/page/G75Y59P874.jpg) 難航したバグ探し • アリエナイTESTINGのバグ • まさかのFortranで配列の添え字に0を使っていた(2か所) • まさかの走ってないtestがあった。 • githubでissue立てた→貢献者リストに名前が載った • https://github.com/Reference-LAPACK/lapack/issues/563 • テストのバグ; xlaenvの初期化忘れ • いくつかのルーチンで、変数をconst受け取るが、呼び出した先のルーチンで変数に代入してる場合 があった(未issue)。 • WindowsはMINGW64でもIP64 ILP64であった。 • 複素数の割り算は難しい: GMPの内部は変わらず(だが関数はlong int受けのみ)→MPFRの fl fl exponentの幅が変わる→MPCの割り算はIP64を前提→over ow, under ow付近での誤差 が大きくなる(未issue) # Page. 65 ![Page Image](https://bcdn.docswell.com/page/9J2GN26VJR.jpg) LAPACKのTESTINGの移植 • 対応すべきxxx.inの数は28個ある。 • feasibleな労力でできた。 • 2年かけて、非対称固有値問題、特異値問題が解けるよ うになった。 # Page. 66 ![Page Image](https://bcdn.docswell.com/page/DEYVPY9QEM.jpg) 昔々:MPACK 0.8.0 (2012-12-25) • 多倍長精度のBLAS, LAPACK • API提供を目的 • MBLASはすべて実装、MPALACPKは100ルーチン実装。 • 実対称、エルミート固有値問題、LU, コレスキー分解、逆行列、条件数推定対応 • Linux/BSD/Mac/Win対応 • GMP, MPFR, binary128, quad-double, double-double, double, long double 対応 • Rgemm (double-double) はCUDAでC2050対応 • Intel Xeon Phi対応 • C++で書かれていて、double/ fl • 2-条項BSDライセンス oatのようにプログラム可能。 # Page. 67 ![Page Image](https://bcdn.docswell.com/page/VJNM5PL278.jpg) MPLAPACK 2.0.1(2022-09-12) • 多倍長精度のBLAS, LAPACK • API提供を目的 • MPLAPACKほぼ実装 実+複素 • 実対称、エルミート固有値問題、LU, コレスキー分解、特異値 • 一般化固有値問題、最小二乗法、逆行列、条件数推定対応 • 未対応: 混合精度 • Linux/Mac/Win対応 • GMP, MPFR, binary128 (̲Float128), quad-double, double-double, double, ̲Float64x 対応 • Rgemm (double-double) はCUDAでA100, V100対応 • C++で書かれていて、double/ fl • 2-条項BSDライセンス oatのようにプログラム可能。 # Page. 68 ![Page Image](https://bcdn.docswell.com/page/YE9NY34DE3.jpg) LAPACK(MPLAPACK)の テストは何をやっているか # Page. 69 ![Page Image](https://bcdn.docswell.com/page/GE8K6MQ57D.jpg) LAPACKのテスト基本方針 • 正解を先に作る。 • Ax=bのサブルーチンをテストするなら、A, xを作って bをまず得ておく • そこから問題を解く。 • ちゃんと解けたらok 解けなかったら失敗 • Aを様々な条件に変化させる。悪条件、良い条件、オー バーフローアンダーフロー付近までスケールなど。 # Page. 70 ![Page Image](https://bcdn.docswell.com/page/LELZ93X3JR.jpg) LAPACKはテストも非常に丁寧 • LAPACKビルド時に→のよう にテストをする。 • xlintstd:binary64版線形方 程式品質保証プログラム • dtest.in : どのルーチンをど うてすとするかが書かれてい て、xlintstdに食わせる • DGE: 一次連立方程式 • DGB: バンド行列DGE • DGT: 三角行列DGE • TESTING/LIN/dchkaa.F # Page. 71 ![Page Image](https://bcdn.docswell.com/page/4JMVZDLMJW.jpg) LAPACKはテストも非常に丁寧 • ./EIG/xeigtstd < svd.in • xeigstd:binary64版固有値問 題、特異値分解、など品質保証 プログラム • svd.in : 特異値分解品質保証の ための入力ファイル • このような.inを23(実) +23(複素)計46個チェック する必要がある # Page. 72 ![Page Image](https://bcdn.docswell.com/page/PJR264KLE9.jpg) 特異値分解のテストは何をしてるか • M=UΣV* • 最初に正解を作る • Σに様々な特異値を代入 • 一様ランダム • 広く分布 • U, V*: ランダムに直交 or ユニタリを発生 • Mを作る。そこからもう一度解きなおす • 求まった特異値Σ*はΣと十分近いか? • |I-U*U| が十分小さいか • ....などをチェックする # Page. 73 ![Page Image](https://bcdn.docswell.com/page/PEXDM2567X.jpg) LAPACKテストの問題点 • テストはどのような入力の場合もテストできてるか • 入力は有限個だが...不可能、または保証する手段はある のか? • テスト行列の次元が50x50程度と比較的小さい • 非対称行列、特異値分解なども。 • コンパイラの最適化、実装の最適化をすると、誤差が 大きくなる、またはテストに通らないことがある。 # Page. 74 ![Page Image](https://bcdn.docswell.com/page/3EKXGM6G7D.jpg) LAPACKのTESTINGの移植 • LAPACKは自分の品質保証のためにTESTINGディレクトリが存在する。 • LIN/xlintstd : binary64の実線形方程式ルーチンのテスター • LIN/xlintstz : binary64の複素線形方程式ルーチンのテスター • EIG/xeigtstd : 実固有値問題、特異値問題のテスター • EIG/xeigtstz : 複素固有値問題、特異値問題のテスター • MATGEN/ テスト用行列生成ライブラリ • これらにxxx.inというファイルを食わせて正常終了したらok # Page. 75 ![Page Image](https://bcdn.docswell.com/page/L735NY65E5.jpg) MPLAPACKを使ってみる # Page. 76 ![Page Image](https://bcdn.docswell.com/page/87DDW5VY7G.jpg) 試し方 • https://github.com/nakatamaho/mplapack を見てね Dockerで試せます Ubuntu, CentOS amd64(x86̲64)/intel one API, aarch64, win64 macOSではmacportsかbrewなどで # Page. 77 ![Page Image](https://bcdn.docswell.com/page/VJPX9GV278.jpg) サンプルプログラム • ビルドが終わると、${pre x}/share/examples/mplapackディレクトリ以下に例 が沢山できる。いくつか例をあげると… • 00̲LinearEquations/*getri̲test* • MPFR/GMP/QD...などで任意精度で実、複素数行列逆行列を求める。 • 00̲LinearEquations/Rgetri̲Hilbert̲XXX • MPFR/GMP/QD...などでヒルベルト行列とその逆行列を求める。 • 04̲NonsymmetricEigenproblems/*geev* • MPFR/GMP/QD...などで非対称行列の固有値問題を解く • 05̲SingularValueDecomposition/*gesvd̲test̲* fi • MPFR/GMP/QD...などで特異値問題を解く # Page. 78 ![Page Image](https://bcdn.docswell.com/page/2EVLM6ZXJQ.jpg) Rgemmのベンチマーク dgemmの多倍長精度版はRgemm 実装、テスト、パフォーマンスをみてよう GMP版512bit(153桁): 600MFlops Binary128(33桁)：2000M ops Double-double(32桁): 18GFlops-550GFlops Ryzen 3970X(32cores), NVIDIA A100 fl double-double計算はA100で実用になる # Page. 79 ![Page Image](https://bcdn.docswell.com/page/57GXGN4REL.jpg) A benchmark of Rgemm # Page. 80 ![Page Image](https://bcdn.docswell.com/page/4EQ8XQGYJP.jpg) A benchmark of Rgemm # Page. 81 ![Page Image](https://bcdn.docswell.com/page/KJ426N5ZE1.jpg) double-doubleでのRgemm on A100 ほぼ四倍精度、まさかの500GFlops超え！ # Page. 82 ![Page Image](https://bcdn.docswell.com/page/LE1MV5WD7G.jpg) 他のベンチマーク結果 • https://github.com/nakatamaho/mplapack/tree/ master/benchmark/results/2022 https://github.com/nakatamaho/mplapack/tree/master/benchmark/results/2022 # Page. 83 ![Page Image](https://bcdn.docswell.com/page/GEW3L568E2.jpg) ソフト開発で学んだこと • 中田からは何もしてない。待っているだけであった。 • Fable: FORTRANからC++へ、可読性のあるコンバータの登場。 • Pythonで書かれているため、パーサへの手入れが簡単。 • 2012年頃はpython知らなかった。2015年ころから使い始めた。 • FEM (Fortran EMulator) の搭載 : write, formatがほぼそのまま使えた。 • Clang-format: C++の整形ツール。sedに渡しやすい。 • Docker: OSやアーキテクチャを手軽に試せた。 • Github: webベースの存在がありがたい。 • デパス: 鈍い頭痛、肩こりに20年位悩まされていた。 • コロナ: 通勤しなくて良くなった。 # Page. 84 ![Page Image](https://bcdn.docswell.com/page/47ZVMNY9E3.jpg) ソフト開発で学んだこと+感謝 • 線形代数は人類永遠の課題。 • どこにでも出てくる。 • いつでも何か課題が潤沢にある。理論的？計算機的？HPC的？ • 人間の頭では非線形を理解出来ないのかもしれない。 • なるべく自分を追い詰め過ぎず、適度に甘やかすこと • 燃え尽き症候群となると回復に時間がかかった。工夫が必要。 • 体調管理。デパスが頭痛、肩こりを随分和らげてくれた。 • 道具がないときは作るじゃなくて待つのも良かろう。 • モチベーションを維持できないときは出てくるまで待てば良い。 • 最大限楽をするために最大限努力せよ(これはよく言われる)。 • 姫野先生に拾ってもらった(多謝) • 失業せず何とかなった。8年の停滞がありながらも続けることができた。 • 下司先生に感謝 • 前回の授業や本執筆でモチベーションもらった。 • 書ききれないほどたくさんの人に暖かい言葉や支援をいただいた。