extra/benchmark/spectral-norm/spectral-norm.factor

   1 ! Copyright (C) 2008, 2010 Slava Pestov.
   2 ! See https://factorcode.org/license.txt for BSD license.
   3 !
   4 ! Factor port of
   5 ! https://shootout.alioth.debian.org/gp4/benchmark.php?test=spectralnorm&lang=all
   6 USING: alien.c-types io kernel math math.functions math.parser
   7 math.vectors sequences sequences.private specialized-arrays
   8 typed locals ;
   9 SPECIALIZED-ARRAY: double
  10 IN: benchmark.spectral-norm
  11
  12 :: inner-loop ( u n quot -- seq )
  13     n <iota> [| i |
  14         n <iota> 0.0 [| j |
  15             u i j quot call +
  16         ] reduce
  17     ] double-array{ } map-as ; inline
  18
  19 : eval-A ( i j -- n )
  20     [ >float ] bi@
  21     [ drop ] [ + [ ] [ 1 + ] bi * 0.5 * ] 2bi
  22     + 1 + recip ; inline
  23
  24 : (eval-A-times-u) ( u i j -- x )
  25     [ swap nth-unsafe ] [ eval-A ] bi-curry bi* * ; inline
  26
  27 : eval-A-times-u ( n u -- seq )
  28     [ (eval-A-times-u) ] inner-loop ; inline
  29
  30 : (eval-At-times-u) ( u i j -- x )
  31     [ swap nth-unsafe ] [ swap eval-A ] bi-curry bi* * ; inline
  32
  33 : eval-At-times-u ( u n -- seq )
  34     [ (eval-At-times-u) ] inner-loop ; inline
  35
  36 : eval-AtA-times-u ( u n -- seq )
  37     [ eval-A-times-u ] [ eval-At-times-u ] bi ; inline
  38
  39 : ones ( n -- seq ) [ 1.0 ] double-array{ } replicate-as ; inline
  40
  41 :: u/v ( n -- u v )
  42     n ones dup
  43     10 [
  44         drop
  45         n eval-AtA-times-u
  46         [ n eval-AtA-times-u ] keep
  47     ] times ; inline
  48
  49 TYPED: spectral-norm ( n: fixnum -- norm )
  50     u/v [ vdot ] [ norm-sq ] bi /f sqrt ;
  51
  52 : spectral-norm-benchmark ( -- )
  53     2000 spectral-norm number>string print ;
  54
  55 MAIN: spectral-norm-benchmark