Previous ToC Next

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. まとめ

  1. 多次元配列を作った。

20.3. 参考資料

特になし。Crystal 言語の Pointer, Slice, Macro 等をみておこう。
Previous ToC Next