extra/math/blas/vectors/vectors.factor

   1 USING: accessors alien alien.c-types alien.complex alien.data
   2 ascii byte-arrays combinators.short-circuit functors kernel math
   3 math.blas.ffi math.order parser prettyprint.custom sequences
   4 specialized-arrays ;
   5 FROM: alien.c-types => float ;
   6 SPECIALIZED-ARRAY: float
   7 SPECIALIZED-ARRAY: double
   8 SPECIALIZED-ARRAY: complex-float
   9 SPECIALIZED-ARRAY: complex-double
  10 IN: math.blas.vectors
  11
  12 TUPLE: blas-vector-base underlying length inc ;
  13
  14 INSTANCE: blas-vector-base virtual-sequence
  15
  16 GENERIC: element-type ( v -- type )
  17
  18 GENERIC: n*V+V! ( alpha x y -- y=alpha*x+y )
  19 GENERIC: n*V!   ( alpha x -- x=alpha*x )
  20 GENERIC: V. ( x y -- x.y )
  21 GENERIC: V.conj ( x y -- xconj.y )
  22 GENERIC: Vnorm ( x -- norm )
  23 GENERIC: Vasum ( x -- sum )
  24 GENERIC: Vswap ( x y -- x=y y=x )
  25 GENERIC: Viamax ( x -- max-i )
  26
  27 <PRIVATE
  28
  29 GENERIC: (blas-vector-like) ( data length inc exemplar -- vector )
  30
  31 GENERIC: (blas-direct-array) ( blas-vector -- direct-array )
  32
  33 : shorter-length ( v1 v2 -- length )
  34     [ length>> ] bi@ min ; inline
  35 : data-and-inc ( v -- data inc )
  36     [ ] [ inc>> ] bi ; inline
  37 : datas-and-incs ( v1 v2 -- v1-data v1-inc v2-data v2-inc )
  38     [ data-and-inc ] bi@ ; inline
  39
  40 :: (prepare-copy)
  41     ( v element-size -- length v-data v-inc v-dest-data v-dest-inc
  42                         copy-data copy-length copy-inc )
  43     v [ length>> ] [ data-and-inc ] bi
  44     v length>> element-size * <byte-array>
  45     1
  46     over v length>> 1 ;
  47
  48 : (prepare-swap)
  49     ( v1 v2 -- length v1-data v1-inc v2-data v2-inc
  50                v1 v2 )
  51     [ shorter-length ] [ datas-and-incs ] [ ] 2tri ;
  52
  53 :: (prepare-axpy)
  54     ( n v1 v2 -- length n v1-data v1-inc v2-data v2-inc
  55                  v2 )
  56     v1 v2 shorter-length
  57     n
  58     v1 v2 datas-and-incs
  59     v2 ;
  60
  61 :: (prepare-scal)
  62     ( n v -- length n v-data v-inc
  63              v )
  64     v length>>
  65     n
  66     v data-and-inc
  67     v ;
  68
  69 : (prepare-dot) ( v1 v2 -- length v1-data v1-inc v2-data v2-inc )
  70     [ shorter-length ] [ datas-and-incs ] 2bi ;
  71
  72 : (prepare-nrm2) ( v -- length data inc )
  73     [ length>> ] [ data-and-inc ] bi ;
  74
  75 PRIVATE>
  76
  77 : n*V+V ( alpha x y -- alpha*x+y ) clone n*V+V! ; inline
  78 : n*V ( alpha x -- alpha*x ) clone n*V! ; inline
  79
  80 :: V+ ( x y -- x+y )
  81     1.0 x y n*V+V ; inline
  82 :: V- ( x y -- x-y )
  83     -1.0 y x n*V+V ; inline
  84
  85 : Vneg ( x -- -x )
  86     -1.0 swap n*V ; inline
  87
  88 : V*n ( x alpha -- x*alpha )
  89     swap n*V ; inline
  90 : V/n ( x alpha -- x/alpha )
  91     recip swap n*V ; inline
  92
  93 : Vamax ( x -- max )
  94     [ Viamax ] keep nth ; inline
  95
  96 :: Vsub ( v start length -- sub )
  97     v inc>> start * v element-type heap-size *
  98     v underlying>> <displaced-alien>
  99     length v inc>> v (blas-vector-like) ;
 100
 101 : <zero-vector> ( exemplar -- zero )
 102     [ element-type heap-size <byte-array> ]
 103     [ length>> 0 ]
 104     [ (blas-vector-like) ] tri ;
 105
 106 : <empty-vector> ( length exemplar -- vector )
 107     [ element-type heap-size * <byte-array> ]
 108     [ 1 swap ] 2bi
 109     (blas-vector-like) ;
 110
 111 M: blas-vector-base equal?
 112     {
 113         [ [ length ] same? ]
 114         [ [ = ] 2all? ]
 115     } 2&& ;
 116
 117 M: blas-vector-base length
 118     length>> ;
 119 M: blas-vector-base virtual-exemplar
 120     (blas-direct-array) ;
 121 M: blas-vector-base virtual@
 122     [ inc>> * ] [ nip (blas-direct-array) ] 2bi ;
 123
 124 : float>arg ( f -- f ) ; inline
 125 : double>arg ( f -- f ) ; inline
 126 : arg>float ( f -- f ) ; inline
 127 : arg>double ( f -- f ) ; inline
 128
 129 <<
 130
 131 <FUNCTOR: (define-blas-vector) ( TYPE T -- )
 132
 133 <DIRECT-ARRAY> IS <direct-${TYPE}-array>
 134 XCOPY          IS ${T}COPY
 135 XSWAP          IS ${T}SWAP
 136 IXAMAX         IS I${T}AMAX
 137
 138 VECTOR         DEFINES-CLASS ${TYPE}-blas-vector
 139 <VECTOR>       DEFINES <${TYPE}-blas-vector>
 140 >VECTOR        DEFINES >${TYPE}-blas-vector
 141
 142 t              [ T >lower ]
 143
 144 XVECTOR{       DEFINES ${t}vector{
 145
 146 XAXPY          IS ${T}AXPY
 147 XSCAL          IS ${T}SCAL
 148
 149 WHERE
 150
 151 TUPLE: VECTOR < blas-vector-base ;
 152 : <VECTOR> ( underlying length inc -- vector ) VECTOR boa ; inline
 153
 154 : >VECTOR ( seq -- v )
 155     [ TYPE >c-array underlying>> ] [ length ] bi 1 <VECTOR> ;
 156
 157 M: VECTOR clone
 158     TYPE heap-size (prepare-copy)
 159     [ XCOPY ] 3dip <VECTOR> ;
 160
 161 M: VECTOR element-type
 162     drop TYPE ;
 163 M: VECTOR Vswap
 164     (prepare-swap) [ XSWAP ] 2dip ;
 165 M: VECTOR Viamax
 166     (prepare-nrm2) IXAMAX 1 - ;
 167
 168 M: VECTOR (blas-vector-like)
 169     drop <VECTOR> ;
 170
 171 M: VECTOR (blas-direct-array)
 172     [ underlying>> ]
 173     [ [ length>> ] [ inc>> ] bi * ] bi
 174     <DIRECT-ARRAY> ;
 175
 176 M: VECTOR n*V+V!
 177     (prepare-axpy) [ XAXPY ] dip ;
 178 M: VECTOR n*V!
 179     (prepare-scal) [ XSCAL ] dip ;
 180
 181 SYNTAX: XVECTOR{ \ } [ >VECTOR ] parse-literal ;
 182
 183 M: VECTOR pprint-delims
 184     drop \ XVECTOR{ \ } ;
 185
 186 ;FUNCTOR>
 187
 188
 189 <FUNCTOR: (define-real-blas-vector) ( TYPE T -- )
 190
 191 VECTOR         IS ${TYPE}-blas-vector
 192 XDOT           IS ${T}DOT
 193 XNRM2          IS ${T}NRM2
 194 XASUM          IS ${T}ASUM
 195
 196 WHERE
 197
 198 M: VECTOR V.
 199     (prepare-dot) XDOT ;
 200 M: VECTOR V.conj
 201     (prepare-dot) XDOT ;
 202 M: VECTOR Vnorm
 203     (prepare-nrm2) XNRM2 ;
 204 M: VECTOR Vasum
 205     (prepare-nrm2) XASUM ;
 206
 207 ;FUNCTOR>
 208
 209
 210 <FUNCTOR: (define-complex-blas-vector) ( TYPE C S -- )
 211
 212 VECTOR         IS ${TYPE}-blas-vector
 213 XDOTU          IS ${C}DOTU
 214 XDOTC          IS ${C}DOTC
 215 XXNRM2         IS ${S}${C}NRM2
 216 XXASUM         IS ${S}${C}ASUM
 217
 218 WHERE
 219
 220 M: VECTOR V.
 221     (prepare-dot) XDOTU ;
 222 M: VECTOR V.conj
 223     (prepare-dot) XDOTC ;
 224 M: VECTOR Vnorm
 225     (prepare-nrm2) XXNRM2 ;
 226 M: VECTOR Vasum
 227     (prepare-nrm2) XXASUM ;
 228
 229 ;FUNCTOR>
 230
 231
 232 : define-real-blas-vector ( TYPE T -- )
 233     [ (define-blas-vector) ]
 234     [ (define-real-blas-vector) ] 2bi ;
 235 : define-complex-blas-vector ( TYPE C S -- )
 236     [ drop (define-blas-vector) ]
 237     [ (define-complex-blas-vector) ] 3bi ;
 238
 239 float  "S" define-real-blas-vector
 240 double "D" define-real-blas-vector
 241 complex-float  "C" "S" define-complex-blas-vector
 242 complex-double "Z" "D" define-complex-blas-vector
 243
 244 >>
 245
 246 M: blas-vector-base >pprint-sequence ;
 247 M: blas-vector-base pprint* pprint-object ;