include/a/linalg.h

   1 /*!
   2  @file linalg.h
   3  @brief linear algebra functions
   4 */
   5
   6 #ifndef LIBA_LINALG_H
   7 #define LIBA_LINALG_H
   8
   9 #include "a.h"
  10
  11 /*!
  12  @ingroup liba
  13  @addtogroup a_linalg linear algebra functions
  14  @{
  15 */
  16
  17 #if defined(__cplusplus)
  18 extern "C" {
  19 #endif /* __cplusplus */
  20
  21 /*!
  22  @brief transpose an n x n square matrix in-place.
  23  @param[in,out] A an n x n square matrix
  24  @param[in] n order of square matrix A
  25  */
  26 A_EXTERN void a_linalg_T1(a_real *A, a_uint n);
  27
  28 /*!
  29  @brief transpose a given m x n matrix A into an n x m matrix T.
  30  @param[out] T the output matrix where the transposed matrix T (n x m) will be stored.
  31  @param[in] A the input matrix A (m x n), stored in row-major order.
  32  @param[in] m rows in the input matrix A.
  33  @param[in] n columns in the input matrix A.
  34 */
  35 A_EXTERN void a_linalg_T2(a_real *__restrict T, a_real const *__restrict A, a_uint m, a_uint n);
  36
  37 /*!
  38  @brief compute the dot product of two vectors.
  39  @param[in] X points to the first vector.
  40  @param[in] Y points to the second vector.
  41  @param[in] n number of elements in each of the vectors X and Y.
  42  @return dot product of vectors X and Y.
  43 */
  44 A_EXTERN a_real a_linalg_dot(a_real const *X, a_real const *Y, a_size n);
  45
  46 /*!
  47  @brief multiply two matrices X and Y, storing the result in Z.
  48  \f{aligned}{
  49   \pmb Z_{rc}&=\pmb X_{rn}\pmb Y_{nc}
  50   \\&=
  51   \begin{bmatrix}
  52   x_{11} & \cdots & x_{1n} \\
  53   \vdots & \ddots & \vdots \\
  54   x_{r1} & \cdots & x_{rn} \\
  55   \end{bmatrix}
  56   \begin{bmatrix}
  57   y_{11} & \cdots & y_{1c} \\
  58   \vdots & \ddots & \vdots \\
  59   y_{n1} & \cdots & y_{nc} \\
  60   \end{bmatrix}
  61   \\&=
  62   \begin{bmatrix}
  63   (x_{11}y_{11}+\ldots+x_{1n}y_{n1}) & \cdots & (x_{11}y_{1c}+\ldots+x_{1n}y_{nc}) \\
  64   \vdots & \ddots & \vdots \\
  65   (x_{r1}y_{11}+\ldots+x_{rn}y_{n1}) & \cdots & (x_{r1}y_{1c}+\ldots+x_{rn}y_{nc}) \\
  66   \end{bmatrix}
  67  \f}
  68  @param[out] Z the output matrix where the result will be stored.
  69  @param[in] X the first input matrix.
  70  @param[in] Y the second input matrix.
  71  @param[in] row rows matrix Z and rows in matrix X.
  72  @param[in] c_r columns in matrix X and rows in matrix Y.
  73  @param[in] col columns in matrix Z and columns in matrix Y.
  74 */
  75 A_EXTERN void a_linalg_mulmm(a_real *__restrict Z, a_real const *__restrict X, a_real const *__restrict Y, a_uint row, a_uint c_r, a_uint col);
  76
  77 /*!
  78  @brief multiply the transpose of matrix X with matrix Y, storing the result in Z.
  79  \f{aligned}{
  80   \pmb Z_{rc}&=\pmb X_{nr}^{T}\pmb Y_{nc}
  81   \\&=
  82   \begin{bmatrix}
  83   x_{11} & \cdots & x_{1r} \\
  84   \vdots & \ddots & \vdots \\
  85   x_{n1} & \cdots & x_{nr} \\
  86   \end{bmatrix}^T
  87   \begin{bmatrix}
  88   y_{11} & \cdots & y_{1c} \\
  89   \vdots & \ddots & \vdots \\
  90   y_{n1} & \cdots & y_{nc} \\
  91   \end{bmatrix}
  92   \\&=
  93   \begin{bmatrix}
  94   (x_{11}y_{11}+\ldots+x_{n1}y_{n1}) & \cdots & (x_{11}y_{1c}+\ldots+x_{n1}y_{nc}) \\
  95   \vdots & \ddots & \vdots \\
  96   (x_{1r}y_{11}+\ldots+x_{nr}y_{n1}) & \cdots & (x_{1r}y_{1c}+\ldots+x_{nr}y_{nc}) \\
  97   \end{bmatrix}
  98   \\&=
  99   \begin{bmatrix}
 100   x_{11}y_{11} & \cdots & x_{11}y_{1c} \\
 101   \vdots & \ddots & \vdots \\
 102   x_{1r}y_{11} & \cdots & x_{1r}y_{1c} \\
 103   \end{bmatrix}+\cdots+
 104   \begin{bmatrix}
 105   x_{n1}y_{n1} & \cdots & x_{n1}y_{nc} \\
 106   \vdots & \ddots & \vdots \\
 107   x_{nr}y_{n1} & \cdots & x_{nr}y_{nc} \\
 108   \end{bmatrix}
 109  \f}
 110  @param[out] Z the output matrix where the result will be stored.
 111  @param[in] X the first input matrix that will be transposed during multiplication.
 112  @param[in] Y the second input matrix.
 113  @param[in] c_r rows in matrix X and rows in matrix Y.
 114  @param[in] row rows in matrix Z and columns in matrix X.
 115  @param[in] col columns in matrix Z and columns in matrix Y.
 116 */
 117 A_EXTERN void a_linalg_mulTm(a_real *__restrict Z, a_real const *__restrict X, a_real const *__restrict Y, a_uint c_r, a_uint row, a_uint col);
 118
 119 /*!
 120  @brief multiply matrix X with the transpose of matrix Y, storing the result in Z.
 121  \f{aligned}{
 122   \pmb Z_{rc}&=\pmb X_{rn}\pmb Y_{cn}^T
 123   \\&=
 124   \begin{bmatrix}
 125   x_{11} & \cdots & x_{1n} \\
 126   \vdots & \ddots & \vdots \\
 127   x_{r1} & \cdots & x_{rn} \\
 128   \end{bmatrix}
 129   \begin{bmatrix}
 130   y_{11} & \cdots & y_{1n} \\
 131   \vdots & \ddots & \vdots \\
 132   y_{c1} & \cdots & y_{cn} \\
 133   \end{bmatrix}^T
 134   \\&=
 135   \begin{bmatrix}
 136   (x_{11}y_{11}+\ldots+x_{1n}y_{1n}) & \cdots & (x_{11}y_{c1}+\ldots+x_{1n}y_{cn}) \\
 137   \vdots & \ddots & \vdots \\
 138   (x_{r1}y_{11}+\ldots+x_{rn}y_{1n}) & \cdots & (x_{r1}y_{c1}+\ldots+x_{rn}y_{cn}) \\
 139   \end{bmatrix}
 140  \f}
 141  @param[out] Z the output matrix where the result will be stored.
 142  @param[in] X the first input matrix.
 143  @param[in] Y the second input matrix that will be transposed during multiplication.
 144  @param[in] row rows matrix Z and rows in matrix X.
 145  @param[in] col columns in matrix Z and rows in matrix Y.
 146  @param[in] c_r columns in matrix X and columns in matrix Y.
 147 */
 148 A_EXTERN void a_linalg_mulmT(a_real *__restrict Z, a_real const *__restrict X, a_real const *__restrict Y, a_uint row, a_uint col, a_uint c_r);
 149
 150 /*!
 151  @brief multiply the transpose of matrix X with the transpose of matrix Y, storing the result in Z.
 152  \f{aligned}{
 153   \pmb Z_{rc}&=\pmb X_{nr}^T\pmb Y_{cn}^T
 154   \\&=
 155   \begin{bmatrix}
 156   x_{11} & \cdots & x_{1r} \\
 157   \vdots & \ddots & \vdots \\
 158   x_{n1} & \cdots & x_{nr} \\
 159   \end{bmatrix}^T
 160   \begin{bmatrix}
 161   y_{11} & \cdots & y_{1n} \\
 162   \vdots & \ddots & \vdots \\
 163   y_{c1} & \cdots & y_{cn} \\
 164   \end{bmatrix}^T
 165   \\&=
 166   \begin{bmatrix}
 167   (x_{11}y_{11}+\ldots+x_{n1}y_{1n}) & \cdots & (x_{11}y_{c1}+\ldots+x_{n1}y_{cn}) \\
 168   \vdots & \ddots & \vdots \\
 169   (x_{1r}y_{11}+\ldots+x_{nr}y_{1n}) & \cdots & (x_{1r}y_{c1}+\ldots+x_{nr}y_{cn}) \\
 170   \end{bmatrix}
 171  \f}
 172  @param[out] Z the output matrix where the result will be stored.
 173  @param[in] X the first input matrix that will be transposed during multiplication.
 174  @param[in] Y the second input matrix that will be transposed during multiplication.
 175  @param[in] row rows matrix Z and columns in matrix X.
 176  @param[in] c_r rows in matrix X and columns in matrix Y.
 177  @param[in] col columns in matrix Z and rows in matrix Y.
 178 */
 179 A_EXTERN void a_linalg_mulTT(a_real *__restrict Z, a_real const *__restrict X, a_real const *__restrict Y, a_uint row, a_uint c_r, a_uint col);
 180
 181 /*!
 182  @brief compute LU decomposition of a square matrix with partial pivoting.
 183  @details This function performs an LU decomposition on the given square matrix A,
 184  where L is a lower triangular matrix, and U is an upper triangular matrix.
 185  Partial pivoting is used to improve numerical stability during the decomposition process.
 186  The result is stored in the original matrix A, with L stored below, and U stored in the diagonal and above.
 187  Additionally, it calculates a permutation matrix P that records the row exchanges made during partial pivoting,
 188  and determines the sign of the permutation (which can be used to find the determinant's sign).
 189  @param[in,out] A an n x n square matrix.
 190   on input, contains the matrix to decompose. on output, contains the L and U matrices.
 191  @param[in] n the order of the square matrix A (number of rows and columns).
 192  @param[out] p the row permutation indices after partial pivoting.
 193  @param[out] sign store the sign of the permutation (+1 or -1).
 194  @return 0 on success, or a non-zero error code if the decomposition fails.
 195   @retval -1 on failure, A is a singular matrix.
 196  */
 197 A_EXTERN int a_linalg_plu(a_real *A, a_uint n, a_uint *p, int *sign);
 198
 199 /*!
 200  @brief construct the permutation matrix P from a permutation vector p.
 201  @param[in] p the row permutation indices after partial pivoting.
 202  @param[in] n the order of the square matrix that was decomposed.
 203  @param[out] P the output matrix where the permutation matrix will be stored.
 204 */
 205 A_EXTERN void a_linalg_plu_get_P(a_uint const *p, a_uint n, a_real *P);
 206
 207 /*!
 208  @brief extract the lower triangular matrix L from matrix A.
 209  @param[in] A the matrix containing L and U in a compact form after LU decomposition.
 210  @param[in] n the order of the square matrix that was decomposed.
 211  @param[out] L the output matrix where the lower triangular matrix will be stored.
 212 */
 213 A_EXTERN void a_linalg_plu_get_L(a_real const *A, a_uint n, a_real *L);
 214
 215 /*!
 216  @brief extract the upper triangular matrix U from matrix A.
 217  @param[in] A the matrix containing L and U in a compact form after LU decomposition.
 218  @param[in] n the order of the square matrix that was decomposed.
 219  @param[out] U the output matrix where the upper triangular matrix will be stored.
 220 */
 221 A_EXTERN void a_linalg_plu_get_U(a_real const *A, a_uint n, a_real *U);
 222
 223 /*!
 224  @brief apply the permutation P to the vector b, producing Pb.
 225  @param[in] p the row permutation indices after partial pivoting.
 226  @param[in] n the order of the square matrix that was decomposed.
 227  @param[in] b the input vector of size n that will be permuted.
 228  @param[out] Pb the output vector where the permuted result will be stored.
 229 */
 230 A_EXTERN void a_linalg_plu_apply(a_uint const *p, a_uint n, a_real const *b, a_real *Pb);
 231
 232 /*!
 233  @brief solve the lower triangular system Ly = Pb for y.
 234  @param[in] L the lower triangular matrix L, stored in row-major order.
 235  @param[in] n the order of the square matrix L (number of rows and columns).
 236  @param[in,out] y on input, contains the permuted vector Pb. on output, contains the solution vector y.
 237 */
 238 A_EXTERN void a_linalg_plu_lower(a_real const *L, a_uint n, a_real *y);
 239
 240 /*!
 241  @brief solve the upper triangular system Ux = y for x.
 242  @param[in] U the upper triangular matrix U, stored in row-major order.
 243  @param[in] n the order of the square matrix U (number of rows and columns).
 244  @param[in,out] x on input, contains the vector y. on output, contains the solution vector x.
 245 */
 246 A_EXTERN void a_linalg_plu_upper(a_real const *U, a_uint n, a_real *x);
 247
 248 /*!
 249  @brief solve the linear system Ax = b using LU decomposition with partial pivoting.
 250  @param[in] A the matrix containing L and U in a compact form after LU decomposition.
 251  @param[in] n the order of the square matrix A (number of rows and columns).
 252  @param[in] p the permutation indices obtained during LU decomposition.
 253  @param[in] b the input vector b of the linear system.
 254  @param[out] x the output vector x where the solution will be stored.
 255 */
 256 A_EXTERN void a_linalg_plu_solve(a_real const *A, a_uint n, a_uint const *p, a_real const *b, a_real *x);
 257
 258 /*!
 259  @brief compute the inverse of a matrix using its LU decomposition and permutation matrix.
 260  @param[in] A the matrix containing L and U in a compact form after LU decomposition.
 261  @param[in] n the order of the square matrix A (number of rows and columns).
 262  @param[in] p the permutation indices obtained during LU decomposition.
 263  @param[in] b a pre-allocated temporary buffer of size n for intermediate computations.
 264  @param[out] I the output matrix where the inverse of A will be stored.
 265 */
 266 A_EXTERN void a_linalg_plu_inv(a_real const *A, a_uint n, a_uint const *p, a_real *b, a_real *I);
 267
 268 /*!
 269  @brief compute the determinant of a matrix using its LU decomposition.
 270  @param[in] A the matrix containing L and U in a compact form after LU decomposition.
 271  @param[in] n the order of the square matrix A (number of rows and columns).
 272  @param[in] sign the sign of the permutation matrix P (+1 or -1).
 273  @return the determinant of matrix A.
 274 */
 275 A_EXTERN a_real a_linalg_plu_det(a_real const *A, a_uint n, int sign);
 276
 277 /*!
 278  @brief compute the natural logarithm of the absolute value of the determinant of a matrix using its LU decomposition.
 279  @param[in] A the matrix containing L and U in a compact form after LU decomposition.
 280  @param[in] n the order of the square matrix A (number of rows and columns).
 281  @return the natural logarithm of the absolute value of the determinant.
 282 */
 283 A_EXTERN a_real a_linalg_plu_lndet(a_real const *A, a_uint n);
 284
 285 /*!
 286  @brief compute the sign of the determinant of a matrix using its LU decomposition.
 287  @param[in] A the matrix containing L and U in a compact form after LU decomposition.
 288  @param[in] n the order of the square matrix A (number of rows and columns).
 289  @param[in] sign the sign of the permutation matrix P (+1 or -1).
 290  @return the sign of the determinant: -1, 0, +1.
 291 */
 292 A_EXTERN int a_linalg_plu_sgndet(a_real const *A, a_uint n, int sign);
 293
 294 /*!
 295  @brief compute Cholesky decomposition of a symmetric positive-definite matrix.
 296  @param[in,out] A an n x n square matrix.
 297   on input, contains the matrix to decompose. on output, contains the L matrix.
 298  @param[in] n the order of the square matrix A (number of rows and columns).
 299  @return 0 on success, or a non-zero error code if the decomposition fails.
 300   @retval -1 on failure, A is a singular matrix.
 301 */
 302 A_EXTERN int a_linalg_cho(a_real *A, a_uint n);
 303
 304 /*!
 305  @brief extract the lower triangular matrix L from matrix A.
 306  @param[in] A the matrix containing L form after Cholesky decomposition.
 307  @param[in] n the order of the square matrix that was decomposed.
 308  @param[out] L the output matrix where the lower triangular matrix will be stored.
 309 */
 310 A_EXTERN void a_linalg_cho_get_L(a_real const *A, a_uint n, a_real *L);
 311
 312 /*!
 313  @brief solve the lower triangular system Ly = b for y.
 314  @param[in] L the lower triangular matrix L, stored in row-major order.
 315  @param[in] n the order of the square matrix L (number of rows and columns).
 316  @param[in,out] y on input, contains the vector b. on output, contains the solution vector y.
 317 */
 318 A_EXTERN void a_linalg_cho_lower(a_real const *L, a_uint n, a_real *y);
 319
 320 /*!
 321  @brief solve the upper triangular system L^T x = y for x.
 322  @param[in] L the lower triangular matrix L, stored in row-major order.
 323  @param[in] n the order of the square matrix L (number of rows and columns).
 324  @param[in,out] x on input, contains the vector y. on output, contains the solution vector x.
 325 */
 326 A_EXTERN void a_linalg_cho_upper(a_real const *L, a_uint n, a_real *x);
 327
 328 /*!
 329  @brief solve the linear system Ax = b using the Cholesky factorization A = LL^T.
 330  @param[in] A the matrix containing L form after Cholesky decomposition.
 331  @param[in] n the order of the square matrix A (number of rows and columns).
 332  @param[in,out] x on input, contains the vector b. on output, contains the solution vector x.
 333 */
 334 A_EXTERN void a_linalg_cho_solve(a_real const *A, a_uint n, a_real *x);
 335
 336 /*!
 337  @brief compute the inverse of a matrix using its Cholesky factorization A = LL^T.
 338  @param[in] A the matrix containing L form after Cholesky decomposition.
 339  @param[in] n the order of the square matrix A (number of rows and columns).
 340  @param[in] b a pre-allocated temporary buffer of size n for intermediate computations.
 341  @param[out] I the output matrix where the inverse of A will be stored.
 342 */
 343 A_EXTERN void a_linalg_cho_inv(a_real const *A, a_uint n, a_real *b, a_real *I);
 344
 345 #if defined(__cplusplus)
 346 } /* extern "C" */
 347 #endif /* __cplusplus */
 348
 349 /*! @} a_linalg */
 350
 351 #endif /* a/linalg.h */