chicken-users
[Top][All Lists]

## [Chicken-users] spectralnorm algorithm from comp. language shootout

 From: Matthias Heiler Subject: [Chicken-users] spectralnorm algorithm from comp. language shootout Date: 05 Jun 2005 17:20:40 +0200 User-agent: Gnus/5.09 (Gnus v5.9.0) Emacs/21.3

```Hi,

I played with the spectral norm algorithm from the computer language

I took the C version and integrated it using the simple FFI.  Also, I
converted it into real scheme (but keeping close to the C style).  The
real scheme version is quite a bit slower than the FFI-Version (about
15x slower).  First I had the temporary vector in AtA-times-u under
suspicion but removing it doesn't improve speed a lot.  Maybe loop is
the culprit?  Or is it the basic arithmetic operators?

Yours,

Matthias

===
;;;; Implement spectralnorm algorithm from Computer Language Shootout

;;;; usage: csc -O3 spectalnorm.scm
;;;;        ./spectralnorm 500

(require-extension srfi-4 loop)

(define (A i j)
(/ 1.0 (+ (* (+ i j)
(/ (+ i j 1.0) 2.0))
i 1.0)))

(define (A-times-u N u Au)
(loop for i from 0 to (- N 1) do
(f64vector-set! Au i
(loop for j from 0 to (- N 1)
sum (* (A i j) (f64vector-ref u j))))))

(define (At-times-u N u Au)
(loop for i from 0 to (- N 1) do
(f64vector-set! Au i
(loop for j from 0 to (- N 1)
sum (* (A j i) (f64vector-ref u j))))))

(define (AtA-times-u N u AtAu)
(let ((v (make-f64vector N)))
(A-times-u N u v)
(At-times-u N v AtAu)))

#>!
double eval_A(int i, int j)
{
return (1.0/((i+j)*(i+j+1)/2+i+1));
}

void eval_A_times_u(int N, const double u[], double Au[])
{
int i,j;
for(i=0;i<N;i++)
{
Au[i]=0;
for(j=0;j<N;j++) Au[i]+=eval_A(i,j)*u[j];
}
}

void eval_At_times_u(int N, const double u[], double Au[])
{
int i,j;
for(i=0;i<N;i++)
{
Au[i]=0;
for(j=0;j<N;j++) Au[i]+=eval_A(j,i)*u[j];
}
}

void eval_AtA_times_u(int N, const double u[], double AtAu[])
{ double v[N]; eval_A_times_u(N,u,v); eval_At_times_u(N,v,AtAu); }

double eval_dot(int N, const double a[], const double b[])
{
int i;
double res = 0.0;
for (i = 0; i < N; ++i) {
res += a[i] * b[i];
}
return (res);
}
<#

(define (dot N u v)
(loop for i from 0 to (- N 1)
sum (* (f64vector-ref u i)
(f64vector-ref v i))))

(define (main-1)
(let* ((N (with-input-from-string (car (command-line-arguments)) read))
(u (make-f64vector N 1.0))
(v (make-f64vector N)))
(loop for i from 0 to 9 do
(eval_AtA_times_u N u v)
(eval_AtA_times_u N v u))
(print (sqrt (/ (dot N v u) (dot N v v))))))

(define (main-2)
(let* ((N (with-input-from-string (car (command-line-arguments)) read))
(u (make-f64vector N 1.0))
(v (make-f64vector N)))
(loop for i from 0 to 9 do
(AtA-times-u N u v)
(AtA-times-u N v u))
(print (sqrt (/ (dot N v u) (dot N v v))))))

(time (main-1))
(time (main-2))
===

```

reply via email to