2013/1/20 Berthold Bäuml <berthold.bae...@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
