[racket] performance problem in math/matrix

Dear Neil,

first I want to thank you for the great work the math package is. We hope to use it heavily in our "Racket in Robotics" environment we are currently developing.

I tried a little test with mildly large Flonum matrices using math/matrix in Typed Racket. Two 1000x1000 dimensional matrices should be multiplied with (array-strictness #t). Running this in DrRacket results in an "out of memory" error (also having set the memory limit to 1GB) and running it with racket on the command line makes the computer start to heavily swap. Where is this high memory consumption coming from? I would have expected that it will not take more than 3*8*10^6 bytes = 24MB to hold the three involved matrices in memory.

I also tried matrix* with two 100x100 dimensional matrices. This finishes but takes about 1s for the multiplication only which is quite long. For comparison, I tried the same with Mathematica, Python/NumPy and C++ using the Eigen-libary (I know, that this is not quite fair as the former two also use a optimized C-library as their backend, but it gives a baseline of what performance is possible)

				100x100	1000x1000				
Typed Racket	1000ms		(out of memory)	
Mathematica		2ms			250ms
Python/NumPy	1ms			240ms
C++/Eigen		<2ms		<300ms

Here the configuration details of the MacBook I used:
Racket 5.3.1.900
OS X 10.8.2
2.53 GHz Intel Core 2 Duo with 4GB memory

And here the code snippets:


Typed Racket:
------------------------------------

#lang typed/racket/base
(require math)

(array-strictness #t)

(define dim 1000)

(define big1 (build-matrix dim dim (lambda (i j) (random))))

(define big2 (build-matrix dim dim (lambda (i j) (random))))

(define res (time (matrix* big1 big2)))


Mathematica:
------------------------------------------
dim = 1000;
big1 = Array[Random[] &, {dim, dim}];
big2 = Array[Random[] &, {dim, dim}];

Timing[res = big1.big2 ;]


Python/NumPy
------------------------------------------
import numpy
import time

dim = 100
big1 = numpy.random.random((dim,dim))
big2 = numpy.random.random((dim,dim))

start = time.time()
res = numpy.dot(big1, big2)
stop = time.time()

print(stop - start)


C++/Eigen (compiled with g++ -O2)
---------------------------------------------

#include <Eigen/Dense>
using Eigen::MatrixXd;
int main()
{
  const int dim = 100;
  MatrixXd big1(dim, dim);
  MatrixXd big2(dim, dim);

  MatrixXd res(dim, dim);

  res = big1 * big2;
}

-- 
-----------------------------------------------------------------------
Berthold Bäuml 
DLR, Robotics and Mechatronics Center (RMC)
Münchner Str. 20, D-82234 Wessling
Phone +49 8153 282489
http://www.robotic.de/Berthold.Baeuml