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

[Matthias Heiler]

Hi,

I played with the spectral norm algorithm from the computer language
shootout http://tinylink.com/?fFObJZJlyQ 
  
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))
===