add another missing benchmark

svn: r5499

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")

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)))