🍝
単精度・倍精度・四倍精度・八倍精度で固有値問題を解く
数値計算 Advent Calendar 2022 19日目です🎄
このノートでは単精度・倍精度・（ついでに）疑似四倍精度・四倍精度・八倍精度でフランク行列・ヒルベルト行列の固有値を求めます. フランク行列の場合, 1000×1000程のサイズで単精度では最小固有値が1桁も正しく計算できなくなることを確かめました. より条件数の大きいヒルベルト行列の場合, さらに小さい5×5～35×35程で単精度・倍精度・疑似四倍精度・四倍精度・八倍精度では最小固有値が1桁も正しく計算できなくなることを確かめました.
パッケージ・環境
Quadmath.jlの四倍精度・Float128はIEEE 754 binary128準拠です. ついでにDoubleFloats.jlの疑似四倍精度・Double64も追加しました. QD.jlの八倍精度・Float256はインストール, ビルド時にエラーが出たので, 今回は代わりにBigFloatを使いました. BigFloatの仮数は自由に変えられるので237ビットにしてIEEE 754 binary256に合わせましたが, 指数部分はデフォルトから変えられないらしいのでIEEE 754 binary256より大きくなっています. なお, BigFloatのバックエンドはMPFRで, Juliaに最初から入っているので特別なパッケージ等は不要です. 後ほど, setprecision(237)のようにして仮数の桁を変えます. 固有値については, Float32とFloat64までは標準ライブラリのLinearAlgebra.jl（固有値部分はLAPACK）で計算できます. Double64, Float128, BigFloatはGenericLinearAlgebra.jlで計算します.
入力
# using Pkg
# Pkg.add("Quadmath")
# Pkg.add("DoubleFloats")
# Pkg.add("GenericLinearAlgebra")
using Printf # 表示
using DoubleFloats # Double64型のサポート
using Quadmath # Float128型のサポート
using LinearAlgebra # Float32, Float64までの固有値計算
using GenericLinearAlgebra # Double64, Float128, BigFloatの固有値計算
versioninfo()
出力
Julia Version 1.8.3
Commit 0434deb161 (2022-11-14 20:14 UTC)
Platform Info:
OS: Windows (x86_64-w64-mingw32)
CPU: 8 × 11th Gen Intel(R) Core(TM) i7-1185G7 @ 3.00GHz
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-13.0.1 (ORCJIT, tigerlake)
Threads: 1 on 8 virtual cores
BigFloatを八倍精度の代わりに使う
まず, BigFloatが単精度・倍精度・疑似四倍精度・四倍精度の代わりに使えるのか試します. それぞれの型の仮数部分のビット数はそれぞれ下記のようになっています.
入力
for t in [Float32, Float64, Double64, Float128, BigFloat]
println(t, " \t", precision(t))
end
出力
Float32 24
Float64 53
Double64 106
Float128 113
BigFloat 256
例えば, setprecision(53)とするとBigFloatとFloat64の仮数部分のビット数が一致しますので, 近い振る舞いをすると考えられます. 有理数・Rational型の
入力
for t in [Float32, Float64, Double64, Float128]
p = precision(t)
setprecision(p) do
println("$t \t", big(convert(t,1//3)))
println("BigFloat($(precision(BigFloat)))\t", BigFloat(1//3))
println()
end
end
setprecision(237) do
println("BigFloat($(precision(BigFloat)))\t", BigFloat(1//3))
println()
end
println("BigFloat($(precision(BigFloat)))\t", BigFloat(1//3))
出力
Float32 0.333333343
BigFloat(24) 0.333333343
Float64 0.33333333333333331
BigFloat(53) 0.33333333333333331
Double64 0.333333333333333333333333333333335
BigFloat(106) 0.333333333333333333333333333333335
Float128 0.333333333333333333333333333333333317
BigFloat(113) 0.333333333333333333333333333333333317
BigFloat(237) 0.3333333333333333333333333333333333333333333333333333333333333333333333326
BigFloat(256) 0.3333333333333333333333333333333333333333333333333333333333333333333333333333348
BigFloatの仮数部分を他の型に合わせると, 少なくとも有理数・Rationalから変換した際の桁数は一致することが確かめられました. 完全には一致しないかもしれませんが, BigFloat(237)で概ねFloat256の目安にはなると思われます. 実際, 後ほどの固有値計算でもFloat128とBigFloat(113)は同じ行列サイズまで計算できるという結果を与えています.
フランク行列
フランク行列の実装例を示します. まず整数型で要素を作ってから, datatypeキーワードで指定した型へ変換するようにしてあります. なお, 速度よりもFortran等への移植性を優先しました.
入力
# 速い例：
# Frank(n; datatype=Int64) = convert.(datatype, [n - max(i,j) + 1 for i=1:n, j=1:n])
# 遅い例：
function Frank(n; datatype=Int64)
A = zeros(datatype, n, n)
for j in 1:n
for i in 1:n
A[i,j] = convert(datatype, n - max(i,j) + 1)
end
end
return A
end
A = Frank(5)
出力
5×5 Matrix{Int64}:
5 4 3 2 1
4 4 3 2 1
3 3 3 2 1
2 2 2 2 1
1 1 1 1 1
幸谷智紀『LAPACK/BLAS入門』(森北出版, 2016) p.76 表4.7を参考に
入力
# using LinearAlgebra
# using GenericLinearAlgebra
for n in [10,100,1000]
println("\nn = $n \tλ[1]")
for datatype in [Float32,Float64]#,Double64,Float128,BigFloat]
A = Frank(n, datatype=datatype)
λ = eigvals(A)
println(datatype, " \t", λ[1])
end
end
出力
n = 10 λ[1]
Float32 0.25567994
Float64 0.2556795627964355
n = 100 λ[1]
Float32 0.2500508
Float64 0.2500610827206968
n = 1000 λ[1]
Float32 0.12658979
Float64 0.2500006162349238
Float32では1000×1000のフランク行列の最小固有値が1桁も正しく計算できないことが確かめられました. 実際には900×900程度から振動し始め, おかしな挙動を見せます. フランク行列だとメモリや計算時間の限界が先に来ますので, より悪条件なヒルベルト行列で単精度・倍精度・疑似四倍精度・四倍精度・八倍精度の限界を見てみます.
ヒルベルト行列
ヒルベルト行列の実装例を示します. まず有理数型で要素を作り, それをdatatypeキーワードで指定した型に変換します. 1 / (i + j - 1)だと勝手にFloat64と解釈されて丸め誤差が発生してしまいますので, 1 // (i + j - 1)とすることでRationalから他の型に直接変換し, 丸め誤差を小さくできます. なお, 速度よりもFortran等への移植性を優先しました.
入力
# 速い例：
# Hilbert(n; datatype=Rational) = convert.(datatype, [1 // (i + j - 1) for i=1:n, j=1:n])
# 遅い例：
function Hilbert(n; datatype=Rational)
A = zeros(datatype, n, n)
for j in 1:n
for i in 1:n
A[i,j] = convert(datatype, 1 // (i + j - 1))
end
end
return A
end
H = Hilbert(5)
出力
5×5 Matrix{Rational}:
1//1 1//2 1//3 1//4 1//5
1//2 1//3 1//4 1//5 1//6
1//3 1//4 1//5 1//6 1//7
1//4 1//5 1//6 1//7 1//8
1//5 1//6 1//7 1//8 1//9
ヒルベルト行列のサイズを変えてFloat32, Float64, Double64, Float128, BigFloat(113), BigFloat(237), BigFloat(256)で固有値を計算します.
入力
# using Printf
# using DoubleFloats
# using Quadmath
# using LinearAlgebra
# using GenericLinearAlgebra
for n in 2:40
println("\nn = $n \tλ[1]")
for datatype in [Float32, Float64, Double64, Float128]
A = Hilbert(n, datatype=datatype)
λ = eigvals(A)
@printf("%s \t%.9e\n", datatype, λ[1])
end
for p in [113,237,256,1024]
setprecision(p) do
A = Hilbert(n, datatype=BigFloat)
λ = eigvals(A)
@printf("%s \t%.9e\n", "BigFloat($(precision(BigFloat)))", λ[1])
end
end
end
# setprecision(BigFloat, 256)
出力
n = 2 λ[1]
Float32 6.574144959e-02
Float64 6.574145409e-02
Double64 6.574145409e-02
Float128 6.574145409e-02
BigFloat(113) 6.574145409e-02
BigFloat(237) 6.574145409e-02
BigFloat(256) 6.574145409e-02
BigFloat(1024) 6.574145409e-02
n = 3 λ[1]
Float32 2.687271917e-03
Float64 2.687340356e-03
Double64 2.687340356e-03
Float128 2.687340356e-03
BigFloat(113) 2.687340356e-03
BigFloat(237) 2.687340356e-03
BigFloat(256) 2.687340356e-03
BigFloat(1024) 2.687340356e-03
n = 4 λ[1]
Float32 9.672460874e-05
Float64 9.670230402e-05
Double64 9.670230402e-05
Float128 9.670230402e-05
BigFloat(113) 9.670230402e-05
BigFloat(237) 9.670230402e-05
BigFloat(256) 9.670230402e-05
BigFloat(1024) 9.670230402e-05
n = 5 λ[1]
Float32 3.361120889e-06
Float64 3.287928772e-06
Double64 3.287928772e-06
Float128 3.287928772e-06
BigFloat(113) 3.287928772e-06
BigFloat(237) 3.287928772e-06
BigFloat(256) 3.287928772e-06
BigFloat(1024) 3.287928772e-06
n = 6 λ[1]
Float32 1.306005259e-07
Float64 1.082799486e-07
Double64 1.082799485e-07
Float128 1.082799485e-07
BigFloat(113) 1.082799485e-07
BigFloat(237) 1.082799485e-07
BigFloat(256) 1.082799485e-07
BigFloat(1024) 1.082799485e-07
n = 7 λ[1]
Float32 7.071950847e-08
Float64 3.493898716e-09
Double64 3.493898606e-09
Float128 3.493898606e-09
BigFloat(113) 3.493898606e-09
BigFloat(237) 3.493898606e-09
BigFloat(256) 3.493898606e-09
BigFloat(1024) 3.493898606e-09
n = 8 λ[1]
Float32 -4.826334887e-08
Float64 1.111539635e-10
Double64 1.111538966e-10
Float128 1.111538966e-10
BigFloat(113) 1.111538966e-10
BigFloat(237) 1.111538966e-10
BigFloat(256) 1.111538966e-10
BigFloat(1024) 1.111538966e-10
n = 9 λ[1]
Float32 2.470161853e-09
Float64 3.499736372e-12
Double64 3.499676403e-12
Float128 3.499676403e-12
BigFloat(113) 3.499676403e-12
BigFloat(237) 3.499676403e-12
BigFloat(256) 3.499676403e-12
BigFloat(1024) 3.499676403e-12
n = 10 λ[1]
Float32 -6.898610394e-08
Float64 1.093880753e-13
Double64 1.093153819e-13
Float128 1.093153819e-13
BigFloat(113) 1.093153819e-13
BigFloat(237) 1.093153819e-13
BigFloat(256) 1.093153819e-13
BigFloat(1024) 1.093153819e-13
n = 11 λ[1]
Float32 -8.925326256e-08
Float64 3.375236254e-15
Double64 3.393218595e-15
Float128 3.393218595e-15
BigFloat(113) 3.393218595e-15
BigFloat(237) 3.393218595e-15
BigFloat(256) 3.393218595e-15
BigFloat(1024) 3.393218595e-15
n = 12 λ[1]
Float32 -3.159367168e-08
Float64 1.630874793e-16
Double64 1.047946398e-16
Float128 1.047946398e-16
BigFloat(113) 1.047946398e-16
BigFloat(237) 1.047946398e-16
BigFloat(256) 1.047946398e-16
BigFloat(1024) 1.047946398e-16
n = 13 λ[1]
Float32 -6.428467003e-08
Float64 -3.723992712e-17
Double64 3.222901015e-18
Float128 3.222901015e-18
BigFloat(113) 3.222901015e-18
BigFloat(237) 3.222901015e-18
BigFloat(256) 3.222901015e-18
BigFloat(1024) 3.222901015e-18
n = 14 λ[1]
Float32 -2.169667468e-08
Float64 6.961897267e-18
Double64 9.877051735e-20
Float128 9.877051735e-20
BigFloat(113) 9.877051735e-20
BigFloat(237) 9.877051735e-20
BigFloat(256) 9.877051735e-20
BigFloat(1024) 9.877051735e-20
n = 15 λ[1]
Float32 -5.260721281e-08
Float64 -1.970813429e-17
Double64 3.017915296e-21
Float128 3.017915296e-21
BigFloat(113) 3.017915296e-21
BigFloat(237) 3.017915296e-21
BigFloat(256) 3.017915296e-21
BigFloat(1024) 3.017915296e-21
n = 16 λ[1]
Float32 -2.375249863e-08
Float64 -7.991212304e-17
Double64 1.971645939e-22
Float128 9.197419821e-23
BigFloat(113) 9.197419821e-23
BigFloat(237) 9.197419821e-23
BigFloat(256) 9.197419821e-23
BigFloat(1024) 9.197419821e-23
n = 17 λ[1]
Float32 -7.643087940e-08
Float64 -3.084917681e-17
Double64 2.796721645e-24
Float128 2.796721632e-24
BigFloat(113) 2.796721632e-24
BigFloat(237) 2.796721632e-24
BigFloat(256) 2.796721632e-24
BigFloat(1024) 2.796721632e-24
n = 18 λ[1]
Float32 -3.352907640e-08
Float64 -2.922712756e-17
Double64 1.804347127e-25
Float128 8.487413205e-26
BigFloat(113) 8.487413251e-26
BigFloat(237) 8.487413230e-26
BigFloat(256) 8.487413230e-26
BigFloat(1024) 8.487413230e-26
n = 19 λ[1]
Float32 -2.656284792e-08
Float64 -1.790409527e-17
Double64 3.747605947e-27
Float128 3.747605889e-27
BigFloat(113) 3.747605889e-27
BigFloat(237) 2.571239705e-27
BigFloat(256) 2.571239705e-27
BigFloat(1024) 2.571239705e-27
n = 20 λ[1]
Float32 -2.125183407e-08
Float64 -2.797264572e-16
Double64 -1.905490628e-25
Float128 1.233279673e-28
BigFloat(113) 1.233279672e-28
BigFloat(237) 7.777377397e-29
BigFloat(256) 7.777377397e-29
BigFloat(1024) 7.777377397e-29
n = 21 λ[1]
Float32 -5.624558241e-08
Float64 -1.348650289e-17
Double64 -2.714301231e-23
Float128 3.678871784e-30
BigFloat(113) 3.678871738e-30
BigFloat(237) 2.349182538e-30
BigFloat(256) 2.349182538e-30
BigFloat(1024) 2.349182538e-30
n = 22 λ[1]
Float32 -9.640318410e-08
Float64 -1.951507152e-16
Double64 -2.703080487e-30
Float128 -2.703634429e-30
BigFloat(113) -2.703493685e-30
BigFloat(237) 7.086841582e-32
BigFloat(256) 7.086841582e-32
BigFloat(1024) 7.086841582e-32
n = 23 λ[1]
Float32 -1.269348502e-07
Float64 -1.840525074e-16
Double64 -1.395136546e-26
Float128 -1.395136676e-26
BigFloat(113) -1.395136676e-26
BigFloat(237) 2.135463438e-33
BigFloat(256) 2.135463438e-33
BigFloat(1024) 2.135463438e-33
n = 24 λ[1]
Float32 -5.662219849e-08
Float64 -1.536692857e-16
Double64 -6.294692986e-25
Float128 -6.294693199e-25
BigFloat(113) -6.294693199e-25
BigFloat(237) 6.428050325e-35
BigFloat(256) 6.428050325e-35
BigFloat(1024) 6.428050325e-35
n = 25 λ[1]
Float32 -9.834495529e-08
Float64 -1.158400702e-16
Double64 -9.786027902e-24
Float128 3.475441588e-36
BigFloat(113) 7.644790285e-36
BigFloat(237) 1.933092811e-36
BigFloat(256) 1.933092811e-36
BigFloat(1024) 1.933092811e-36
n = 26 λ[1]
Float32 -8.256382245e-08
Float64 -1.745151220e-16
Double64 -9.530727276e-23
Float128 -3.164126445e-29
BigFloat(113) -3.164126445e-29
BigFloat(237) 5.808263531e-38
BigFloat(256) 5.808263531e-38
BigFloat(1024) 5.808263531e-38
n = 27 λ[1]
Float32 -1.064803996e-07
Float64 -3.982034472e-17
Double64 -2.741452206e-34
Float128 -1.690538599e-27
BigFloat(113) -1.690538566e-27
BigFloat(237) 1.743773036e-39
BigFloat(256) 1.743773036e-39
BigFloat(1024) 1.743773036e-39
n = 28 λ[1]
Float32 -8.385376304e-08
Float64 -3.548348719e-16
Double64 -1.942851263e-34
Float128 -2.793018499e-26
BigFloat(113) -2.793018543e-26
BigFloat(237) 5.231305960e-41
BigFloat(256) 5.231305960e-41
BigFloat(1024) 5.231305960e-41
n = 29 λ[1]
Float32 -1.160353804e-07
Float64 -6.199549616e-17
Double64 2.031105710e-35
Float128 -2.822973889e-25
BigFloat(113) -2.822973889e-25
BigFloat(237) 1.568304495e-42
BigFloat(256) 1.568304495e-42
BigFloat(1024) 1.568304495e-42
n = 30 λ[1]
Float32 -5.065797382e-08
Float64 -5.289389957e-17
Double64 -1.102984167e-34
Float128 -2.146399554e-24
BigFloat(113) -2.146399554e-24
BigFloat(237) 4.698636565e-44
BigFloat(256) 4.698636565e-44
BigFloat(1024) 4.698636565e-44
n = 31 λ[1]
Float32 -7.601806118e-08
Float64 -7.329944390e-17
Double64 -2.916407324e-34
Float128 -1.255399562e-26
BigFloat(113) -1.255399570e-26
BigFloat(237) 1.406868380e-45
BigFloat(256) 1.406868380e-45
BigFloat(1024) 1.406868380e-45
n = 32 λ[1]
Float32 -6.675949749e-08
Float64 -2.374236517e-17
Double64 -1.025210430e-34
Float128 -3.257260908e-36
BigFloat(113) -3.257273966e-36
BigFloat(237) 4.210099032e-47
BigFloat(256) 4.210099032e-47
BigFloat(1024) 4.210099032e-47
n = 33 λ[1]
Float32 -3.973698597e-07
Float64 -6.226694532e-16
Double64 7.189674368e-35
Float128 -2.193596836e-36
BigFloat(113) -2.196233777e-36
BigFloat(237) 1.259226436e-48
BigFloat(256) 1.259226436e-48
BigFloat(1024) 1.259226436e-48
n = 34 λ[1]
Float32 -2.554655794e-07
Float64 -3.335481520e-16
Double64 -3.129737789e-34
Float128 -1.976669371e-25
BigFloat(113) -1.976656135e-25
BigFloat(237) 3.764454505e-50
BigFloat(256) 3.764454505e-50
BigFloat(1024) 3.764454505e-50
n = 35 λ[1]
Float32 -2.991216945e-07
Float64 -5.137067885e-16
Double64 -2.859278165e-34
Float128 -2.330807894e-36
BigFloat(113) -4.558199889e-36
BigFloat(237) 2.318729279e-51
BigFloat(256) 1.124863268e-51
BigFloat(1024) 1.124863268e-51
n = 36 λ[1]
Float32 -1.958887168e-07
Float64 -4.311007920e-16
Double64 -2.255124275e-34
Float128 -3.238481422e-36
BigFloat(113) -3.284703834e-36
BigFloat(237) 6.919587088e-53
BigFloat(256) 3.359762135e-53
BigFloat(1024) 3.359762135e-53
n = 37 λ[1]
Float32 -1.831468239e-07
Float64 -6.535979384e-16
Double64 -2.116415176e-34
Float128 -4.635128397e-36
BigFloat(113) -4.608945882e-36
BigFloat(237) 2.064207412e-54
BigFloat(256) 1.003088162e-54
BigFloat(1024) 1.003088162e-54
n = 38 λ[1]
Float32 -3.336613190e-07
Float64 -4.512357590e-16
Double64 -4.432253018e-34
Float128 -7.250642462e-36
BigFloat(113) -4.329796608e-36
BigFloat(237) 4.617382896e-56
BigFloat(256) 6.155701539e-56
BigFloat(1024) 2.993651511e-56
n = 39 λ[1]
Float32 -2.744432948e-07
Float64 -5.247024262e-16
Double64 -1.290923053e-34
Float128 -3.133177866e-36
BigFloat(113) -3.137829276e-36
BigFloat(237) 1.376487504e-57
BigFloat(256) 1.835102354e-57
BigFloat(1024) 8.931079340e-58
n = 40 λ[1]
Float32 -2.591677912e-07
Float64 -4.441535386e-16
Double64 -4.783821138e-34
Float128 -1.174660257e-36
BigFloat(113) -1.174234562e-36
BigFloat(237) 3.420000952e-59
BigFloat(256) 4.102179409e-59
BigFloat(1024) 2.663517273e-59
Float32は
まとめ
単精度・倍精度・四倍精度(ついでに疑似四倍精度)・八倍精度でフランク行列・ヒルベルト行列の固有値を求めました. フランク行列の場合, 1000×1000程のサイズで単精度では固有値が1桁も正しく計算できなくなることを確かめました. より条件数の大きいヒルベルト行列の場合, 下表のような結果になりました.
| 型 | 固有値が1桁も正しく計算できなくなるサイズ |
|---|---|
| Float32 | |
| Float64 | |
| Double64 | |
| Float128 | |
| BigFloat(113) ≒ Float128 | |
| BigFloat(237) ≒ Float256 | |
| BigFloat(256) |
数値計算では速度至上主義のような価値観に陥りがちですが, 正しい結果が得られないことには速度を追求しても仕方がありません. 四倍精度, 八倍精度は倍精度よりも格段に遅いとされますが, この記事の執筆を通して, 高精度演算のライブラリ開発のモチベーションが理解できました. Juliaの高精度計算のポテンシャルも感じられたと思います.
付録：メモリ消費量
各要素のメモリ消費量をそのまま要素数にかけているだけですので, 実際にはもう少し大きい値になります. 正確な値はvarinfo()で確認できます.
入力
# using Printf
function KMG(byte)
if byte > 10^9
return @sprintf("%5.1f GB", byte/10^9)
elseif byte > 10^6
return @sprintf("%5.0f MB", byte/10^6)
elseif byte > 10^3
return @sprintf("%5.0f KB", byte/10^3)
else
return @sprintf("%5.0f B", byte)
end
end
println("size\t\t", [@sprintf("%9s", "Float$j \t") for j in [32,64,128,256]]...)
for i in union(1,10,100,1000:1000:20000, 30000:10000:100000)
println("$(i)²\t\t", [KMG(i^2*j/8)*"\t" for j in [32,64,128,256]]...)
end
出力
size Float32 Float64 Float128 Float256
1² 4 B 8 B 16 B 32 B
10² 400 B 800 B 2 KB 3 KB
100² 40 KB 80 KB 160 KB 320 KB
1000² 4 MB 8 MB 16 MB 32 MB
2000² 16 MB 32 MB 64 MB 128 MB
3000² 36 MB 72 MB 144 MB 288 MB
4000² 64 MB 128 MB 256 MB 512 MB
5000² 100 MB 200 MB 400 MB 800 MB
6000² 144 MB 288 MB 576 MB 1.2 GB
7000² 196 MB 392 MB 784 MB 1.6 GB
8000² 256 MB 512 MB 1.0 GB 2.0 GB
9000² 324 MB 648 MB 1.3 GB 2.6 GB
10000² 400 MB 800 MB 1.6 GB 3.2 GB
11000² 484 MB 968 MB 1.9 GB 3.9 GB
12000² 576 MB 1.2 GB 2.3 GB 4.6 GB
13000² 676 MB 1.4 GB 2.7 GB 5.4 GB
14000² 784 MB 1.6 GB 3.1 GB 6.3 GB
15000² 900 MB 1.8 GB 3.6 GB 7.2 GB
16000² 1.0 GB 2.0 GB 4.1 GB 8.2 GB
17000² 1.2 GB 2.3 GB 4.6 GB 9.2 GB
18000² 1.3 GB 2.6 GB 5.2 GB 10.4 GB
19000² 1.4 GB 2.9 GB 5.8 GB 11.6 GB
20000² 1.6 GB 3.2 GB 6.4 GB 12.8 GB
30000² 3.6 GB 7.2 GB 14.4 GB 28.8 GB
40000² 6.4 GB 12.8 GB 25.6 GB 51.2 GB
50000² 10.0 GB 20.0 GB 40.0 GB 80.0 GB
60000² 14.4 GB 28.8 GB 57.6 GB 115.2 GB
70000² 19.6 GB 39.2 GB 78.4 GB 156.8 GB
80000² 25.6 GB 51.2 GB 102.4 GB 204.8 GB
90000² 32.4 GB 64.8 GB 129.6 GB 259.2 GB
100000² 40.0 GB 80.0 GB 160.0 GB 320.0 GB
参考文献
- W. L. Frank, Jpn J Ind Appl Math, 6, 378 (1958) https://doi.org/10.1137/0106026
- D. Hilbert, Acta Math., 18, 155 (1894) https://doi.org/10.1007/BF02418278
- 幸谷智紀『LAPACK/BLAS入門』(森北出版, 2016)
追伸
Fortran/LAPACKでの実装例を作りました.
Discussion