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
10 :: inner-loop ( u n quot -- seq )
12 n iota [| j | u i j quot call ] [ v+ ] map-reduce
13 ] double-4-array{ } map-as ; inline
17 [ drop ] [ + [ ] [ 1 + ] bi * 0.5 * ] 2bi
20 : vrecip ( u -- v ) double-4{ 1.0 1.0 1.0 1.0 } swap v/ ; inline
22 :: eval4-A ( i j -- n )
27 double-4-boa vrecip ; inline
29 : (eval-A-times-u) ( u i j -- x )
30 [ swap nth-unsafe ] [ eval4-A ] bi-curry bi* n*v ; inline
32 : eval-A-times-u ( n u -- seq )
33 [ (eval-A-times-u) ] inner-loop ; inline
35 :: eval4-A' ( i j -- n )
40 double-4-boa vrecip ; inline
42 : (eval-At-times-u) ( u i j -- x )
43 [ swap nth-unsafe ] [ eval4-A' ] bi-curry bi* n*v ; inline
45 : eval-At-times-u ( u n -- seq )
46 [ double cast-array ] dip [ (eval-At-times-u) ] inner-loop ; inline
48 : eval-AtA-times-u ( u n -- seq )
49 [ double cast-array ] dip [ eval-A-times-u ] [ eval-At-times-u ] bi ; inline
52 4 /i [ double-4{ 1.0 1.0 1.0 1.0 } ] double-4-array{ } replicate-as ; inline
59 [ n eval-AtA-times-u ] keep
62 TYPED: spectral-norm ( n: fixnum -- norm )
63 u/v [ double cast-array ] bi@ [ v. ] [ norm-sq ] bi /f sqrt ;
65 : spectral-norm-simd-benchmark ( -- )
66 2000 spectral-norm . ;
68 MAIN: spectral-norm-simd-benchmark