spectral-norm Racket #2 program
source code
#lang racket/base
;;; The Computer Language Benchmarks Game
;;; https://salsa.debian.org/benchmarksgame-team/benchmarksgame/
;;; Translated directly from the C# version by Isaac Gouy
;;; contributed by Matthew Flatt
;; This version uses unsafe operations
(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))
notes, command-line, and program output
NOTES:
64-bit Ubuntu quad core
Racket v8.14 [cs].
Fri, 06 Sep 2024 19:13:36 GMT
MAKE:
make: *** No rule to make target 'spectralnorm.racket-2.racket_run'. Stop.
0.07s to complete and log all make actions
COMMAND LINE:
/opt/src/./racket-8.14/bin/racket spectralnorm.racket-2.racket 5500
PROGRAM OUTPUT:
1.274224153