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

On Sat, Nov 06, 2004 at 07:00:23AM +0900, Michael Neumann scribed:
> > 
> > If we had a superfast matrix class, it would give (a) along
> > with decreasing the runtime by a factor of ten or so.  Sounds
> 
> We have a superfast matrix class: NArray (not in standard lib).
> It claimed to be 3x faster than NumPy for some tests and a bit faster
> than Octave for for matrix operations.
> 
> http://www.ir.isas.ac.jp/~masa/ruby/index-e.html

This looks great, btw.  Not having looked at the code
size, I don't know if it's worth having it replace the library
version of 'Matrix' or not.  If we actually just want some
quick speedups to our current matrix class, I offer the 
following bit of code that performs most of the matrix multiplication
in C.  It could be faster if it assumed anything about the contents
of the matrix;  as it is, it's pretty much a straightforward
translation of the current matrix code into C.  It speeds matrix
multiplication up by about an order of magnitude.  Division is
implemented as multiplication by the inverse, so it's also accelerated
somewhat.

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
it to see how well it would work to just make a minor tweak with a
little bit of C code, while preserving the majority of the 1,200
lines of ruby.  The answer seems to be "pretty well!"

It takes 100 iterations of the computer shootout matrix
test from 7.5 seconds (their version) / 13 seconds (ruby 'matrix')
down to about 1.3 seconds.  The perl version takes 3.4 seconds.
The modified python version takes under .2 seconds, in comparison,
as would a Ruby version using NArray...

The happy part about this is that it's all of 103 lines of
C, since all of the wrapper code for matrix handling is still done
in ruby.  To use it, modify matrix.rb to include fastmath after
the class definition, and dispatch to the mult version instead.
Bonus points if you want to make it depend on having fastmath
installed;  should be easy. :)

***************
*** 107,112 ****
--- 107,113 ----
  class Matrix
    @RCS_ID='-$Id: matrix.rb,v 1.11 1999/10/06 11:01:53 keiju Exp keiju $-'
    
+   require 'fastmath'
  #  extend Exception2MessageMapper
    include ExceptionForMatrix
    
***************
*** 465,470 ****
--- 466,472 ----
        return r.column(0)
      when Matrix
        Matrix.Raise ErrDimensionMismatch if column_size != m.row_size
+       return Matrix[ *self.mult(m) ]
      
        rows = (0 .. row_size - 1).collect {
          |i|

  -Dave

=========== ext/matrix/fastmath/extconf.rb =======
require 'mkmf'
create_makefile 'matrix/fastmath'

=========== ext/matrix/fastmath/fastmath.c =======

/*
    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.
    (I hereby assign copyright to matz if he wants it).

*/

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


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

#define FASTMATH_VERSION "0.0.1"

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

/*
=begin
--- mult
    Multiplies two matrices.  Returns an array that must be
    coerced back into a matrix by the caller.
=end
*/

static VALUE
fastmath_mult(VALUE self, VALUE other)
{
	VALUE a, b, c, item, tmp, val, row, m1i, at, bt;
	int arows, acols;
	int brows, bcols;
	int crows, ccols;
	int i, j, k;
	/* Note:  Might be better to call to_a .. this is not speed critical */
	a = rb_iv_get(self, "@rows");
	b = rb_iv_get(other, "@rows");
	if (rb_obj_is_instance_of(other, matclass) != Qtrue) {
		rb_raise(rb_eTypeError, "type of arg must be Matrix");
		return self;
	}
	tmp = rb_funcall(a, id_length, 0);
	arows = NUM2INT(tmp);
	tmp = rb_funcall(b, id_length, 0);
	brows = NUM2INT(tmp);
	if (!arows || !brows) return self;

	item = rb_ary_entry(a, 0);
	if (NIL_P(item) || TYPE(item) != T_ARRAY) {
		raise(rb_eIndexError, "Invalid self matrix");
		return self;
	}
	tmp = rb_funcall(item, id_length, 0);
	acols = NUM2INT(tmp);

	item = rb_ary_entry(b, 0);
	if (NIL_P(item) || TYPE(item) != T_ARRAY) {
		raise(rb_eIndexError, "Invalid self matrix");
		return self;
	}
	tmp = rb_funcall(item, id_length, 0);
	bcols = NUM2INT(tmp);

	crows = arows;
	ccols = bcols;

	c = rb_ary_new2(crows);
	for (i = 0; i < crows; i++) {
		row = rb_ary_new2(ccols);
		m1i = rb_ary_entry(a, i);
		for (j = 0; j < ccols; j++) {
			val = INT2FIX(0);
			for (k = 0; k < ccols; k++) {
				at = rb_ary_entry(m1i, k);
				bt = rb_ary_entry(rb_ary_entry(b, k), j);
				val = rb_funcall(val, id_plus, 1, 
					rb_funcall(at, id_mul, 1, bt));
			}
			rb_ary_store(row, j, val);
		}
		rb_ary_store(c, i, row);
	}

	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("*");
}