Linear Algebra

This chapter describes the real matrix, complex cmatrix and integer imatrix objects as supported by MAD-NG. The module for Vector and Matrix is not exposed, only the contructors are visible from the MAD environment and thus, matrices are handled directly by their methods or by the generic functions of the same name from the module MAD.gmath. The imatrix, i.e. matrix of integers, are mainly used for indexing other types of matrix and therefore supports only a limited subset of the features. Column and row vectors are shortcuts for \([n\times 1]\) and \([1\times n]\) matrices respectively. Note that matrix, cmatrix and imatrix are all defined as C structures containing their elements in row-major order for direct compliance with the C API.

Types promotion

The matrix operations may involve other data types like real and complex numbers leading to many combinations of types. In order to simplify the descriptions, the generic names num, cpx and idx (indexes) are used for real, complex and integer numbers respectively, vec, cvec and ivec for real, complex and integer vectors respectively, and mat, cmat and imat for real, complex and integer matrices respectively. For example, the sum of a complex number cpx and a real matrix mat gives a complex matrix cmat. The case of idx means that a number will be interpreted as an index and automatically rounded if it does not hold an integer value. The following table summarizes all valid combinations of types for binary operations involving at least one matrix type:

Left Operand Type	Right Operand Type	Result Type
number	imatrix	imatrix
imatrix	number	imatrix
imatrix	imatrix	imatrix
number	matrix	matrix
matrix	number	matrix
matrix	matrix	matrix
number	cmatrix	cmatrix
complex	matrix	cmatrix
complex	cmatrix	cmatrix
matrix	complex	cmatrix
matrix	cmatrix	cmatrix
cmatrix	number	cmatrix
cmatrix	complex	cmatrix
cmatrix	matrix	cmatrix
cmatrix	cmatrix	cmatrix

Constructors

The constructors for vectors and matrices are directly available from the MAD environment. Note that real, complex or integer matrix with zero size are not allowed, i.e. the smallest allowed matrix has sizes of \([1\times 1]\).

vector(nrow)

cvector(nrow)

ivector(nrow)

