1.4.1 A dataset on the GPU🔗ℹ

The datasets package loads Fisher’s iris flowers as a Polars dataframe:

> (require datasets            (only-in polars in-dataframe-columns series->list))

> (define iris (load-iris))

> iris

shape: (150, 5) ┌──────────────┬─────────────┬──────────────┬─────────────┬────────────────┐ │ sepal-length ┆ sepal-width ┆ petal-length ┆ petal-width ┆ species        │ │ ---          ┆ ---         ┆ ---          ┆ ---         ┆ ---            │ │ f64          ┆ f64         ┆ f64          ┆ f64         ┆ str            │ ╞══════════════╪═════════════╪══════════════╪═════════════╪════════════════╡ │ 5.1          ┆ 3.5         ┆ 1.4          ┆ 0.2         ┆ Iris-setosa    │ │ 4.9          ┆ 3.0         ┆ 1.4          ┆ 0.2         ┆ Iris-setosa    │ │ 4.7          ┆ 3.2         ┆ 1.3          ┆ 0.2         ┆ Iris-setosa    │ │ 4.6          ┆ 3.1         ┆ 1.5          ┆ 0.2         ┆ Iris-setosa    │ │ 5.0          ┆ 3.6         ┆ 1.4          ┆ 0.2         ┆ Iris-setosa    │ │ …            ┆ …           ┆ …            ┆ …           ┆ …              │ │ 6.7          ┆ 3.0         ┆ 5.2          ┆ 2.3         ┆ Iris-virginica │ │ 6.3          ┆ 2.5         ┆ 5.0          ┆ 1.9         ┆ Iris-virginica │ │ 6.5          ┆ 3.0         ┆ 5.2          ┆ 2.0         ┆ Iris-virginica │ │ 6.2          ┆ 3.4         ┆ 5.4          ┆ 2.3         ┆ Iris-virginica │ │ 5.9          ┆ 3.0         ┆ 5.1          ┆ 1.8         ┆ Iris-virginica │ └──────────────┴─────────────┴──────────────┴─────────────┴────────────────┘

The four measurement columns, as a list of rows:

