| title |
|---|
Multidimensional arrays |
Common Lisp has native support for multidimensional arrays, with some
special treatment for 1-D arrays, called vectors. Arrays can be
generalised and contain any type (element-type t), or they
can be specialised to contain specific types such as single-float
or integer. A good place to start is
Practical Common Lisp Chapter 11, Collections by
Peter Seibel.
A quick reference to some common operations on arrays is given in the section on Arrays and vectors.
Some libraries available on Quicklisp for manipulating arrays:
- array-operations defines
functions
generate,permute,displace,flatten,split,combine,reshape. It also defineseach, for element-wise operations. This library is not maintained by the original author, but there is an actively maintained fork. - cmu-infix includes array indexing syntax for multidimensional arrays.
- lla is a library for linear algebra, calling BLAS and LAPACK libraries. It differs from most CL linear algebra packages in using intuitive function names, and can operate on native arrays as well as CLOS objects.
This page covers what can be done with the built-in multidimensional arrays, but there are limitations. In particular:
- Interoperability with foreign language arrays, for example when calling libraries such as BLAS, LAPACK or GSL.
- Extending arithmetic and other mathematical operators to handle
arrays, for example so that
(+ a b)works whenaand/orbare arrays.
Both of these problems can be solved by using CLOS to define an extended array class, with native arrays as a special case. Some libraries available through quicklisp which take this approach are:
- matlisp, some of which is described in sections below.
- MGL-MAT, which has a manual and provides bindings to BLAS and CUDA. This is used in a machine learning library MGL.
- cl-ana, a data analysis package with a manual, which includes operations on arrays.
- Antik, used in GSLL, a binding to the GNU Scientific Library.
A relatively new but actively developed package is
MAGICL, which provides
wrappers around BLAS and LAPACK libraries. At the time of writing this
package is not on Quicklisp, and only works under SBCL and CCL. It
seems to be particularly focused on complex arrays, but not
exclusively.
To install, clone the repository in your quicklisp local-projects
directory e.g. under Linux/Unix:
$ cd ~/quicklisp/local-projects
$ git clone https://github.com/rigetticomputing/magicl.gitInstructions for installing dependencies (BLAS, LAPACK and Expokit)
are given on the github web pages.
Low-level routines wrap foreign functions, so have the Fortran names
e.g magicl.lapack-cffi::%zgetrf. Higher-level interfaces to some of
these functions also exist, see the
source directory and documentation.
Taking this further, domain specific languages have been built on Common Lisp, which can be used for numerical calculations with arrays. At the time of writing the most widely used and supported of these are:
CLASP is a project which aims to ease interoperability of Common Lisp with other languages (particularly C++), by using LLVM. One of the main applications of this project is to numerical/scientific computing.
The function CLHS: make-array can create arrays filled with a single value
* (defparameter *my-array* (make-array '(3 2) :initial-element 1.0))
*MY-ARRAY*
* *my-array*
#2A((1.0 1.0) (1.0 1.0) (1.0 1.0))More complicated array values can be generated by first making an array, and then iterating over the elements to fill in the values (see section below on element access).
The array-operations
library provides generate, a convenient
function for creating arrays which wraps this iteration.
* (ql:quickload :array-operations)
To load "array-operations":
Load 1 ASDF system:
array-operations
; Loading "array-operations"
(:ARRAY-OPERATIONS)
* (aops:generate #'identity 7 :position)
#(0 1 2 3 4 5 6)Note that the nickname for array-operations is aops. The
generate function can also iterate over the array subscripts by
passing the key :subscripts. See the
github repository for
more examples.
To create an 3x3 array containing random numbers drawn from a uniform
distribution, generate can be used to call the CL
random function:
* (aops:generate (lambda () (random 1.0)) '(3 3))
#2A((0.99292254 0.929777 0.93538976)
(0.31522608 0.45167792 0.9411855)
(0.96221936 0.9143338 0.21972346))An array of Gaussian (normal) random numbers with mean of zero and standard deviation of one, using the alexandria package:
* (ql:quickload :alexandria)
To load "alexandria":
Load 1 ASDF system:
alexandria
; Loading "alexandria"
(:ALEXANDRIA)
* (aops:generate #'alexandria:gaussian-random 4)
#(0.5522547885338768d0 -1.2564808468164517d0 0.9488161476129733d0
-0.10372852118266523d0)Note that this is not particularly efficient: It requires a function
call for each element, and although gaussian-random returns
two random numbers, only one of them is used.
For more efficient implementations, and a wider range of probability distributions, there are packages available on Quicklisp. See CLiki for a list.
To access the individual elements of an array there are the aref and row-major-aref functions.
The aref function takes the same number of index arguments as the array has dimensions. Indexing is from 0 and row-major as in C, but not Fortran.
* (defparameter *a* #(1 2 3 4))
*A*
* (aref *a* 0)
1
* (aref *a* 3)
4
* (defparameter *b* #2A((1 2 3) (4 5 6)))
*B*
* (aref *b* 1 0)
4
* (aref *b* 0 2)
3The range of these indices can be found using array-dimensions:
* (array-dimensions *a*)
(4)
* (array-dimensions *b*)
(2 3)
or the rank of the array can be found, and then the size of each dimension queried:
* (array-rank *a*)
1
* (array-dimension *a* 0)
4
* (array-rank *b*)
2
* (array-dimension *b* 0)
2
* (array-dimension *b* 1)
3To loop over an array nested loops can be used, such as
* (defparameter a #2A((1 2 3) (4 5 6)))
A
* (destructuring-bind (n m) (array-dimensions a)
(loop for i from 0 below n do
(loop for j from 0 below m do
(format t "a[~a ~a] = ~a~%" i j (aref a i j)))))
a[0 0] = 1
a[0 1] = 2
a[0 2] = 3
a[1 0] = 4
a[1 1] = 5
a[1 2] = 6
NILA utility macro which does this for multiple dimensions is nested-loop:
(defmacro nested-loop (syms dimensions &body body)
"Iterates over a multidimensional range of indices.
SYMS must be a list of symbols, with the first symbol
corresponding to the outermost loop.
DIMENSIONS will be evaluated, and must be a list of
dimension sizes, of the same length as SYMS.
Example:
(nested-loop (i j) '(10 20) (format t '~a ~a~%' i j))
"
(unless syms (return-from nested-loop `(progn ,@body))) ; No symbols
;; Generate gensyms for dimension sizes
(let* ((rank (length syms))
(syms-rev (reverse syms)) ; Reverse, since starting with innermost
(dims-rev (loop for i from 0 below rank collecting (gensym))) ; innermost dimension first
(result `(progn ,@body))) ; Start with innermost expression
;; Wrap previous result inside a loop for each dimension
(loop for sym in syms-rev for dim in dims-rev do
(unless (symbolp sym) (error "~S is not a symbol. First argument to nested-loop must be a list of symbols" sym))
(setf result
`(loop for ,sym from 0 below ,dim do
,result)))
;; Add checking of rank and dimension types, and get dimensions into gensym list
(let ((dims (gensym)))
`(let ((,dims ,dimensions))
(unless (= (length ,dims) ,rank) (error "Incorrect number of dimensions: Expected ~a but got ~a" ,rank (length ,dims)))
(dolist (dim ,dims)
(unless (integerp dim) (error "Dimensions must be integers: ~S" dim)))
(destructuring-bind ,(reverse dims-rev) ,dims ; Dimensions reversed so that innermost is last
,result)))))
so that the contents of a 2D array can be printed using
* (defparameter a #2A((1 2 3) (4 5 6)))
A
* (nested-loop (i j) (array-dimensions a)
(format t "a[~a ~a] = ~a~%" i j (aref a i j)))
a[0 0] = 1
a[0 1] = 2
a[0 2] = 3
a[1 0] = 4
a[1 1] = 5
a[1 2] = 6
NIL[Note: This macro is available in this fork of array-operations, but not Quicklisp]
In some cases, particularly element-wise operations, the number of dimensions does not matter. To write code which is independent of the number of dimensions, array element access can be done using a single flattened index via row-major-aref. The array size is given by array-total-size, with the flattened index starting at 0.
* (defparameter a #2A((1 2 3) (4 5 6)))
A
* (array-total-size a)
6
* (loop for i from 0 below (array-total-size a) do
(setf (row-major-aref a i) (+ 2.0 (row-major-aref a i))))
NIL
* a
#2A((3.0 4.0 5.0) (6.0 7.0 8.0))The cmu-infix library provides some different syntax which can make mathematical expressions easier to read:
* (ql:quickload :cmu-infix)
To load "cmu-infix":
Load 1 ASDF system:
cmu-infix
; Loading "cmu-infix"
(:CMU-INFIX)
* (named-readtables:in-readtable cmu-infix:syntax)
(("COMMON-LISP-USER" . #<NAMED-READTABLE CMU-INFIX:SYNTAX {10030158B3}>)
...)
* (defparameter arr (make-array '(3 2) :initial-element 1.0))
ARR
* #i(arr[0 1] = 2.0)
2.0
* arr
#2A((1.0 2.0) (1.0 1.0) (1.0 1.0))A matrix-matrix multiply operation can be implemented as:
(let ((A #2A((1 2) (3 4)))
(B #2A((5 6) (7 8)))
(result (make-array '(2 2) :initial-element 0.0)))
(loop for i from 0 to 1 do
(loop for j from 0 to 1 do
(loop for k from 0 to 1 do
#i(result[i j] += A[i k] * B[k j]))))
result)See the section below on linear algebra, for alternative matrix-multiply implementations.
To multiply two arrays of numbers of the same size, pass a function
to each in the array-operations library:
* (aops:each #'* #(1 2 3) #(2 3 4))
#(2 6 12)For improved efficiency there is the aops:each* function, which
takes a type as first argument to specialise the result array.
To add a constant to all elements of an array:
* (defparameter *a* #(1 2 3 4))
*A*
* (aops:each (lambda (it) (+ 42 it)) *a*)
#(43 44 45 46)
* *a*
#(1 2 3 4)Note that each is not destructive, but makes a new array.
All arguments to each must be arrays of the same size,
so (aops:each #'+ 42 *a*) is not valid.
An alternative approach to the each function above, is to use a
macro to iterate over all elements of an array:
(defmacro vectorize (variables &body body)
;; Check that variables is a list of only symbols
(dolist (var variables)
(if (not (symbolp var))
(error "~S is not a symbol" var)))
;; Get the size of the first variable, and create a new array
;; of the same type for the result
`(let ((size (array-total-size ,(first variables))) ; Total array size (same for all variables)
(result (make-array (array-dimensions ,(first variables)) ; Returned array
:element-type (array-element-type ,(first variables)))))
;; Check that all variables have the same sizeo
,@(mapcar (lambda (var) `(if (not (equal (array-dimensions ,(first variables))
(array-dimensions ,var)))
(error "~S and ~S have different dimensions" ',(first variables) ',var)))
(rest variables))
(dotimes (indx size)
;; Locally redefine variables to be scalars at a given index
(let ,(mapcar (lambda (var) (list var `(row-major-aref ,var indx))) variables)
;; User-supplied function body now evaluated for each index in turn
(setf (row-major-aref result indx) (progn ,@body))))
result))[Note: Expanded versions of this macro are available in this fork of array-operations, but not Quicklisp]
This can be used as:
* (defparameter *a* #(1 2 3 4))
*A*
* (vectorize (*a*) (* 2 *a*))
#(2 4 6 8)Inside the body of the expression (second form in vectorize
expression) the symbol *a* is bound to a single element.
This means that the built-in mathematical functions can be used:
* (defparameter a #(1 2 3 4))
A
* (defparameter b #(2 3 4 5))
B
* (vectorize (a b) (* a (sin b)))
#(0.9092974 0.28224 -2.2704074 -3.8356972)and combined with cmu-infix
* (vectorize (a b) #i(a * sin(b)) )
#(0.9092974 0.28224 -2.2704074 -3.8356972)Several packages provide wrappers around BLAS, for fast matrix manipulation.
The lla package in quicklisp includes calls to some functions:
scaling by a constant factor:
* (defparameter a #(1 2 3))
* (lla:scal! 2.0 a)
* a
#(2.0d0 4.0d0 6.0d0)This calculates a * x + y where a is a constant, x and y are arrays.
The lla:axpy! function is destructive, modifying the last argument (y).
* (defparameter x #(1 2 3))
A
* (defparameter y #(2 3 4))
B
* (lla:axpy! 0.5 x y)
#(2.5d0 4.0d0 5.5d0)
* x
#(1.0d0 2.0d0 3.0d0)
* y
#(2.5d0 4.0d0 5.5d0)If the y array is complex, then this operation calls the complex
number versions of these operators:
* (defparameter x #(1 2 3))
* (defparameter y (make-array 3 :element-type '(complex double-float)
:initial-element #C(1d0 1d0)))
* y
#(#C(1.0d0 1.0d0) #C(1.0d0 1.0d0) #C(1.0d0 1.0d0))
* (lla:axpy! #C(0.5 0.5) a b)
#(#C(1.5d0 1.5d0) #C(2.0d0 2.0d0) #C(2.5d0 2.5d0))The dot product of two vectors:
* (defparameter x #(1 2 3))
* (defparameter y #(2 3 4))
* (lla:dot x y)
20.0d0The reduce function operates on sequences, including vectors (1D arrays), but not on multidimensional arrays. To get around this, multidimensional arrays can be displaced to create a 1D vector. Displaced arrays share storage with the original array, so this is a fast operation which does not require copying data:
* (defparameter a #2A((1 2) (3 4)))
A
* (reduce #'max (make-array (array-total-size a) :displaced-to a))
4The array-operations package contains flatten, which returns a
displaced array i.e doesn't copy data:
* (reduce #'max (aops:flatten a))An SBCL extension, array-storage-vector provides an efficient but not portable way to achieve the same thing:
* (reduce #'max (array-storage-vector a))
4More complex reductions are sometimes needed, for example finding the maximum absolute difference between two arrays. Using the above methods we could do:
* (defparameter a #2A((1 2) (3 4)))
A
* (defparameter b #2A((1 3) (5 4)))
B
* (reduce #'max (aops:flatten
(aops:each
(lambda (a b) (abs (- a b))) a b)))
2This involves allocating an array to hold the intermediate result,
which for large arrays could be inefficient. Similarly to vectorize
defined above, a macro which does not allocate can be defined as:
(defmacro vectorize-reduce (fn variables &body body)
"Performs a reduction using FN over all elements in a vectorized expression
on array VARIABLES.
VARIABLES must be a list of symbols bound to arrays.
Each array must have the same dimensions. These are
checked at compile and run-time respectively.
"
;; Check that variables is a list of only symbols
(dolist (var variables)
(if (not (symbolp var))
(error "~S is not a symbol" var)))
(let ((size (gensym)) ; Total array size (same for all variables)
(result (gensym)) ; Returned value
(indx (gensym))) ; Index inside loop from 0 to size
;; Get the size of the first variable
`(let ((,size (array-total-size ,(first variables))))
;; Check that all variables have the same size
,@(mapcar (lambda (var) `(if (not (equal (array-dimensions ,(first variables))
(array-dimensions ,var)))
(error "~S and ~S have different dimensions" ',(first variables) ',var)))
(rest variables))
;; Apply FN with the first two elements (or fewer if size < 2)
(let ((,result (apply ,fn (loop for ,indx below (min ,size 2) collecting
(let ,(map 'list (lambda (var) (list var `(row-major-aref ,var ,indx))) variables)
(progn ,@body))))))
;; Loop over the remaining indices
(loop for ,indx from 2 below ,size do
;; Locally redefine variables to be scalars at a given index
(let ,(mapcar (lambda (var) (list var `(row-major-aref ,var ,indx))) variables)
;; User-supplied function body now evaluated for each index in turn
(setf ,result (funcall ,fn ,result (progn ,@body)))))
,result))))
[Note: This macro is available in this fork of array-operations, but not Quicklisp]
Using this macro, the maximum value in an array A (of any shape) is:
* (vectorize-reduce #'max (a) a)The maximum absolute difference between two arrays A and B, of any shape as long as they have the same shape, is:
* (vectorize-reduce #'max (a b) (abs (- a b)))Several packages provide bindings to BLAS and LAPACK libraries, including:
A longer list of available packages is on CLiki's linear algebra page.
In the examples below the lla package is loaded:
* (ql:quickload :lla)
To load "lla":
Load 1 ASDF system:
lla
; Loading "lla"
.
(:LLA)The lla function mm performs
vector-vector, matrix-vector and matrix-matrix
multiplication.
Note that one vector is treated as a row vector, and the other as column:
* (lla:mm #(1 2 3) #(2 3 4))
20* (lla:mm #2A((1 1 1) (2 2 2) (3 3 3)) #(2 3 4))
#(9.0d0 18.0d0 27.0d0)which has performed the sum over j of A[i j] * x[j]
* (lla:mm #2A((1 2 3) (1 2 3) (1 2 3)) #2A((2 3 4) (2 3 4) (2 3 4)))
#2A((12.0d0 18.0d0 24.0d0) (12.0d0 18.0d0 24.0d0) (12.0d0 18.0d0 24.0d0))which summed over j in A[i j] * B[j k]
Note that the type of the returned arrays are simple arrays,
specialised to element type double-float
* (type-of (lla:mm #2A((1 0 0) (0 1 0) (0 0 1)) #(1 2 3)))
(SIMPLE-ARRAY DOUBLE-FLOAT (3))The array-operations package contains a generalised outer product function:
* (ql:quickload :array-operations)
To load "array-operations":
Load 1 ASDF system:
array-operations
; Loading "array-operations"
(:ARRAY-OPERATIONS)
* (aops:outer #'* #(1 2 3) #(2 3 4))
#2A((2 3 4) (4 6 8) (6 9 12))which has created a new 2D array A[i j] = B[i] * C[j]. This outer
function can take an arbitrary number of inputs, and inputs with
multiple dimensions.
The direct inverse of a dense matrix can be calculated with invert
* (lla:invert #2A((1 0 0) (0 1 0) (0 0 1)))
#2A((1.0d0 0.0d0 -0.0d0) (0.0d0 1.0d0 -0.0d0) (0.0d0 0.0d0 1.0d0))e.g
* (defparameter a #2A((1 2 3) (0 2 1) (1 3 2)))
A
* (defparameter b (lla:invert a))
B
* (lla:mm a b)
#2A((1.0d0 2.220446049250313d-16 0.0d0)
(0.0d0 1.0d0 0.0d0)
(0.0d0 1.1102230246251565d-16 0.9999999999999998d0))Calculating the direct inverse is generally not advisable, particularly for large matrices. Instead the LU decomposition can be calculated and used for multiple inversions.
* (defparameter a #2A((1 2 3) (0 2 1) (1 3 2)))
A
* (defparameter b (lla:mm a #(1 2 3)))
B
* (lla:solve (lla:lu a) b)
#(1.0d0 2.0d0 3.0d0)The svd function calculates the singular value decomposition
of a given matrix, returning an object with slots for the three returned
matrices:
* (defparameter a #2A((1 2 3) (0 2 1) (1 3 2)))
A
* (defparameter a-svd (lla:svd a))
A-SVD
* a-svd
#S(LLA:SVD
:U #2A((-0.6494608633564334d0 0.7205486773948702d0 0.24292013188045855d0)
(-0.3744175632000917d0 -0.5810891192666799d0 0.7225973455785591d0)
(-0.6618248071322363d0 -0.3783451320875919d0 -0.6471807210432038d0))
:D #S(CL-NUM-UTILS.MATRIX:DIAGONAL-MATRIX
:ELEMENTS #(5.593122609997059d0 1.2364443401235103d0
0.43380279311714376d0))
:VT #2A((-0.2344460799312531d0 -0.7211054639318696d0 -0.6519524104506949d0)
(0.2767642134809678d0 -0.6924017945853318d0 0.6663192365460215d0)
(-0.9318994611765425d0 -0.02422116311440764d0 0.3619070730398283d0)))The diagonal matrix (singular values) and vectors can be accessed with functions:
(lla:svd-u a-svd)
#2A((-0.6494608633564334d0 0.7205486773948702d0 0.24292013188045855d0)
(-0.3744175632000917d0 -0.5810891192666799d0 0.7225973455785591d0)
(-0.6618248071322363d0 -0.3783451320875919d0 -0.6471807210432038d0))
* (lla:svd-d a-svd)
#S(CL-NUM-UTILS.MATRIX:DIAGONAL-MATRIX
:ELEMENTS #(5.593122609997059d0 1.2364443401235103d0 0.43380279311714376d0))
* (lla:svd-vt a-svd)
#2A((-0.2344460799312531d0 -0.7211054639318696d0 -0.6519524104506949d0)
(0.2767642134809678d0 -0.6924017945853318d0 0.6663192365460215d0)
(-0.9318994611765425d0 -0.02422116311440764d0 0.3619070730398283d0))The Matlisp scientific computation library provides high performance operations on arrays, including wrappers around BLAS and LAPACK functions. It can be loaded using quicklisp:
* (ql:quickload :matlisp)The nickname for matlisp is m. To avoid typing matlisp: or
m: in front of each symbol, you can define your own package which
uses matlisp
(See the PCL section on packages):
* (defpackage :my-new-code
(:use :common-lisp :matlisp))
#<PACKAGE "MY-NEW-CODE">
* (in-package :my-new-code)and to use the #i infix reader (note the same name as for
cmu-infix), run:
* (named-readtables:in-readtable :infix-dispatch-table)* (matlisp:zeros '(2 2))
#<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: DOUBLE-FLOAT>| #(2 2)
0.000 0.000
0.000 0.000
>Note that by default matrix storage types are double-float.
To create a complex array using zeros, ones and eye, specify the type:
* (matlisp:zeros '(2 2) '((complex double-float)))
#<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: (COMPLEX DOUBLE-FLOAT)>| #(2 2)
0.000 0.000
0.000 0.000
>As well as zeros and ones there is eye which creates an identity
matrix:
* (matlisp:eye '(3 3) '((complex double-float)))
#<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: (COMPLEX DOUBLE-FLOAT)>| #(3 3)
1.000 0.000 0.000
0.000 1.000 0.000
0.000 0.000 1.000
>To generate 1D arrays there are the range and linspace functions:
* (matlisp:range 1 10)
#<|<SIMPLE-DENSE-TENSOR: (INTEGER 0 4611686018427387903)>| #(9)
1 2 3 4 5 6 7 8 9
>The range function rounds down it's final argument to an integer:
* (matlisp:range 1 -3.5)
#<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: SINGLE-FLOAT>| #(5)
1.000 0.000 -1.000 -2.000 -3.000
>
* (matlisp:range 1 3.3)
#<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: SINGLE-FLOAT>| #(3)
1.000 2.000 3.000
>Linspace is a bit more general, and the values returned include the end point.
* (matlisp:linspace 1 10)
#<|<SIMPLE-DENSE-TENSOR: (INTEGER 0 4611686018427387903)>| #(10)
1 2 3 4 5 6 7 8 9 10
>* (matlisp:linspace 0 (* 2 pi) 5)
#<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: DOUBLE-FLOAT>| #(5)
0.000 1.571 3.142 4.712 6.283
>Currently linspace requires real inputs, and doesn't work with complex numbers.
* (matlisp:random-uniform '(2 2))
#<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: DOUBLE-FLOAT>| #(2 2)
0.7287 0.9480
2.6703E-2 0.1834
>(matlisp:random-normal '(2 2))
#<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: DOUBLE-FLOAT>| #(2 2)
0.3536 -1.291
-0.3877 -1.371
>There are functions for other distributions, including
random-exponential, random-beta, random-gamma and
random-pareto.
The #d and #e reader macros provide a way to create double-float
and single-float tensors:
* #d[1,2,3]
#<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: DOUBLE-FLOAT>| #(3)
1.000 2.000 3.000
>
* #d[[1,2,3],[4,5,6]]
#<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: DOUBLE-FLOAT>| #(2 3)
1.000 2.000 3.000
4.000 5.000 6.000
>Note that the comma separators are needed.
Common lisp arrays can be converted to Matlisp tensors by copying:
* (copy #2A((1 2 3)
(4 5 6))
'#.(tensor 'double-float))
#<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: DOUBLE-FLOAT>| #(2 3)
1.000 2.000 3.000
4.000 5.000 6.000
>Instances of the tensor class can also be created, specifying the
dimensions. The internal storage of tensor objects is a 1D array
(simple-vector) in a slot store.
For example, to create a double-float type tensor:
(make-instance (tensor 'double-float)
:dimensions (coerce '(2) '(simple-array index-type (*)))
:store (make-array 2 :element-type 'double-float))The array store can be accessed using slots:
* (defparameter vec (m:range 0 5))
* vec
#<|<SIMPLE-DENSE-TENSOR: (INTEGER 0 4611686018427387903)>| #(5)
0 1 2 3 4
>
* (slot-value vec 'm:store)
#(0 1 2 3 4)Multidimensional tensors are also stored in 1D arrays, and are stored in column-major order rather than the row-major ordering used for common lisp arrays. A displaced array will therefore be transposed.
The contents of a tensor can be copied into an array
* (let ((tens (m:ones '(2 3))))
(m:copy tens 'array))
#2A((1.0d0 1.0d0 1.0d0) (1.0d0 1.0d0 1.0d0))or a list:
* (m:copy (m:ones '(2 3)) 'cons)
((1.0d0 1.0d0 1.0d0) (1.0d0 1.0d0 1.0d0))The ref function is the equivalent of aref for standard CL
arrays, and is also setf-able:
* (defparameter a (matlisp:ones '(2 3)))
* (setf (ref a 1 1) 2.0)
2.0d0
* a
#<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: DOUBLE-FLOAT>| #(2 3)
1.000 1.000 1.000
1.000 2.000 1.000
>The matlisp-user package, loaded when matlisp is loaded, contains
functions for operating element-wise on tensors.
* (matlisp-user:* 2 (ones '(2 3)))
#<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: DOUBLE-FLOAT>| #(2 3)
2.000 2.000 2.000
2.000 2.000 2.000
>This includes arithmetic operators '+', '-', '*', '/' and 'expt', but
also sqrt,sin,cos,tan, hyperbolic functions, and their inverses.
The #i reader macro recognises many of these, and uses the
matlisp-user functions:
* (let ((a (ones '(2 2)))
(b (random-normal '(2 2))))
#i( 2 * a + b ))
#<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: DOUBLE-FLOAT>| #(2 2)
0.9684 3.250
1.593 1.508
>
* (let ((a (ones '(2 2)))
(b (random-normal '(2 2))))
(macroexpand-1 '#i( 2 * a + b )))
(MATLISP-USER:+ (MATLISP-USER:* 2 A) B)