.gitignore

Makefile
*.o
*.so
mkmf.log

extconf.rb

require "mkmf"
RbConfig::CONFIG["CPP"] = "g++ -E"
have_library("stdc++")
have_header("boost/multi_array.hpp")
create_makefile("plain")

plain.cpp

#include <ruby.h>
#include <boost/multi_array.hpp>

typedef boost::multi_array<double, 2> plain_t;
typedef boost::multi_array_types::index_range range;

VALUE Plain = Qnil;

static VALUE
constructor(VALUE klass, VALUE x, VALUE y) {
    VALUE self = Data_Wrap_Struct(
        klass, 0, free,
        new plain_t(boost::extents[NUM2INT(x)][NUM2INT(y)])
    );
    VALUE argv[2] = { x, y };
    rb_obj_call_init(self, 2, argv);
    return self;
}

static VALUE
init(VALUE self, VALUE x, VALUE y) {
    return self;
}

static VALUE
shape(VALUE self) {
    plain_t *p;
    Data_Get_Struct(self, plain_t, p);
    VALUE result = rb_ary_new3(
        2,
        INT2NUM(p->shape()[0]),
        INT2NUM(p->shape()[1])
    );
    return result;
}

static VALUE
at(VALUE self, VALUE x, VALUE y) {
    plain_t *p;
    Data_Get_Struct(self, plain_t, p);
    int i = NUM2INT(x);
    int j = NUM2INT(y);
    if (!(i < p->shape()[0])) rb_raise(rb_eException, "over x range");
    if (!(j < p->shape()[1])) rb_raise(rb_eException, "over y range");
    return rb_float_new((*p)[i][j]);
}

static VALUE
set(VALUE self, VALUE x, VALUE y, VALUE v) {
    plain_t *p;
    Data_Get_Struct(self, plain_t, p);
    int i = NUM2INT(x);
    int j = NUM2INT(y);
    if (!(i < p->shape()[0])) rb_raise(rb_eException, "over x range");
    if (!(j < p->shape()[1])) rb_raise(rb_eException, "over y range");
    (*p)[i][j] = NUM2DBL(v);
    return v;
}

static VALUE
elems(VALUE self) {
    plain_t *p;
    Data_Get_Struct(self, plain_t, p);
    return INT2NUM(p->num_elements());
}

static VALUE
msize(VALUE self) {
    plain_t *p;
    Data_Get_Struct(self, plain_t, p);
    return INT2NUM(sizeof(double)*p->num_elements());
}

#define RANGE_BEG(r) (RSTRUCT(r)->as.ary[0])
#define RANGE_END(r) (RSTRUCT(r)->as.ary[1])
static VALUE
slice(VALUE self, VALUE xrange, VALUE yrange) {
    plain_t *p;
    Data_Get_Struct(self, plain_t, p);
    int i = FIX2INT(RANGE_BEG(xrange));
    int I = FIX2INT(RANGE_END(xrange));
    int j = FIX2INT(RANGE_BEG(yrange));
    int J = FIX2INT(RANGE_END(yrange));
    VALUE result = Data_Wrap_Struct(
        rb_singleton_class(self), 0, free,
        new plain_t( (*p)[ boost::indices[range(i,I)][range(j,J)] ] )
    );
    return result;
}

static VALUE
each_with_index(int argc, VALUE *argv, VALUE self) {
    VALUE block = Qnil;
    rb_scan_args(argc, argv, "0&", &block);
    if (!RTEST(block)) {
        rb_raise(rb_eArgError, "a block is required");
    }

    plain_t *p;
    Data_Get_Struct(self, plain_t, p);
    for (int i = 0; i < p->shape()[0]; i++) {
        for (int j = 0; j < p->shape()[1]; j++) {
            rb_funcall(block, rb_intern("call"), 3, DBL2NUM((*p)[i][j]), INT2NUM(i), INT2NUM(j));
        }
    }

    return self;
}

#define RMETHOD VALUE(*)(...)
extern "C" void Init_plain() {
    Plain = rb_define_class("Plain", rb_cObject);
    rb_define_singleton_method(Plain,   "new",              (RMETHOD)constructor,        2);
    rb_define_method(Plain,             "initialize",       (RMETHOD)init,               2);
    rb_define_method(Plain,             "shape",            (RMETHOD)shape,              0);
    rb_define_method(Plain,             "[]",               (RMETHOD)at,                 2);
    rb_define_method(Plain,             "at",               (RMETHOD)at,                 2);
    rb_define_method(Plain,             "[]=",              (RMETHOD)set,                3);
    rb_define_method(Plain,             "set",              (RMETHOD)set,                3);
    rb_define_method(Plain,             "elems",            (RMETHOD)elems,              0);
    rb_define_method(Plain,             "msize",            (RMETHOD)msize,              0);
    rb_define_method(Plain,             "slice",            (RMETHOD)slice,              2);
    rb_define_method(Plain,             "each_with_index",  (RMETHOD)each_with_index,   -1);
}

test.rb

$:.unshift(".")
require "plain"

p = Plain.new(4,4)
x, y = p.shape

puts p.elems
puts p.msize

p.instance_eval do
    each_with_index do |v, i, j|
        self[i,j] = i+j
    end
end

p.each_with_index do |v, i, j|
    puts "#{i},#{j}: #{v}"
end

s = p.slice(0..2,0..2)

s.each_with_index do |v, i, j|
    puts "#{i},#{j}: #{v}"
end