diff --git a/collects/tests/mzscheme/benchmarks/shootout/run.ss b/collects/tests/mzscheme/benchmarks/shootout/run.ss index 5c3d2cc2a5..2277a77123 100644 --- a/collects/tests/mzscheme/benchmarks/shootout/run.ss +++ b/collects/tests/mzscheme/benchmarks/shootout/run.ss @@ -34,6 +34,7 @@ ("reversefile.ss") ("sieve.ss" "1200") ("spellcheck.ss") + ("spectralnorm.ss" "5500") ("strcat.ss" "40000") ("sumcol.ss" #f ,(lambda () (mk-sumcol-input))) ("wc.ss") diff --git a/collects/tests/mzscheme/benchmarks/shootout/spectralnorm.ss b/collects/tests/mzscheme/benchmarks/shootout/spectralnorm.ss new file mode 100644 index 0000000000..1eea79e1c3 --- /dev/null +++ b/collects/tests/mzscheme/benchmarks/shootout/spectralnorm.ss @@ -0,0 +1,62 @@ + +(module spectralnorm mzscheme + (require (lib "string.ss")) + + (define (Approximate n) + (let ([u (make-vector n 1.0)] + [v (make-vector n 0.0)]) + ;; 20 steps of the power method + (let loop ([i 0]) + (unless (= i 10) + (MultiplyAtAv n u v) + (MultiplyAtAv n v u) + (loop (add1 i)))) + + ;; B=AtA A multiplied by A transposed + ;; v.Bv /(v.v) eigenvalue of v + (let loop ([i 0][vBv 0][vv 0]) + (if (= i n) + (sqrt (/ vBv vv)) + (let ([vi (vector-ref v i)]) + (loop (add1 i) + (+ vBv (* (vector-ref u i) vi)) + (+ vv (* vi vi)))))))) + + ;; return element i,j of infinite matrix A + (define (A i j) + (/ 1.0 (+ (* (+ i j) (/ (+ i j 1) 2)) i 1))) + + ;; multiply vector v by matrix A + (define (MultiplyAv n v Av) + (let loop ([i 0]) + (unless (= i n) + (let jloop ([j 0][r 0]) + (if (= j n) + (vector-set! Av i r) + (jloop (add1 j) + (+ r (* (A i j) (vector-ref v j)))))) + (loop (add1 i))))) + + ;; multiply vector v by matrix A transposed + (define (MultiplyAtv n v Atv) + (let loop ([i 0]) + (unless (= i n) + (let jloop ([j 0][r 0]) + (if (= j n) + (vector-set! Atv i r) + (jloop (add1 j) + (+ r (* (A j i) (vector-ref v j)))))) + (loop (add1 i))))) + + ;; multiply vector v by matrix A and then by matrix A transposed + (define (MultiplyAtAv n v AtAv) + (let ([u (make-vector n 0.0)]) + (MultiplyAv n v u) + (MultiplyAtv n u AtAv))) + + (printf "~a\n" + (real->decimal-string + (Approximate (string->number (vector-ref + (current-command-line-arguments) + 0))) + 9)))