-
Notifications
You must be signed in to change notification settings - Fork 14
/
Copy pathREADME
1964 lines (1389 loc) · 71.1 KB
/
README
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
# MAT Manual
###### \[in package MGL-MAT\]
## The MGL-MAT ASDF System
- Version: 0.1.0
- Description: MAT is library for working with multi-dimensional
arrays which supports efficient interfacing to foreign and CUDA
code. BLAS and CUBLAS bindings are available.
- Licence: MIT, see COPYING.
- Author: Gábor Melis <[email protected]>
- Mailto: [[email protected]](mailto:[email protected])
- Homepage: [http://melisgl.github.io/mgl-mat](http://melisgl.github.io/mgl-mat)
- Bug tracker: [https://github.com/melisgl/mgl-mat/issues](https://github.com/melisgl/mgl-mat/issues)
- Source control: [GIT](https://github.com/melisgl/mgl-mat.git)
## Links
Here is the [official
repository](https://github.com/melisgl/mgl-mat) and the [HTML
documentation](http://melisgl.github.io/mgl-mat/mat-manual.html)
for the latest version.
## Introduction
### What's MGL-MAT?
MGL-MAT is library for working with multi-dimensional arrays
which supports efficient interfacing to foreign and CUDA code with
automatic translations between cuda, foreign and lisp storage. BLAS
and CUBLAS bindings are available.
### What kind of matrices are supported?
Currently only row-major single and double float matrices are
supported, but it would be easy to add single and double precision
complex types too. Other numeric types, such as byte and native
integer, can be added too, but they are not supported by CUBLAS.
There are no restrictions on the number of dimensions, and reshaping
is possible. All functions operate on the visible portion of the
matrix (which is subject to displacement and shaping), invisible
elements are not affected.
### Where to Get it?
All dependencies are in quicklisp except for
[CL-CUDA](https://github.com/takagi/cl-cuda) that needs to be
fetched from github. Just clone CL-CUDA and MGL-MAT into
`quicklisp/local-projects/` and you are set. MGL-MAT itself lives
[at github](https://github.com/melisgl/mgl-mat), too.
Prettier-than-markdown HTML documentation cross-linked with other
libraries is
[available](http://melisgl.github.io/mgl-pax-world/mat-manual.html)
as part of [PAX World](http://melisgl.github.io/mgl-pax-world/).
## Tutorial
We are going to see how to create matrices, access their contents.
Creating matrices is just like creating lisp arrays:
```commonlisp
(make-mat '6)
==> #<MAT 6 A #(0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0)>
(make-mat '(2 3) :ctype :float :initial-contents '((1 2 3) (4 5 6)))
==> #<MAT 2x3 AB #2A((1.0 2.0 3.0) (4.0 5.0 6.0))>
(make-mat '(2 3 4) :initial-element 1)
==> #<MAT 2x3x4 A #3A(((1.0d0 1.0d0 1.0d0 1.0d0)
--> (1.0d0 1.0d0 1.0d0 1.0d0)
--> (1.0d0 1.0d0 1.0d0 1.0d0))
--> ((1.0d0 1.0d0 1.0d0 1.0d0)
--> (1.0d0 1.0d0 1.0d0 1.0d0)
--> (1.0d0 1.0d0 1.0d0 1.0d0)))>
```
The most prominent difference from lisp arrays is that `MAT`s are
always numeric and their element type (called CTYPE here) defaults
to :DOUBLE.
Individual elements can be accessed or set with MREF:
```commonlisp
(let ((m (make-mat '(2 3))))
(setf (mref m 0 0) 1)
(setf (mref m 0 1) (* 2 (mref m 0 0)))
(incf (mref m 0 2) 4)
m)
==> #<MAT 2x3 AB #2A((1.0d0 2.0d0 4.0d0) (0.0d0 0.0d0 0.0d0))>
```
Compared to AREF MREF is a very expensive operation and it's best
used sparingly. Instead, typical code relies much more on matrix
level operations:
```commonlisp
(princ (scal! 2 (fill! 3 (make-mat 4))))
.. #<MAT 4 BF #(6.0d0 6.0d0 6.0d0 6.0d0)>
==> #<MAT 4 ABF #(6.0d0 6.0d0 6.0d0 6.0d0)>
```
The content of a matrix can be accessed in various representations
called *facets*. MGL-MAT takes care of synchronizing these facets
by copying data around lazily, but doing its best to share storage
for facets that allow it.
Notice the `ABF` in the printed results. It illustrates that behind
the scenes FILL! worked on the BACKING-ARRAY
facet (hence the `B`) that's basically a 1d lisp array. SCAL! on the
other hand made a foreign call to the BLAS `dscal` function for
which it needed the FOREIGN-ARRAY facet (`F`).
Finally, the `A` stands for the ARRAY facet that was
created when the array was printed. All facets are up-to-date (else
some of the characters would be lowercase). This is possible because
these three facets actually share storage which is never the case
for the CUDA-ARRAY facet. Now if we have a
CUDA-capable GPU, CUDA can be enabled with WITH-CUDA\*:
```commonlisp
(with-cuda* ()
(princ (scal! 2 (fill! 3 (make-mat 4)))))
.. #<MAT 4 C #(6.0d0 6.0d0 6.0d0 6.0d0)>
==> #<MAT 4 A #(6.0d0 6.0d0 6.0d0 6.0d0)>
```
Note the lonely `C` showing that only the CUDA-ARRAY
facet was used for both FILL! and SCAL!. When WITH-CUDA\* exits and
destroys the CUDA context, it destroys all CUDA facets, moving their
data to the ARRAY facet, so the returned MAT only has
that facet.
When there is no high-level operation that does what we want, we may
need to add new operations. This is usually best accomplished by
accessing one of the facets directly, as in the following example:
```commonlisp
(defun logdet (mat)
"Logarithm of the determinant of MAT. Return -1, 1 or 0 (or
equivalent) to correct for the sign, as the second value."
(with-facets ((array (mat 'array :direction :input)))
(lla:logdet array)))
```
Notice that LOGDET doesn't know about CUDA at all. WITH-FACETS
gives it the content of the matrix as a normal multidimensional lisp
array, copying the data from the GPU or elsewhere if necessary. This
allows new representations (`FACET`s) to be added easily and it also
avoids copying if the facet is already up-to-date. Of course, adding
CUDA support to LOGDET could make it more efficient.
Adding support for matrices that, for instance, live on a remote
machine is thus possible with a new facet type and existing code
would continue to work (albeit possibly slowly). Then one could
optimize the bottleneck operations by sending commands over the
network instead of copying data.
It is a bad idea to conflate resource management policy and
algorithms. MGL-MAT does its best to keep them separate.
## Basics
- [class] MAT CUBE
A MAT is a data CUBE that is much like a lisp
array, it supports DISPLACEMENT, arbitrary DIMENSIONS and
INITIAL-ELEMENT with the usual semantics. However, a MAT supports
different representations of the same data. See @MAT-TUTORIAL for
an introduction.
- [reader] MAT-CTYPE MAT (:CTYPE = \*DEFAULT-MAT-CTYPE\*)
One of *SUPPORTED-CTYPES*. The matrix can hold
only values of this type.
- [reader] MAT-DISPLACEMENT MAT (:DISPLACEMENT = 0)
A value in the `[0,MAX-SIZE]` interval. This is
like the DISPLACED-INDEX-OFFSET of a lisp array, but displacement
is relative to the start of the underlying storage vector.
- [reader] MAT-DIMENSIONS MAT (:DIMENSIONS)
Like ARRAY-DIMENSIONS. It holds a list of
dimensions, but it is allowed to pass in scalars too.
- [function] MAT-DIMENSION MAT AXIS-NUMBER
Return the dimension along AXIS-NUMBER. Similar to
ARRAY-DIMENSION.
- [reader] MAT-INITIAL-ELEMENT MAT (:INITIAL-ELEMENT = 0)
If non-nil, then when a facet is created, it is
filled with INITIAL-ELEMENT coerced to the appropriate numeric
type. If NIL, then no initialization is performed.
- [reader] MAT-SIZE MAT
The number of elements in the visible portion of
the array. This is always the product of the elements
MAT-DIMENSIONS and is similar to ARRAY-TOTAL-SIZE.
- [reader] MAT-MAX-SIZE MAT (:MAX-SIZE)
The number of elements for which storage may be
allocated. This is DISPLACEMENT + MAT-SIZE + `SLACK` where `SLACK`
is the number of trailing invisible elements.
- [function] MAKE-MAT DIMENSIONS &REST ARGS &KEY (CTYPE \*DEFAULT-MAT-CTYPE\*) (DISPLACEMENT 0) MAX-SIZE INITIAL-ELEMENT INITIAL-CONTENTS (SYNCHRONIZATION \*DEFAULT-SYNCHRONIZATION\*) DISPLACED-TO (CUDA-ENABLED \*DEFAULT-MAT-CUDA-ENABLED\*)
Return a new MAT object. If INITIAL-CONTENTS is given then the
matrix contents are initialized with REPLACE!. See class MAT for the
description of the rest of the parameters. This is exactly
what (MAKE-INSTANCE 'MAT ...) does except DIMENSIONS is not a
keyword argument so that MAKE-MAT looks more like MAKE-ARRAY. The
semantics of SYNCHRONIZATION are desribed in the
MGL-CUBE::@CUBE-SYNCHRONIZATION section.
If specified, DISPLACED-TO must be a MAT object large enough (in the
sense of its MAT-SIZE), to hold DISPLACEMENT plus `(REDUCE #'*
DIMENSIONS)` elements. Just like with MAKE-ARRAY, INITIAL-ELEMENT
and INITIAL-CONTENTS must not be supplied together with
DISPLACED-TO. See @MAT-SHAPING for more.
- [function] ARRAY-TO-MAT ARRAY &KEY CTYPE (SYNCHRONIZATION \*DEFAULT-SYNCHRONIZATION\*)
Create a MAT that's equivalent to ARRAY. Displacement of the
created array will be 0 and the size will be equal to
ARRAY-TOTAL-SIZE. If CTYPE is non-nil, then it will be the ctype of
the new matrix. Else ARRAY's type is converted to a ctype. If there
is no corresponding ctype, then *DEFAULT-MAT-CTYPE* is used.
Elements of ARRAY are coerced to CTYPE.
Also see MGL-CUBE::@CUBE-SYNCHRONIZATION.
- [function] MAT-TO-ARRAY MAT
- [function] REPLACE! MAT SEQ-OF-SEQS
Replace the contents of MAT with the elements of SEQ-OF-SEQS.
SEQ-OF-SEQS is a nested sequence of sequences similar to the
INITIAL-CONTENTS argument of MAKE-ARRAY. The total number of
elements must match the size of MAT. Returns MAT.
SEQ-OF-SEQS may contain multi-dimensional arrays as *leafs*, so the
following is legal:
```common-lisp
(replace! (make-mat '(1 2 3)) '(#2A((1 2 3) (4 5 6))))
==> #<MAT 1x2x3 AB #3A(((1.0d0 2.0d0 3.0d0) (4.0d0 5.0d0 6.0d0)))>
```
- [function] MREF MAT &REST INDICES
Like AREF for arrays. Don't use this if you care about performance
at all. SETFable. When set, the value is coerced to the ctype of MAT
with COERCE-TO-CTYPE. Note that currently MREF always operates on
the BACKING-ARRAY facet so it can trigger copying of facets. When
it's SETF'ed, however, it will update the CUDA-ARRAY if cuda is
enabled and it is up-to-date or there are no facets at all.
- [function] ROW-MAJOR-MREF MAT INDEX
Like ROW-MAJOR-AREF for arrays. Don't use this if you care about
performance at all. SETFable. When set, the value is coerced to the
ctype of MAT with COERCE-TO-CTYPE. Note that currently
ROW-MAJOR-MREF always operates on the BACKING-ARRAY facet so it can
trigger copying of facets. When it's SETF'ed, however, it will
update the CUDA-ARRAY if cuda is enabled and it is up-to-date or
there are no facets at all.
- [function] MAT-ROW-MAJOR-INDEX MAT &REST SUBSCRIPTS
Like ARRAY-ROW-MAJOR-INDEX for arrays.
## Element types
- [variable] *SUPPORTED-CTYPES* (:FLOAT :DOUBLE)
- [type] CTYPE
This is basically `(MEMBER :FLOAT :DOUBLE)`.
- [variable] *DEFAULT-MAT-CTYPE* :DOUBLE
By default MATs are created with this ctype. One of :FLOAT
or :DOUBLE.
- [function] COERCE-TO-CTYPE X &KEY (CTYPE \*DEFAULT-MAT-CTYPE\*)
Coerce the scalar X to the lisp type corresponding to CTYPE.
## Printing
- [variable] *PRINT-MAT* T
Controls whether the contents of a MAT object are printed as an
array (subject to the standard printer control variables).
- [variable] *PRINT-MAT-FACETS* T
Controls whether a summary of existing and up-to-date facets is
printed when a MAT object is printed. The summary that looks like
`ABcfh` indicates that all five facets (ARRAY,
BACKING-ARRAY, CUDA-ARRAY,
FOREIGN-ARRAY, CUDA-HOST-ARRAY) are
present and the first two are up-to-date. A summary of a single #-
indicates that there are no facets.
## Shaping
We are going to discuss various ways to change the visible portion
and dimensions of matrices. Conceptually a matrix has an *underlying
non-displaced storage vector*. For `(MAKE-MAT 10 :DISPLACEMENT
7 :MAX-SIZE 21)` this underlying vector looks like this:
displacement | visible elements | slack
. . . . . . . 0 0 0 0 0 0 0 0 0 0 . . . .
Whenever a matrix is reshaped (or *displaced to* in lisp
terminology), its displacement and dimensions change but the
underlying vector does not.
The rules for accessing displaced matrices is the same as always:
multiple readers can run in parallel, but attempts to write will
result in an error if there are either readers or writers on any of
the matrices that share the same underlying vector.
### Comparison to Lisp Arrays
One way to reshape and displace MAT objects is with MAKE-MAT and
its DISPLACED-TO argument whose semantics are similar to that of
MAKE-ARRAY in that the displacement is *relative* to the
displacement of DISPLACED-TO.
```commonlisp
(let* ((base (make-mat 10 :initial-element 5 :displacement 1))
(mat (make-mat 6 :displaced-to base :displacement 2)))
(fill! 1 mat)
(values base mat))
==> #<MAT 1+10+0 A #(5.0d0 5.0d0 1.0d0 1.0d0 1.0d0 1.0d0 1.0d0 1.0d0 5.0d0
--> 5.0d0)>
==> #<MAT 3+6+2 AB #(1.0d0 1.0d0 1.0d0 1.0d0 1.0d0 1.0d0)>
```
There are important semantic differences compared to lisp arrays all
which follow from the fact that displacement operates on the
underlying conceptual non-displaced vector.
- Matrices can be displaced and have slack even without DISPLACED-TO
just like `BASE` in the above example.
- It's legal to alias invisible elements of DISPLACED-TO as long as
the new matrix fits into the underlying storage.
- Negative displacements are allowed with DISPLACED-TO as long as
the adjusted displacement is non-negative.
- Further shaping operations can make invisible portions of the
DISPLACED-TO matrix visible by changing the displacement.
- In contrast to ARRAY-DISPLACEMENT, MAT-DISPLACEMENT only returns
an offset into the underlying storage vector.
### Functional Shaping
The following functions are collectively called the functional
shaping operations, since they don't alter their arguments in any
way. Still, since storage is aliased modification to the returned
matrix will affect the original.
- [function] RESHAPE-AND-DISPLACE MAT DIMENSIONS DISPLACEMENT
Return a new matrix of DIMENSIONS that aliases MAT's storage at
offset DISPLACEMENT. DISPLACEMENT 0 is equivalent to the start of
the storage of MAT regardless of MAT's displacement.
- [function] RESHAPE MAT DIMENSIONS
Return a new matrix of DIMENSIONS whose displacement is the same as
the displacement of MAT.
- [function] DISPLACE MAT DISPLACEMENT
Return a new matrix that aliases MAT's storage at offset
DISPLACEMENT. DISPLACEMENT 0 is equivalent to the start of the
storage of MAT regardless of MAT's displacement. The returned matrix
has the same dimensions as MAT.
### Destructive Shaping
The following destructive operations don't alter the contents of
the matrix, but change what is visible. ADJUST! is the odd one out,
it may create a new MAT.
- [function] RESHAPE-AND-DISPLACE! MAT DIMENSIONS DISPLACEMENT
Change the visible (or active) portion of MAT by altering its
displacement offset and dimensions. Future operations will only
affect this visible portion as if the rest of the elements were not
there. Return MAT.
DISPLACEMENT + the new size must not exceed MAT-MAX-SIZE.
Furthermore, there must be no facets being viewed (with WITH-FACETS)
when calling this function as the identity of the facets is not
stable.
- [function] RESHAPE! MAT DIMENSIONS
Like RESHAPE-AND-DISPLACE! but only alters the dimensions.
- [function] DISPLACE! MAT DISPLACEMENT
Like RESHAPE-AND-DISPLACE! but only alters the displacement.
- [function] RESHAPE-TO-ROW-MATRIX! MAT ROW
Reshape the 2d MAT to make only a single ROW visible. This is made
possible by the row-major layout, hence no column counterpart.
Return MAT.
- [macro] WITH-SHAPE-AND-DISPLACEMENT (MAT &OPTIONAL (DIMENSIONS NIL) (DISPLACEMENT NIL)) &BODY BODY
Reshape and displace MAT if DIMENSIONS and/or DISPLACEMENT is given
and restore the original shape and displacement after BODY is
executed. If neither is specificed, then nothing will be changed,
but BODY is still allowed to alter the shape and displacement.
- [function] ADJUST! MAT DIMENSIONS DISPLACEMENT &KEY (DESTROY-OLD-P T)
Like RESHAPE-AND-DISPLACE! but creates a new matrix if MAT isn't
large enough. If a new matrix is created, the contents are not
copied over and the old matrix is destroyed with DESTROY-CUBE if
DESTROY-OLD-P.
## Assembling
The functions here assemble a single MAT from a number of
MATs.
- [function] STACK! AXIS MATS MAT
Stack MATS along AXIS into MAT and return MAT. If AXIS is 0, place
MATS into MAT below each other starting from the top. If AXIS is 1,
place MATS side by side starting from the left. Higher AXIS are also
supported. All dimensions except for AXIS must be the same for all
MATS.
- [function] STACK AXIS MATS &KEY (CTYPE \*DEFAULT-MAT-CTYPE\*)
Like STACK! but return a new MAT of CTYPE.
```commonlisp
(stack 1 (list (make-mat '(3 2) :initial-element 0)
(make-mat '(3 1) :initial-element 1)))
==> #<MAT 3x3 B #2A((0.0d0 0.0d0 1.0d0)
--> (0.0d0 0.0d0 1.0d0)
--> (0.0d0 0.0d0 1.0d0))>
```
## Caching
Allocating and initializing a MAT object and its necessary facets
can be expensive. The following macros remember the previous value
of a binding in the same thread and /place/. Only weak references
are constructed so the cached objects can be garbage collected.
While the cache is global, thread safety is guaranteed by having
separate subcaches per thread. Each subcache is keyed by a /place/
object that's either explicitly specified or else is unique to each
invocation of the caching macro, so different occurrences of caching
macros in the source never share data. Still, recursion could lead
to data sharing between different invocations of the same function.
To prevent this, the cached object is removed from the cache while
it is used so other invocations will create a fresh one which isn't
particularly efficient but at least it's safe.
- [macro] WITH-THREAD-CACHED-MAT (VAR DIMENSIONS &REST ARGS &KEY (PLACE :SCRATCH) (CTYPE '\*DEFAULT-MAT-CTYPE\*) (DISPLACEMENT 0) MAX-SIZE (INITIAL-ELEMENT 0) INITIAL-CONTENTS) &BODY BODY
Bind VAR to a matrix of DIMENSIONS, CTYPE, etc. Cache this matrix,
and possibly reuse it later by reshaping it. When BODY exits the
cached object is updated with the binding of VAR which BODY may
change.
There is a separate cache for each thread and each `PLACE` (under
EQ). Since every cache holds exactly one MAT per CTYPE, nested
WITH-THREAD-CACHED-MAT often want to use different `PLACE`s. By
convention, these places are called `:SCRATCH-1`, `:SCRATCH-2`,
etc.
- [macro] WITH-THREAD-CACHED-MATS SPECS &BODY BODY
A shorthand for writing nested WITH-THREAD-CACHED-MAT calls.
```
(with-thread-cached-mat (a ...)
(with-thread-cached-mat (b ...)
...))
```
is equivalent to:
```
(with-thread-cached-mat ((a ...)
(b ...))
...)
```
- [macro] WITH-ONES (VAR DIMENSIONS &KEY (CTYPE '\*DEFAULT-MAT-CTYPE\*) (PLACE :ONES)) &BODY BODY
Bind VAR to a matrix of DIMENSIONS whose every element is 1. The
matrix is cached for efficiency.
## BLAS Operations
Only some BLAS functions are implemented, but it should be easy to
add more as needed. All of them default to using CUDA, if it is
initialized and enabled (see USE-CUDA-P).
Level 1 BLAS operations
- [function] ASUM X &KEY (N (MAT-SIZE X)) (INCX 1)
Return the l1 norm of X, that is, sum of the absolute values of its
elements.
- [function] AXPY! ALPHA X Y &KEY (N (MAT-SIZE X)) (INCX 1) (INCY 1)
Set Y to ALPHA \* X + Y. Return Y.
- [function] COPY! X Y &KEY (N (MAT-SIZE X)) (INCX 1) (INCY 1)
Copy X into Y. Return Y.
- [function] DOT X Y &KEY (N (MAT-SIZE X)) (INCX 1) (INCY 1)
Return the dot product of X and Y.
- [function] NRM2 X &KEY (N (MAT-SIZE X)) (INCX 1)
Return the l2 norm of X, which is the square root of the sum of the
squares of its elements.
- [function] SCAL! ALPHA X &KEY (N (MAT-SIZE X)) (INCX 1)
Set X to ALPHA \* X. Return X.
Level 3 BLAS operations
- [function] GEMM! ALPHA A B BETA C &KEY TRANSPOSE-A? TRANSPOSE-B? M N K LDA LDB LDC
Basically C = ALPHA \* A' \* B' + BETA \* C. A' is A or its transpose
depending on TRANSPOSE-A?. B' is B or its transpose depending on
TRANSPOSE-B?. Returns C.
A' is an MxK matrix. B' is a KxN matrix. C is an MxN matrix.
LDA is the width of the matrix A (not of A'). If A is not transposed,
then K <= LDA, if it's transposed then M <= LDA.
LDB is the width of the matrix B (not of B'). If B is not transposed,
then N <= LDB, if it's transposed then K <= LDB.
In the example below M=3, N=2, K=5, LDA=6, LDB=3, LDC=4. The cells
marked with + do not feature in the calculation.
N
--+
--+
K -B+
--+
--+
+++
K
-----+ --++
M --A--+ -C++
-----+ --++
++++++ ++++
## Destructive API
- [function] .SQUARE! X &KEY (N (MAT-SIZE X))
Set X to its elementwise square. Return X.
- [function] .SQRT! X &KEY (N (MAT-SIZE X))
Set X to its elementwise square root. Return X.
- [function] .LOG! X &KEY (N (MAT-SIZE X))
Set X to its elementwise natural logarithm. Return X.
- [function] .EXP! X &KEY (N (MAT-SIZE X))
Apply EXP elementwise to X in a destructive manner. Return X.
- [function] .EXPT! X POWER
Raise matrix X to POWER in an elementwise manner. Return X. Note
that CUDA and non-CUDA implementations may disagree on the treatment
of NaNs, infinities and complex results. In particular, the lisp
implementation always computes the REALPART of the results while
CUDA's pow() returns NaNs instead of complex numbers.
- [function] .INV! X &KEY (N (MAT-SIZE X))
Set X to its elementwise inverse `(/ 1 X)`. Return X.
- [function] .LOGISTIC! X &KEY (N (MAT-SIZE X))
Destructively apply the logistic function to X in an elementwise
manner. Return X.
- [function] .+! ALPHA X
Add the scalar ALPHA to each element of X destructively modifying
X. Return X.
- [function] .*! X Y
- [function] GEEM! ALPHA A B BETA C
Like GEMM!, but multiplication is elementwise. This is not a
standard BLAS routine.
- [function] GEERV! ALPHA A X BETA B
GEneric Elementwise Row - Vector multiplication. `B = beta * B +
alpha a .* X*` where `X*` is a matrix of the same shape as A whose
every row is X. Perform elementwise multiplication on each row of A
with the vector X and add the scaled result to the corresponding row
of B. Return B. This is not a standard BLAS routine.
- [function] .<! X Y
For each element of X and Y set Y to 1 if the element in Y is
greater than the element in X, and to 0 otherwise. Return Y.
- [function] .MIN! ALPHA X
Set each element of X to ALPHA if it's greater than ALPHA. Return
X.
- [function] .MAX! ALPHA X
Set each element of X to ALPHA if it's less than ALPHA. Return X.
- [function] ADD-SIGN! ALPHA A BETA B
Add the elementwise sign (-1, 0 or 1 for negative, zero and
positive numbers respectively) of A times ALPHA to BETA \* B. Return
B.
- [function] FILL! ALPHA X &KEY (N (MAT-SIZE X))
Fill matrix X with ALPHA. Return X.
- [function] SUM! X Y &KEY AXIS (ALPHA 1) (BETA 0)
Sum matrix X along AXIS and add ALPHA \* SUMS to BETA \* Y
destructively modifying Y. Return Y. On a 2d matrix (nothing else is
supported currently), if AXIS is 0, then columns are summed, if AXIS
is 1 then rows are summed.
- [function] SCALE-ROWS! SCALES A &KEY (RESULT A)
Set RESULT to `DIAG(SCALES)*A` and return it. `A` is an `MxN`
matrix, SCALES is treated as a length `M` vector.
- [function] SCALE-COLUMNS! SCALES A &KEY (RESULT A)
Set RESULT to `A*DIAG(SCALES)` and return it. `A` is an `MxN`
matrix, SCALES is treated as a length `N` vector.
- [function] .SIN! X &KEY (N (MAT-SIZE X))
Apply SIN elementwise to X in a destructive manner. Return X.
- [function] .COS! X &KEY (N (MAT-SIZE X))
Apply COS elementwise to X in a destructive manner. Return X.
- [function] .TAN! X &KEY (N (MAT-SIZE X))
Apply TAN elementwise to X in a destructive manner. Return X.
- [function] .SINH! X &KEY (N (MAT-SIZE X))
Apply SINH elementwise to X in a destructive manner. Return X.
- [function] .COSH! X &KEY (N (MAT-SIZE X))
Apply COSH elementwise to X in a destructive manner. Return X.
- [function] .TANH! X &KEY (N (MAT-SIZE X))
Apply TANH elementwise to X in a destructive manner. Return X.
Finally, some neural network operations.
- [function] CONVOLVE! X W Y &KEY START STRIDE ANCHOR BATCHED
Y = Y + conv(X, W) and return Y. If BATCHED, then the first
dimension of X and Y is the number of elements in the batch (B),
else B is assumed to be 1. The rest of the dimensions encode the
input (X) and output (Y\} N dimensional feature maps. START, STRIDE
and ANCHOR are lists of length N. START is the multi-dimensional
index of the first element of the input feature map (for each
element in the batch) for which the convolution must be computed.
Then (ELT STRIDE (- N 1)) is added to the last element of START and
so on until (ARRAY-DIMENSION X 1) is reached. Then the last element
of START is reset, (ELT STRIDE (- N 2)) is added to the first but
last element of START and we scan the last dimension again. Take a
2d example, START is (0 0), STRIDE is (1 2), and X is a B\*2x7
matrix.
W is:
1 2 1
2 4 2
1 2 1
and ANCHOR is (1 1) which refers to the element of W whose value is
4. This anchor point of W is placed over elements of X whose multi
dimensional index is in numbers in this figure (only one element in
the batch is shown):
0,0 . 0,2 . 0,4 . 0,6
1,0 . 1,2 . 1,4 . 1,6
When applying W at position P of X, the convolution is the sum of
the products of overlapping elements of X and W when W's ANCHOR is
placed at P. Elements of W over the edges of X are multiplied with 0
so are effectively ignored. The order of application of W to
positions defined by START, STRIDE and ANCHOR is undefined.
Y must be a B\*2x4 (or 2x4 if not BATCHED) matrix in this example,
just large enough to hold the results of the convolutions.
- [function] DERIVE-CONVOLVE! X XD W WD YD &KEY START STRIDE ANCHOR BATCHED
Add the dF/dX to XD and and dF/dW to WD where YD is dF/dY for some
function F where Y is the result of convolution with the same
arguments.
- [function] MAX-POOL! X Y &KEY START STRIDE ANCHOR BATCHED POOL-DIMENSIONS
- [function] DERIVE-MAX-POOL! X XD Y YD &KEY START STRIDE ANCHOR BATCHED POOL-DIMENSIONS
Add the dF/dX to XD and and dF/dW to WD where YD is dF/dY for some
function F where Y is the result of MAX-POOL! with the same
arguments.
## Non-destructive API
- [function] COPY-MAT A
Return a copy of the active portion with regards to displacement
and shape of A.
- [function] COPY-ROW A ROW
Return ROW of A as a new 1d matrix.
- [function] COPY-COLUMN A COLUMN
Return COLUMN of A as a new 1d matrix.
- [function] MAT-AS-SCALAR A
Return the first element of A. A must be of size 1.
- [function] SCALAR-AS-MAT X &KEY (CTYPE (LISP-\>CTYPE (TYPE-OF X)))
Return a matrix of one dimension and one element: X. CTYPE, the
type of the matrix, defaults to the ctype corresponding to the type
of X.
- [function] M= A B
Check whether A and B, which must be matrices of the same size, are
elementwise equal.
- [function] TRANSPOSE A
Return the transpose of A.
- [function] M* A B &KEY TRANSPOSE-A? TRANSPOSE-B?
Compute op(A) \* op(B). Where op is either the identity or the
transpose operation depending on TRANSPOSE-A? and TRANSPOSE-B?.
- [function] MM* M &REST ARGS
Convenience function to multiply several matrices.
(mm\* a b c) => a \* b \* c
- [function] M- A B
Return A - B.
- [function] M+ A B
Return A + B.
- [function] INVERT A
Return the inverse of A.
- [function] LOGDET MAT
Logarithm of the determinant of MAT. Return -1, 1 or 0 (or
equivalent) to correct for the sign, as the second value.
## Mappings
- [function] MAP-CONCAT FN MATS MAT &KEY KEY PASS-RAW-P
Call FN with each element of MATS and MAT temporarily reshaped to
the dimensions of the current element of MATS and return MAT. For
the next element the displacement is increased so that there is no
overlap.
MATS is keyed by KEY just like the CL sequence functions. Normally,
FN is called with the matrix returned by KEY. However, if
PASS-RAW-P, then the matrix returned by KEY is only used to
calculate dimensions and the element of MATS that was passed to KEY
is passed to FN, too.
```
(map-concat #'copy! (list (make-mat 2) (make-mat 4 :initial-element 1))
(make-mat '(2 3)))
==> #<MAT 2x3 AB #2A((0.0d0 0.0d0 1.0d0) (1.0d0 1.0d0 1.0d0))>
```
- [function] MAP-DISPLACEMENTS FN MAT DIMENSIONS &KEY (DISPLACEMENT-START 0) DISPLACEMENT-STEP
Call FN with MAT reshaped to DIMENSIONS, first displaced by
DISPLACEMENT-START that's incremented by DISPLACEMENT-STEP each
iteration while there are enough elements left for DIMENSIONS at the
current displacement. Returns MAT.
```commonlisp
(let ((mat (make-mat 14 :initial-contents '(-1 0 1 2 3
4 5 6 7
8 9 10 11 12))))
(reshape-and-displace! mat '(4 3) 1)
(map-displacements #'print mat 4))
..
.. #<MAT 1+4+9 B #(0.0d0 1.0d0 2.0d0 3.0d0)>
.. #<MAT 5+4+5 B #(4.0d0 5.0d0 6.0d0 7.0d0)>
.. #<MAT 9+4+1 B #(8.0d0 9.0d0 10.0d0 11.0d0)>
```
- [function] MAP-MATS-INTO RESULT-MAT FN &REST MATS
Like CL:MAP-INTO but for MAT objects. Destructively modifies
RESULT-MAT to contain the results of applying FN to each element in
the argument MATS in turn.
## Random numbers
Unless noted these work efficiently with CUDA.
- [generic-function] COPY-RANDOM-STATE STATE
Return a copy of STATE be it a lisp or cuda random
state.
- [function] UNIFORM-RANDOM! MAT &KEY (LIMIT 1)
Fill MAT with random numbers sampled uniformly from the \[0,LIMIT)
interval of MAT's type.
- [function] GAUSSIAN-RANDOM! MAT &KEY (MEAN 0) (STDDEV 1)
Fill MAT with independent normally distributed random numbers with
MEAN and STDDEV.
- [function] MV-GAUSSIAN-RANDOM &KEY MEANS COVARIANCES
Return a column vector of samples from the multivariate normal
distribution defined by MEANS (Nx1) and COVARIANCES (NxN). No CUDA
implementation.
- [function] ORTHOGONAL-RANDOM! M &KEY (SCALE 1)
Fill the matrix M with random values in such a way that `M^T * M`
is the identity matrix (or something close if M is wide). Return M.
## I/O
- [variable] *MAT-HEADERS* T
If true, a header with MAT-CTYPE and MAT-SIZE is written by
WRITE-MAT before the contents and READ-MAT checks that these match
the matrix into which it is reading.
- [generic-function] WRITE-MAT MAT STREAM
Write MAT to binary STREAM in portable binary
format. Return MAT. Displacement and size are taken into account,
only visible elements are written. Also see *MAT-HEADERS*.
- [generic-function] READ-MAT MAT STREAM
Destructively modify the visible portion (with
regards to displacement and shape) of MAT by reading MAT-SIZE number
of elements from binary STREAM. Return MAT. Also see
*MAT-HEADERS*.
## Debugging
The largest class of bugs has to do with synchronization of facets
being broken. This is almost always caused by an operation that
mispecifies the DIRECTION argument of WITH-FACET. For example, the
matrix argument of SCAL! should be accessed with direciton :IO. But
if it's :INPUT instead, then subsequent access to the ARRAY facet
will not see the changes made by AXPY!, and if it's :OUTPUT, then
any changes made to the ARRAY facet since the last update of the
CUDA-ARRAY facet will not be copied and from the wrong input SCAL!
will compute the wrong result.
Using the SLIME inspector or trying to access the
CUDA-ARRAY facet from threads other than the one in
which the corresponding CUDA context was initialized will fail. For
now, the easy way out is to debug the code with CUDA disabled (see
*CUDA-ENABLED*).
Another thing that tends to come up is figuring out where memory is
used.
- [function] MAT-ROOM &KEY (STREAM \*STANDARD-OUTPUT\*) (VERBOSE T)
Calls FOREIGN-ROOM and CUDA-ROOM.
- [macro] WITH-MAT-COUNTERS (&KEY COUNT N-BYTES) &BODY BODY
Count all MAT allocations and also the number of bytes they may
require. *May require* here really means an upper bound,
because `(MAKE-MAT (EXPT 2 60))` doesn't actually uses memory until
one of its facets is accessed (don't simply evaluate it though,
printing the result will access the ARRAY facet if
*PRINT-MAT*). Also, while facets today all require the same number
of bytes, this may change in the future. This is a debugging tool,
don't use it in production.
```common-lisp
(with-mat-counters (:count count :n-bytes n-bytes)
(assert (= count 0))
(assert (= n-bytes 0))
(make-mat '(2 3) :ctype :double)
(assert (= count 1))
(assert (= n-bytes (* 2 3 8)))
(with-mat-counters (:n-bytes n-bytes-1 :count count-1)
(make-mat '7 :ctype :float)
(assert (= count-1 1))
(assert (= n-bytes-1 (* 7 4))))
(assert (= n-bytes (+ (* 2 3 8) (* 7 4))))
(assert (= count 2)))
```
## Facet API
### Facets
A MAT is a CUBE (see MGL-CUBE::@CUBE-MANUAL) whose facets are
different representations of numeric arrays. These facets can be
accessed with WITH-FACETS with one of the following
FACET-NAME locatives:
- [facet-name] BACKING-ARRAY
The corresponding facet's value is a one dimensional lisp array or
a static vector that also looks exactly like a lisp array but is
allocated in foreign memory. See *FOREIGN-ARRAY-STRATEGY*.
- [facet-name] ARRAY
Same as BACKING-ARRAY if the matrix is one-dimensional, all
elements are visible (see @MAT-SHAPING), else it's a lisp array
displaced to the backing array.
- [facet-name] FOREIGN-ARRAY