basis/math/blas/vectors/vectors.factor

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