20. 多次元配列
赤木 ：色々忙しいと思わず逃避行動にでているんだけど(だれが？)、Crystal
で多次元配列ってどうなのかしら？
学生C：どうなのかしら？といわれましても、、、ないですよね？ Ruby にな
いし。もちろん、配列の配列つくれば多次元配列だから、速度とか気にしなけ
ればそれでいいし、Ruby はそのレベルで速度とか気にするなら使えないし。
赤木 ：でもほら、田中さんの narray とか、あと Python の numpy とか、ラ
イブラリ使えば速いみたいなのがあるじゃない?そういう感じのが Crystal
にもあってもいいと思わない?
学生C：あればいいんでしょうけど、えーと、具体的に何が欲しいと?
赤木 ：要するに、普通に多次元配列宣言して、その要素にアクセスしたいわ
け。例えば Fortran なら
subroutine matmul(a,b,c,n)
integer*4 n
real*8 a(n,n), b(n,n), c(n,n)
...
end
みたいな感じで配列かけて、多次元配列のサイズを関数に渡せるし、 C も C99 ならできるじゃない? C++ は未来永劫これができそうにないけど。
学生C：C++ なら、サイズと、実際のデータエリアへのポインタもった構造体とか作りますよね？
赤木 ：まあそうね。
学生C：で、Crystal でどうしたいと?
赤木 ：どっちいうとC++ でやるみたいにサイズと、実際のデータエリアへの
ポインタもったクラス作って、でも、なるべくそれを自然に使うみたいな。例
えば Ruby の NArray だと
x =Narray.float(10)
xy =Narray.float(3,3)
とかで、長さ10の1次元配列とか、3x3の2次元配列ができて、2次元配列のほう
は
xy[1,2] = 3
pp! xy[1,2]
みたいな感じで、 [1][2] みたいな書き方じゃなくてコンマ区切りでいいのと、
あと、Fortran や Cのプログラムに配列のポインタを渡せる、言い換えると、
内部で Ruby のArrayを使ってなくて、連続アドレス領域をとってデータがお
いてあるわけ。だから、これ使うと色々なライブラリが簡単に使えるのね。
学生C：えー、でも、Array じゃないんだとどうするんですか？全部 C で書くとか？そんなのできませんよ。
赤木 ：あ、Crystal その辺便利で、Pointer ってのがあるの。それで C みた
いに malloc とかできるの。で、面白いのは、 free はないのね。プログラム
で使わなくなったらガーベージコレクタが自動的にケアするから。なので、
こういうの使っても C/C++ みたいなややこしいバグとかは起こらないの。
学生C：まあそれはよさそうですが、、、多次元配列を a[2,3] みたいに書くっ
て、 [] を2つ引数とる関数として定義するみたいなのがいりますよね？それ
できるんですか？
赤木 ： def [](i,j) とか def[]=(i,j,val) でできるみたいよ。
学生C： それならなんとかなりそうですね。Ruby の narray とか Python の
numpy と違って、ちゃんとコンパイラだから演算自体は Crystal で書いて別
に遅くないはずだし。ちょっとかいてみます。
学生C： えーと、Pointer というのがあるんですね。で、赤木さんがいってたみた
いに malloc できると。あれ、Slice というのもあると。
While a pointer is unsafe because no bound checks are performed
when reading from and writing to it, reading from and writing to a
slice involve bound checks. In this way, a slice is a safe
alternative to Pointer.
えーと、範囲チェックが必要ならこっちを使えということですね。でも、いつもチェックして
たら遅い気がします。コンパイラオプションで切換えるとかですかね？まずはそういうの考
えないで最低限だと、2次元までなら
1: class Narray_F64
2: property :data
3: def initialize(nx : Int32,ny : Int32 = 1)
4: @data = Pointer(Float64).malloc(nx*ny)
5: @nx=nx
6: @ny=ny
7: end
8: def [](i) @data[i] end
9: def [](i,j) @data[i*@ny+j] end
10: def []=(i,x) @data[i]=x end
11: def []=(i,j,x) @data[i*@ny+j]=x end
12: end
3次元以上にしたければ @nz とかもつければいいと。initialize 最初は配列
の次元数毎に書いたんだけど、中身は要するに使わない次元を1にするだけだっ
たのでデフォルト引数ですむ。
一応これ呼ぶプログラム作ると、
1:require "./narray0.cr"
2:
3:x =Narray_F64.new(10)
4:x[1]=2.0
5:pp! x
6:xy =Narray_F64.new(4,4)
7:16.times{|i| xy[i]=i.to_f64}
8:xy[2,3]=-1.0
9:4.times{|i|
10: 4.times{|j|print " ",xy[i,j]}
11: print "\n"
12:}
13:pp! xy
これ実行してみると
gravity> crystal narray0test.cr
x # => #<Narray_F64:0x7ff2a772ee60
@data=Pointer(Float64)@0x7ff2a7739f60,
@nx=10,
@ny=1>
0.0 1.0 2.0 3.0
4.0 5.0 6.0 7.0
8.0 9.0 10.0 -1.0
12.0 13.0 14.0 15.0
xy # => #<Narray_F64:0x7ff2a772e940
@data=Pointer(Float64)@0x7ff2a773dea0,
@nx=4,
@ny=4>
値もはいっているしそれっぽいですね。
赤木 ：そうね。でもこれだと Float64 しか使えないから、Float32 も使える
ようにできる？
学生C：2つ書けばよくないですか？
赤木 ：まあ、そういうことをしなくていいための機能が C++ とかでいうとこ
ろのテンプレートで、 Crystal だと Generics というのね。多項式のところ
で、 Array から新しいクラス作るとか、あと MathVector でも使ってるけど、
class Foo
の代わりに
class Foo(T)
として、別にTでなくてもいいはずだけど、クラスが引数として型を持つよう
にできるのね。Array も
class Array(T)
で、 new は
Array(Float64).new(...)
とかするのと同じわけ。
学生C：そうすると、
class Narray(T)
で、
def initialize(nx : Int32,ny : Int32 = 1)
@data = Pointer(T).malloc(nx*ny)
@nx=nx
@ny=ny
end
ですね。これを narray1.cr にして
1:require "./narray1.cr"
2:
3:x =Narray(Float64).new(10)
4:x[1]=2.0
5:pp! x
6:xy =Narray(Float32).new(4,4)
7:16.times{|i| xy[i]=i.to_f32}
8:xy[2,3]=-1.0_f32
9:4.times{|i|
10: 4.times{|j|print " ",xy[i,j]}
11: print "\n"
12:}
13:pp! xy
これ実行してみると
gravity> crystal narray1test.cr
x # => #<Narray(Float64):0x7f60e559fe60
@data=Pointer(Float64)@0x7f60e55aaf60,
@nx=10,
@ny=1>
0.0 1.0 2.0 3.0
4.0 5.0 6.0 7.0
8.0 9.0 10.0 -1.0
12.0 13.0 14.0 15.0
xy # => #<Narray(Float32):0x7f60e559f9a0
@data=Pointer(Float32)@0x7f60e55b0f50,
@nx=4,
@ny=4>
なるほど、大丈夫ですね。ちゃんと x と xy で違う型ですね。
あとなんでしたっけ？範囲チェックですね？
赤木 ：そう。Slice 使ってみて。
学生C：単に使うだけなら
def initialize(nx : Int32,ny : Int32 = 1)
@data = Slice(T).new(Pointer(T).malloc(nx*ny), nx*ny)
@nx=nx
@ny=ny
end
ですかね。これで @data は Pointer じゃなくて Slice になるんですよね？
赤木 ：そのはず。なんかテスト書いて範囲外アクセスとかしてみて。
学生C：ではこんなので
1:require "./narray2.cr"
2:
3:x =Narray(Float64).new(10)
4:x[1]=2.0
5:pp! x
6:xy =Narray(Float64).new(4,4)
7:16.times{|i| xy[i]=i.to_f64}
8:xy[2,3]=-1.0
9:4.times{|i|
10: 4.times{|j|print " ",xy[i,j]}
11: print "\n"
12:}
13:pp! x[11]
gravity> crystal narray2test.cr
x # => #<Narray(Float64):0x7f794d130f00
@data=Slice[0.0, 2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
@nx=10,
@ny=1>
0.0 1.0 2.0 3.0
4.0 5.0 6.0 7.0
8.0 9.0 10.0 -1.0
12.0 13.0 14.0 15.0
Unhandled exception: Index out of bounds (IndexError)
from /usr/share/crystal/src/indexable.cr:903:8 in '[]'
from narray2.cr:8:17 in '[]'
from narray2test.cr:13:1 in '__crystal_main'
from /usr/share/crystal/src/crystal/main.cr:115:5 in 'main_user_code'
from /usr/share/crystal/src/crystal/main.cr:101:7 in 'main'
from /usr/share/crystal/src/crystal/main.cr:127:3 in 'main'
from /lib/x86_64-linux-gnu/libc.so.6 in '__libc_start_main'
from /home/makino/.cache/crystal/crystal-run-narray2test.tmp in '_start'
from ???
x[11] # =>
学生C：おお、ちゃんと x[11]は範囲外だとでますね。narray1.cr を使うと
1:require "./narray1.cr"
2:
3:x =Narray(Float64).new(10)
4:x[1]=2.0
5:pp! x
6:xy =Narray(Float64).new(4,4)
7:16.times{|i| xy[i]=i.to_f64}
8:xy[2,3]=-1.0
9:4.times{|i|
10: 4.times{|j|print " ",xy[i,j]}
11: print "\n"
12:}
13:pp! x[11]
gravity> crystal narray2test1.cr
x # => #<Narray(Float64):0x7f6165ad9e60
@data=Pointer(Float64)@0x7f6165ae4f60,
@nx=10,
@ny=1>
0.0 1.0 2.0 3.0
4.0 5.0 6.0 7.0
8.0 9.0 10.0 -1.0
12.0 13.0 14.0 15.0
x[11] # => 0.0
なんか知らないけど0が戻りますね。
赤木 ：これ、Slice にしておくと、範囲チェックができる他に色々いいこと
があるのね。ドキュメントみてもらうとわかるけど、Slice にはArray にある
ような便利な関数がほとんどあるから、each とか map! とか使えるようにな
るわけ。map もあるけどこれは配列になっちゃうからそのままではちょっとね。
学生C：sort! とかもできるわけですか？
赤木 ：できるわよ。ただ、今のだと、 x.@data.sort! とかでちょっと見苦し
いわね。method_missing 使って
macro method_missing(call)
@data.{{call}}
end
としておけば、 Narray に対して知らない関数がきたら @data に投げる、と
できるから、これで sort! とか map とかなんでも使えるわね。
学生C：ちょっと試しに、、、
1:require "./narray2.cr"
2:
3:xy =Narray(Float64).new(4,4)
4:16.times{|i| xy[i]=i.to_f64}
5:xy[2,3]=-1.0
6:4.times{|i|
7: 4.times{|j|print " ",xy[i,j]}
8: print "\n"
9:}
10:xy.sort!
11:4.times{|i|
12: 4.times{|j|print " ",xy[i,j]}
13: print "\n"
14:}
15:
を実行すると
gravity> crystal narray2test2.cr
0.0 1.0 2.0 3.0
4.0 5.0 6.0 7.0
8.0 9.0 10.0 -1.0
12.0 13.0 14.0 15.0
-1.0 0.0 1.0 2.0
3.0 4.0 5.0 6.0
7.0 8.0 9.0 10.0
12.0 13.0 14.0 15.0
おおお、確かに実行できてソートもされますね。2次元配列こんなふうにソー
トすることはあんまりないとは思いますが、map! とかはいいかもですね。
赤木 ：で、これでいいんだけど、「コンパイルオプションによって Pointer
か Slice か切換える」「でも、 Slice の機能は使えるようにする」ってでき
る？
学生C：まずポインタで領域とって、それに対してSlice作成だから、添字での
アクセスでどっち使うかですね。コンパイルオプションで選ぶのは、少なくと
もクラス全体をマクロにすればよかったような、、、なので、こんなのですね。
1:# narray.cr
2:#
3:# multi-dimensional array similar to NUMO:Narray
4:#
5:# Copyright 2020- J. Makino
6:#
7:macro use_narray(bound_check = false)
8:class Narray(T)
9: property :data
10: def Narray.range_check?
11: {% if bound_check %}
12: true
13: {% else %}
14: false
15: {% end %}
16: end
17: def initialize(@nx : Int32, @ny : Int32 = 1, @nz : Int32 = 1)
18: n=@nx*@ny*@nz
19: @data = Slice(T).new(Pointer(T).malloc(n), n)
20: @datasl = @data
21: {% if bound_check %}
22: @data = @datasl
23: {% end %}
24: end
25: def [](i) @data[i] end
26: def [](i,j) @data[i*@ny+j] end
27: def [](i,j,k) @data[(i*@ny+j)*@nz+k] end
28: def []=(i,x) @data[i]=x end
29: def []=(i,j,x) @data[i*@ny+j]=x end
30: def []=(i,j,k, x) @data[(i*@ny+j)*@nz+k]=x end
31: macro method_missing(call)
32: @datasl.\{{call}}
33: end
34: def each_index
35: @nx.times{|i|
36: @ny.times{|j|
37: @nz.times{|k|
38: yield i,j,k
39: }
40: }
41: }
42: end
43: def each_index(idim)
44: if idim == 0
45: @nx.times{|i| yield i}
46: elsif idim == 1
47: @ny.times{|j|yield j}
48: elsif idim == 2
49: @nz.times{|k|yield k}
50: end
51: end
52:end
53:end
54:
55:{% if flag?(:range_check) %}
56: use_narray(true)
57: {% else %}
58: use_narray(false)
59:{% end %}
まず、クラス定義全体を
macro use_narray(bound_check = false)
...
end
で囲って、
use_narray(true)
で範囲チェックあり、 false か引数なしならチェックなしのコードがコンパ
イル時に生成されるようにします。そのなかでは
{% if bound_check %}
foo
{% else %}
bar
{% end %}
でどっちのコードになるか選べるわけです。チェックするかどうかは、
initializeのほうで
{% if bound_check %}
@data = Slice(T).new(Pointer(T).malloc(nx*ny*nz), nx*ny*nz)
@datasl = @data
{% else %}
@data = Pointer(T).malloc(nx*ny*nz)
@datasl = Slice(T).new(@data, nx*ny*nz)
{% end %}
として、チェックするなら data も datasl も Slice、しない時は data は
Pointer のままで、としました。
@data = Pointer(T).malloc(nx*ny*nz)
@datasl = Slice(T).new(@data, nx*ny*nz)
{% if bound_check %}
@data = @datasl
{% end %}
のほうが格好よいですがまあ。あと、range_check? というクラス関数も作っ
てみました。で、最後に
{% if flag?(:range_check) %}
use_narray(true)
{% else %}
use_narray(false)
{% end %}
で、コンパイル時に -Drange_check されるかどうかでどっちかのマクロが生
成されます。
1:require "narray"
2:x =Narray(Float64).new(10)
3:x[1]=2.0
4:xy =Narray(Float64).new(4,4)
5:16.times{|i| xy[i]=i.to_f64}
6:xy[2,3]=-1.0
7:4.times{|i|
8: 4.times{|j|print " ",xy[i,j]}
9: print "\n"
10:}
11:xy.sort!
12:pp! xy
13:pp! Narray.range_check?
14:pp! x[11]
を範囲チェックありで実行すると
gravity> crystal run -Drange_check narray_test.cr
0.0 1.0 2.0 3.0
4.0 5.0 6.0 7.0
8.0 9.0 10.0 -1.0
12.0 13.0 14.0 15.0
xy # => #<Narray(Float64):0x7f0fdfae8680
@data=
Slice[-1.0,
0.0,
1.0,
2.0,
3.0,
4.0,
5.0,
6.0,
7.0,
8.0,
9.0,
10.0,
12.0,
13.0,
14.0,
15.0],
@datasl=
Slice[-1.0,
0.0,
1.0,
2.0,
3.0,
4.0,
5.0,
6.0,
7.0,
8.0,
9.0,
10.0,
12.0,
13.0,
14.0,
15.0],
@nx=4,
@ny=4,
@nz=1>
Narray.range_check? # => true
Unhandled exception: Index out of bounds (IndexError)
from /usr/share/crystal/src/indexable.cr:903:8 in '[]'
from lib/narray/src/narray.cr:56:1 in '[]'
from narray_test.cr:14:1 in '__crystal_main'
from /usr/share/crystal/src/crystal/main.cr:115:5 in 'main_user_code'
from /usr/share/crystal/src/crystal/main.cr:101:7 in 'main'
from /usr/share/crystal/src/crystal/main.cr:127:3 in 'main'
from /lib/x86_64-linux-gnu/libc.so.6 in '__libc_start_main'
from /home/makino/.cache/crystal/crystal-run-narray_test.tmp in '_start'
from ???
x[11] # =>
なしだと
gravity> crystal run narray_test.cr
0.0 1.0 2.0 3.0
4.0 5.0 6.0 7.0
8.0 9.0 10.0 -1.0
12.0 13.0 14.0 15.0
xy # => #<Narray(Float64):0x7f4ac75dced0
@data=Pointer(Float64)@0x7f4ac75e0f30,
@datasl=
Slice[-1.0,
0.0,
1.0,
2.0,
3.0,
4.0,
5.0,
6.0,
7.0,
8.0,
9.0,
10.0,
12.0,
13.0,
14.0,
15.0],
@nx=4,
@ny=4,
@nz=1>
Narray.range_check? # => false
x[11] # => 0.0
で、ちゃんとコンパイルオプションで @data の型が変わって、範囲チェック
の結果も変わってるのがわかります。
赤木 ：なるほど。最低限のことはこれでできる感じね。あとは、C/C++ に比
べて速度は、とか、 SIMD 化はどうするの？とかまだ色々あるけど、それはお
いおいね。
20.1. 課題
この章は息抜きなので課題はありません。
20.2. まとめ
-
多次元配列を作った。
20.3. 参考資料
特になし。Crystal 言語の Pointer, Slice, Macro 等をみておこう。