[racket] performance problem in math/matrix

2013/1/20 Berthold Bäuml <berthold.baeuml at dlr.de>:
>> Ah! The Numpy example uses floats and not doubles.
>
> On my machine Numpy says it is using float64, hence, double.
>  big1.dtype.name -> float64

Sorry for the confusion. I am used to the Racket documentation
were float means single precision, so I misinterpreted the Numpy
documentation.

Also

> numpy.random.random((5,5))
array([[ 0.04748764,  0.15242073,  0.2366255 ,  0.86926654,  0.16586114],
       [ 0.33117124,  0.7504899 ,  0.43179013,  0.1992083 ,  0.83266119],
       [ 0.25741902,  0.36912007,  0.38273193,  0.29622522,  0.75817261],
       [ 0.49313726,  0.20514653,  0.45329344,  0.964105  ,  0.77459153],
       [ 0.5267172 ,  0.87269101,  0.34530438,  0.38218051,  0.05107181]])

suggested, that the numbers were single precision. It turns out
that random only produces 53 bit random floats.

Anyways, here is a standard matrix multiplication algorithm in
Typed Racket.

#lang typed/racket
(require racket/unsafe/ops
         racket/flonum)

(define-type NF Nonnegative-Fixnum)

(: matrix* :
   (Vectorof Flonum) (Vectorof Flonum) Index Index Index
   -> (Vectorof Flonum))
(define (matrix* A B m p n)
  (define size (* m n))
  (define: v : (Vectorof Flonum) (make-vector size 0.0))
  (define index 0)
  (: loop-i : NF -> (Vectorof Flonum))
  (define (loop-i i)
    (cond [(< i m)
           (loop-j i 0)
           (loop-i (+ i 1))]
          [else v]))
  (: loop-j : NF NF -> Void)
  (define (loop-j i j)
    (when (< j n)
      (loop-k i j 0 0.0)
      (set! index (+ index 1))
      (loop-j i (+ j 1))))
  (: loop-k : NF NF NF Flonum -> Void)
  (define (loop-k i j k sum)
    (define i*p (* i p))
    (define k*n (* k n))
    (cond
      [(< k p)
       (loop-k i j (+ k 1)
               (+ sum
                  (* (unsafe-vector-ref A (+ i*p k))
                     (unsafe-vector-ref B (+ k*n j)))))]
      [else
       (vector-set! v index sum)]))
  (loop-i 0))

(define dim 200)
(define A (build-vector (* dim dim) (λ (_) (random))))
(define B (build-vector (* dim dim) (λ (_) (random))))
(collect-garbage)
(collect-garbage)
(define A*B (time (matrix* A B dim dim dim)))

/Jens Axel