[racket] performance problem in math/matrix

Recall that last.  Always measure at the command line <sigh>.  I'm seeing
FlVector at 30% faster.

#lang typed/racket

(provide main)

(require racket/unsafe/ops
 racket/fixnum
         racket/flonum)

(define-type NF Nonnegative-Fixnum)

(: matrix* :
   FlVector FlVector Index Index Index
   -> FlVector)
(define (matrix* A B m p n)
  (define size (* m n))
  (define: v : FlVector (make-flvector size 0.0))
  (define index 0)
  (: loop-i : NF -> FlVector)
  (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-flvector-ref A (+ i*p k))
(unsafe-flvector-ref B (+ k*n j)))))]
     [else
      (unsafe-flvector-set! v index sum)]))
  (loop-i 0))

(: mkFlVector (Fixnum (-> Flonum) -> FlVector))
(define (mkFlVector sz init)
  (let ((v (make-flvector sz)))
    (do ((i #{0 : Fixnum} (+ i 1)))
((>= i sz) v)
      (unsafe-flvector-set! v i (init)))))

(define: dim : Index 300)
(define: A : FlVector (mkFlVector (fx* dim dim) random))
(define: B : FlVector (mkFlVector (fx* dim dim) random))

(define (main)
  (collect-garbage)
  (collect-garbage)
  (let ((m (time (matrix* A B dim dim dim))))
    (void)))



On Sun, Jan 20, 2013 at 4:25 PM, Ray Racine <ray.racine at gmail.com> wrote:

> Modification of the above to use FlVector did not show a difference.
>
>
> On Sun, Jan 20, 2013 at 2:35 PM, Jens Axel Søgaard <jensaxel at soegaard.net>wrote:
>
>> 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
>>
>> ____________________
>>   Racket Users list:
>>   http://lists.racket-lang.org/users
>>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.racket-lang.org/users/archive/attachments/20130120/517a83d1/attachment-0001.html>