Return a real, complex or integer column vector (i.e. a matrix of size \([n_{\text{row}}\times 1]\)) filled with zeros. If nrow is a table, it is equivalent to vector(#nrow):fill(nrow).

matrix(nrow, ncol_)

cmatrix(nrow, ncol_)

imatrix(nrow, ncol_)

Return a real, complex or integer matrix of size \([n_{\text{row}}\times n_{\text{col}}]\) filled with zeros. If nrow is a table, it is equivalent to matrix(#nrow, #nrow[1] or 1):fill(nrow), and ignoring ncol. Default: ncol_ = rnow.

linspace([start_, ] stop, size_)

Return a real or complex column vector of length size filled with values equally spaced between start and stop on a linear scale. Note that numerical range can generate the same real sequence of values in a more compact form. Default: start_ = 0, size_ = 100.

logspace([start_, ] stop, size_)

Return a real or complex column vector of length size filled with values equally spaced between start and stop on a logarithmic scale. Note that numerical logrange can generate the same real sequence of values in a more compact form. Default: start_ = 1, size_ = 100.

Attributes

mat.nrow

The number of rows of the real, complex or integer matrix mat.

mat.ncol

The number of columns of the real, complex or integer matrix mat.

Functions

is_vector(a)

is_cvector(a)

is_ivector(a)

Return true if a is respectively a real, complex or integer matrix of size \([n_{\text{row}}\times 1]\) or \([1\times n_{\text{row}}]\), false otherwise. These functions are only available from the module MAD.typeid.

isa_vector(a)

Return true if a is a real or complex vector (i.e. is-a vector), false otherwise. This function is only available from the module MAD.typeid.

isy_vector(a)

Return true if a is a real, complex or integer vector (i.e. is-any vector), false otherwise. This function is only available from the module MAD.typeid.

is_matrix(a)

is_cmatrix(a)

is_imatrix(a)

Return true if a is respectively a real, complex or integer matrix, false otherwise. These functions are only available from the module MAD.typeid.

isa_matrix(a)

Return true if a is a real or complex matrix (i.e. is-a matrix), false otherwise. This function is only available from the module MAD.typeid.

isy_matrix(a)

Return true if a is a real, complex or integer matrix (i.e. is-any matrix), false otherwise. This function is only available from the module MAD.typeid.

Methods

Special Constructors

mat:vec()

Return a vector of the same type as mat filled with the values of the elements of the vectorized real, complex or integer matrix mat equivalent to mat:t():reshape(#mat,1).

mat:vech()

Return a vector of the same type as mat filled with the values of the elements of the half vectorized real, complex or integer symmetric matrix mat. The symmetric property can be pre-checked by the user with mat:is_symm().

mat:diag(k_)

If mat is a matrix, return a column vector containing its \(k\)-th diagonal equivalent to mat:getdiag(k). If mat is a vector, return a square matrix with its \(k\)-th diagonal set to the values of the elements of mat equivalent to mat:same(n,n):setdiag(mat,k) where n = #mat+abs(k). Default: k_ = 0.

Sizes and Indexing

mat:size()

Return the number of elements nrow * ncol of the real, complex or integer matrix mat equivalent to #mat.

mat:bytesize()

Return the number of bytes used by the data storage of the real, complex or integer matrix mat equivalent to #mat * sizeof(mat[1]).

mat:sizes()

Return the number of rows nrow and columns ncol of the real, complex or integer matrix mat. Note that #mat returns the full size nrow * ncol of the matrix.

mat:tsizes()

Return the number of columns ncol and rows nrow (i.e. transposed sizes) of the real, complex or integer matrix mat equivalent to swap(mat:sizes()).

mat:getij(ij_, ri_, rj_)

Return two ivector or ri and rj containing the indexes (i,j) extracted from the iterable ij for the real, complex or integer matrix mat. If ij is a number, the two returned items are also numbers. This method is the reverse method of mat:getidx() to convert 1D indexes into 2D indexes for the given matrix sizes. Default: ij_ = 1..#mat.

mat:getidx(ir_, jc_, rij_)

Return an ivector or rij containing #ir * #jc vector indexes in row-major order given by the iterable ir and jc of the real, complex or integer matrix mat, followed by ir and jc potentially set from defaults. If both ir and jc are numbers, it returns a single number. This method is the reverse method of mat:getij() to convert 2D indexes into 1D indexes for the given matrix sizes. Default: ir_ = 1..nrow, jc_ = 1..ncol.

mat:getdidx(k_)

Return an iterable describing the indexes of the \(k\)-th diagonal of the real, complex or integer matrix mat where -nrow <= k <= ncol. This method is useful to build the 1D indexes of the \(k\)-th diagonal for the given matrix sizes. Default k_ = 0

Getters and Setters

mat:get(i, j)

Return the value of the element at the indexes (i,j) of the real, complex or integer matrix mat for 1 <= i <= nrow and 1 <= j <= ncol, nil otherwise.

mat:set(i, j, v)

Assign the value v to the element at the indexes (i,j) of the real, complex or integer matrix mat for 1 <= i <= nrow and 1 <= j <= ncol and return the matrix, otherwise raise an “index out of bounds” error.

mat:geti(n)

Return the value of the element at the vector index n of the real, complex or integer matrix mat for 1 <= n <= #mat, i.e. interpreting the matrix as a vector, nil otherwise.

mat:seti(n, v)

Assign the value v to the element at the vector index n of the real, complex or integer matrix mat for 1 <= n <= #mat and return the matrix, i.e. interpreting the matrix as a vector, otherwise raise an “index out of bounds” error.

mat:getvec(ij_, r_)

Return a column vector or r containing the values of the elements at the vector indexes given by the iterable ij of the real, complex or integer matrix mat, i.e. interpreting the matrix as a vector. Default: ij_ = 1..#mat.

mat:setvec(ij_, a, p_, s_)

Return the real, complex or integer matrix mat after filling the elements at the vector indexes given by the iterable ij, i.e. interpreting the matrix as a vector, with the values given by a depending of its kind:

• if a is a scalar, it is will be used repetitively.

• if a is an iterable then the matrix will be filled with values from a[n] for 1 <= n <= #a and recycled repetitively if #a < #ij.

• if a is a callable, then a is considered as a stateless iterator, and the matrix will be filled with the values v returned by iterating s, v = a(p, s).

Default: ij_ = 1..#mat.

mat:insvec(ij, a)

Return the real, complex or integer matrix mat after inserting the elements at the vector indexes given by the iterable ij, i.e. interpreting the matrix as a vector, with the values given by a depending of its kind:

• if a is a scalar, it is will be used repetitively.

• if a is an iterable then the matrix will be filled with values from a[n] for 1 <= n <= #a.

The elements after the inserted indexes are shifted toward the end of the matrix in row-major order and discarded if they go beyond the last index.

mat:remvec(ij)

Return the real, complex or integer matrix mat after removing the elements at the vector indexes given by the iterable ij, i.e. interpreting the matrix as a shrinking vector, and reshaped as a column vector of size #mat - #ij.

mat:swpvec(ij, ij2)

Return the real, complex or integer matrix mat after swapping the values of the elements at the vector indexes given by the iterable ij and ij2, i.e. interpreting the matrix as a vector.

mat:getsub(ir_, jc_, r_)

Return a \([\) #ir \(\times\) #jc \(]\) matrix or r containing the values of the elements at the indexes given by the iterable ir and jc of the real, complex or integer matrix mat. If ir = nil, jc ~= nil and r is a 1D iterable, then the latter is filled using column-major indexes. Default: as mat:getidx().

mat:setsub(ir_, jc_, a, p_, s_)

Return the real, complex or integer matrix mat after filling the elements at the indexes given by the iterable ir and jc with the values given by a depending of its kind:

• if a is a scalar, it is will be used repetitively.

• if a is an iterable then the rows and columns will be filled with values from a[n] for 1 <= n <= #a and recycled repetitively if #a < #ir * #ic.

• if a is a callable, then a is considered as a stateless iterator, and the columns will be filled with the values v returned by iterating s, v = a(p, s).

If ir = nil, jc ~= nil and a is a 1D iterable, then the latter is used to filled the matrix in the column-major order. Default: as mat:getidx().

mat:inssub(ir_, jc_, a)

Return the real, complex or integer matrix mat after inserting elements at the indexes (i,j) given by the iterable ir and jc with the values given by a depending of its kind:

• if a is a scalar, it is will be used repetitively.

• if a is an iterable then the rows and columns will be filled with values from a[n] for 1 <= n <= #a and recycled repetitively if #a < #ir * #ic.

The values after the inserted indexes are pushed toward the end of the matrix, i.e. interpreting the matrix as a vector, and discarded if they go beyond the last index. If ir = nil, jc ~= nil and a is a 1D iterable, then the latter is used to filled the matrix in the column-major order. Default: as mat:getidx().

mat:remsub(ir_, jc_)

Return the real, complex or integer matrix mat after removing the rows and columns at the indexes given by the iterable ir and jc and reshaping the matrix accordingly. Default: as mat:getidx().

mat:swpsub(ir_, jc_, ir2_, jc2_)

Return the real, complex or integer matrix mat after swapping the elements at indexes given by the iterable iterable ir and jc with the elements at indexes given by iterable ir2 and jc2. Default: as mat:getidx().

mat:getrow(ir, r_)

Equivalent to mat:getsub() with jc = nil.

mat:setrow(ir, a, p_, s_)

Equivalent to mat:setsub() with jc = nil.

mat:insrow(ir, a)

Equivalent to mat:inssub() with jc = nil.

mat:remrow(ir)

Equivalent to mat:remsub() with jc = nil.

mat:swprow(ir, ir2)

Equivalent to mat:swpsub() with jc = nil and jc2 = nil.

mat:getcol(jc, r_)

Equivalent to mat:getsub() with ir = nil.

mat:setcol(jc, a, p_, s_)

Equivalent to mat:setsub() with ir = nil.

mat:inscol(jc, a)

Equivalent to mat:inssub() with ir = nil. If a is a matrix with ncol > 1 then a = 0 and it is followed by mat:setsub() with ir = nil to obtain the expected result.

mat:remcol(jc)

Equivalent to mat:remsub() with ir = nil.

mat:swpcol(jc, jc2)

Equivalent to mat:swpsub() with ir = nil and ir2 = nil.

mat:getdiag([k_, ] r_)

Return a column vector of length min(nrow, ncol)-abs(k) or r containing the values of the elements on the \(k\)-th diagonal of the real, complex or integer matrix mat using mat:getvec(). Default: as mat:getdidx().

mat:setdiag(a, [k_, ] p_, s_)

Return the real, complex or integer matrix mat after filling the elements on its \(k\)-th diagonal with the values given by a using mat:setvec(). Default: as mat:getdidx().

Cloning and Reshaping

mat:copy(r_)

Return a matrix or r filled with a copy of the real, complex or integer matrix mat.

mat:same([nr_, nc_, ] v_)

Return a matrix with elements of the type of v and with nr rows and nc columns. Default: v_ = mat[1], nr_ = nrow, nc_ = ncol.

mat:reshape(nr_, nc_)

Return the real, complex or integer matrix mat reshaped to the new sizes nr and nc that must result into an equal or smaller number of elements, or it will raise an invalid new sizes error. Default: nr_ = #mat, nc_ = 1.

mat:_reshapeto(nr, nc_)

Same as mat:reshape() except that nr must be explicitly provided as this method allows for a larger size than mat current size. A typical use of this method is to expand a vector after an explicit shrinkage, while keeping track of its original size, e.g. similar to vector(100) :reshape(1):seti(1,1) :_reshapeto(2):seti(2,1) that would raise an “index out of bounds” error without the call to _reshapeto(). Default nc_ = 1.

WARNING: This method is unsafe and may crash MAD-NG with, e.g. a Segmentation fault , if wrongly used. It is the responsibility of the user to ensure that mat was created with a size greater than or equal to the new size.

Matrix Properties

mat:is_const(tol_)

Return true if all elements are equal within the tolerance tol, false otherwise. Default: tol_ = 0.

mat:is_real(tol_)

Return true if the imaginary part of all elements are equal to zero within the tolerance tol, false otherwise. Default: tol_ = 0.

mat:is_imag(tol_)

Return true if the real part of all elements are equal to zero within the tolerance tol, false otherwise. Default: tol_ = 0.

mat:is_diag(tol_)

Return true if all elements off the diagonal are equal to zero within the tolerance tol, false otherwise. Default: tol_ = 0.

mat:is_symm([tol_, ] [sk_, ] c_)

Return true if mat is a symmetric matrix, i.e. \(M = M^*\) within the tolerance tol, false otherwise. It checks for a skew-symmetric matrix if sk = true, and for a Hermitian matrix if c ~= false, and a skew-Hermitian matrix if both are combined. Default: tol_ = 0.

mat:is_symp(tol_)

Return true if mat is a symplectic matrix, i.e. \(M^* S_{2n} M = S_{2n}\) within the tolerance tol, false otherwise. Default: tol_ = eps.

Filling and Moving

mat:zeros()

Return the real, complex or integer matrix mat filled with zeros.

mat:ones(v_)

Return the real, complex or integer matrix mat filled with the value of v. Default: v_ = 1.

mat:eye(v_)

Return the real, complex or integer matrix mat filled with the value of v on the diagonal and zeros elsewhere. The name of this method comes from the spelling of the Identity matrix \(I\). Default: v_ = 1.

mat:seq([v_, ] d_)

Return the real, complex or integer matrix mat filled with the indexes of the elements (i.e. starting at 1) and shifted by the value of v. The matrix is filled in the row-major order unless d = 'col'. Default: v_ = 0.

mat:random(f_, ...)

Return the real, complex or integer matrix mat filled with random values generated by f(...), and called twice for each element for a cmatrix. Default: f_ = math.random.

mat:shuffle()

Return the real, complex or integer matrix mat with its elements randomly swapped using the Fisher–Yates or Knuth shuffle algorithm and math.random as the PRNG.

mat:symp()

Return the real or complex matrix mat filled with the block diagonal unitary Symplectic matrix sometimes named \(J_{2n}\) or \(S_{2n}\). The matrix mat must be square with even number of rows and columns otherwise a “2n square matrix expected” error is raised.

mat:circ(v)

Return the real or complex matrix mat filled as a Circulant matrix using the values from the iterable v, and rotating elements for each row or column depending on the shape of v.

mat:fill(a, p_, s_)

Return the real, complex or integer matrix mat filled with values provided by a depending of its kind:

• if a is a scalar, it is equivalent to mat:ones(a).

• if a is a callable, then:

  • if p and s are provided, then a is considered as a stateless iterator, and the matrix will be filled with the values v returned by iterating s, v = a(p, s).

  • otherwise a is considered as a generator, and the matrix will be filled with values returned by calling a(mat:get(i,j), i, j).

• if a is an iterable then:

  • if a[1] is also an iterable, the matrix will be filled with the values from a[i][j] for 1 <= i <= nrow and 1 <= j <= ncol, i.e. treated as a 2D container.

  • otherwise the matrix will be filled with values from a[n] for 1 <= n <= #mat, i.e. treated as a 1D container.

mat:rev(d_)

Reverse the elements of the matrix mat according to the direction d:

• If d = 'vec', it reverses the entire matrix.

• If d = 'row', it reverses each row.

• If d = 'col', it reverses each column.

• If d = 'diag', it reverse the only the diagonal.

Default: d_ = 'vec'.

mat:roll(nr_, nc_)

Return the real, complex or integer matrix mat after rolling its rows by nr \(\in \mathbb{Z}\) and then columns by nc \(\in \mathbb{Z}\). Default: nr_ = 0, nc_ = 0.

mat:movev(i, j, k, r_)

Return the real, complex or integer matrix r after moving the elements in mat[i..j] to r[k..k+j-i] with 1 <= i <= j <= #mat and 1 <= k <= k+j-i <= #r. Default: r_ = mat.

mat:shiftv(i, n_)

Return the real, complex or integer matrix mat after shifting the elements in mat[i..] to mat[i+n..] if n > 0 and in the opposite direction if n < 0, i.e. it is equivalent to mat:movev(i, #mat-n, i+n) for n > 0 and to mat:movev(i, #mat, i+n) for n < 0. Default: n_ = 1.

Mapping and Folding

This section lists the high-order functions map, fold and their variants useful in functional programming [1], followed by sections that list their direct application.

mat:foreach([ij_, ] f)

Return the real, complex or integer matrix mat after applying the callable f to the elements at the indexes given by the iterable ij using f(mat[n], n), i.e. interpreting the matrix as a vector. Default: ij_ = 1..#mat.

mat:filter([ij_, ] p, r_)

Return a matrix or r filled with the values of the elements of the real, complex or integer matrix mat at the indexes given by the iterable ij if they are selected by the callable predicate p using p(mat[n], n) = true, i.e. interpreting the matrix as a vector. This method returns next to the matrix, a table if r is a table or a ivector otherwise, containing the indexes of the selected elements returned. Default: ij_ = 1..#mat.

mat:filter_out([ij_, ] p, r_)

Equivalent to map:filter(ij_, compose(lnot,p), r_), where the functions compose() and lnot() are provided by the module MAD.gfunc.

mat:map([ij_, ] f, r_)

Return a matrix or r filled with the values returned by the callable (or the operator string) f applied to the elements of the real, complex or integer matrix mat at the indexes given by the iterable ij using f(mat[n], n), i.e. interpreting the matrix as a vector. If r = 'in' or r = nil and ij ~= nil then it is assigned mat, i.e. map in place. If r = nil still, then the type of the returned matrix is determined by the type of the value returned by f() called once before mapping. Default: ij_ = 1..#mat.

mat:map2(y, [ij_, ] f, r_)

Equivalent to mat:map() but with two arguments passed to f, i.e. using f(mat[n], y[n], n).

mat:map3(y, z, [ij_, ] f, r_)

Equivalent to mat:map() but with three arguments passed to f, i.e. using f(mat[n], y[n], z[n], n). Note that f cannot be an operator string, as only unary and binary operators are avalaible in such form.

mat:foldl(f, [x0_, ] [d_, ] r_)

Return a scalar, a vector or r filled with the values returned by the callable (or the operator string) f applied iteratively to the elements of the real, complex or integer matrix mat using the folding left (forward with increasing indexes) expression v = f(v, mat[n]) starting at x0 and running in the direction depending on the string d:

• If d = 'vec', the folding left iteration runs on the entire matrix mat interpreted as a vector and a scalar is returned.

• If d = 'row', the folding left iteration runs on the rows of the matrix mat and a column vector is returned.

• If d = 'col', the folding left iteration runs on the columns of the matrix mat and a row vector is returned.

• If d = 'diag', the folding left iteration runs on the diagonal of the matrix mat and a scalar is returned.

Note that ommitting both x0 and d implies to not specify r as well, otherwise the latter will be interpreted as x0. If r = nil and d = 'row' or d = 'col', then the type of the returned vector is determined by the type of the value returned by f() called once before folding. Default: x0 = mat[1] (or first row or column element), d = 'vec'.

mat:foldr(f, [x0_, ] [d_, ] r_)

Same as mat:foldl() but the callable (or the operator string) f is applied iteratively using the folding right (backward with decreasing indexes) expression v = f(mat[n], v). Default: x0 = mat[#mat] (or last row or column element), d = 'vec'.

mat:scanl(f, [x0_, ] [d_, ] r_)

Return a vector, a matrix or r filled with the values returned by the callable (or the operator string) f applied iteratively to the elements of the real, complex or integer matrix mat using the scanning left (forward with increasing indexes) expression v = f(v, mat[n]) starting at x0 and running in the direction depending on the string d:

• If d = 'vec', the scanning left iteration runs on the entire matrix mat interpreted as a vector and a vector is returned.

• If d = 'row', the scanning left iteration runs on the rows of the matrix mat and a matrix is returned.

• If d = 'col', the scanning left iteration runs on the columns of the matrix mat and a matrix is returned.

• If d = 'diag', the scanning left iteration runs on the diagonal of the matrix mat and a vector is returned.

Note that ommitting both x0 and d implies to not specify r as well, otherwise the latter will be interpreted as x0. If r = nil, then the type of the returned matrix is determined by the type of the value returned by f() called once before scanning. Default: x0 = mat[1] (or first row or column element), d = 'vec'.

mat:scanr(f, [x0_, ] [d_, ] r_)

Same as mat:scanl() but the callable (or the operator string) f is applied iteratively using the scanning right (backward with decreasing indexes) expression v = f(mat[n], v). Default: x0 = mat[#mat] (or last row or column element), d = 'vec'.

Mapping Real-like Methods

The following table lists the methods built from the application of mat:map() and variants to the real-like functions from the module MAD.gmath for matrix and cmatrix. The methods mat:sign(), mat:sign1() and mat:atan2() are not available for cmatrix, and only the methods mat:abs(), mat:sqr() and mat:sign() are available for imatrix.

Functions	Equivalent Mapping
mat:abs(r_)	mat:map(abs,r_)
mat:acos(r_)	mat:map(acos,r_)
mat:acosh(r_)	mat:map(acosh,r_)
mat:acot(r_)	mat:map(acot,r_)
mat:acoth(r_)	mat:map(acoth,r_)
mat:asin(r_)	mat:map(asin,r_)
mat:asinh(r_)	mat:map(asinh,r_)
mat:asinc(r_)	mat:map(asinc,r_)
mat:asinhc(r_)	mat:map(asinhc,r_)
mat:atan(r_)	mat:map(atan,r_)
mat:atan2(y,r_)	mat:map2(y,atan2,r_)
mat:atanh(r_)	mat:map(atanh,r_)
mat:ceil(r_)	mat:map(ceil,r_)
mat:cos(r_)	mat:map(cos,r_)
mat:cosh(r_)	mat:map(cosh,r_)
mat:cot(r_)	mat:map(cot,r_)
mat:coth(r_)	mat:map(coth,r_)
mat:exp(r_)	mat:map(exp,r_)
mat:floor(r_)	mat:map(floor,r_)
mat:frac(r_)	mat:map(frac,r_)
mat:hypot(y,r_)	mat:map2(y,hypot,r_)
mat:hypot3(y,z,r_)	mat:map3(y,z,hypot3,r_)
mat:invsqrt([v_,]r_)	mat:map2(v_ or 1,invsqrt,r_)
mat:log(r_)	mat:map(log,r_)
mat:log10(r_)	mat:map(log10,r_)
mat:round(r_)	mat:map(round,r_)
mat:sign(r_)	mat:map(sign,r_)
mat:sign1(r_)	mat:map(sign1,r_)
mat:sin(r_)	mat:map(sin,r_)
mat:sinc(r_)	mat:map(sinc,r_)
mat:sinh(r_)	mat:map(sinh,r_)
mat:sinhc(r_)	mat:map(sinhc,r_)
mat:sqr(r_)	mat:map(sqr,r_)
mat:sqrt(r_)	mat:map(sqrt,r_)
mat:tan(r_)	mat:map(tan,r_)
mat:tanh(r_)	mat:map(tanh,r_)
mat:trunc(r_)	mat:map(trunc,r_)

Mapping Complex-like Methods

The following table lists the methods built from the application of mat:map() to the the complex-like functions from the module MAD.gmath for matrix and cmatrix.

Functions	Equivalent Mapping
mat:cabs(r_)	mat:map(cabs,r_)
mat:carg(r_)	mat:map(carg,r_)
mat:conj(r_)	mat:map(conj,r_)
mat:cplx(im_,r_)	mat:map2(im_, cplx, r_)
mat:fabs(r_)	mat:map(fabs,r_)
mat:imag(r_)	mat:map(imag,r_)
mat:polar(r_)	mat:map(polar,r_)
mat:proj(r_)	mat:map(proj,r_)
mat:real(r_)	mat:map(real,r_)
mat:rect(r_)	mat:map(rect,r_)
mat:reim(re_, im_)	mat:real(re_), mat:imag(im_)

The method mat:cplx() has a special implementation that allows to used it without a real part, e.g. im.cplx(nil, im, r_).

The method mat:conjugate() is also available as an alias for mat:conj().

Mapping Error-like Methods

The following table lists the methods built from the application of mat:map() to the error-like functions from the module MAD.gmath for matrix and cmatrix.

Functions	Equivalent Mapping
mat:erf([rtol_,]r_)	mat:map2(rtol_,erf,r_)
mat:erfc([rtol_,]r_)	mat:map2(rtol_,erfc,r_)
mat:erfcx([rtol_,]r_)	mat:map2(rtol_,erfcx,r_)
mat:erfi([rtol_,]r_)	mat:map2(rtol_,erfi,r_)
mat:wf([rtol_,]r_)	mat:map2(rtol_,wf,r_)

Mapping Vector-like Methods

The following table lists the methods built from the application of mat:map2() to the vector-like functions from the module MAD.gfunc for matrix, cmatrix, and imatrix.

Functions	Equivalent Mapping
mat:emul(mat2,r_)	mat:map2(mat2,mul,r_)
mat:ediv(mat2,r_)	mat:map2(mat2,div,r_)
mat:emod(mat2,r_)	mat:map2(mat2,mod,r_)
mat:epow(mat2,r_)	mat:map2(mat2,pow,r_)

Folding Methods

The following table lists the methods built from the application of mat:foldl() to the functions from the module MAD.gmath for matrix, cmatrix, and imatrix. The methods mat:min() and mat:max() are not available for cmatrix.

Functions	Equivalent Folding
mat:all(p,d_,r_)	mat:foldl(all(p),false,d_,r_)
mat:any(p,d_,r_)	mat:foldl(any(p),true,d_,r_)
mat:min(d_,r_)	mat:foldl(min,nil,d_,r_)
mat:max(d_,r_)	mat:foldl(max,nil,d_,r_)
mat:sum(d_,r_)	mat:foldl(add,nil,d_,r_)
mat:prod(d_,r_)	mat:foldl(mul,nil,d_,r_)
mat:sumsqr(d_,r_)	mat:foldl(sumsqrl,0,d_,r_)
mat:sumabs(d_,r_)	mat:foldl(sumabsl,0,d_,r_)
mat:minabs(d_,r_)	mat:foldl(minabsl,inf,d_,r_)
mat:maxabs(d_,r_)	mat:foldl(maxabsl,0,d_,r_)

Where any() and all() are functions that bind the predicate p to the propagation of the logical AND and the logical OR respectively, that can be implemented like:

• all = \p -> \r,x -> lbool(land(r, p(x)))

• any = \p -> \r,x -> lbool(lor (r, p(x)))

Scanning Methods

The following table lists the methods built from the application of mat:scanl() and mat:scanr() to the functions from the module MAD.gmath for matrix and cmatrix. The methods mat:accmin(), mat:raccmin(), mat:accmax() and mat:raccmax() are not available for cmatrix.

Functions	Equivalent Scanning
mat:accmin(d_,r_)	mat:scanl(min,nil,d_,r_)
mat:accmax(d_,r_)	mat:scanl(max,nil,d_,r_)
mat:accsum(d_,r_)	mat:scanl(add,nil,d_,r_)
mat:accprod(d_,r_)	mat:scanl(mul,nil,d_,r_)
mat:accsumsqr(d_,r_)	mat:scanl(sumsqrl,0,d_,r_)
mat:accsumabs(d_,r_)	mat:scanl(sumabsl,0,d_,r_)
mat:accminabs(d_,r_)	mat:scanl(minabsl,inf,d_,r_)
mat:accmaxabs(d_,r_)	mat:scanl(maxabsl,0,d_,r_)
mat:raccmin(d_,r_)	mat:scanr(min,nil,d_,r_)
mat:raccmax(d_,r_)	mat:scanr(max,nil,d_,r_)
mat:raccsum(d_,r_)	mat:scanr(add,nil,d_,r_)
mat:raccprod(d_,r_)	mat:scanr(mul,nil,d_,r_)
mat:raccsumsqr(d_,r_)	mat:scanr(sumsqrr,0,d_,r_)
mat:raccsumabs(d_,r_)	mat:scanr(sumabsr,0,d_,r_)
mat:raccminabs(d_,r_)	mat:scanr(minabsr,inf,d_,r_)
mat:raccmaxabs(d_,r_)	mat:scanr(maxabsr,0,d_,r_)

The method mat:accumulate() is also available as an alias for mat:accsum().

Matrix Functions

The following table lists the methods built from the application of mat:mfun() to the real-like functions from the module MAD.gmath for matrix and cmatrix.

Functions	Equivalent Matrix Function
mat:macos()	mat:mfun(acos)
mat:macosh()	mat:mfun(acosh)
mat:macot()	mat:mfun(acot)
mat:macoth()	mat:mfun(acoth)
mat:masin()	mat:mfun(asin)
mat:masinh()	mat:mfun(asinh)
mat:masinc()	mat:mfun(asinc)
mat:masinhc()	mat:mfun(asinhc)
mat:matan()	mat:mfun(atan)
mat:matanh()	mat:mfun(atanh)
mat:mcos()	mat:mfun(cos)
mat:mcosh()	mat:mfun(cosh)
mat:mcot()	mat:mfun(cot)
mat:mcoth()	mat:mfun(coth)
mat:mexp()	mat:mfun(exp)
mat:mlog()	mat:mfun(log)
mat:mlog10()	mat:mfun(log10)
mat:msin()	mat:mfun(sin)
mat:msinc()	mat:mfun(sinc)
mat:msinh()	mat:mfun(sinh)
mat:msinhc()	mat:mfun(sinhc)
mat:msqrt()	mat:mfun(sqrt)
mat:mtan()	mat:mfun(tan)
mat:mtanh()	mat:mfun(tanh)

Operator-like Methods

mat:unm(r_)

Equivalent to -mat with the possibility to place the result in r.

mat:add(a, r_)

Equivalent to mat + a with the possibility to place the result in r.

mat:sub(a, r_)

Equivalent to mat - a with the possibility to place the result in r.

mat:mul(a, r_)

Equivalent to mat * a with the possibility to place the result in r.

mat:tmul(mat2, r_)

Equivalent to mat:t() * mat2 with the possibility to place the result in r.

mat:mult(mat2, r_)

Equivalent to mat * mat2:t() with the possibility to place the result in r.

mat:dmul(mat2, r_)

Equivalent to mat:getdiag():diag() * mat2 with the possibility to place the result in r. If mat is a vector, it will be interpreted as the diagonal of a square matrix like in mat:diag(), i.e. omitting mat:getdiag() is the previous expression.

mat:muld(mat2, r_)

Equivalent to mat * mat2:getdiag():diag() with the possibility to place the result in r. If mat2 is a vector, it will be interpreted as the diagonal of a square matrix like in mat2:diag(), i.e. omitting mat2:getdiag() is the previous expression.

mat:div(a, [r_, ] rcond_)

Equivalent to mat / a with the possibility to place the result in r, and to specify the conditional number rcond used by the solver to determine the effective rank of non-square systems. Default: rcond = eps.

mat:inv([r_, ] rcond_)

Equivalent to mat.div(1, mat, r_, rcond_).

mat:pinv([r_, ] [rcond_, ] ncond_)

Return a real or complex matrix or r resulting from the pseudo-inverse of the matrix mat using its SVD factorisation returned by mat:svd(rcond_), followed by the rank of the system. The positive (resp. negative) integer ncond specifies the number of smallest (resp. largest) singular values to be discarded.

mat:pow(n, r_)

Equivalent to mat ^ n with the possibility to place the result in r.

mat:eq(a, tol_)

Equivalent to mat == a with the possibility to specify the tolerance tol of the comparison. Default: tol_ = 0.

mat:concat(mat2, [d_, ] r_)

Equivalent to mat .. mat2 with the possibility to place the result in r and to specify the direction of the concatenation:

• If d = 'vec', it concatenates the matrices (appended as vectors)

• If d = 'row', it concatenates the rows (horizontal)

• If d = 'col', it concatenates the columns (vectical)

Default: d_ = 'row'.

Special Methods

mat:transpose([c_, ] r_)

mat:t([c_, ] r_)

Return a real, complex or integer matrix or r resulting from the conjugate transpose \(M^*\) of the matrix mat unless c = false which disables the conjugate to get \(M^\tau\). If r = 'in' then it is assigned mat.

mat:trace()

mat:tr()

Return the Trace of the real or complex mat equivalent to mat:sum('diag').

mat:inner(mat2)

mat:dot(mat2)

Return the Inner Product of the two real or complex matrices mat and mat2 with compatible sizes, i.e. return \(M^* . M_2\) interpreting matrices as vectors. Note that multiple dot products, i.e. not interpreting matrices as vectors, can be achieved with mat:tmul().

mat:outer(mat2, r_)

Return the real or complex matrix resulting from the Outer Product of the two real or complex matrices mat and mat2, i.e. return \(M . M_2^*\) interpreting matrices as vectors.

mat:cross(mat2, r_)

Return the real or complex matrix resulting from the Cross Product of the two real or complex matrices mat and mat2 with compatible sizes, i.e. return \(M \times M_2\) interpreting matrices as a list of \([3 \times 1]\) column vectors.

mat:mixed(mat2, mat3, r_)

Return the real or complex matrix resulting from the Mixed Product of the three real or complex matrices mat, mat2 and mat3 with compatible sizes, i.e. return \(M^* . (M_2 \times M_3)\) interpreting matrices as a list of \([3 \times 1]\) column vectors.

mat:norm()

Return the Frobenius norm of the matrix \(\| M \|_2\). Other \(L_p\) matrix norms and variants can be easily calculated using already provided methods, e.g. \(L_1\) = mat:sumabs('col'):max(), \(L_{\infty}\) = mat:sumabs('row'):max(), and \(L_2\) = mat:svd():max().

mat:dist(mat2)

Equivalent to (mat - mat2):norm().

mat:unit()

Return the scaled matrix mat to the unit norm equivalent to mat:div(mat:norm(), mat).

mat:center(d_)

Return the centered matrix mat to have zero mean equivalent to mat:sub(mat:mean(),mat). The direction d indicates how the centering must be performed:

• If d = 'vec', it centers the entire matrix by substracting its mean.

• If d = 'row', it centers each row by substracting their mean.

• If d = 'col' , it centers each column by substracting their mean.

• If d = 'diag', it centers the diagonal by substracting its mean.

Default: d_ = 'vec'.

mat:angle(mat2, n_)

Return the angle between the two real or complex vectors mat and mat2 using the method mat:inner(). If n is provided, the sign of mat:mixed(mat2, n) is used to define the angle in \([-\pi,\pi]\), otherwise it is defined in \([0,\pi]\).

mat:minmax(abs_)

Return the minimum and maximum values of the elements of the real, complex or integer matrix mat. If abs = true, it returns the minimum and maximum absolute values of the elements. Default: abs_ = false.

mat:iminmax(abs_)

Return the two vector-like indexes of the minimum and maximum values of the elements of the real, complex or integer matrix mat. If abs = true, it returns the indexes of the minimum and maximum absolute values of the elements. Default: abs_ = false.

mat:mean()

Equivalent to mat:sum()/#mat, i.e. interpreting the matrix as a vector.

mat:variance()

Equivalent to (mat - mat:mean()):sumsqr()/(#mat-1), i.e. return the unbiased estimator of the variance with second order Bessel’s correction, interpreting the matrix as a vector.

mat:ksum()

mat:kdot(mat2)

Same as mat:sum() and mat:dot() respectively, except that they use the more accurate Kahan Babushka Neumaier algorithm for the summation, e.g. the sum of the elements of the vector \([1,10^{100},1,-10^{100}]\) should return \(0\) with sum() and the correct answer \(2\) with ksum().

mat:kadd(a, x)

Return the real or complex matrix mat filled with the linear combination of the compatible matrices stored in x times the scalars stored in a, i.e. mat = a[1]*x[1] + a[2]*x[2] ...

mat:eval(x0)

Return the evaluation of the real or complex matrix mat at the value x0, i.e. interpreting the matrix as a vector of polynomial coefficients of increasing orders in x evaluated at x = x0 using Horner’s method.

mat:sympconj(r_)

mat:bar(r_)

Return a real or complex matrix or r resulting from the symplectic conjugate of the matrix mat, with \(\bar{M} = -S_{2n} M^* S_{2n}\), and \(M^{-1} = \bar{M}\) if \(M\) is symplectic. If r = 'in' then it is assigned mat.

mat:symperr(r_)

Return the norm of the symplectic deviation matrix given by \(M^* S_{2n} M - S_{2n}\) of the real or complex matrix mat. If r is provided, it is filled with the symplectic deviation matrix.

mat:dif(mat2, r_)

Return a real or complex matrix or r resulting from the term-by-term difference between the matrices mat and mat2 using the absolute difference for values with magnitude below 1 and the relative difference otherwise, i.e. \(r_i = (x_i - y_i) / \max(|x_i|, 1)\).

Solvers and Decompositions

Except for nsolve(), the solvers hereafter are wrappers around the library Lapack [2].

mat:solve(b, rcond_)

Return the real or complex \([ n \times p ]\) matrix \(x\) as the minimum-norm solution of the linear least square problem \(\min \| A x - B \|\) where \(A\) is the real or complex \([ m \times n ]\) matrix mat and \(B\) is a \([ m \times p ]\) matrix b of the same type as mat, using LU, QR or LQ factorisation depending on the shape of the system. The conditional number rcond is used by the solver to determine the effective rank of non-square system. This method also returns the rank of the system. Default: rcond_ = eps.

mat:ssolve(b, rcond_, ncond_)

Return the real or complex \([ n \times p ]\) matrix \(x\) as the minimum \(L_2\)-norm solution of the linear least square problem \(\min \| A x - B \|\) where \(A\) is the real or complex \([ m \times n ]\) matrix mat and \(B\) is a \([ m \times p ]\) matrix b of the same type as mat, using SVD factorisation. The conditional number rcond is used by the solver to determine the effective rank of the system. If ncond is non-zero then mat:pinv() is used to solve the system. This method also returns the rank of the system followed by the real \([ \min(m,n) \times 1 ]\) vector of singluar values. Default: rcond_ = eps, ncond_ = 0.

mat:gsolve(b, c, d)

Return the real or complex \([ n \times 1 ]\) vector x as the minimum-norm solution of the linear least square problem \(\min \| A x - C \|\) under the constraint \(B x = D\) where \(A\) is the real or complex \([ m \times n ]\) matrix mat, \(B\) is a \([ p \times n ]\) matrix b, \(C\) is a \([ m \times 1 ]\) vector c and \(D\) is a \([ p \times 1 ]\) vector d, all of the same type as mat, using QR or LQ factorisation depending on the shape of the system. This method also returns the norm of the residues and the status info.

mat:gmsolve(b, d)

Return the real or complex \([ n \times 1 ]\) vector x and \([ p \times 1 ]\) matrix y as the minimum-norm solution of the linear Gauss-Markov problem \(\min_x \| y \|\) under the constraint \(A x + B y = D\) where \(A\) is the \([ m \times n ]\) real or complex matrix mat, \(B\) is a \([ m \times p ]\) matrix b, and \(D\) is a \([ m \times 1 ]\) vector d, both of the same type as mat, using QR or LQ factorisation depending on the shape of the system. This method also returns the status info.

mat:nsolve(b, nc_, tol_)

Return the real \([ n \times 1 ]\) vector x (of correctors kicks) as the minimum-norm solution of the linear (best-kick) least square problem \(\min \| A x - B \|\) where \(A\) is the real \([ m \times n ]\) (response) matrix mat and \(B\) is a real \([ m \times 1 ]\) vector b (of monitors readings), using the MICADO [3] algorithm based on the Householder-Golub method [MICADO]. The argument nc is the maximum number of correctors to use with \(0 < n_c \leq n\) and the argument tol is a convergence threshold (on the residues) to stop the (orbit) correction if \(\| A x - B \| \leq m \times\) tol. This method also returns the updated number of correctors \(n_c\) effectively used during the correction followed by the real \([ m \times 1 ]\) vector of residues. Default: nc_ = ncol, tol_ = eps.

mat:pcacnd(ns_, rcond_)

Return the integer column vector ic containing the indexes of the columns to remove from the real or complex \([ m \times n ]\) matrix mat using the Principal Component Analysis. The argument ns is the maximum number of singular values to consider and rcond is the conditioning number used to select the singular values versus the largest one, i.e. consider the ns larger singular values \(\sigma_i\) such that \(\sigma_i > \sigma_{\max}\times\)rcond. This method also returns the real \([ \min(m,n) \times 1 ]\) vector of singluar values. Default: ns_ = ncol, rcond_ = eps.

mat:svdcnd(ns_, rcond_, tol_)

Return the integer column vector ic containing the indexes of the columns to remove from the real or complex \([ m \times n ]\) matrix mat based on the analysis of the right matrix \(V\) from the SVD decomposition \(U S V\). The argument ns is the maximum number of singular values to consider and rcond is the conditioning number used to select the singular values versus the largest one, i.e. consider the ns larger singular values \(\sigma_i\) such that \(\sigma_i > \sigma_{\max}\times\)rcond. The argument tol is a threshold similar to rcond used to reject components in \(V\) that have similar or opposite effect than components already encountered. This method also returns the real \([ \min(m,n) \times 1 ]\) vector of singluar values. Default: ns_ = min(nrow,ncol), rcond_ = eps.

mat:svd()

Return the real vector of the singular values and the two real or complex matrices resulting from the SVD factorisation of the real or complex matrix mat, followed the status info. The singular values are positive and sorted in decreasing order of values, i.e. largest first.

mat:eigen(vr_, vl_)

Return the complex vector filled with the eigenvalues followed by the by the status info and the two optional real or complex matrices vr and vl containing the right and the transposed left eigenvectors resulting from the Eigen Decomposition of the real or complex square matrix mat. The eigenvectors are normalized to have unit Euclidean norm and their largest component real, and satisfy \(A v_r = \lambda v_r\) and \(v_l^\tau A = \lambda v_l^\tau\).

mat:det()

Return the Determinant of the real or complex square matrix mat using LU factorisation for better numerical stability, followed by the status info.

mat:mfun(fun)

Return the real or complex matrix resulting from the matrix function fun applied to the real or complex matrix mat. So far, mat:mfun() uses the eigen decomposition of the matrix mat, which must be Diagonalizable. See the section Matrix Functions for the list of matrix functions already provided. Future versions of this method may be extended to use the more general Schur-Parlett algorithm [MATFUN], and other specialized versions for msqrt(), mpow, mexp, and mlog may be implemented too.

Fourier Transforms and Convolutions

The methods described is this section are based on the FFTW and NFFT libraries.

mat:fft([d_, ] r_)

Return the complex \([n_r \times n_c]\) vector, matrix or r resulting from the 1D or 2D Fourier Transform of the real or complex \([n_r \times n_c]\) vector or matrix mat in the direction given by d:

• If d = 'vec', it returns a 1D vector FFT of length \(n_r n_c\).

• If d = 'row', it returns \(n_r\) 1D row FFTs of length \(n_c\).

• If d = 'col', it returns \(n_c\) 1D column FFTs of length \(n_r\).

• otherwise, it returns a 2D FFT of sizes \([n_r \times n_c]\).

mat:ifft([d_, ] r_)

Return the complex \([n_r \times n_c]\) vector, matrix or r resulting from the 1D or 2D inverse Fourier Transform of the complex \([n_r \times n_c]\) vector or matrix mat. See mat:fft() for the direction d.

mat:rfft([d_, ] r_)

Return the complex \([n_r \times \lfloor n_c/2+1\rfloor]\) vector, matrix or r resulting from the 1D or 2D Fourier Transform of the real \([n_r \times n_c]\) vector or matrix mat. This method used an optimized version of the FFT for real data, which is about twice as fast and compact as the method mat:fft(). See mat:fft() for the direction d.

mat:irfft([d_, ] r)

Return the real \([n_r \times n_c]\) vector, matrix or r resulting from the 1D or 2D inverse Fourier Transform of the complex \([n_r \times \lfloor n_c/2+1\rfloor]\) vector or matrix mat as computed by the method mat:rfft(). See mat:fft() for the direction d. Note that r must be provided to specify the correct \(n_c\) of the result.

mat:nfft(p_, r_)

Return the complex vector, matrix or r resulting from the 1D or 2D Nonequispaced Fourier Transform of the real or complex vector or matrix mat respectively at p time nodes.

mat:infft(p_, r_)

Return the complex vector, matrix or r resulting from the 1D or 2D Nonequispaced inverse Fourier Transform of the real or complex vector or matrix mat respectively at p frequency nodes.

mat:conv([y_, ] [d_], r_)

Return the real or complex vector, matrix or r resulting from the 1D or 2D Convolution between the compatible real or complex vectors or matrices mat and y respectively. See mat:fft() for the direction d. Default: y = mat.

mat:corr([y_, ] [d_], r_)

Return the real or complex vector, matrix or r resulting from the 1D or 2D Correlation between the compatible real or complex vectors or matrices mat and y respectively. See mat:fft() for the direction d. Default: y = mat.

mat:covar([y_, ] [d_, ] r_)

Return the real or complex vector, matrix or r resulting from the 1D or 2D Covariance between the compatible real or complex vectors or matrices mat and y respectively. See mat:fft() for the direction d. Default: y = mat.

mat:zpad(nr, nc, d_)

Return the real or complex vector or matrix resulting from the zero padding of the matrix mat extended to the sizes nr and nc, following the direction d:

• If d = 'vec', it pads the zeros at the end of the matrix equivalent x:same(nr,nc) :setvec(1..#x,x), i.e. interpreting the matrix as a vector.

• If d = 'row', it pads the zeros at the end of the rows equivalent x:same(x.nrow,nc) :setsub(1..x.nrow,1..x.ncol,x), i.e. ignoring nr.

• If d = 'col', it pads the zeros at the end of the columns equivalent x:same(nr,x.ncol) :setsub(1..x.nrow,1..x.ncol,x), i.e. ignoring nc.

• otherwise, it pads the zeros at the end of the rows and the columns equivalent to x:same(nr,nc) :setsub(1..x.nrow,1..x.ncol,x).

If the zero padding does not change the size of mat, the orignal mat is returned unchanged.

Rotations

This section describe methods dealing with 2D and 3D rotations (see Rotation Matrix) with angles in radians and trigonometric (counter-clockwise) direction for a right-handed frame, and where the following convention hold: ax = -phi is the elevation angle, ay =  theta is the azimuthal angle and az =  psi is the roll/tilt angle.

mat:rot(a)

Return the \([2\times 2]\) real matrix mat filled with a 2D rotation of angle a.

mat:rotx(a)

mat:roty(a)

mat:rotz(a)

Return the \([3\times 3]\) real matrix mat filled with a 3D rotation of angle a around the x-axis, y-axis and z-axis respectively.

mat:rotxy(ax, ay, inv_)

mat:rotxz(ax, az, inv_)

mat:rotyx(ay, ax, inv_)

mat:rotyz(ay, az, inv_)

mat:rotzx(az, ax, inv_)

mat:rotzy(az, ay, inv_)

Return the \([3\times 3]\) real matrix mat filled with a 3D rotation of the first angle argument ax, ay or az around the x-axis, y-axis or z-axis respectively followed by another 3D rotation of the second angle argument ax, ay or az around the x-axis, y-axis or z-axis respectively of the frame rotated by the first rotation. If inv is true, the returned matrix is the inverse rotation, i.e. the transposed matrix.

mat:rotxyz(ax, ay, az, inv_)

mat:rotxzy(ax, az, ay, inv_)

mat:rotyxz(ay, ax, az, inv_)

mat:rotyzx(ay, az, ax, inv_)

mat:rotzxy(az, ax, ay, inv_)

mat:rotzyx(az, ay, ax, inv_)

Return the \([3\times 3]\) real matrix mat filled with a 3D rotation of the first angle argument ax, ay or az around the x-axis, y-axis or z-axis respectively followed by another 3D rotation of the second angle argument ax, ay or az around the x-axis, y-axis or z-axis respectively of the frame rotated by the first rotation, and followed by a last 3D rotation of the third angle argument ax, ay or az around the x-axis, y-axis or z-axis respectively of the frame already rotated by the two first rotations. If inv is true, the returned matrix is the inverse rotations, i.e. the transposed matrix.

mat:torotxyz(inv_)

mat:torotxzy(inv_)

mat:torotyxz(inv_)

mat:torotyzx(inv_)

mat:torotzxy(inv_)

mat:torotzyx(inv_)

Return three real number representing the three angles ax, ay and az (always in this order) of the 3D rotations stored in the \([3\times 3]\) real matrix mat by the methods with corresponding names. If inv is true, the inverse rotations are returned, i.e. extracted from the transposed matrix.

mat:rotv(v, av, inv_)

Return the \([3\times 3]\) real matrix mat filled with a 3D rotation of angle av around the axis defined by the 3D vector-like v (see Axis-Angle representation). If inv is true, the returned matrix is the inverse rotation, i.e. the transposed matrix.

mat:torotv(v_, inv_)

Return a real number representing the angle of the 3D rotation around the axis defined by a 3D vector as stored in the \([3\times 3]\) real matrix mat by the method mat:rotv(). If the iterable v is provided, it is filled with the components of the unit vector that defines the axis of the rotation. If inv is true, the inverse rotation is returned, i.e. extracted from the transposed matrix.

mat:rotq(q, inv_)

Return the \([3\times 3]\) real matrix mat filled with a 3D rotation defined by the quaternion q (see Axis-Angle representation). If inv is true, the returned matrix is the inverse rotation, i.e. the transposed matrix.

mat:torotq(q_, inv_)

Return a quaternion representing the 3D rotation as stored in the \([3\times 3]\) real matrix mat by the method mat:rotq(). If the iterable q is provided, it is filled with the components of the quaternion otherwise the quaternion is returned in a list of length 4. If inv is true, the inverse rotation is returned, i.e. extracted from the transposed matrix.

Conversions

mat:tostring(sep_, lsep_)

Return the string containing the real, complex or integer matrix converted to string. The argument sep and lsep are used as separator for columns and rows respectively. The elements values are formated using tostring() that follows the option.numfmt string format for real numbers. Default: sep = " ", lsep = "\n".

mat:totable([d_, ] r_)

Return the table or r containing the real, complex or integer matrix converted to tables, i.e. one per row unless mat is a vector or the direction d = 'vec'.

Input and Output

mat:write(filename_, name_, eps_, line_, nl_)

Return the real, complex or integer matrix after writing it to the file filename opened with MAD.utility.openfile(). The content of the matrix mat is preceded by a header containing enough information to read it back. If name is provided, it is part of the header. If line = 'line', the matrix is displayed on a single line with rows separated by a semicolumn, otherwise it is displayed on multiple lines separated by nl. Elements with absolute value below eps are displayed as zeros. The formats defined by MAD.option.numfmt and MAD.option.intfmt are used to format numbers of matrix, cmatrix and imatrix respectively. Default: filename_ = io.stdout, name_ = '', eps_ = 0, line_ = nil, nl_ = '\n'.

mat:print(name_, eps_, line_, nl_)

Equivalent to mat:write(nil, name_, eps_, line_, nl_).

mat:read(filename_)

Return the real, complex or integer matrix read from the file filename opened with MAD.utility.openfile(). Note that the matrix mat is only used to call the method :read() and has no impact on the type and sizes of the returned matrix fully characterized by the content of the file. Default: filename_ = io.stdin.

Operators

#mat

Return the size of the real, complex or integer matrix mat, i.e. the number of elements interpreting the matrix as a vector.

mat[n]

Return the value of the element at index n of the real, complex or integer matrix mat for 1 <= n <= #mat, i.e. interpreting the matrix as a vector, nil otherwise.

mat[n] = v

Assign the value v to the element at index n of the real, complex or integer matrix mat for 1 <= n <= #mat, i.e. interpreting the matrix as a vector, otherwise raise an “out of bounds” error.

-mat

Return a real, complex or integer matrix resulting from the unary minus applied individually to all elements of the matrix mat.

num + mat

mat + num

mat + mat2

Return a matrix resulting from the sum of the left and right operands that must have compatible sizes. If one of the operand is a scalar, the operator will be applied individually to all elements of the matrix.

num + cmat

cpx + mat

cpx + cmat

mat + cpx

mat + cmat

cmat + num

cmat + cpx

cmat + mat

cmat + cmat2

Return a cmatrix resulting from the sum of the left and right operands that must have compatible sizes. If one of the operand is a scalar, the operator will be applied individually to all elements of the matrix.

idx + imat

imat + idx

imat + imat

Return a imatrix resulting from the sum of the left and right operands that must have compatible sizes. If one of the operand is a scalar, the operator will be applied individually to all elements of the matrix.

num - mat

mat - num

mat - mat2

Return a matrix resulting from the difference of the left and right operands that must have compatible sizes. If one of the operand is a scalar, the operator will be applied individually to all elements of the matrix.

num - cmat

cpx - mat

cpx - cmat

mat - cpx

mat - cmat

cmat - num

cmat - cpx

cmat - mat

cmat - cmat2

Return a cmatrix resulting from the difference of the left and right operands that must have compatible sizes. If one of the operand is a scalar, the operator will be applied individually to all elements of the matrix.

idx - imat

imat - idx

imat - imat

Return a imatrix resulting from the difference of the left and right operands that must have compatible sizes. If one of the operand is a scalar, the operator will be applied individually to all elements of the matrix.

num * mat

mat * num

mat * mat2

Return a matrix resulting from the product of the left and right operands that must have compatible sizes. If one of the operand is a scalar, the operator will be applied individually to all elements of the matrix. If the two operands are matrices, the mathematical matrix multiplication is performed.

num * cmat

cpx * mat

cpx * cmat

mat * cpx

mat * cmat

cmat * num

cmat * cpx

cmat * mat

cmat * cmat2

Return a cmatrix resulting from the product of the left and right operands that must have compatible sizes. If one of the operand is a scalar, the operator will be applied individually to all elements of the matrix. If the two operands are matrices, the mathematical matrix multiplication is performed.

idx * imat

imat * idx

Return a imatrix resulting from the product of the left and right operands that must have compatible sizes. If one of the operand is a scalar, the operator will be applied individually to all elements of the matrix.

num / mat

mat / num

mat / mat2

Return a matrix resulting from the division of the left and right operands that must have compatible sizes. If the right operand is a scalar, the operator will be applied individually to all elements of the matrix. If the left operand is a scalar the operation x/Y is converted to x (I/Y) where I is the identity matrix with compatible sizes. If the right operand is a matrix, the operation X/Y is performed using a system solver based on LU, QR or LQ factorisation depending on the shape of the system.

num / cmat

cpx / mat

cpx / cmat

mat / cpx

mat / cmat

cmat / num

cmat / cpx

cmat / mat

cmat / cmat2

Return a cmatrix resulting from the division of the left and right operands that must have compatible sizes. If the right operand is a scalar, the operator will be applied individually to all elements of the matrix. If the left operand is a scalar the operation x/Y is converted to x (I/Y) where I is the identity matrix with compatible sizes. If the right operand is a matrix, the operation X/Y is performed using a system solver based on LU, QR or LQ factorisation depending on the shape of the system.

imat / idx

Return a imatrix resulting from the division of the left and right operands, where the operator will be applied individually to all elements of the matrix.

mat % num

mat % mat

Return a matrix resulting from the modulo between the elements of the left and right operands that must have compatible sizes. If the right operand is a scalar, the operator will be applied individually to all elements of the matrix.

cmat % num

cmat % cpx

cmat % mat

cmat % cmat

Return a cmatrix resulting from the modulo between the elements of the left and right operands that must have compatible sizes. If the right operand is a scalar, the operator will be applied individually to all elements of the matrix.

imat % idx

imat % imat

Return a imatrix resulting from the modulo between the elements of the left and right operands that must have compatible sizes. If the right operand is a scalar, the operator will be applied individually to all elements of the matrix.

mat ^ n

cmat ^ n

Return a matrix or cmatrix resulting from n products of the square input matrix by itself. If n is negative, the inverse of the matrix is used for the product.

num == mat

num == cmat

num == imat

cpx == mat

cpx == cmat

mat == num

mat == cpx

mat == mat2

mat == cmat

cmat == num

cmat == cpx

cmat == mat

cmat == cmat2

imat == num

imat == imat2

Return false if the left and right operands have incompatible sizes or if any element differ in a one-to-one comparison, true otherwise. If one of the operand is a scalar, the operator will be applied individually to all elements of the matrix.

mat .. mat2

mat .. imat

imat .. mat

Return a matrix resulting from the row-oriented (horizontal) concatenation of the left and right operands. If the first element of the right operand mat (third case) is an integer, the resulting matrix will be a imatrix instead.

mat .. cmat

imat .. cmat

cmat .. mat

cmat .. imat

cmat .. cmat2

Return a cmatrix resulting from the row-oriented (horizontal) concatenation of the left and right operands.

imat .. imat2

Return a imatrix resulting from the row-oriented (horizontal) concatenation of the left and right operands.

Iterators

ipairs(mat)

Return an ipairs iterator suitable for generic for loops. The returned values are those given by mat[i].

C API

This C Application Programming Interface describes only the C functions declared in the scripting language and used by the higher level functions and methods presented before in this chapter. For more functions and details, see the C headers. The const vectors and matrices are inputs, while the non-const vectors and matrices are outpouts or are modified inplace.

Vector

void mad_vec_fill(num_t x, num_t r[], ssz_t n)

void mad_cvec_fill(cpx_t x, cpx_t r[], ssz_t n)

void mad_cvec_fill_r(num_t x_re, num_t x_im, cpx_t r[], ssz_t n)

void mad_ivec_fill(idx_t x, idx_t r[], ssz_t n)

Return the vector r of size n filled with the value of x.

void mad_vec_roll(num_t x[], ssz_t n, int nroll)

void mad_cvec_roll(cpx_t x[], ssz_t n, int nroll)

void mad_ivec_roll(idx_t x[], ssz_t n, int nroll)

Roll in place the values of the elements of the vector x of size n by nroll.

void mad_vec_copy(const num_t x[], num_t r[], ssz_t n)

void mad_vec_copyv(const num_t x[], cpx_t r[], ssz_t n)

void mad_cvec_copy(const cpx_t x[], cpx_t r[], ssz_t n)

void mad_ivec_copy(const idx_t x[], idx_t r[], ssz_t n)

Fill the vector r of size n with the content of the vector x.

void mad_vec_minmax(const num_t x[], log_t absf, idx_t r[2], ssz_t n)

void mad_cvec_minmax(const cpx_t x[], idx_t r[2], ssz_t n)

void mad_ivec_minmax(const idx_t x[], log_t absf, idx_t r[2], ssz_t n)

Return in r the indexes of the minimum and maximum values of the elements of the vector x of size n. If absf = TRUE, the function abs() is applied to the elements before comparison.

num_t mad_vec_eval(const num_t x[], num_t x0, ssz_t n)

cpx_t mad_cvec_eval(const cpx_t x[], cpx_t x0, ssz_t n)

void mad_cvec_eval_r(const cpx_t x[], num_t x0_re, num_t x0_im, cpx_t *r, ssz_t n)

Return in r or directly the evaluation of the vector x of size n at the point x0 using Honer’s scheme.

num_t mad_vec_sum(const num_t x[], ssz_t n)

cpx_t mad_cvec_sum(const cpx_t x[], ssz_t n)

void mad_cvec_sum_r(const cpx_t x[], cpx_t *r, ssz_t n)

num_t mad_vec_ksum(const num_t x[], ssz_t n)

cpx_t mad_cvec_ksum(const cpx_t x[], ssz_t n)

void mad_cvec_ksum_r(const cpx_t x[], cpx_t *r, ssz_t n)

Return in r or directly the sum of the values of the elements of the vector x of size n. The k versions use the Neumaier variants of the Kahan sum.

num_t mad_vec_mean(const num_t x[], ssz_t n)

cpx_t mad_cvec_mean(const cpx_t x[], ssz_t n)

void mad_cvec_mean_r(const cpx_t x[], cpx_t *r, ssz_t n)

Return in r or directly the mean of the vector x of size n.

num_t mad_vec_var(const num_t x[], ssz_t n)

cpx_t mad_cvec_var(const cpx_t x[], ssz_t n)

void mad_cvec_var_r(const cpx_t x[], cpx_t *r, ssz_t n)

Return in r or directly the unbiased variance with 2nd order correction of the vector x of size n.

void mad_vec_center(const num_t x[], num_t r[], ssz_t n)

void mad_cvec_center(const cpx_t x[], cpx_t r[], ssz_t n)

Return in r the centered, vector x of size n equivalent to x[i] - mean(x).

num_t mad_vec_norm(const num_t x[], ssz_t n)

num_t mad_cvec_norm(const cpx_t x[], ssz_t n)

Return the norm of the vector x of size n.

num_t mad_vec_dist(const num_t x[], const num_t y[], ssz_t n)

num_t mad_vec_distv(const num_t x[], const cpx_t y[], ssz_t n)

num_t mad_cvec_dist(const cpx_t x[], const cpx_t y[], ssz_t n)

num_t mad_cvec_distv(const cpx_t x[], const num_t y[], ssz_t n)

Return the distance between the vectors x and y of size n equivalent to norm(x - y).

num_t mad_vec_dot(const num_t x[], const num_t y[], ssz_t n)

cpx_t mad_cvec_dot(const cpx_t x[], const cpx_t y[], ssz_t n)

void mad_cvec_dot_r(const cpx_t x[], const cpx_t y[], cpx_t *r, ssz_t n)

cpx_t mad_cvec_dotv(const cpx_t x[], const num_t y[], ssz_t n)

void mad_cvec_dotv_r(const cpx_t x[], const num_t y[], cpx_t *r, ssz_t n)

num_t mad_vec_kdot(const num_t x[], const num_t y[], ssz_t n)

cpx_t mad_cvec_kdot(const cpx_t x[], const cpx_t y[], ssz_t n)

void mad_cvec_kdot_r(const cpx_t x[], const cpx_t y[], cpx_t *r, ssz_t n)

cpx_t mad_cvec_kdotv(const cpx_t x[], const num_t y[], ssz_t n)

void mad_cvec_kdotv_r(const cpx_t x[], const num_t y[], cpx_t *r, ssz_t n)

Return in r or directly the dot product between the vectors x and y of size n. The k versions use the Neumaier variants of the Kahan sum.

void mad_vec_cplx(const num_t re_[], const num_t im_[], cpx_t r[], ssz_t n)

Convert the real and imaginary vectors re and im of size n into the complex vector r.

void mad_cvec_reim(const cpx_t x[], num_t re_[], num_t ri_[], ssz_t n)

Split the complex vector x of size n into the real vector re and the imaginary vector ri.

void mad_cvec_conj(const cpx_t x[], cpx_t r[], ssz_t n)

Return in r the conjugate of the complex vector x of size n.

void mad_vec_abs(const num_t x[], num_t r[], ssz_t n)

void mad_cvec_abs(const cpx_t x[], num_t r[], ssz_t n)

Return in r the absolute value of the vector x of size n.

void mad_vec_add(const num_t x[], const num_t y[], num_t r[], ssz_t n)

void mad_vec_addn(const num_t x[], num_t y, num_t r[], ssz_t n)

void mad_vec_addc(const num_t x[], cpx_t y, cpx_t r[], ssz_t n)

void mad_vec_addc_r(const num_t x[], num_t y_re, num_t y_im, cpx_t r[], ssz_t n)

void mad_cvec_add(const cpx_t x[], const cpx_t y[], cpx_t r[], ssz_t n)

void mad_cvec_addv(const cpx_t x[], const num_t y[], cpx_t r[], ssz_t n)

void mad_cvec_addn(const cpx_t x[], num_t y, cpx_t r[], ssz_t n)

void mad_cvec_addc(const cpx_t x[], cpx_t y, cpx_t r[], ssz_t n)

void mad_cvec_addc_r(const cpx_t x[], num_t y_re, num_t y_im, cpx_t r[], ssz_t n)

void mad_ivec_add(const idx_t x[], const idx_t y[], idx_t r[], ssz_t n)

void mad_ivec_addn(const idx_t x[], idx_t y, idx_t r[], ssz_t n)

Return in r the sum of the scalar or vectors x and y of size n.

void mad_vec_sub(const num_t x[], const num_t y[], num_t r[], ssz_t n)

void mad_vec_subv(const num_t x[], const cpx_t y[], cpx_t r[], ssz_t n)

void mad_vec_subn(const num_t y[], num_t x, num_t r[], ssz_t n)

void mad_vec_subc(const num_t y[], cpx_t x, cpx_t r[], ssz_t n)

void mad_vec_subc_r(const num_t y[], num_t x_re, num_t x_im, cpx_t r[], ssz_t n)

void mad_cvec_sub(const cpx_t x[], const cpx_t y[], cpx_t r[], ssz_t n)

void mad_cvec_subv(const cpx_t x[], const num_t y[], cpx_t r[], ssz_t n)

void mad_cvec_subn(const cpx_t y[], num_t x, cpx_t r[], ssz_t n)

void mad_cvec_subc(const cpx_t y[], cpx_t x, cpx_t r[], ssz_t n)

void mad_cvec_subc_r(const cpx_t y[], num_t x_re, num_t x_im, cpx_t r[], ssz_t n)

void mad_ivec_sub(const idx_t x[], const idx_t y[], idx_t r[], ssz_t n)

void mad_ivec_subn(const idx_t y[], idx_t x, idx_t r[], ssz_t n)

Return in r the difference between the scalar or vectors x and y of size n.

void mad_vec_mul(const num_t x[], const num_t y[], num_t r[], ssz_t n)

void mad_vec_muln(const num_t x[], num_t y, num_t r[], ssz_t n)

void mad_vec_mulc(const num_t x[], cpx_t y, cpx_t r[], ssz_t n)

void mad_vec_mulc_r(const num_t x[], num_t y_re, num_t y_im, cpx_t r[], ssz_t n)

void mad_cvec_mul(const cpx_t x[], const cpx_t y[], cpx_t r[], ssz_t n)

void mad_cvec_mulv(const cpx_t x[], const num_t y[], cpx_t r[], ssz_t n)

void mad_cvec_muln(const cpx_t x[], num_t y, cpx_t r[], ssz_t n)

void mad_cvec_mulc(const cpx_t x[], cpx_t y, cpx_t r[], ssz_t n)

void mad_cvec_mulc_r(const cpx_t x[], num_t y_re, num_t y_im, cpx_t r[], ssz_t n)

void mad_ivec_mul(const idx_t x[], const idx_t y[], idx_t r[], ssz_t n)

void mad_ivec_muln(const idx_t x[], idx_t y, idx_t r[], ssz_t n)

Return in r the product of the scalar or vectors x and y of size n.

void mad_vec_div(const num_t x[], const num_t y[], num_t r[], ssz_t n)

void mad_vec_divv(const num_t x[], const cpx_t y[], cpx_t r[], ssz_t n)

void mad_vec_divn(const num_t y[], num_t x, num_t r[], ssz_t n)

void mad_vec_divc(const num_t y[], cpx_t x, cpx_t r[], ssz_t n)

void mad_vec_divc_r(const num_t y[], num_t x_re, num_t x_im, cpx_t r[], ssz_t n)

void mad_cvec_div(const cpx_t x[], const cpx_t y[], cpx_t r[], ssz_t n)

void mad_cvec_divv(const cpx_t x[], const num_t y[], cpx_t r[], ssz_t n)

void mad_cvec_divn(const cpx_t y[], num_t x, cpx_t r[], ssz_t n)

void mad_cvec_divc(const cpx_t y[], cpx_t x, cpx_t r[], ssz_t n)

void mad_cvec_divc_r(const cpx_t y[], num_t x_re, num_t x_im, cpx_t r[], ssz_t n)

void mad_ivec_divn(const idx_t x[], idx_t y, idx_t r[], ssz_t n)

Return in r the division of the scalar or vectors x and y of size n.

void mad_ivec_modn(const idx_t x[], idx_t y, idx_t r[], ssz_t n)

Return in r the modulo of the integer vector x of size n by the integer y.

void mad_vec_kadd(int k, const num_t a[], const num_t *x[], num_t r[], ssz_t n)

void mad_cvec_kadd(int k, const cpx_t a[], const cpx_t *x[], cpx_t r[], ssz_t n)

Return in r the linear combination of the k vectors in x of size n scaled by the k scalars in a.

void mad_vec_fft(const num_t x[], cpx_t r[], ssz_t n)

void mad_cvec_fft(const cpx_t x[], cpx_t r[], ssz_t n)

void mad_cvec_ifft(const cpx_t x[], cpx_t r[], ssz_t n)

Return in the vector r the 1D FFT and inverse of the vector x of size n.

void mad_vec_rfft(const num_t x[], cpx_t r[], ssz_t n)

Return in the vector r of size n/2+1 the 1D real FFT of the vector x of size n.

void mad_cvec_irfft(const cpx_t x[], num_t r[], ssz_t n)

Return in the vector r of size n the 1D real FFT inverse of the vector x of size n/2+1.

void mad_vec_nfft(const num_t x[], const num_t x_node[], cpx_t r[], ssz_t n, ssz_t nr)

void mad_cvec_nfft(const cpx_t x[], const num_t x_node[], cpx_t r[], ssz_t n, ssz_t nr)

Return in the vector r of size nr the 1D non-equispaced FFT of the vectors x and x_node of size n.

void mad_cvec_infft(const cpx_t x[], const num_t r_node[], cpx_t r[], ssz_t n, ssz_t nx)

Return in the vector r of size n the 1D non-equispaced FFT inverse of the vector x of size nx and the vector r_node of size n. Note that r_node here is the same vector as x_node in the 1D non-equispaced forward FFT.

Matrix

void mad_mat_rev(num_t x[], ssz_t m, ssz_t n, int d)

void mad_cmat_rev(cpx_t x[], ssz_t m, ssz_t n, int d)

void mad_imat_rev(idx_t x[], ssz_t m, ssz_t n, int d)

Reverse in place the matrix x following the direction d in {0,1,2,3} for respectively the entire matrix, each row, each column and the diagonal.

void mad_mat_center(num_t x[], ssz_t m, ssz_t n, int d)

void mad_cmat_center(cpx_t x[], ssz_t m, ssz_t n, int d)

Center in place the matrix x following the direction d in {0,1,2,3} for respectively the entire matrix, each row, each column and the diagonal.

void mad_mat_roll(num_t x[], ssz_t m, ssz_t n, int mroll, int nroll)

void mad_cmat_roll(cpx_t x[], ssz_t m, ssz_t n, int mroll, int nroll)

void mad_imat_roll(idx_t x[], ssz_t m, ssz_t n, int mroll, int nroll)

Roll in place the values of the elements of the matrix x of sizes [m, n] by mroll and nroll.

void mad_mat_eye(num_t x[], num_t v, ssz_t m, ssz_t n, ssz_t ldr)

void mad_cmat_eye(cpx_t x[], cpx_t v, ssz_t m, ssz_t n, ssz_t ldx)

void mad_cmat_eye_r(cpx_t x[], num_t v_re, num_t v_im, ssz_t m, ssz_t n, ssz_t ldr)

void mad_imat_eye(idx_t x[], idx_t v, ssz_t m, ssz_t n, ssz_t ldr)

Fill in place the matrix x of sizes [m, n] with zeros and v on the diagonal.

void mad_mat_copy(const num_t x[], num_t r[], ssz_t m, ssz_t n, ssz_t ldx, ssz_t ldr)

void mad_mat_copym(const num_t x[], cpx_t r[], ssz_t m, ssz_t n, ssz_t ldx, ssz_t ldr)

void mad_cmat_copy(const cpx_t x[], cpx_t r[], ssz_t m, ssz_t n, ssz_t ldx, ssz_t ldr)

void mad_imat_copy(const idx_t x[], idx_t r[], ssz_t m, ssz_t n, ssz_t ldx, ssz_t ldr)

void mad_imat_copym(const idx_t x[], num_t r[], ssz_t m, ssz_t n, ssz_t ldx, ssz_t ldr)

Fill the matrix r of sizes [m, n] and leading dimension ldr with the content of the matrix x of sizes [m, n] and leading dimension ldx.

void mad_mat_trans(const num_t x[], num_t r[], ssz_t m, ssz_t n)

void mad_cmat_trans(const cpx_t x[], cpx_t r[], ssz_t m, ssz_t n)

void mad_cmat_ctrans(const cpx_t x[], cpx_t r[], ssz_t m, ssz_t n)

void mad_imat_trans(const idx_t x[], idx_t r[], ssz_t m, ssz_t n)

Fill the matrix r of sizes [n, m] with the (conjugate) transpose of the matrix x of sizes [m, n].

void mad_mat_mul(const num_t x[], const num_t y[], num_t r[], ssz_t m, ssz_t n, ssz_t p)

void mad_mat_mulm(const num_t x[], const cpx_t y[], cpx_t r[], ssz_t m, ssz_t n, ssz_t p)

void mad_cmat_mul(const cpx_t x[], const cpx_t y[], cpx_t r[], ssz_t m, ssz_t n, ssz_t p)

void mad_cmat_mulm(const cpx_t x[], const num_t y[], cpx_t r[], ssz_t m, ssz_t n, ssz_t p)

Fill the matrix r of sizes [m, n] with the product of the matrix x of sizes [m, p] by the matrix y of sizes [p, n].

void mad_mat_tmul(const num_t x[], const num_t y[], num_t r[], ssz_t m, ssz_t n, ssz_t p)

void mad_mat_tmulm(const num_t x[], const cpx_t y[], cpx_t r[], ssz_t m, ssz_t n, ssz_t p)

void mad_cmat_tmul(const cpx_t x[], const cpx_t y[], cpx_t r[], ssz_t m, ssz_t n, ssz_t p)

void mad_cmat_tmulm(const cpx_t x[], const num_t y[], cpx_t r[], ssz_t m, ssz_t n, ssz_t p)

Fill the matrix r of sizes [m, n] with the product of the transposed matrix x of sizes [p, m] by the matrix y of sizes [p, n].

void mad_mat_mult(const num_t x[], const num_t y[], num_t r[], ssz_t m, ssz_t n, ssz_t p)

void mad_mat_multm(const num_t x[], const cpx_t y[], cpx_t r[], ssz_t m, ssz_t n, ssz_t p)

void mad_cmat_mult(const cpx_t x[], const cpx_t y[], cpx_t r[], ssz_t m, ssz_t n, ssz_t p)

void mad_cmat_multm(const cpx_t x[], const num_t y[], cpx_t r[], ssz_t m, ssz_t n, ssz_t p)

Fill the matrix r of sizes [m, n] with the product of the matrix x of sizes [m, p] and the transposed matrix y of sizes [n, p].

void mad_mat_dmul(const num_t x[], const num_t y[], num_t r[], ssz_t m, ssz_t n, ssz_t p)

void mad_mat_dmulm(const num_t x[], const cpx_t y[], cpx_t r[], ssz_t m, ssz_t n, ssz_t p)

void mad_cmat_dmul(const cpx_t x[], const cpx_t y[], cpx_t r[], ssz_t m, ssz_t n, ssz_t p)

void mad_cmat_dmulm(const cpx_t x[], const num_t y[], cpx_t r[], ssz_t m, ssz_t n, ssz_t p)

Fill the matrix r of size [m, n] with the product of the diagonal of the matrix x of sizes [m, p] by the matrix y of sizes [p, n]. If p = 1 then x will be interpreted as the diagonal of a square matrix.

void mad_mat_muld(const num_t x[], const num_t y[], num_t r[], ssz_t m, ssz_t n, ssz_t p)

void mad_mat_muldm(const num_t x[], const cpx_t y[], cpx_t r[], ssz_t m, ssz_t n, ssz_t p)

void mad_cmat_muld(const cpx_t x[], const cpx_t y[], cpx_t r[], ssz_t m, ssz_t n, ssz_t p)

void mad_cmat_muldm(const cpx_t x[], const num_t y[], cpx_t r[], ssz_t m, ssz_t n, ssz_t p)

Fill the matrix r of sizes [m, n] with the product of the matrix x of sizes [m, p] by the diagonal of the matrix y of sizes [p, n]. If p = 1 then y will be interpreted as the diagonal of a square matrix.

int mad_mat_div(const num_t x[], const num_t y[], num_t r[], ssz_t m, ssz_t n, ssz_t p, num_t rcond)

int mad_mat_divm(const num_t x[], const cpx_t y[], cpx_t r[], ssz_t m, ssz_t n, ssz_t p, num_t rcond)

int mad_cmat_div(const cpx_t x[], const cpx_t y[], cpx_t r[], ssz_t m, ssz_t n, ssz_t p, num_t rcond)

int mad_cmat_divm(const cpx_t x[], const num_t y[], cpx_t r[], ssz_t m, ssz_t n, ssz_t p, num_t rcond)

Fill the matrix r of sizes [m, n] with the division of the matrx x of sizes [m, p] by the matrix y of sizes [n, p]. The conditional number rcond is used by the solver to determine the effective rank of non-square systems. It returns the rank of the system.

int mad_mat_invn(const num_t y[], num_t x, num_t r[], ssz_t m, ssz_t n, num_t rcond)

int mad_mat_invc(const num_t y[], cpx_t x, cpx_t r[], ssz_t m, ssz_t n, num_t rcond)

int mad_mat_invc_r(const num_t y[], num_t x_re, num_t x_im, cpx_t r[], ssz_t m, ssz_t n, num_t rcond)

int mad_cmat_invn(const cpx_t y[], num_t x, cpx_t r[], ssz_t m, ssz_t n, num_t rcond)

int mad_cmat_invc(const cpx_t y[], cpx_t x, cpx_t r[], ssz_t m, ssz_t n, num_t rcond)

int mad_cmat_invc_r(const cpx_t y[], num_t x_re, num_t x_im, cpx_t r[], ssz_t m, ssz_t n, num_t rcond)

Fill the matrix r of sizes [n, m] with the inverse of the matrix y of sizes [m, n] scaled by the scalar x. The conditional number rcond is used by the solver to determine the effective rank of non-square systems. It returns the rank of the system.

int mad_mat_solve(const num_t a[], const num_t b[], num_t x[], ssz_t m, ssz_t n, ssz_t p, num_t rcond)

int mad_cmat_solve(const cpx_t a[], const cpx_t b[], cpx_t x[], ssz_t m, ssz_t n, ssz_t p, num_t rcond)

Fill the matrix x of sizes [n, p] with the minimum-norm solution of the linear least square problem \(\min \| A x - B \|\) where \(A\) is the matrix a of sizes [m, n] and \(B\) is the matrix b of sizes [m, p], using LU, QR or LQ factorisation depending on the shape of the system. The conditional number rcond is used by the solver to determine the effective rank of non-square system. It returns the rank of the system.

int mad_mat_ssolve(const num_t a[], const num_t b[], num_t x[], ssz_t m, ssz_t n, ssz_t p, num_t rcond, num_t s_[])

int mad_cmat_ssolve(const cpx_t a[], const cpx_t b[], cpx_t x[], ssz_t m, ssz_t n, ssz_t p, num_t rcond, num_t s_[])

Fill in the matrix x of sizes [n, p] with the minimum-norm solution of the linear least square problem \(\min \| A x - B \|\) where \(A\) is the matrix a of sizes [m, n] and \(B\) is the matrix b of sizes [m, p], using SVD factorisation. The conditional number rcond is used by the solver to determine the effective rank of non-square system. It returns the rank of the system and fill the optional column vector s of size min(m,n) with the singular values.

int mad_mat_gsolve(const num_t a[], const num_t b[], const num_t c[], const num_t d[], num_t x[], ssz_t m, ssz_t n, ssz_t p, num_t *nrm_)

int mad_cmat_gsolve(const cpx_t a[], const cpx_t b[], const cpx_t c[], const cpx_t d[], cpx_t x[], ssz_t m, ssz_t n, ssz_t p, num_t *nrm_)

Fill the column vector x of size n with the minimum-norm solution of the linear least square problem \(\min \| A x - C \|\) under the constraint \(B x = D\) where \(A\) is a matrix a of sizes [m, n], \(B\) is a matrix b of sizes [p, n], \(C\) is a column vector of size m and \(D\) is a column vector of size p, using QR or LQ factorisation depending on the shape of the system. This function also returns the status info and optionally the norm of the residues in the nrm.

int mad_mat_gmsolve(const num_t a[], const num_t b[], const num_t d[], num_t x[], num_t y[], ssz_t m, ssz_t n, ssz_t p)

int mad_cmat_gmsolve(const cpx_t a[], const cpx_t b[], const cpx_t d[], cpx_t x[], cpx_t y[], ssz_t m, ssz_t n, ssz_t p)

Fill the column vector x of size n and column vector y of size p with the minimum-norm solution of the linear Gauss-Markov problem \(\min_x \| y \|\) under the constraint \(A x + B y = D\) where \(A\) is a matrix a of sizes [m, n], \(B\) is a matrix b of sizes [m, p], and \(D\) is a column vector of size m, using QR or LQ factorisation depending on the shape of the system. This function also returns the status info.

int mad_mat_nsolve(const num_t a[], const num_t b[], num_t x[], ssz_t m, ssz_t n, ssz_t nc, num_t rcond, num_t r_[])

Fill the column vector x (of correctors kicks) of size n with the minimum-norm solution of the linear (best-kick) least square problem \(\min \| A x - B \|\) where \(A\) is the (response) matrix a of sizes [m, n] and \(B\) is a column vector (of monitors readings) of size m, using the MICADO [3] algorithm based on the Householder-Golub method [MICADO]. The argument nc is the maximum number of correctors to use with \(0 < n_c \leq n\) and the argument tol is a convergence threshold (on the residues) to stop the (orbit) correction if \(\| A x - B \| \leq m \times\) tol. This function also returns the updated number of correctors nc effectively used during the correction and the residues in the optional column vector r of size m.

int mad_mat_pcacnd(const num_t a[], idx_t ic[], ssz_t m, ssz_t n, ssz_t ns, num_t cut, num_t s_[])

int mad_cmat_pcacnd(const cpx_t a[], idx_t ic[], ssz_t m, ssz_t n, ssz_t ns, num_t cut, num_t s_[])

Fill the column vector ic of size n with the indexes of the columns to remove from the matrix a of sizes [m, n] using the Principal Component Analysis. The argument ns is the maximum number of singular values to consider and rcond is the conditioning number used to select the singular values versus the largest one, i.e. consider the ns larger singular values \(\sigma_i\) such that \(\sigma_i > \sigma_{\max}\times\)rcond. This function also returns the column vector of size min(m,n) filled with the singluar values. Default: ns_ = ncol, rcond_ = eps.

int mad_mat_svdcnd(const num_t a[], idx_t ic[], ssz_t m, ssz_t n, ssz_t ns, num_t cut, num_t s_[], num_t tol)

int mad_cmat_svdcnd(const cpx_t a[], idx_t ic[], ssz_t m, ssz_t n, ssz_t ns, num_t cut, num_t s_[], num_t tol)

Fill the column vector ic of size n with the indexes of the columns to remove from the matrix a of sizes [m, n] based on the analysis of the right matrix \(V\) from the SVD decomposition \(U S V\). The argument ns is the maximum number of singular values to consider and rcond is the conditioning number used to select the singular values versus the largest one, i.e. consider the ns larger singular values \(\sigma_i\) such that \(\sigma_i > \sigma_{\max}\times\)rcond. The argument tol is a threshold similar to rcond used to reject components in \(V\) that have similar or opposite effect than components already encountered. This function also returns the real column vector of size min(m,n) filled with the singluar values. Default: ns_ = min(m,n), rcond_ = eps.

int mad_mat_svd(const num_t x[], num_t u[], num_t s[], num_t v[], ssz_t m, ssz_t n)

int mad_cmat_svd(const cpx_t x[], cpx_t u[], num_t s[], cpx_t v[], ssz_t m, ssz_t n)

Fill the column vector s of size min(m,n) with the singular values, and the two matrices u of sizes [m, m] and v of sizes [n, n] with the SVD factorisation of the matrix x of sizes [m,n], and returns the status info. The singular values are positive and sorted in decreasing order of values, i.e. largest first.

int mad_mat_eigen(const num_t x[], cpx_t w[], num_t vl[], num_t vr[], ssz_t n)

int mad_cmat_eigen(const cpx_t x[], cpx_t w[], cpx_t vl[], cpx_t vr[], ssz_t n)

Fill the column vector w of size n with the eigenvalues followed by the status info and the two optional matrices vr and vl of sizes [n, n] containing the left and right eigenvectors resulting from the Eigen Decomposition of the square matrix x of sizes [n, n]. The eigenvectors are normalized to have unit Euclidean norm and their largest component real, and satisfy \(X v_r = \lambda v_r\) and \(v_l X = \lambda v_l\).

int mad_mat_det(const num_t x[], num_t *r, ssz_t n)

int mad_cmat_det(const cpx_t x[], cpx_t *r, ssz_t n)

Return in r, the Determinant of the square matrix mat of sizes [n, n] using LU factorisation for better numerical stability, and return the status info.

void mad_mat_fft(const num_t x[], cpx_t r[], ssz_t m, ssz_t n)

void mad_cmat_fft(const cpx_t x[], cpx_t r[], ssz_t m, ssz_t n)

void mad_cmat_ifft(const cpx_t x[], cpx_t r[], ssz_t m, ssz_t n)

Fill the matrix r with the 2D FFT and inverse of the matrix x of sizes [m, n].

void mad_mat_rfft(const num_t x[], cpx_t r[], ssz_t m, ssz_t n)

Fill the matrix r of sizes [m, n/2+1] with the 2D real FFT of the matrix x of sizes [m, n].

void mad_cmat_irfft(const cpx_t x[], num_t r[], ssz_t m, ssz_t n)

Fill the matrix r of sizes [m, n] with the 1D real FFT inverse of the matrix x of sizes [m, n/2+1].

void mad_mat_nfft(const num_t x[], const num_t x_node[], cpx_t r[], ssz_t m, ssz_t n, ssz_t nr)

void mad_cmat_nfft(const cpx_t x[], const num_t x_node[], cpx_t r[], ssz_t m, ssz_t n, ssz_t nr)

Fill the matrix r of sizes [m, nr] with the 2D non-equispaced FFT of the matrices x and x_node of sizes [m, n].

void mad_cmat_infft(const cpx_t x[], const num_t r_node[], cpx_t r[], ssz_t m, ssz_t n, ssz_t nx)

Fill the matrix r of sizes [m, n] with the 2D non-equispaced FFT inverse of the matrix x of sizes [m, nx] and the matrix r_node of sizes [m, n]. Note that r_node here is the same matrix as x_node in the 2D non-equispaced forward FFT.

void mad_mat_sympconj(const num_t x[], num_t r[], ssz_t n)

void mad_cmat_sympconj(const cpx_t x[], cpx_t r[], ssz_t n)

Return in r the symplectic ‘conjugate’ of the vector x of size n.

num_t mad_mat_symperr(const num_t x[], num_t r_[], ssz_t n, num_t *tol_)

num_t mad_cmat_symperr(const cpx_t x[], cpx_t r_[], ssz_t n, num_t *tol_)

Return the norm of the symplectic error and fill the optional matrix r with the symplectic deviation of the matrix x. The optional argument tol is used as the tolerance to check if the matrix x is symplectic or not, and saves the result as 0 (non-symplectic) or 1 (symplectic) within tol for output.

void mad_vec_dif(const num_t x[], const num_t y[], num_t r[], ssz_t n)

void mad_vec_difv(const num_t x[], const cpx_t y[], cpx_t r[], ssz_t n)

void mad_cvec_dif(const cpx_t x[], const cpx_t y[], cpx_t r[], ssz_t n)

void mad_cvec_difv(const cpx_t x[], const num_t y[], cpx_t r[], ssz_t n)

Fill the matrix r of sizes [m, n] with the absolute or relative differences between the elements of the matrix x and y with compatible sizes. The relative difference is taken for the values with magnitude greater than 1, otherwise it takes the absolute difference.

Rotations

void mad_mat_rot(num_t x[2 * 2], num_t a)

Fill the matrix x with a 2D rotation of angle a.

void mad_mat_rotx(num_t x[3 * 3], num_t ax)

void mad_mat_roty(num_t x[3 * 3], num_t ay)

void mad_mat_rotz(num_t x[3 * 3], num_t az)

Fill the matrix x with the 3D rotation of angle a? around the axis given by the suffix ? in {x,y,z}.

void mad_mat_rotxy(num_t x[3 * 3], num_t ax, num_t ay, log_t inv)

void mad_mat_rotxz(num_t x[3 * 3], num_t ax, num_t az, log_t inv)

void mad_mat_rotyz(num_t x[3 * 3], num_t ay, num_t az, log_t inv)

Fill the matrix x with the two successive 3D rotations of angles a? around the two axis given by the suffixes ? in {x,y,z}. If inv = 1 returns the inverse rotations, i.e. the transpose of the matrix x. Note that the first rotation changes the axis orientation of the second rotation.

void mad_mat_rotxyz(num_t x[3 * 3], num_t ax, num_t ay, num_t az, log_t inv)

void mad_mat_rotxzy(num_t x[3 * 3], num_t ax, num_t ay, num_t az, log_t inv)

void mad_mat_rotyxz(num_t x[3 * 3], num_t ax, num_t ay, num_t az, log_t inv)

Fill the matrix x with the three successive 3D rotations of angles a? around the three axis given by the suffixes ? in {x,y,z}. If inv = 1 returns the inverse rotations, i.e. the transpose of the matrix x. Note that the first rotation changes the axis orientation of the second rotation, which changes the axis orientation of the third rotation.

void mad_mat_torotxyz(const num_t x[3 * 3], num_t r[3], log_t inv)

void mad_mat_torotxzy(const num_t x[3 * 3], num_t r[3], log_t inv)

void mad_mat_torotyxz(const num_t x[3 * 3], num_t r[3], log_t inv)

Fill the vector of the three angles r around the axis {x,y,z}, {x,z,y} and {y,x,z} from the matrix x. If inv = 1, it takes the inverse rotations, i.e. the transpose of the matrix x.

void mad_mat_rotv(num_t x[3 * 3], const num_t v[3], num_t a, log_t inv)

Fill the matrix x with the 3D rotation of angle a around the vector v. If inv = 1 returns the inverse rotations, i.e. the transpose of the matrix.

num_t mad_mat_torotv(const num_t x[3 * 3], num_t v_[3], log_t inv)

Return the angle and fill the optional vector v with the 3D rotations in x. If inv = 1, it takes the inverse rotations, i.e. the transpose of the matrix x.

void mad_mat_rotq(num_t x[3 * 3], const num_t q[4], log_t inv)

Fill the matrix x with the 3D rotation given by the quaternion q. If inv = 1 returns the inverse rotations, i.e. the transpose of the matrix.

void mad_mat_torotq(const num_t x[3 * 3], num_t q[4], log_t inv)

Fill the quaternion q with the 3D rotations in x. If inv = 1, it takes the inverse rotations, i.e. the transpose of the matrix x.

Misalignments

void mad_mat_rtbar(num_t Rb[3 * 3], num_t Tb[3], num_t el, num_t ang, num_t tlt, const num_t R_[3 * 3], const num_t T[3])

Compute as output the rotation matrix Rb, i.e. \(\bar{R}\), and the translation vector Tb, i.e. \(\bar{T}\), used to restore the global frame at exit of a misaligned element in survey, given as input the element length el, angle ang, tilt tlt, and the rotation matrix R and the translation vector T at entry.

Miscellaneous

void mad_fft_cleanup(void)

Cleanup data allocated by the FFTW library.

References

[MICADO] (1,2)

2. Autin, and Y. Marti, “Closed Orbit Correction of Alternating Gradient Machines using a Small Number of Magnets”, CERN ISR-MA/73-17, Mar. 1973.

[MATFUN]

N.J. Higham, and X. Liu, “A Multiprecision Derivative-Free Schur–Parlett Algorithm for Computing Matrix Functions”, SIAM J. Matrix Anal. Appl., Vol. 42, No. 3, pp. 1401-1422, 2021.

Footnotes

[1]

For true Functional Programming, see the module MAD.lfun, a binding of the LuaFun library adapted to the ecosystem of MAD-NG.

[2]

The solvers are based, among others, on the following Lapack drivers:

• dgesv() and zgesv() for LU factorization.

• dgelsy() and zgelsy() for QR or LQ factorization.

• dgelsd() and zgelsd() for SVD factorisation.

• dgees() and zgees() for Schur factorisation.

• dgglse() and zgglse() for equality-constrained linear Least Squares problems.

• dggglm() and zggglm() for general Gauss-Markov linear model problems.

[3] (1,2)

MICADO stands for “Minimisation des CArrés des Distortions d’Orbite” in french.