On 05 Jun 2005 17:20:40 +0200, Matthias Heiler <[EMAIL PROTECTED]> wrote:
> Hi,
>
> I played with the spectral norm algorithm from the computer language
> shootout http://tinylink.com/?fFObJZJlyQ
>
Here a few points:
- loop is (relatively) slow - it expands into at least one call/cc for example
(this isn't as bad as it sounds, continuations are fast with chicken, but if
you want to squeeze out perfornance, a simple 'do' loop is faster).
- I would add `(declare (block) (disable-interrupts))'.
- The real whopper: f64vector. Extracting doubles from a f64vector is
really costly: we have to cons flonums all the time. Putting the flonums
into a vector instead reduces the need to create flonums (there are already
consed and pointers to them are stored in the vector). This doesn't
work with the FFI-wrapped code, of course, which wants f64vectors.
Attached an altenative version, here the Scheme code needs about
7 seconds (compared to 13 seconds for the original version).
cheers,
felix
;;;; Implement spectralnorm algorithm from Computer Language Shootout
;;;; usage: csc -O3 spectalnorm.scm
;;;; ./spectralnorm 500
;(require-extension srfi-4)
(declare (block) (disable-interrupts))
(define (A i j)
(/ 1.0 (+ (* (+ i j)
(/ (+ i j 1.0) 2.0))
i 1.0)))
(define (A-times-u N u Au)
(do ((i 0 (fx+ i 1)))
((fx>= i N))
(vector-set! Au i
(let loop ((j 0) (sum 0))
(if (fx>= j N)
sum
(loop (fx+ j 1)
(+ sum (* (A i j) (vector-ref u j)))))))))
(define (At-times-u N u Au)
(do ((i 0 (fx+ i 1)))
((fx>= i N))
(vector-set! Au i
(let loop ((j 0) (sum 0))
(if (fx>= j N)
sum
(loop (fx+ j 1)
(+ sum (* (A j i) (vector-ref u j)))))))))
(define (AtA-times-u N u AtAu)
(let ((v (make-vector 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)
(do ((i 0 (fx+ i 1))
(sum 0 (+ sum (* (vector-ref u i)
(vector-ref v i)))))
((fx>= i N) sum) ))
(define (main-1)
(let* ((N (with-input-from-string (car (command-line-arguments)) read))
(u (make-vector N 1.0))
(v (make-vector N)))
(do ((i 0 (fx+ i 1)))
((fx>= i 10))
(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-vector N 1.0))
(v (make-vector N)))
(do ((i 0 (fx+ i 1)))
((fx>= i 10))
(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))
_______________________________________________
Chicken-users mailing list
Chicken-users@nongnu.org
http://lists.nongnu.org/mailman/listinfo/chicken-users
