basis/math/blas/vectors/vectors.factor

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