# The Common Lisp Cookbook – Multidimensional arrays

📹 Discover Alexander's videos on Common Lisp development (also on the Fediverse).

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 maintained by @Symbolics defines functions `generate`, `permute`, `displace`, `flatten`, `split`, `combine`, `reshape`. It also defines `each`, for element-wise operations. This is a fork of bendudson/array-operations which is a fork of tpapp/array-operations, the original author.
• 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 when `a` and/or `b` are 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.git
``````

Instructions 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.

## Creating

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)
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 Array Operations manual on generate for more examples.

### Random numbers

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)
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.

## Accessing elements

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)
3
``````

The 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)
3
``````

To 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
NIL
``````

A 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
;; 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]

### Row major indexing

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))
``````

### Infix syntax

The cmu-infix library provides some different syntax which can make mathematical expressions easier to read:

``````* (ql:quickload :cmu-infix)
cmu-infix

(:CMU-INFIX)

...)

* (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.

## Element-wise operations

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.

### Vectorising expressions

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)
``````

### Calling BLAS

Several packages provide wrappers around BLAS, for fast matrix manipulation.

The lla package in quicklisp includes calls to some functions:

#### Scale an array

scaling by a constant factor:

``````* (defparameter a #(1 2 3))
* (lla:scal! 2.0 a)
* a
#(2.0d0 4.0d0 6.0d0)
``````

#### AXPY

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))
``````

#### Dot product

The dot product of two vectors:

``````* (defparameter x #(1 2 3))
* (defparameter y #(2 3 4))
* (lla:dot x y)
20.0d0
``````

### Reductions

The 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))
4
``````

The `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))
4
``````

More 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)))
2
``````

This 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)))
``````

## Linear algebra

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)

lla
.
(:LLA)
``````

### Matrix multiplication

The lla function `mm` performs vector-vector, matrix-vector and matrix-matrix multiplication.

#### Vector dot product

Note that one vector is treated as a row vector, and the other as column:

``````* (lla:mm #(1 2 3) #(2 3 4))
20
``````

#### Matrix-vector product

``````* (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]`

#### Matrix-matrix multiply

``````* (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))
``````

#### Outer product

The array-operations package contains a generalised outer product function:

``````* (ql:quickload :array-operations)
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.

### Matrix inverse

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)
``````

### Singular value decomposition

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))
``````

## Matlisp

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)
``````

### Creating tensors

``````* (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
>
``````

#### Ranges

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.

#### Random 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.

#### Tensors from arrays

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))
``````

#### Arrays from tensors

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))
``````

### Element access

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
>
``````

### Element-wise operations

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)
``````

Page source: arrays.md

T
O
C