.\"
.\" Copyright (c) 2006-2022 Julien Nadeau Carriere <vedge@csoft.net>
.\"
.\" Redistribution and use in source and binary forms, with or without
.\" modification, are permitted provided that the following conditions
.\" are met:
.\" 1. Redistributions of source code must retain the above copyright
.\"    notice, this list of conditions and the following disclaimer.
.\" 2. Redistributions in binary form must reproduce the above copyright
.\"    notice, this list of conditions and the following disclaimer in the
.\"    documentation and/or other materials provided with the distribution.
.\" 
.\" THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
.\" IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
.\" WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
.\" ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT,
.\" INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
.\" (INCLUDING BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
.\" SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
.\" HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
.\" STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
.\" IN ANY WAY OUT OF THE USE OF THIS SOFTWARE EVEN IF ADVISED OF THE
.\" POSSIBILITY OF SUCH DAMAGE.
.\"
.Dd December 21, 2022
.Dt M_MATRIX 3
.Os Agar 1.7
.Sh NAME
.Nm M_Matrix
.Nd Agar-Math matrix-related functions
.Sh SYNOPSIS
.Bd -literal
#include <agar/core.h>
#include <agar/gui.h>
#include <agar/math/m.h>
.Ed
.Sh DESCRIPTION
.\" IMAGE(/widgets/M_Matview.png, "Graphical representation of a matrix")
.\" BEGIN DUP in M_Vector(3)
The
.Xr M_Vector 3
and
.Nm
interfaces implement linear algebra operations on (real or complex valued)
.Va n
dimensional vectors, and
.Va m
by
.Va n
matrices.
Optimized interfaces are provided for fixed-dimensional types (which have
entries directly accessible as
.Va x ,
.Va y ,
.Va z
and
.Va w ) .
Arbitrary-dimensional types may or may not use fixed arrays in memory.
For example, the "sparse" backend uses a sparse matrix representation,
and the "db" backend stores vector entries in a database.
.Pp
Backends can be selected at run-time, or Agar-Math can be compiled to
provide inline expansions of all operations of a specific backend.
Vector extensions (such as SSE and AltiVec) are used by default, if
a runtime cpuinfo check determines that they are available (the build
remains compatible with non-vector platforms, at the cost of extra
function calls).
For best performance, Agar should be compiled with "--with-sse=inline",
or "--with-altivec=inline".
.\" END DUP
.Sh M-BY-N MATRICES
The following routines operate on dynamically-allocated
.Va m
by
.Va n
matrices:
.Pp
.Bl -tag -width "sparse " -compact
.It fpu
Native scalar floating point methods.
.It sparse
Methods optimized for large, sparse matrices.
Based on the excellent Sparse 1.4 package by Kenneth Kundert.
.El
.Sh M-BY-N MATRICES: INITIALIZATION
.nr nS 1
.Ft "M_Matrix *"
.Fn M_New "Uint m" "Uint n"
.Pp
.Ft "void"
.Fn M_Free "M_Matrix *M"
.Pp
.Ft "int"
.Fn M_Resize "M_Matrix *M" "Uint m" "Uint n"
.Pp
.Ft "void"
.Fn M_SetIdentity "M_Matrix *M"
.Pp
.Ft "void"
.Fn M_SetZero "M_Matrix *M"
.Pp
.Ft "int"
.Fn M_Copy "M_Matrix *D" "const M_Matrix *A"
.Pp
.Ft "M_Matrix *"
.Fn M_Dup "const M_Matrix *M"
.Pp
.Ft "M_Matrix *"
.Fn M_ReadMatrix "AG_DataSource *ds"
.Pp
.Ft "void"
.Fn M_WriteMatrix "AG_DataSource *ds" "const M_Matrix *A"
.Pp
.nr nS 0
The
.Fn M_New
function allocates a new
.Fa m
by
.Fa n
matrix.
.Fn M_Free
releases all resources allocated for the specified matrix.
.Fn M_Resize
resizes
.Fa M
to
.Fa m
by
.Fa n .
Existing entries are preserved, but new entries are left uninitialized.
If insufficient memory is available, -1 is returned and an error message
is set.
On success, the function returns 0.
.Pp
.Fn M_SetIdentity
initializes
.Fa M
to the identity matrix.
.Fn M_SetZero
initializes
.Fa M
to all zeros.
.Pp
.Fn M_Copy
copies the contents of matrix
.Fa A
into
.Fa D ,
which is assumed to have the same dimensions (otherwise, -1 is returned).
.Fn M_Dup
returns a duplicate of
.Fa M .
.Pp
The
.Fn M_ReadMatrix
and
.Fn M_WriteMatrix
functions are used to (de)serialize the contents of matrix
.Fa A
from/to the specified
.Xr AG_DataSource 3 .
.Sh M-BY-N MATRICES: ACCESSING ELEMENTS
.nr nS 1
.Ft "M_Real"
.Fn M_Get "M_Matrix *M" "Uint i" "Uint j"
.Pp
.Ft "void"
.Fn M_Set "M_Matrix *M" "Uint i" "Uint j" "M_Real val"
.Pp
.Ft "M_Real *"
.Fn M_GetElement "M_Matrix *M" "Uint i" "Uint j"
.Pp
.Ft "void"
.Fn M_ToFloats "float *values" "const M_Matrix *A"
.Pp
.Ft "void"
.Fn M_ToDoubles "double *values" "const M_Matrix *A"
.Pp
.Ft "void"
.Fn M_FromFloats "M_Matrix *A" "const float *values"
.Pp
.Ft "void"
.Fn M_FromDoubles "M_Matrix *A" "const double *values"
.Pp
.Ft "void"
.Fn M_Print "const M_Matrix *A"
.Pp
.nr nS 0
The
.Fn M_Get
and
.Fn M_Set
routines respectively retrieve and set the element
.Fa i ,
.Fa j .
.Pp
.Fn M_GetElement
returns a pointer to the element
.Fa i ,
.Fa j .
As long as the entry exists, it is safe to read and write the element.
.Pp
The
.Fn M_ToFloats
and
.Fn M_ToDoubles
functions return a representation of matrix
.Fa A
as an array of
.Ft float
or
.Ft double
values in row-major order.
The
.Fn M_FromFloats
and
.Fn M_FromDoubles
functions initialize matrix
.Fa A
from an array of
.Ft float
or
.Ft double
values in row-major order.
In both cases, it is assumed that the arrays are of the correct size for
the given matrix dimensions.
.Pp
.Fn M_Print
dumps the individual matrix entries to the standard error output.
It is only for debugging purposes. Agar GUI applications can use the provided
.Xr M_Matview 3
widget to display matrix contents.
.Sh M-BY-N MATRICES: OPERATIONS
.nr nS 1
.Ft "M_Matrix *"
.Fn M_Transpose "M_Matrix *M"
.Pp
.Ft "M_Matrix *"
.Fn M_Add "const M_Matrix *A" "const M_Matrix *B"
.Pp
.Ft "int"
.Fn M_Addv "M_Matrix *A" "const M_Matrix *B"
.Pp
.Ft "void"
.Fn M_AddToDiag "M_Matrix *A" "M_Real value"
.Pp
.Ft "M_Matrix *"
.Fn M_DirectSum "const M_Matrix *A" "const M_Matrix *B"
.Pp
.Ft "M_Matrix *"
.Fn M_Mul "const M_Matrix *A" "const M_Matrix *B"
.Pp
.Ft "int"
.Fn M_Mulv "const M_Matrix *A" "const M_Matrix *B" "M_Matrix *AB"
.Pp
.Ft "M_Matrix *"
.Fn M_EntMul "const M_Matrix *A" "const M_Matrix *B"
.Pp
.Ft "int"
.Fn M_EntMulv "const M_Matrix *A" "const M_Matrix *B" "M_Matrix *AB"
.Pp
.Ft "void"
.Fn M_Compare "const M_Matrix *A" "const M_Matrix *B" "M_Real *diff"
.Pp
.Ft "int"
.Fn M_Trace "M_Real *trace" "const M_Matrix *A"
.Pp
.Ft "void"
.Fn M_IsSquare "M_Matrix *A"
.Pp
.Ft "M_Matrix *"
.Fn M_GaussJordan "const M_Matrix *A" "M_Matrix *b"
.Pp
.Ft "int"
.Fn M_GaussJordanv "M_Matrix *A" "M_Matrix *b"
.Pp
.Ft "int"
.Fn M_FactorizeLU "M_Matrix *A"
.Pp
.Ft "void"
.Fn M_BacksubstLU "M_Matrix *LU" "M_Vector *b"
.Pp
.Ft "void"
.Fn M_MNAPreorder "M_Matrix *A"
.Pp
.nr nS 0
The
.Fn M_Transpose
function returns the transpose of
.Fa M
(i.e., all
.Fa i ,
.Fa j
elements are swapped against
.Fa j ,
.Fa i
elements).
.Pp
.Fn M_Add
returns the sum of the matrices
.Fa A
and
.Fa B .
The
.Fn M_Addv
variant returns the sum into an existing matrix, returning -1 if the
dimensions are incorrect.
.Pp
The
.Fn M_AddToDiag
routine adds
.Va value
to each diagonal entry
.Fa i ,
.Fa i
of matrix
.Fa A .
.Pp
.Fn M_DirectSum
returns the direct sum of
.Fa A
and
.Fa B .
.Pp
.Fn M_Mul
returns the product of matrices
.Fa A
and
.Fa B .
The
.Fn M_Mulv
variant returns the product into an existing matrix, returning -1 if the
dimensions are incorrect.
.Fn M_EntMul
and
.Fn M_EntMulv
perform entrywise multiplication as opposed to matrix multiplication.
.Pp
The
.Fn M_Compare
function compares each entry of
.Fa A
and
.Fa B ,
returning the largest difference into
.Fa diff .
.Pp
.Fn M_Trace
returns the trace (the sum of elements on the diagonal) of a square matrix
.Fa A
into
.Fa trace .
.Pp
The
.Fn M_IsSquare
function returns 1 if
.Fa A
is a square (n-by-n) matrix.
.Pp
The
.Fn M_GaussJordan
function solves for
.Va x
in
.Va Ax = b .
The solver replaces the contents of
.Va A
by its inverse, and returns the solution vector
into
.Fa b .
.Pp
The
.Fn M_FactorizeLU
routine computes the
.Va LU
factorization of a square matrix
.Fa A .
If successful, the original contents of
.Fa A
are destroyed and replaced by the
.Va LU
factorization.
On error, -1 is returned.
Partial pivoting information is recorded in the
.Nm
structure for subsequent backsubstitution.
.Pp
The
.Fn M_BacksubstLU
routine solves a system of linear equations represented by a LU factorization
.Fa LU
(previously computed by
.Fn M_FactorizeLU )
and a right-hand side
.Fa b .
The solution vector is returned into
.Fa b .
.Pp
The
.Fn M_MNAPreorder
routine attempts to remove zeros from the diagonal, by taking into
account the structure of modified node admittance matrices (found in
applications such as electronic simulators).
.Sh 4-BY-4 MATRICES
The following routines are optimized for 4x4 matrices, as frequently
encountered in computer graphics.
Entries are directly accessible as structure members.
Available backends include:
.Pp
.Bl -tag -width "fpu " -compact
.It fpu
Native scalar floating point methods.
.It sse
Accelerate operations using Streaming SIMD Extensions (SSE).
.El
.\" MANLINK(M_Matrix44)
.Sh 4-BY-4 MATRICES: INITIALIZATION
.nr nS 1
.Ft "M_Matrix44"
.Fn M_MatZero44 "void"
.Pp
.Ft "void"
.Fn M_MatZero44v "M_Matrix44 *Z"
.Pp
.Ft "M_Matrix44"
.Fn M_MatIdentity44 "void"
.Pp
.Ft "void"
.Fn M_MatIdentity44v "M_Matrix44 *I"
.Pp
.Ft "void"
.Fn M_MatCopy44 "M_Matrix44 *Mdst" "const M_Matrix44 *Msrc"
.Pp
.nr nS 0
The
.Fn M_MatZero44
and
.Fn M_MatZero44v
functions initializes the target matrix
.Fa Z
to the zero matrix.
.Pp
.Fn M_MatIdentity44
and
.Fn M_MatIdentity44v
initializes the target matrix
.Fa I
to the identity matrix.
.Pp
The
.Fn M_MatCopy44
routine copies the contents of matrix
.Fa Msrc
into
.Fa Mdst.
The original contents of
.Fa Mdst
are overwritten.
.Sh 4-BY-4 MATRICES: ACCESSING ELEMENTS
The elements of
.Ft M_Matrix44
are directly accessible via the
.Va m[4][4]
member of the structure.
Elements of the matrix are stored in row-major format.
The structure is defined as:
.Bd -literal
.\" SYNTAX(c)
#ifdef HAVE_SSE
typedef union m_matrix44 {
	struct { __m128 m1, m2, m3, m4; };
	float m[4][4];
} M_Matrix44;
#else
typedef struct m_matrix44 {
	M_Real m[4][4];
} M_Matrix44;
#endif
.Ed
.Pp
Notice that SIMD extensions force single-precision floats, regardless of
the precision for which Agar-Math was built (if a 4x4 matrix of higher
precision is required,
the general
.Ft M_Matrix
type may be used).
.Pp
The following functions convert between
.Ft M_Matrix44
and numerical arrays:
.Pp
.nr nS 1
.Ft "void"
.Fn M_MatToFloats44 "float *flts" "const M_Matrix44 *A"
.Pp
.Ft "void"
.Fn M_MatToDoubles44 "double *dbls" "const M_Matrix44 *A"
.Pp
.Ft "void"
.Fn M_MatFromFloats44 "M_Matrix44 *M" "const float *flts"
.Pp
.Ft "void"
.Fn M_MatFromDoubles44 "M_Matrix44 *M" "const double *dbls"
.Pp
.nr nS 0
.Fn M_MatToFloats44
converts matrix
.Fa A
to a 4x4 array of floats
.Fa flts .
.Fn M_MatToDoubles44
converts matrix
.Fa A
to a 4x4 array of doubles
.Fa dbls .
.Fn M_MatFromFloats44
initializes matrix
.Fa M
from the contents of a 4x4 array of floats
.Fa flts .
.Fn M_MatFromDoubles44
initializes matrix
.Fa M
from the contents of a 4x4 array of doubles
.Fa dbls .
.Sh 4-BY-4 MATRICES: OPERATIONS
.nr nS 1
.Ft "M_Matrix44"
.Fn M_MatTranspose44 "M_Matrix44 A"
.Pp
.Ft "M_Matrix44"
.Fn M_MatTranspose44p "const M_Matrix44 *A"
.Pp
.Ft "void"
.Fn M_MatTranspose44v "M_Matrix44 *A"
.Pp
.Ft "M_Matrix44"
.Fn M_MatInvert44 "M_Matrix44 A"
.Pp
.Ft "int"
.Fn M_MatInvertElim44 "M_Matrix44 A" "M_Matrix44 *Ainv"
.Pp
.Ft "M_Matrix44"
.Fn M_MatMult44 "M_Matrix44 A" "M_Matrix44 B"
.Pp
.Ft "void"
.Fn M_MatMult44v "M_Matrix44 *A" "const M_Matrix44 *B"
.Pp
.Ft "void"
.Fn M_MatMult44pv "M_Matrix44 *AB" "const M_Matrix44 *A" "const M_Matrix44 *B"
.Pp
.Ft "M_Vector4"
.Fn M_MatMultVector44 "M_Matrix44 A" "M_Vector4 x"
.Pp
.Ft "M_Vector4"
.Fn M_MatMultVector44p "const M_Matrix44 *A" "const M_Vector4 *x"
.Pp
.Ft "void"
.Fn M_MatMultVector44v "M_Vector4 *x" "const M_Matrix44 *A"
.Pp
.Ft "void"
.Fn M_MatRotateAxis44 "M_Matrix44 *T" "M_Real theta" "M_Vector3 axis"
.Pp
.Ft "void"
.Fn M_MatOrbitAxis44 "M_Matrix44 *T" "M_Vector3 center" "M_Vector3 axis" "M_Real theta"
.Pp
.Ft "void"
.Fn M_MatRotateEul44 "M_Matrix44 *T" "M_Real pitch" "M_Real roll" "M_Real yaw"
.Pp
.Ft "void"
.Fn M_MatRotate44I "M_Matrix44 *T" "M_Real theta"
.Pp
.Ft "void"
.Fn M_MatRotate44J "M_Matrix44 *T" "M_Real theta"
.Pp
.Ft "void"
.Fn M_MatRotate44K "M_Matrix44 *T" "M_Real theta"
.Pp
.Ft "void"
.Fn M_MatTranslate44v "M_Matrix44 *T" "M_Vector3 v"
.Pp
.Ft "void"
.Fn M_MatTranslate44 "M_Matrix44 *T" "M_Real x" "M_Real y" "M_Real z"
.Pp
.Ft "void"
.Fn M_MatTranslate44X "M_Matrix44 *T" "M_Real c"
.Pp
.Ft "void"
.Fn M_MatTranslate44Y "M_Matrix44 *T" "M_Real c"
.Pp
.Ft "void"
.Fn M_MatTranslate44Z "M_Matrix44 *T" "M_Real c"
.Pp
.Ft "void"
.Fn M_MatScale44 "M_Matrix44 *T" "M_Real x" "M_Real y" "M_Real z" "M_Real w"
.Pp
.Ft "void"
.Fn M_MatUniScale44 "M_Matrix44 *T" "M_Real c"
.Pp
.nr nS 0
The
.Fn M_MatTranspose44 ,
.Fn M_MatTranspose44p
and
.Fn M_MatTranspose44v
function compute and return the transpose of matrix
.Fa A
(i.e., all elements
.Va i,j
are swapped for elements
.Va j,i ) .
.Pp
The function
.Fn M_MatInvert44
computes the inverse of
.Fa A
using Cramer's rule and cofactors.
If the matrix is not invertible, the return value is undefined.
.Pp
The
.Fn M_MatInvertElim44
function computes the inverse of
.Fa A
by systematic Gaussian elimination.
If the matrix is not invertible (singular up to
.Dv M_MACHEP
precision), the function fails.
.Pp
.Fn M_MatMult44 ,
.Fn M_MatMult44v
and
.Fn M_MatMult44pv
compute the product of matrices
.Fa A
and
.Fa B .
.Pp
The
.Fn M_MatMultVector44 ,
.Fn M_MatMultVector44p
and
.Fn M_MatMultVector44v
functions perform matrix-vector multiplication
.Va Ax ,
and returns
.Va x .
.Pp
.Fn M_MatRotateAxis44
multiplies matrix
.Fa T
against a rotation matrix describing a rotation of
.Fa theta
radians about
.Fa axis
(relative to the origin).
The
.Fn M_MatOrbitAxis44
variant takes
.Fa axis
to be relative to the specified
.Fa center
point as opposed to the origin.
.Pp
.Fn M_MatRotateEul44
multiplies
.Fa T
against a matrix describing a rotation about the origin in terms of Euler
angles
.Fa pitch ,
.Fa roll
and
.Fa yaw
(given in radians).
.Pp
.Fn M_MatRotate44I ,
.Fn M_MatRotate44J
and
.Fn M_MatRotate44K
multiply
.Fa T
with a matrix describing a rotation of
.Fa theta
radians about the basis vector
.Va i ,
.Va j
or
.Va k ,
respectively.
.Pp
.Fn M_MatTranslate44v
multiplies
.Fa T
against a matrix describing a translation by vector
.Fa v .
.Fn M_MatTranslate44 ,
.Fn M_MatTranslate44X ,
.Fn M_MatTranslate44Y
and
.Fn M_MatTranslate44Z
accept individual coordinate arguments.
.Pp
.Fn M_MatScale44
multiplies
.Fa T
against a matrix describing uniform/non-uniform scaling by [x,y,z,w].
.Fn M_MatUniScale44
performs uniform scaling by
.Fa c .
.Sh SEE ALSO
.Xr AG_Intro 3 ,
.Xr M_Complex 3 ,
.Xr M_Quaternion 3 ,
.Xr M_Real 3 ,
.Xr M_Vector 3
.Pp
.Bl -tag -compact
.It Lk https://sparse.sourceforge.net/ Sparse 1.4
.El
.Sh HISTORY
The
.Nm
interface first appeared in Agar 1.3.3.