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

   1 ! Copyright (C) 2010 Marc Fauconneau.
   2 ! See http://factorcode.org/license.txt for BSD license.
   3 USING: alien.c-types alien.data specialized-arrays kernel math
   4 math.functions math.vectors sequences sequences.private
   5 prettyprint words typed locals math.vectors.simd
   6 math.vectors.simd.cords ;
   7 SPECIALIZED-ARRAYS: double double-4 ;
   8 IN: benchmark.spectral-norm-simd
   9
  10 :: inner-loop ( u n quot -- seq )
  11     n 4 /i iota [| i |
  12         n iota [| j | u i j quot call ] [ v+ ] map-reduce
  13     ] double-4-array{ } map-as ; inline
  14
  15 : eval-A ( i j -- n )
  16     [ >float ] bi@
  17     [ drop ] [ + [ ] [ 1 + ] bi * 0.5 * ] 2bi
  18     + 1 + ; inline
  19
  20 : vrecip ( u -- v ) double-4{ 1.0 1.0 1.0 1.0 } swap v/ ; inline
  21
  22 :: eval4-A ( i j -- n )
  23     i 4 * 0 + j eval-A
  24     i 4 * 1 + j eval-A
  25     i 4 * 2 + j eval-A
  26     i 4 * 3 + j eval-A
  27     double-4-boa vrecip ; inline
  28
  29 : (eval-A-times-u) ( u i j -- x )
  30     [ swap nth-unsafe ] [ eval4-A ] bi-curry bi* n*v ; inline
  31
  32 : eval-A-times-u ( n u -- seq )
  33     [ (eval-A-times-u) ] inner-loop ; inline
  34
  35 :: eval4-A' ( i j -- n )
  36     j i 4 * 0 + eval-A
  37     j i 4 * 1 + eval-A
  38     j i 4 * 2 + eval-A
  39     j i 4 * 3 + eval-A
  40     double-4-boa vrecip ; inline
  41
  42 : (eval-At-times-u) ( u i j -- x )
  43     [ swap nth-unsafe ] [ eval4-A' ] bi-curry bi* n*v ; inline
  44
  45 : eval-At-times-u ( u n -- seq )
  46     [ double cast-array ] dip [ (eval-At-times-u) ] inner-loop ; inline
  47
  48 : eval-AtA-times-u ( u n -- seq )
  49     [ double cast-array ] dip [ eval-A-times-u ] [ eval-At-times-u ] bi ; inline
  50
  51 : ones ( n -- seq )
  52     4 /i [ double-4{ 1.0 1.0 1.0 1.0 } ] double-4-array{ } replicate-as ; inline
  53
  54 :: u/v ( n -- u v )
  55     n ones dup
  56     10 [
  57         drop
  58         n eval-AtA-times-u
  59         [ n eval-AtA-times-u ] keep
  60     ] times ; inline
  61
  62 TYPED: spectral-norm ( n: fixnum -- norm )
  63     u/v [ double cast-array ] bi@ [ v. ] [ norm-sq ] bi /f sqrt ;
  64
  65 : spectral-norm-main ( -- )
  66     2000 spectral-norm . ;
  67
  68 MAIN: spectral-norm-main