1 USING: accessors alien alien.c-types arrays byte-arrays combinators
2 combinators.short-circuit fry kernel math math.blas.cblas
3 math.complex math.functions math.order sequences.complex
4 sequences.complex-components sequences sequences.private
6 specialized-arrays.float specialized-arrays.double
7 specialized-arrays.direct.float specialized-arrays.direct.double ;
10 TUPLE: blas-vector-base underlying length inc ;
12 INSTANCE: blas-vector-base virtual-sequence
14 GENERIC: element-type ( v -- type )
16 GENERIC: n*V+V! ( alpha x y -- y=alpha*x+y )
17 GENERIC: n*V! ( alpha x -- x=alpha*x )
18 GENERIC: V. ( x y -- x.y )
19 GENERIC: V.conj ( x y -- xconj.y )
20 GENERIC: Vnorm ( x -- norm )
21 GENERIC: Vasum ( x -- sum )
22 GENERIC: Vswap ( x y -- x=y y=x )
23 GENERIC: Viamax ( x -- max-i )
27 GENERIC: (blas-vector-like) ( data length inc exemplar -- vector )
29 GENERIC: (blas-direct-array) ( blas-vector -- direct-array )
31 : shorter-length ( v1 v2 -- length )
32 [ length>> ] bi@ min ; inline
33 : data-and-inc ( v -- data inc )
34 [ underlying>> ] [ inc>> ] bi ; inline
35 : datas-and-incs ( v1 v2 -- v1-data v1-inc v2-data v2-inc )
36 [ data-and-inc ] bi@ ; inline
39 ( v element-size -- length v-data v-inc v-dest-data v-dest-inc
40 copy-data copy-length copy-inc )
41 v [ length>> ] [ data-and-inc ] bi
42 v length>> element-size * <byte-array>
47 ( v1 v2 -- length v1-data v1-inc v2-data v2-inc
49 [ shorter-length ] [ datas-and-incs ] [ ] 2tri ;
52 ( n v1 v2 -- length n v1-data v1-inc v2-data v2-inc
60 ( n v -- length n v-data v-inc
67 : (prepare-dot) ( v1 v2 -- length v1-data v1-inc v2-data v2-inc )
68 [ shorter-length ] [ datas-and-incs ] 2bi ;
70 : (prepare-nrm2) ( v -- length data inc )
71 [ length>> ] [ data-and-inc ] bi ;
75 : n*V+V ( alpha x y -- alpha*x+y ) clone n*V+V! ; inline
76 : n*V ( alpha x -- alpha*x ) clone n*V! ; inline
79 1.0 -rot n*V+V ; inline
81 -1.0 spin n*V+V ; inline
84 -1.0 swap n*V ; inline
86 : V*n ( x alpha -- x*alpha )
88 : V/n ( x alpha -- x/alpha )
89 recip swap n*V ; inline
92 [ Viamax ] keep nth ; inline
94 :: Vsub ( v start length -- sub )
95 v inc>> start * v element-type heap-size *
96 v underlying>> <displaced-alien>
97 length v inc>> v (blas-vector-like) ;
99 : <zero-vector> ( exemplar -- zero )
100 [ element-type <c-object> ]
102 [ (blas-vector-like) ] tri ;
104 : <empty-vector> ( length exemplar -- vector )
105 [ element-type <c-array> ]
109 M: blas-vector-base equal?
115 M: blas-vector-base length
117 M: blas-vector-base virtual-seq
118 (blas-direct-array) ;
119 M: blas-vector-base virtual@
120 [ inc>> * ] [ nip (blas-direct-array) ] 2bi ;
122 : float>arg ( f -- f ) ; inline
123 : double>arg ( f -- f ) ; inline
124 : arg>float ( f -- f ) ; inline
125 : arg>double ( f -- f ) ; inline
129 FUNCTOR: (define-blas-vector) ( TYPE T -- )
131 <DIRECT-ARRAY> IS <direct-${TYPE}-array>
132 >ARRAY IS >${TYPE}-array
133 XCOPY IS cblas_${T}copy
134 XSWAP IS cblas_${T}swap
135 IXAMAX IS cblas_i${T}amax
137 VECTOR DEFINES ${TYPE}-blas-vector
138 <VECTOR> DEFINES <${TYPE}-blas-vector>
139 >VECTOR DEFINES >${TYPE}-blas-vector
143 TUPLE: VECTOR < blas-vector-base ;
144 : <VECTOR> ( underlying length inc -- vector ) VECTOR boa ; inline
146 : >VECTOR ( seq -- v )
147 [ >ARRAY execute underlying>> ] [ length ] bi 1 <VECTOR> execute ;
150 TYPE heap-size (prepare-copy)
151 [ XCOPY execute ] 3dip <VECTOR> execute ;
153 M: VECTOR element-type
156 (prepare-swap) [ XSWAP execute ] 2dip ;
158 (prepare-nrm2) IXAMAX execute ;
160 M: VECTOR (blas-vector-like)
161 drop <VECTOR> execute ;
163 M: VECTOR (blas-direct-array)
165 [ [ length>> ] [ inc>> ] bi * ] bi
166 <DIRECT-ARRAY> execute ;
171 FUNCTOR: (define-real-blas-vector) ( TYPE T -- )
173 VECTOR IS ${TYPE}-blas-vector
174 XDOT IS cblas_${T}dot
175 XNRM2 IS cblas_${T}nrm2
176 XASUM IS cblas_${T}asum
177 XAXPY IS cblas_${T}axpy
178 XSCAL IS cblas_${T}scal
183 (prepare-dot) XDOT execute ;
185 (prepare-dot) XDOT execute ;
187 (prepare-nrm2) XNRM2 execute ;
189 (prepare-nrm2) XASUM execute ;
191 (prepare-axpy) [ XAXPY execute ] dip ;
193 (prepare-scal) [ XSCAL execute ] dip ;
198 FUNCTOR: (define-complex-helpers) ( TYPE -- )
200 <DIRECT-COMPLEX-ARRAY> DEFINES <direct-${TYPE}-complex-array>
201 >COMPLEX-ARRAY DEFINES >${TYPE}-complex-array
202 ARG>COMPLEX DEFINES arg>${TYPE}-complex
203 COMPLEX>ARG DEFINES ${TYPE}-complex>arg
204 <DIRECT-ARRAY> IS <direct-${TYPE}-array>
205 >ARRAY IS >${TYPE}-array
209 : <DIRECT-COMPLEX-ARRAY> ( alien len -- sequence )
210 1 shift <DIRECT-ARRAY> execute <complex-sequence> ;
211 : >COMPLEX-ARRAY ( sequence -- sequence )
212 <complex-components> >ARRAY execute ;
213 : COMPLEX>ARG ( complex -- alien )
214 >rect 2array >ARRAY execute underlying>> ;
215 : ARG>COMPLEX ( alien -- complex )
216 2 <DIRECT-ARRAY> execute first2 rect> ;
221 FUNCTOR: (define-complex-blas-vector) ( TYPE C S -- )
223 VECTOR IS ${TYPE}-blas-vector
224 XDOTU_SUB IS cblas_${C}dotu_sub
225 XDOTC_SUB IS cblas_${C}dotc_sub
226 XXNRM2 IS cblas_${S}${C}nrm2
227 XXASUM IS cblas_${S}${C}asum
228 XAXPY IS cblas_${C}axpy
229 XSCAL IS cblas_${C}scal
230 TYPE>ARG IS ${TYPE}>arg
231 ARG>TYPE IS arg>${TYPE}
236 (prepare-dot) TYPE <c-object>
237 [ XDOTU_SUB execute ] keep
240 (prepare-dot) TYPE <c-object>
241 [ XDOTC_SUB execute ] keep
244 (prepare-nrm2) XXNRM2 execute ;
246 (prepare-nrm2) XXASUM execute ;
248 [ TYPE>ARG execute ] 2dip
249 (prepare-axpy) [ XAXPY execute ] dip ;
251 [ TYPE>ARG execute ] dip
252 (prepare-scal) [ XSCAL execute ] dip ;
257 : define-real-blas-vector ( TYPE T -- )
258 [ (define-blas-vector) ]
259 [ (define-real-blas-vector) ] 2bi ;
260 :: define-complex-blas-vector ( TYPE C S -- )
261 TYPE (define-complex-helpers)
262 TYPE "-complex" append
263 [ C (define-blas-vector) ]
264 [ C S (define-complex-blas-vector) ] bi ;
266 "float" "s" define-real-blas-vector
267 "double" "d" define-real-blas-vector
268 "float" "c" "s" define-complex-blas-vector
269 "double" "z" "d" define-complex-blas-vector