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