From: "David G. Andersen" <dga@...>
Date: 2004-11-06T10:25:39+09:00
Subject: Re: Ruby is a slow performer

On Sat, Nov 06, 2004 at 09:06:38AM +0900, David G. Andersen scribed:

> I offer the following bit of code that performs most of the matrix
> multiplication in C.

> The code is attached at the end of this message.  It should probably
> be extended with a few more sanity checks by someone who knows what
> they're doing.  I haven't tested it extensively.  I mostly wrote
> [...]

  I've added the requisite sanity checks to the code and
converted it to directly access the array elements.  It's about
10% faster, a bit more readable, and decidedly safer. :).  As
before, it preserves the semantics of matrix.rb, so it should
be drop-in compatible.

  (Sorry for the 100 line email)

  -Dave

/*
    fastmath.c -- matrix manipulation core
  
    Copyright (c) 2004 David Andersen <dga@pobox.com>
  
    This library is free software.
    You can distribute/modify this program under the same terms of ruby.
*/

#include "ruby.h"
#include <stdio.h>

VALUE matclass;  /* Set in Init */
static int id_length, id_plus, id_mul;

#define FASTMATH_VERSION "0.0.1"

/*
=begin
= Fastmath methods.
= end
 */

/*
=begin
--- mult
    Multiplies two matrices
=end
*/

static VALUE
fastmath_mult(VALUE self, VALUE other)
{
	VALUE a, b, c, val, row, m1i, at, bt;
	int arows, brows, bcols, i, j, k;
	
	if (rb_obj_is_instance_of(other, matclass) != Qtrue) {
		rb_raise(rb_eTypeError, "type of arg must be Matrix");
	}

	a = rb_iv_get(self, "@rows");
	b = rb_iv_get(other, "@rows");
	if (NIL_P(a) || NIL_P(b) || TYPE(a) != T_ARRAY || TYPE(b) != T_ARRAY) {
		rb_raise(rb_eTypeError, "invalid matrix arguments");
	}
	arows = RARRAY(a)->len;
	brows = RARRAY(b)->len;
	if (!arows || !brows) {
		rb_raise(rb_eIndexError, "Zero-length matrix");
	}

	row = RARRAY(b)->ptr[0];
	if (NIL_P(row) || TYPE(row) != T_ARRAY) {
		rb_raise(rb_eIndexError, "Invalid other matrix");
	}
	bcols = RARRAY(row)->len;

	for (i = 0; i < brows; i++) {
		row = RARRAY(b)->ptr[i];
		if (TYPE(row) != T_ARRAY || RARRAY(row)->len != arows) {
			rb_raise(rb_eIndexError, "Invalid other row");
		}
	}

	c = rb_ary_new2(arows);
	for (i = 0; i < arows; i++) {
		row = rb_ary_new2(bcols);
		m1i = RARRAY(a)->ptr[i];
		if (TYPE(m1i) != T_ARRAY || RARRAY(m1i)->len != brows) {
			rb_raise(rb_eIndexError, "Invalid self row len");
		}
		for (j = 0; j < bcols; j++) {
			val = INT2FIX(0);
			for (k = 0; k < bcols; k++) {
				at = RARRAY(m1i)->ptr[k];
				bt = RARRAY(RARRAY(b)->ptr[k])->ptr[j];
				val = rb_funcall(val, id_plus, 1, 
					rb_funcall(at, id_mul, 1, bt));
			}
			RARRAY(row)->ptr[j] = val;
		}
		RARRAY(row)->len = bcols;
		rb_ary_store(c, i, row);
	}

	if (OBJ_TAINTED(self) || OBJ_TAINTED(other)) OBJ_TAINT(c);
	return c; /* remember to coerce to matrix in caller */
}

void
Init_fastmath()
{
	matclass = rb_const_get(rb_cObject, rb_intern("Matrix"));
	rb_define_method(matclass, "mult", fastmath_mult, 1);
	id_length = rb_intern("length");
	id_plus = rb_intern("+");
	id_mul = rb_intern("*");
}