Fixed size matrix decompositions

Gandalf supports several of the standard matrix decompositions. Cholesky factorisation applies to any positive definite symmetric matrix $S$, producing the lower triangular matrix $L$ so that

$$S = L L^{\top}$$

It can be computed for $3\times 3$ symmetric matrices using the routines

      Gan_SquMatrix33 sm33S, sm33L;

      /* ... set up sm33S using e.g. gan_symmat33_fill_q() ... */
      gan_symmat33_cholesky_q ( &sm33S, &sm33L ); /* L = chol(S) */
      sm33L = gan_symmat33_cholesky_s ( &sm33S ); /* L = chol(S) */

There is also a routine for computing the Cholesky factorisation in-place in the input matrix, converting an input symmetric matrix $A$ into a lower triangular matrix:

      Gan_SquMatrix33 sm33A;

      /* ... set up sm33A as symmetric using e.g. gan_symmat33_fill_q() ... */
      gan_symmat33_cholesky_i ( &sm33A ); /* A = chol(A) */

The routines gan_symmat33_cholesky_[qi]() return NULL and invoke the Gandalf error handler
gan_err_register() if the matrix is not positive definite (gan_symmat33_cholesky_s() aborts the program on error, so don't use it unless you're SURE your matrix is OK!). If you don't want to invoke the error package when factorisation is attempted on a non-positive-definite matrix, there is a set of routines which allows to instead to return the error condition as part of the result. For instance the code fragment

      Gan_SquMatrix33 sm33S, sm33L;
      int error_code;

      /* ... set up sm33S using e.g. gan_symmat33_fill_q() ... */
      if ( gan_symmat33_choleksy ( &sm33S, &sm33L, &error_code ) == NULL )
      {
         /* error found, act on it ... */
      }

      /* no error found */

attempts to factorise matrix $S$, and if an error is found, returns NULL, with the error condition returned in the error_code variable. For non-positive-definite matrices the error condition is GAN_ERROR_NOT_POSITIVE_DEFINITE.

Other factorisations are available in Gandalf. Singular value decomposition (SVD) can be used to decompose almost any matrix $A$ into factors as

$$A = UWV^{\top}$$

where $U$ and $V$ are orthogonal matrices and $W$ is diagonal. Currently Gandalf supports SVD for $3\times 3$ and $4\times 4$ matrices. To use the functions for $3\times 3$ matrices, include the header file

      #include <gandalf/linalg/3x3matrix_svd.h>

There are routines for SVD of a matrix or its transpose, as follows

      Gan_Matrix33 m33A, m33U, m33VT;
      Gan_Vector3 v3W;

      /* ... set up m33A using e.g. gan_mat33_fill_q() ... */
      gan_mat33_svd ( &m33A, &m33U, &v3W, &m33VT ); /* A = U*W*V^T, OR */
      gan_mat33T_svd ( &m33A, &m33U, &v3W, &m33VT ); /* A^T = U*W*V^T */

These routines return a Gan_Bool result, which is GAN_TRUE on success and GAN_FALSE on failure.

There are also a number of routines for computing the eigenvalues and eigenvectors of fixed size matrices. For $3\times 3$ matrices only there is a routine to compute the eigenvectors and complex eigenvalues of a $3\times 3$ matrix. To use the routine include the header file

      #include <gandalf/linalg/3x3matrix_eigen.h>

The matrix $A$ has ``left'' and ``right'' eigenvectors associated with the same eigenvalues $\lambda_i$, so that the equation

$$A {\bf v}_i = \lambda_i {\bf v}_i$$

defines the right eigenvectors, and

$${\bf u}_i^H A = \lambda_i {\bf u}_i^H$$

defines the left eigenvectors, where ${\bf u}_i^H$ denotes the conjugate transpose of vector ${\bf u}_i$. The computed eigenvectors are normalized to have Euclidean norm equal to one and largest component real. The Gandalf routine that implements this is built on the equivalent CLAPACK routine dgeev():

      Gan_Matrix33 m33A; /* matrix to be decomposed */
      Gan_Matrix33 m33UL, m33VR; /* left and right eigenvectors */
      Gan_Vector3 v3lr, v3li; /* real and imaginary parts of eigenvalues */

      /* ... set up m33A using e.g. gan_mat33_fill_q() ... */
      gan_mat33_eigen ( &m33A, &v3lr, &v3li, &m33UL, &m33VR );

The eigenvalues of symmetric matrices are guaranteed to be real. Routines are available for computing the eigenvalues and eigenvectors of $3\times 3$ and $4\times 4$ symmetric matrices, based on either the CLAPACK routine dspev() or the CCMath routine eigval(). For $3\times 3$ matrices the routine is declared in the header file

      #include <gandalf/linalg/3x3matrix_eigsym.h>

Here is an example using the routine

      Gan_SquMatrix33 sm33S; /* symmetric matrix to be decomposed */
      Gan_Matrix33 m33Z; /* (right) eigenvectors of A */
      Gan_Vector3 v3W; /* eigenvalues */

      /* ... set up sm33S using e.g. gan_symmat33_fill_q() ... */
      gan_symmat33_eigen ( &sm33S, &v3W, &m33Z );