PC・WS 用の多倍長計算環境の開発と その利用法 (抜粋

Motivation
PC・WS 用の多倍長計算環境の開発と
その利用法
(抜粋)
藤原 宏志
京都大学大学院・情報学研究科
多倍長精度計算フォーラム
2009 年 2 月 24 日
工学院大学
Motivation
逆問題などの非適切問題 ill-posed problems の数値的取り扱い
数値解析理論の構築，数値計算の実現，計算力学への応用
磯祐介 (京都大)，今井仁司 (徳島大) らとの共同研究
逆問題 : Inverse Problem
トモグラフィ，リモートセンシングによる非破壊検査，地中探査など
EIT
: 楕円型方程式の Cauchy 問題
X 線 CT : 第一種積分方程式
FUJIWARA Hiroshi, Feb. 24th, 2009
Tutorial of exflib in 多倍長精度計算フォーラム@工学院大学
Well-Posedness and Ill-posedness in the sense of Hadamard
問題が Hadamard の意味で適切 well-posed とは，
FUJIWARA Hiroshi, Feb. 24th, 2009
Tutorial of exflib in 多倍長精度計算フォーラム@工学院大学
Difficulty : Instability
Difficulty
非適切性と，それに起因する数値的不安定性による誤差の増大
解が存在し，
解が唯一つで，
考えているノルムで，解がデータに連続に依存する (安定性)
本講演では，ある入力に対して計算値が “増大” し，それが計算
非適切性 ill-posedness は，適切性の否定概念である．
数値計算においては，安定性の欠如が問題となることが多い．
結果に “顕著な悪影響” を与える入力値を入力し得るとき，その
計算過程は「数値的に不安定」であるという．
適切性は，考える枠組 (函数空間とノルム) に依存する．
浮動小数点演算におけるスキームの数値的不安定性は，計算 (演
数学的に適切でも，離散化が数値的に不安定となる場合もある．
算) の精度，入力値の集合など，枠組に応じて論じるのが自然．
FUJIWARA Hiroshi, Feb. 24th, 2009
Tutorial of exflib in 多倍長精度計算フォーラム@工学院大学
FUJIWARA Hiroshi, Feb. 24th, 2009
Tutorial of exflib in 多倍長精度計算フォーラム@工学院大学
Example of Numerical Instability
34
an+2 =
11
an+1 −
a0 = 1,
a1 =
3
11
1
数値的に不安定な問題に対し，計算の大規模・高速化だけでは
an
an =
1
11
n
11
計算環境によって，同じプログラムの計算結果が異なる．
倍精度での数値計算では，どちらも厳密解を近似しない．
2500
2000
Xeon
Alpha
n
1500
an
1000
500
0
-500
-1000
-1500
-2000
36
38
40
42
44
Xeon(Linux)
Alpha
2
5
10
0.00826446
6.20921 × 10−6
3.86014 × 10−11
0.00826446
6.20921 × 10−6
3.85159 × 10−11
40
50
9.68369
571812
−7.91327
−467270
高精度かつ高信頼の数値計算は期待されない．
多倍長計算により，充分桁数を確保し，丸め誤差の影響が計算結
果に現れないという意味で，仮想的に丸め誤差のない数値計算を
実現する．
非適切問題に対する数値解析の枠組を構築する．
科学技術計算，大規模計算
100 桁–1000 桁程度
64 ビット環境での高速な演算
46
n
FORTRAN, C 言語で利用
計算誤差 (本例では丸め誤差) が急激に増大する．
40 項，50 項の計算に倍精度計算では不充分である．
FUJIWARA Hiroshi, Feb. 24th, 2009
を目的とした多倍長計算環境の設計と実装
Tutorial of exflib in 多倍長精度計算フォーラム@工学院大学
浮動小数点数
FUJIWARA Hiroshi, Feb. 24th, 2009
Tutorial of exflib in 多倍長精度計算フォーラム@工学院大学
Floating-Point Numbers
10 進 4 桁の浮動小数点表記の例:
2723 = +2.723 × 103
0.00037 = +3.700 × 10−4
符号部，仮数部，指数部で表現する．
初期の頃について，講演者の調べた範囲では，
1910 年代，L. Torres y Quevedo (スペイン) が
10 進法浮動小数点演算装置の設計図？, 10 進？桁
1930 年代，Konrad Zuse (ドイツ) が 2 進法浮動小数点演算
(指数 7，仮数 16 ビット) を実装した Z1 を稼働
type
size
bit
significand∗
digits
bit
IEEE754 single
IEEE754 double
32
64
24
53
bits in
exponent
magnitude
8
11
10±38
10±308
IEEE P754 binary128 128
113
34.0
15
IA32 extend
80
65
19.6
15
SPARC extend
128
113
34.0
15
†
POWER extend
128 ≥ 106
≥ 31.9
11
(*) including the hidden one
(†) POWER extended is sum of two IEEE754 doubles
10±4932
10±4932
10±4932
10±308
7.2
16.0
1985 年，IEEE754 規格の策定
FUJIWARA Hiroshi, Feb. 24th, 2009
Tutorial of exflib in 多倍長精度計算フォーラム@工学院大学
FUJIWARA Hiroshi, Feb. 24th, 2009
Tutorial of exflib in 多倍長精度計算フォーラム@工学院大学
Data Structure of the Proposed Multiple-Precision Number
多倍長数の表現に，64 ビット整数の配列をもちいる．
Proposed Multiple-Precision Arithmetic Software : exflib
60–20,000 桁 (10 進法)
四則演算，組み込み関数，丸め操作，入出力
eb
s
1
f1
-
f2
-
63
-
64
···
f3
fn
-
64
-
アセンブリ言語での設計と実装
64
(−1)s × 2eb −BIAS ×
算術式の記述は，組込みの型と同じ
64
f0 +
n
fi
i=1
AMD64 (Opteron, Athlon64), Intel 64 (Xeon, Core, Pentium),
IA-32, SPARC V9, Alpha
264i
Solaris, Linux, MacOS X (64bit), Windows (32bit)
64 ビットの配列ひとつは，10 進法では 19.6 桁に相当
60 桁 : n = 4, 80 桁 : n = 5, 100 桁 : n = 6, 500 桁 : n = 26.
18
10−10
から
18
1010
の範囲
FUJIWARA Hiroshi, Feb. 24th, 2009
Tutorial of exflib in 多倍長精度計算フォーラム@工学院大学
Opteron150 2.4GHz
100, mul
div
1000, mul
div
10000, mul
div
FUJIWARA Hiroshi, Feb. 24th, 2009
Tutorial of exflib in 多倍長精度計算フォーラム@工学院大学
LU decomp of the Hilbert matrix with exflib
unit : µ sec.
64-bit, Linux, GCC 3
digits, op
MPI, OpenMP
http://www-an.acs.i.kyoto-u.ac.jp/~fujiwara/exflib
Arithmetic Benchmark (2007 Jan)
mpfun90
C++, FORTRAN90
FMLIB
unit : sec.
Opteron246 2GHz
Maple
Mathematica
FORTRAN90
MPFR
C
Pari
exflib
C
C
3.0
1.1
22
1.5
0.24
0.29
0.14
3.5
1.4
26
5.9
0.52
1.2
0.33
60
17
100
17
6.2
15
3.7
63
20
93
50
16
23
8.3
–
1105
3400
580
300
1400
280
–
1251
3000
1800
730
1500
630
digits, size
MB
Cond.
Num.
·2
Athlon64
3200+
2.2GHz
Core i7
920
2.66GHz
Core2 Duo
Pentium4 EM64T
PentiumM
E6600
2.4GHz
(Prescott)
3.6(F)GHz
(Banias)
1.0GHz
500, 100
200
2
9
0.57
4.1
0.53
3.99
0.70
4.89
1.16
8.5
5.5
38.6
1000, 100
200
400
5
17
68
1.65
11.5
88
1.51
11.5
91.5
2.06
14.1
108
3.81
26.9
207
18.7
131
1005
2000, 100
200
400
800
9
34
135
538
5.8
38
284
2282
5.1
39.7
311
2355
7.4
48
360
2821
14.6
97
732
5789
74
489
3645
28565
10150
10303
10609
101221
cf. http://www.medicis.polytechnique.fr/~pphd/mpfr/timings-220.html
FUJIWARA Hiroshi, Feb. 24th, 2009
Tutorial of exflib in 多倍長精度計算フォーラム@工学院大学
FUJIWARA Hiroshi, Feb. 24th, 2009
Tutorial of exflib in 多倍長精度計算フォーラム@工学院大学
LU decomp. on T2K@kudpc, Jun 2008
exflib for FORTRAN90 — Required Files
Opteron 8350 (Quad-Core, 2.3GHz) x 4 CPU x 8 node
RHE Linux V4, Fujitsu MPI and FCC
decimal digits
size
500
8000
1000
8000
1971 16612
1.00
1.00
2 node 32 proc
1013
8369
0.97
0.99
8 node 128 proc
289
2202
0.85
0.94
upper : computational time (sec.)
lower : paralell efficiency by node
47175
1.00
23559
1.00
6144
0.96
1 node
100
8000
16 proc
http://www-an.acs.i.kyoto-u.ac.jp/~fujiwara/exflib
配布物
libexfloat.a, exfloat.h, exflib.F90
サンプル・プログラム (C++言語, FORTRAN90, 並列計算)
ドキュメント
FORTRAN90 でのコンパイルに必要なファイル
ユーザ・プログラム
libexfloat.a
exflib.F90
8000 × 8000 : 1.706347 × 1011 MAC ops
128 proc, 100 digits : 590 MFLOPS, 500 digits : 77 MFLOPS, 1000 digits : 28 MFLOPS
FUJIWARA Hiroshi, Feb. 24th, 2009
Tutorial of exflib in 多倍長精度計算フォーラム@工学院大学
exflib for FORTRAN90 — Compiling
G95
Intel Fortran (V9,10,11)
PGI Compiler (Ver. 5,6,7,8)
% pgf90 -DEXFLIB EXFLOAT PRECISION10=100 \
-O2 -Mfree exflib.F90 user.f90 libexfloat.a
代入
x = y
x = 1
x = i
! 多倍長数の代入
! リテラル整数の代入
! 整数型の変数の代入
x = 0.1
! 単精度数，倍精度数の代入は禁止
x = ’0.1’
x = ’1/10’
x = ’#PI/2’
! 文字列による代入
! 文字列中の四則演算も可
! 組み込みの定数 (円周率)
GCC gfortran (ver.4) にも対応
FUJIWARA Hiroshi, Feb. 24th, 2009
型宣言
USE exflib
TYPE(exfloat) :: x, y, v(N), a(N, M)
% ifort -DEXFLIB EXFLOAT PRECISION10=100 \
-O2 -free exflib.F90 user.f90 libexfloat.a
Tutorial of exflib in 多倍長精度計算フォーラム@工学院大学
exflib for FORTRAN90 – Type Declaration
% g95 -DEXFLIB EXFLOAT PRECISION10=100 \
-m64 -O2 -ffree-form \
exflib.F90 user.f90 libexfloat.a
FUJIWARA Hiroshi, Feb. 24th, 2009
Tutorial of exflib in 多倍長精度計算フォーラム@工学院大学
FUJIWARA Hiroshi, Feb. 24th, 2009
Tutorial of exflib in 多倍長精度計算フォーラム@工学院大学
exflib for FORTRAN90 – Assignment
代入 (続き)
exflib for FORTRAN90 — Arithmetic, Built-in Functions
演算，組み込み関数
! リテラル/組み込み定数の演算も可
x = ’12.34*(5/6)*(7e-8)*#E’
!
!
x
x
! 整数との四則演算は可
x = x + 1
単 (倍) 精度数の代入は，明示的に型変換する
代入された値の精度は単 (倍) 精度になる
= exflib_cast(0.1)
= exflib_cast(double_x)
! 単 (倍) 精度数との四則演算は禁止
h = h + 0.1
TYPE(exfloat) :: u(N), v(N)
u = v
! 配列の代入には対応しない
double_x = x
! 倍精度数への代入は可
FUJIWARA Hiroshi, Feb. 24th, 2009
x = a * b
y = exp(x)
FUJIWARA Hiroshi, Feb. 24th, 2009
Tutorial of exflib in 多倍長精度計算フォーラム@工学院大学
exflib for FORTRAN90 — Output
出力
比較
IF ( x < y ) THEN
...
! 通常の表記，小数点以下 10 進 100 桁
CHARACTER(120) :: str
str = exflib_format(’F.100’, x)
WRITE(*,*) TRIM(str)
! 単 (倍) 精度数との比較は可
FUJIWARA Hiroshi, Feb. 24th, 2009
! 倍精度に型変換して出力
WRITE(*,*) exflib_double(x)
IF ( x < 0.1 ) THEN
...
組み込み関数は C 言語に準拠．
Tutorial of exflib in 多倍長精度計算フォーラム@工学院大学
exflib for FORTRAN90 — Comparison
! 配列演算には対応しない
u(:) = u(:) + v(:)
Tutorial of exflib in 多倍長精度計算フォーラム@工学院大学
! 浮動小数点表記，10 進 100 桁
str = exflib_format(’E.100’, x)
WRITE(*,*) TRIM(str)
FUJIWARA Hiroshi, Feb. 24th, 2009
Tutorial of exflib in 多倍長精度計算フォーラム@工学院大学
exflib for FORTRAN90 – Parallelization
Translation from IEEE754 double to exfloat
MPI では，メンバ num(0:exflib exfloat size-1) を操作
CALL MPI_SEND(x%num(0), exflib_exfloat_size, &
MPI REAL8, dst, TAG, MPI_COMM_WORLD, ierr)
exfloat.h のインクルード (C++言語)
use exflib (FORTRAN90)
型宣言
倍精度のリテラル定数の代入・演算
CALL MPI_RECV(x%num(0), exflib_exfloat_size, &
MPI REAL8, src, TAG, MPI_COMM_WORLD, ierr, ierr)
CALL MPI_BCAST(x%num(0), exflib_exfloat_size, &
MPI REAL8, src, MPI_COMM_WORLD)
出力 (FORTRAN90)
特殊函数，配列演算，並列演算，その他ライブラリ
コンパイル・コマンド
machine epsilon, ULP などの設定
OpenMP は，縮約 (reduction) 演算以外には，おおむね対応
FUJIWARA Hiroshi, Feb. 24th, 2009
Tutorial of exflib in 多倍長精度計算フォーラム@工学院大学
計算精度についての注意
log 264 ≈ 19.6 桁単位での設定となる．
指定桁数
decimal
計算精度
decimal
メモリ
bytes
仮数部の情報量
bytes
58 – 77
78 – 96
97 – 115
..
.
77.06
96.33
115.60
40
48
56
24.08–31.97
32.39–39.86
40.28–47.75
200
500
211.93
500.91
96
216
83.05
207.62
精度に関して以下の変数をもつ．sample-f90/version.f90 を参照．
exflib_exfloat_precision10 指定桁数
exflib_exfloat_precision
実際の精度，10 進桁数
exflib_exfloat_size
多倍長数 1 つのメモリ量
FUJIWARA Hiroshi, Feb. 24th, 2009
Tutorial of exflib in 多倍長精度計算フォーラム@工学院大学
FUJIWARA Hiroshi, Feb. 24th, 2009
Tutorial of exflib in 多倍長精度計算フォーラム@工学院大学