suzanne.soy/racket
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Rewrote spectral-norm benchmark based on Lua version.

Browse Source
This commit is contained in:
Sam Tobin-Hochstadt 2010-07-01 11:34:36 -04:00
parent 75bd26d326
commit c0a9704ebc
1 changed files with 44 additions and 63 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

107
collects/tests/racket/benchmarks/shootout/spectralnorm-unsafe.rkt
View File

@ -1,68 +1,49 @@
#lang racket/base
;; The Computer Language Benchmarks Game
;; http://shootout.alioth.debian.org/
;; Translated from Mike Pall's Lua version.
;; Translated directly from the C# version, which was:
;; contributed by Isaac Gouy
(require racket/cmdline racket/trace racket/contract
racket/unsafe/ops racket/flonum)
(require racket/cmdline
racket/require (for-syntax racket/base)
(rename-in
(filtered-in
(lambda (name) (regexp-replace #rx"unsafe-" name ""))
racket/unsafe/ops)
[fx->fl ->fl])
(only-in racket/flonum make-flvector))
(define (Approximate n)
(let ([u (make-flvector n 1.0)]
[v (make-flvector n 0.0)])
;; 20 steps of the power method
(for ([i (in-range 10)])
(MultiplyAtAv n u v)
(MultiplyAtAv n v u))
;; B=AtA A multiplied by A transposed
;; v.Bv /(v.v) eigenvalue of v
(let loop ([i 0][vBv 0.0][vv 0.0])
(if (= i n)
(flsqrt (fl/ vBv vv))
(let ([vi (flvector-ref v i)])
(loop (add1 i)
(fl+ vBv (fl* (flvector-ref u i) vi))
(fl+ vv (fl* vi vi))))))))
;; return element i,j of infinite matrix A
(define (A i j)
(fl/ 1.0 (fl+ (fl* (->fl (+ i j))
(fl/ (->fl (+ i (+ j 1))) 2.0))
(->fl (+ i 1)))))
;; multiply vector v by matrix A
(define (MultiplyAv n v Av)
(for ([i (in-range n)])
(flvector-set! Av i
(for/fold ([r 0.0])
([j (in-range n)])
(fl+ r (fl* (A i j) (flvector-ref v j)))))))
;; multiply vector v by matrix A transposed
(define (MultiplyAtv n v Atv)
(for ([i (in-range n)])
(flvector-set! Atv i
(for/fold ([r 0.0])
([j (in-range n)])
(fl+ r (fl* (A j i) (flvector-ref v j)))))))
;; multiply vector v by matrix A and then by matrix A transposed
(define (MultiplyAtAv n v AtAv)
(let ([u (make-flvector n 0.0)])
(MultiplyAv n v u)
(MultiplyAtv n u AtAv)))
(printf "~a\n"
(real->decimal-string
(Approximate (command-line #:args (n) (string->number n)))
9))
(let* ([A (lambda (i j)
(let ([ij (unsafe-fx+ i j)])
(unsafe-fl/ 1.0 (unsafe-fl+ (unsafe-fl* (unsafe-fl* (unsafe-fx->fl ij)
(unsafe-fx->fl (unsafe-fx+ ij 1)))
0.5)
(unsafe-fx->fl (unsafe-fx+ i 1))))))]
[Av
(lambda (x y N)
(for ([i (in-range N)])
(unsafe-flvector-set!
y i
(let L ([a 0.0] [j 0])
(if (unsafe-fx= j N) a
(L (unsafe-fl+ a (unsafe-fl* (unsafe-flvector-ref x j) (A i j)))
(unsafe-fx+ j 1)))))))]
[Atv
(lambda (x y N)
(for ([i (in-range N)])
(unsafe-flvector-set!
y i
(let L ([a 0.0] [j 0])
(if (unsafe-fx= j N) a
(L (unsafe-fl+ a (unsafe-fl* (unsafe-flvector-ref x j) (A j i)))
(unsafe-fx+ j 1)))))))]
[AtAv (lambda (x y t N) (Av x t N) (Atv t y N))]
[N (command-line #:args (n) (string->number n))]
[u (make-flvector N 1.0)]
[v (make-flvector N)]
[t (make-flvector N)])
(for ([i (in-range 10)])
(AtAv u v t N)
(AtAv v u t N))
(displayln (real->decimal-string
(unsafe-flsqrt
(let L ([vBv 0.0] [vv 0.0] [i 0])
(if (unsafe-fx= i N) (unsafe-fl/ vBv vv)
(let ([ui (unsafe-flvector-ref u i)] [vi (unsafe-flvector-ref v i)])
(L (unsafe-fl+ vBv (unsafe-fl* ui vi))
(unsafe-fl+ vv (unsafe-fl* vi vi))
(unsafe-fx+ i 1))))))
9)))