Hi! I would like to propose a function that could be used as alias for scipy.matrixmultiply (code attached in 'matrixmultiply.py'). The function delegate the work to scipy.dot if the arguments it receives aren't 2D, otherwise it uses the BLAS routine 'gemm'. If BLAS fails, an appropriate error message is raised. The function should now be backward-compatible and is sensibly faster than the original one for matrices larger than 40x40. Here are the benchmarks (the code to perform the benchmarks is included in 'matmult_bench.py'): Multiplying matrices of type f ================================== | contiguous | non-contiguous ---------------------------------------------- size | scipy | matmult | scipy | matmult 20 | 0.11 | 0.24 | 0.11 | 0.24 (secs for 2000 calls) 50 | 0.99 | 0.60 | 1.01 | 0.60 (secs for 2000 calls) 75 | 3.10 | 1.44 | 3.09 | 1.46 (secs for 2000 calls) 100 | 3.64 | 0.80 | 3.62 | 0.83 (secs for 1000 calls) 500 | 5.06 | 0.25 | 5.08 | 0.26 (secs for 4 calls) 1000 | 28.05 | 0.94 | 27.96 | 0.94 (secs for 2 calls) Multiplying matrices of type d ================================== | contiguous | non-contiguous ---------------------------------------------- size | scipy | matmult | scipy | matmult 20 | 0.10 | 0.26 | 0.10 | 0.26 (secs for 2000 calls) 50 | 1.01 | 0.68 | 1.02 | 0.69 (secs for 2000 calls) 75 | 3.44 | 1.87 | 3.48 | 1.91 (secs for 2000 calls) 100 | 3.59 | 1.55 | 3.59 | 1.49 (secs for 1000 calls) 500 | 8.02 | 0.50 | 8.01 | 0.52 (secs for 4 calls) 1000 | 36.31 | 1.86 | 36.15 | 1.83 (secs for 2 calls) Regards, Pietro.