> (define measurements '("sepal-length" "sepal-width" "petal-length" "petal-width"))

> (define rows     (apply map            list            (for/list ([column (in-dataframe-columns iris #:columns measurements)])              (series->list column))))

> (length rows)

150

> (take rows 3)

'((5.1 3.5 1.4 0.2) (4.9 3.0 1.4 0.2) (4.7 3.2 1.3 0.2))

On the GPU, as a 'float32 matrix:

> (define X (list*->device-matrix rows #:dtype 'float32))

> X

#<device-matrix float32[150×4] row-major cuda:0  [[5.1 3.5 1.4 0.2]   [4.9 3.0 1.4 0.2]   [4.7 3.2 1.3 0.2]   ...   [6.5 3.0 5.2 2.0]   [6.2 3.4 5.4 2.3]   [5.9 3.0 5.1 1.8]]>

Printing copies the values back, so it waits for the GPU. The matrix knows its shape, element type, layout and strides:

> (shape X)

'(150 4)

> (dtype X)

'float32

> (layout X)

'row-major

> (strides X)

'(4 1)

> (numel X)

600

Each row is contiguous: the next element along a row is 1 away, the next row 4 away. Strides count elements, not bytes.

1.4.2 Changing the layout🔗ℹ

contiguous copies a matrix into the other layout on the GPU:

> (define F (contiguous X #:layout 'col-major))

> (layout F)

'col-major

> (strides F)

'(1 150)

> (contiguous? F #:layout 'col-major)

#t

> (equal? (device-matrix->list* F) (device-matrix->list* X))

#t

Same values, different memory order. contiguous copies only when it has to:

> (eq? (contiguous F #:layout 'col-major) F)

#t

> (eq? (contiguous X) X)

#t

Data from Racket can be packed column-major on the host instead:

> (define F* (list*->device-matrix rows #:dtype 'float32 #:layout 'col-major))

> (strides F*)

'(1 150)

> (equal? (device-matrix->list* F*) (device-matrix->list* F))

#t

1.4.3 Other forms of Racket data🔗ℹ

raft/compat converts the other forms Racket keeps numbers in, by the same dtype and layout rules:

> (require raft/compat            ffi/vector            math/array            math/matrix            racket/file)

list*->device-array reads the rank from the nesting:

> (shape (list*->device-array rows))

'(150 4)

> (shape (list*->device-array (first rows)))

'(4)

A vector: the species, numbered in order of appearance, as 'int32 labels:

> (define species     (for/first ([column (in-dataframe-columns iris #:columns '("species"))])       (series->list column)))

> (define names (remove-duplicates species))

> (define labels     (vector->device-vector (for/vector ([s (in-list species)])                              (index-of names s))                            #:dtype 'int32))

> labels

#<device-vector int32[150] cuda:0 [0 0 0 ... 2 2 2]>

A math/matrix matrix: the measurements’ covariance, computed on the host:

> (define M (list*->matrix rows))

> (define n (matrix-num-rows M))

> (define means (matrix-scale (matrix* (make-matrix 1 n 1) M) (/ 1.0 n)))

> (define centered (matrix- M (matrix* (make-matrix n 1 1) means)))

> (define covariance     (matrix-scale (matrix* (matrix-transpose centered) centered) (/ 1.0 (sub1 n))))

> (define C (matrix->device-matrix covariance #:dtype 'float32))

> C

#<device-matrix float32[4×4] row-major cuda:0  [[   0.6856935 -0.042434003    1.2743155    0.5162707]   [-0.042434003   0.18997942  -0.32965636  -0.12163937]   [   1.2743155  -0.32965636     3.116278    1.2956094]   [   0.5162707  -0.12163937    1.2956094    0.5810063]]>

An f64vector, as foreign code fills one; a loop stands in for it:

> (define weights (make-f64vector n 1.0))

> (f64vector-set! weights 0 2.0)

> (define w (f64vector->device-vector weights #:dtype 'float32))

> w

#<device-vector float32[150] cuda:0 [2.0 1.0 1.0 ... 1.0 1.0 1.0]>

Narrowing to 'float32 raises on a value too large rather than giving an infinity.

1.4.4 Bringing values back🔗ℹ

The petal widths, from an flvector:

> (define widths     (for/flvector ([row (in-list rows)])       (last row)))

> (define y (flvector->device-vector widths #:dtype 'float32))

> y

#<device-vector float32[150] cuda:0 [0.2 0.2 0.2 ... 2.0 2.3 1.8]>

> (take (device-vector->list y) 4)

'(0.20000000298023224   0.20000000298023224   0.20000000298023224   0.20000000298023224)

Float elements come back as flonums, integer elements as exact integers. A 'float32 cannot hold 0.2 exactly, so what comes back is the nearest 'float32, widened; keep exact values in 'float64, the type inferred for rows:

> (first (device-matrix->list* X))

'(5.099999904632568 3.5 1.399999976158142 0.20000000298023224)

> (first (device-matrix->list* (list*->device-matrix rows)))

'(5.1 3.5 1.4 0.2)

Each conversion names the Racket form it makes. The covariance comes back as a math/array flonum array, off only by the 'float32 rounding; column means through math/array; species counts from a vector; the weights as raw bytes, through a file:

> (< (array-all-max (array-abs (array- (device-matrix->matrix C) covariance))) 1e-6)

#t

> (array-scale (array-axis-sum (device-array->array X) 0) (/ 1.0 n))

(array  #[5.843333326975505 3.057333335876465 3.757999988396963 1.1993333247800668])

> (define counts (make-vector (length names) 0))

> (for ([label (in-vector (device-vector->vector labels))])     (vector-set! counts label (add1 (vector-ref counts label))))

> counts

'#(50 50 50)

> (define saved (make-temporary-file))

> (display-to-file (device-vector->bytes w) saved #:exists 'truncate)

> (file-size saved)

600

> (take (device-vector->list (bytes->device-vector (file->bytes saved) #:dtype 'float32)) 3)

'(2.0 1.0 1.0)

1.4.5 Letting go🔗ℹ

The program never frees an array: when one becomes unreachable, its finalizer returns the memory. Each buffer counts its size as phantom bytes, so GPU memory held brings collections sooner:

> (collect-garbage)

> (define before (current-memory-use))

> (define scratch (device-matrix 1024 1024))

> (quotient (- (current-memory-use) before) (* 1024 1024))

4

> (set! scratch #f)

> (collect-garbage)

> (< (- (current-memory-use) before) (* 1024 1024))

#t

The memory comes back at a collection after the last use.