高精度線形代数演算ライブラリ MPACKの開発 中田真秀

ペタ・エクサへ向けた高精度線形代数演算の必要性
.
.
.
.
..
.
高精度線形代数演算ライブラリ MPACK の開発
中田真秀
理化学研究所 情報基盤センター
第 2 回戦略的高性能計算システム開発に関するワークショップ
∼2010/11/27@秋葉原コンベンションホール
riken-logo
.
中田 真秀
.
.
.
高精度線形代数演算ライブラリ MPACK の開発
.
.
ペタ・エクサへ向けた高精度線形代数演算の必要性
. MPACK 0.6.7: Multiple precision version of BLAS and LAPACK
http://mplapack.sourceforge.net/
中田真秀 @ 理化学研究所
MPACK (MBLAS and MLAPACK) BLAS/LAPACK の任意精度版。
参照実装 (API) の提供が第一目標。
最新バージョン 0.6.7 (2010/8/20); 666 個中、100 個の
MLAPACK, 全ての MBLAS をサポート。
高い移植性: Linux/Windows/MacOSX/FreeBSD
多くの多倍長精度演算ライブラリに対応:
MPFR/MPC/GMP/DD/QD/double
C++; よりわかりやすく、より高速、よりプログラマにやさしく。 riken-logo
2-条項 BSD ライセンス; 商用利用、改変、再配付自由。
.
中田 真秀
.
.
.
高精度線形代数演算ライブラリ MPACK の開発
.
.
ペタ・エクサへ向けた高精度線形代数演算の必要性
. 高精度線形代数計算は重要
線形代数は人類の歴史とともにある重要なテーマ。逆にいう
とそれしかできない。
Top500 での Linpack の計算結果の精度が 1 桁くらいしか
ない。
規模が小さくても難しい問題もある。
精度の保証がないと不安。本当に正しい答えか。
精度はどこまで必要かは... 必ずしもよくわからない。
高精度計算線形代数ライブラリがほとんど存在しない。
過去の資産活用をしたい。プログラムは大きく変えたくない。
riken-logo
BLAS/LAPACK に似た高精度計算線形代数ライブラリを作ってみよう
.
中田 真秀
.
.
.
高精度線形代数演算ライブラリ MPACK の開発
.
.
ペタ・エクサへ向けた高精度線形代数演算の必要性
. BLAS、LAPACK とは何?
BLAS: ベクトル-ベクトル, 行列-ベクトル, 行列-行列 演算
ルーチン群で、参照実装をしてる。API を決めたい。高速な
実装もある: GotoBLAS, Intel MKL, ATLAS など。
LAPACK: 線形連立方程式、固有値問題、最小二乗フィッティ
ング、特異値分解などを BLAS を用いて行う。
標準的なライブラリでみんなが使っている (気づかないで
使っている場合が多い)
Netlib へのアクセスは、8400 万回以上 (2010/11/27)
人類にとって BLAS/LAPACK は非常に重要 riken-logo
.
中田 真秀
.
.
.
高精度線形代数演算ライブラリ MPACK の開発
.
.
ペタ・エクサへ向けた高精度線形代数演算の必要性
. BLAS、LAPACK とは何?
BLAS: ベクトル-ベクトル, 行列-ベクトル, 行列-行列 演算
ルーチン群で、参照実装をしてる。API を決めたい。高速な
実装もある: GotoBLAS, Intel MKL, ATLAS など。
LAPACK: 線形連立方程式、固有値問題、最小二乗フィッティ
ング、特異値分解などを BLAS を用いて行う。
標準的なライブラリでみんなが使っている (気づかないで
使っている場合が多い)
Netlib へのアクセスは、8400 万回以上 (2010/11/27)
人類にとって BLAS/LAPACK は非常に重要 riken-logo
.
中田 真秀
.
.
.
高精度線形代数演算ライブラリ MPACK の開発
.
.
ペタ・エクサへ向けた高精度線形代数演算の必要性
. MPACK 0.6.7: Multiple precision version of BLAS and LAPACK
http://mplapack.sourceforge.net/
中田真秀 @ 理化学研究所
MPACK (MBLAS and MLAPACK) BLAS/LAPACK の任意精度版。
参照実装 (API) の提供が第一目標。
最新バージョン 0.6.7 (2010/8/20); 666 個中、100 個の
MLAPACK, 全ての MBLAS をサポート。
高い移植性: Linux/Windows/MacOSX/FreeBSD
多くの多倍長精度演算ライブラリに対応:
MPFR/MPC/GMP/DD/QD/double
C++; よりわかりやすく、より高速、よりプログラマにやさしく。 riken-logo
2-条項 BSD ライセンス; 商用利用、改変、再配付自由。
.
中田 真秀
.
.
.
高精度線形代数演算ライブラリ MPACK の開発
.
.
ペタ・エクサへ向けた高精度線形代数演算の必要性
. MPACK でつかっている高精度計算の手法
約 4 倍,8 倍精度: 高速なもの: QD, DD
任意精度:とりあえず高速な高精度計算: GMP :
より精密な高精度計算 (丸めモード制御など): MPFR/MPC
参照用: double
riken-logo
.
中田 真秀
.
.
.
高精度線形代数演算ライブラリ MPACK の開発
.
.
ペタ・エクサへ向けた高精度線形代数演算の必要性
. 参照実装の提供; 資産活用; 命名規則
Prefix の変化
float, double → “R”eal,
complex, double complex → “C”omplex.
daxpy, zaxpy → Raxpy, Caxpy
dgemm, zgemm → Rgemm, Cgemm
dsterf, dsyev → Rsterf, Rsyev
dzabs1, dzasum → RCabs1, RCasum
riken-logo
.
中田 真秀
.
.
.
高精度線形代数演算ライブラリ MPACK の開発
.
.
ペタ・エクサへ向けた高精度線形代数演算の必要性
. MBLAS 0.6.7 全ルーチン
Crotg
Rswap
Rdot
iRamax
Cgemv
Rsbmv
Chpmv
Rtpmv
Cgerc
Rspr2
Cgemm
Csyr2k
Cscal
CRscal
Cdotc
RCabs1
Rgemv
Ctrmv
Rsymv
Ctrsv
Cher
Rgemm
Rsyr2k
LEVEL1 MBLAS
Rrotg
Rrot
Rrotm
Rscal
Ccopy
Rcopy
Cdotu
RCnrm2 Rnrm2
Mlsame Mxerbla
Cgbmv
Cgemv
Rsbmv
Rtrsv
Chpr
LEVEL2 MBLAS
Rgbmv Chemv
Rgemv Cgbmv
Rspmv
Ctrmv
Ctbsv
Rtbsv
Cher2
Chpr2
Chbmv
Rgemv
Rtrmv
Ctpsv
Rsyr
LEVEL3 MBLAS
Csymm Rsymm Chemm
Cher2k
Ctrmm
Rtrmm
.
中田 真秀
CRrot
Caxpy
Rasum
Chpmv
Chemv
Ctbmv
Rger
Rspr
Csyrk
Ctrsm
.
Cswap
Raxpy
iCasum
Rsymv
Chbmv
Ctpmv
Cgeru
Rsyr2
Rsyrk
Rtrsm
.
.
高精度線形代数演算ライブラリ MPACK の開発
Cherk
riken-logo
.
.
ペタ・エクサへ向けた高精度線形代数演算の必要性
. MLAPACK 0.6.7 全ルーチン
Mutils
Rlanst
Rlapy3
Claset
Clascl
Rlarf
Clarft
Rlatrd
Cheev
Ctrtri
Rgetri
Ctrtrs
Cspmv
Rposv
Rpotri
Rlamch
Clanht
Rladiv
Rlasr
Rlasrt
Clarf
Rlarfb
Clatrd
Rpotrf
Rgetf2
Cgetri
Rlasyf
Cspr
Rgeequ
Rpocon
Rlae2
Rlansy
Cladiv
Clasr
Rsytd2
Rorg2l
Clarfb
Rsytrd
Cpotrf
Cgetf2
Rgetrs
Clasyf
Csymv
Rlatrs
Rlaev2
Clansy
Clarfg
Rpotf2
Chetd2
Cung2l
Rorgqr
Chetrd
Clacrm
Rlaswp
Cgetrs
Clahef
Csyr
Rlange
Claev2
Clanhe
Rlartg
Clacgv
Rsteqr
Rorg2r
Cungqr
Rorgtr
Rtrti2
Claswp
Rgesv
Clacrt
iCmax1
Rgecon
Rlassq
Rlapy2
Clartg
Cpotf2
Csteqr
Cung2r
Rorgql
Cungtr
Ctrti2
Rgetrf
Cgesv
Claesy
RCsum1
Rlauu2
Classq
Rlarfg
Rlaset
Rlascl
Rsterf
Rlarft
Cungql
Rsyev
Rtrtri
Cgetrf
Rtrtrs
Crot
Rpotrs
Rlauum
riken-logo
赤字は 0.7.6 で新規さポートされた 10 ルーチン
.
中田 真秀
.
.
.
高精度線形代数演算ライブラリ MPACK の開発
.
.
ペタ・エクサへ向けた高精度線形代数演算の必要性
. 参照実装の提供; 資産活用; 関数のコール例
違いは call by value or call by reference
MBLAS/MLAPACK:
Rgemm("n", "n",
Rgetrf(n, n, A,
Rgetri(n, A, n,
Rsyev("V", "U",
n, n, n, alpha, A, n, B, n, beta, C, n);
n, ipiv, &info);
ipiv, work, lwork, &info);
n, A, n, w, work, &lwork, &info);
BLAS/LAPACK:
dgemm_f77("N", "N", &n, &n, &n, &One, A, &n, A, &n, &Zero, C, &
dgetri_f77(&n, A, &n, ipiv, work, &lwork, &info);
(C++/C から BLAS/LAPACK を呼べる lapack.h, blas.h 提供中!)
riken-logo
.
中田 真秀
.
.
.
高精度線形代数演算ライブラリ MPACK の開発
.
.
ペタ・エクサへ向けた高精度線形代数演算の必要性
. MBLAS ソースコードから抜粋
コードはすべて共通:メタ・プログラミング Caxpy; axpy の複素
数版
void Caxpy(INTEGER n, COMPLEX ca, COMPLEX * cx, INTEGER incx, COMPLEX
{
REAL Zero = 0.0;
if (n <= 0)
return;
if (RCabs1(ca) == Zero)
return;
INTEGER ix = 0;
INTEGER iy = 0;
if (incx < 0)
ix = (-n + 1) * incx;
if (incy < 0)
iy = (-n + 1) * incy;
for (INTEGER i = 0; i < n; i++) {
cy[iy] = cy[iy] + ca * cx[ix];
ix = ix + incx;
riken-logo
.
中田 真秀
.
.
.
高精度線形代数演算ライブラリ MPACK の開発
.
.
ペタ・エクサへ向けた高精度線形代数演算の必要性
. MLAPACK ソースコードから抜粋
Rsyev; 対称行列対角化ルーチン
Rlascl(uplo, 0, 0, One, sigma, n, n, A, lda, info);
}
//Call DSYTRD to reduce symmetric matrix to tridiagonal form.
inde = 1;
indtau = inde + n;
indwrk = indtau + n;
llwork = *lwork - indwrk + 1;
Rsytrd(uplo, n, &A[0], lda, &w[0], &work[inde - 1], &work[indtau &work[indwrk - 1], llwork, &iinfo);
//For eigenvalues only, call DSTERF. For eigenvectors, first call
//DORGTR to generate the orthogonal matrix, then call DSTEQR.
if (!wantz) {
Rsterf(n, &w[0], &work[inde - 1], info);
} else {
Rorgtr(uplo, n, A, lda, &work[indtau - 1], &work[indwrk - 1],
&iinfo);
Rsteqr(jobz, n, w, &work[inde - 1], A, lda, &work[indtau - 1],
}
//If matrix was scaled, then rescale eigenvalues appropriately. riken-logo
if (iscale == 1) {
if (*info == 0) {
.
中田 真秀
.
.
.
高精度線形代数演算ライブラリ MPACK の開発
.
.
ペタ・エクサへ向けた高精度線形代数演算の必要性
. MPACK (MBLAS/MLAPACK) に関する事実
研究で使われつつある:半正定値計画法の高精度
版:SDPA-GMP, -QD, -DD (中田が主に; 論文もいくつか)
Google で”Multiple precision BLAS” で検索するとほぼ
MPACK しか出てこない。
Project web hits 20000+, download 1100+ (2010/11/24)
riken-logo
.
中田 真秀
.
.
.
高精度線形代数演算ライブラリ MPACK の開発
.
.
ペタ・エクサへ向けた高精度線形代数演算の必要性
. MPACK (MBLAS/MLAPACK) のパフォーマンス
Core i7 920/Ubuntu 10.04/9.04 amd64 上で計測、
DD(double-double) 型
Raxpy 600MFlops (mpack 0.6.6; gotoblas 1.5GFlops)
Rgemv 140MFlops (gotoblas 3.8GFlops)
Rgemm 140MFlops (gotoblas 42.5GFlops; 中里ら 2010:
30GFlops on RADEON HD5870、椋木ら 10GFlops)
Rsyev 112s (1000 次元; gotoblas 1.83s)
riken-logo
.
中田 真秀
.
.
.
高精度線形代数演算ライブラリ MPACK の開発
.
.
ペタ・エクサへ向けた高精度線形代数演算の必要性
. アプリケーション
半正定値計画法ソルバー: SDPA-GMP, SDPA-QD,
SDPA-DD: 非常に正確に解けるようになった。Powered by
MPACK
riken-logo
.
中田 真秀
.
.
.
高精度線形代数演算ライブラリ MPACK の開発
.
.
ペタ・エクサへ向けた高精度線形代数演算の必要性
. MPACK 0.6.7: Multiple precision version of BLAS and LAPACK
http://mplapack.sourceforge.net/
中田真秀 @ 理化学研究所
MPACK (MBLAS and MLAPACK) BLAS/LAPACK の任意精度版。
参照実装 (API) の提供が第一目標。
最新バージョン 0.6.7 (2010/8/20); 666 個中、100 個の
MLAPACK, 全ての MBLAS をサポート。
高い移植性: Linux/Windows/MacOSX/FreeBSD
多くの多倍長精度演算ライブラリに対応:
MPFR/MPC/GMP/DD/QD/double
C++; よりわかりやすく、より高速、よりプログラマにやさしく。 riken-logo
2-条項 BSD ライセンス; 商用利用、改変、再配付自由。
.
中田 真秀
.
.
.
高精度線形代数演算ライブラリ MPACK の開発
.
.