David G. Andersen
11/6/2004 1:26:00 AM
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("*");
}