arb_mat.h – matrices over the real numbers¶
An arb_mat_t represents a dense matrix over the real numbers,
implemented as an array of entries of type arb_struct.
The dimension (number of rows and columns) of a matrix is fixed at initialization, and the user must ensure that inputs and outputs to an operation have compatible dimensions. The number of rows or columns in a matrix can be zero.
Types, macros and constants¶
-
arb_mat_struct¶
-
arb_mat_t¶ Contains a pointer to a flat array of the entries (entries), an array of pointers to the start of each row (rows), and the number of rows (r) and columns (c).
An arb_mat_t is defined as an array of length one of type arb_mat_struct, permitting an arb_mat_t to be passed by reference.
-
arb_mat_entry(mat, i, j)¶ Macro giving a pointer to the entry at row i and column j.
-
arb_mat_nrows(mat)¶ Returns the number of rows of the matrix.
-
arb_mat_ncols(mat)¶ Returns the number of columns of the matrix.
Memory management¶
Conversions¶
-
void
arb_mat_set_fmpz_mat(arb_mat_t dest, const fmpz_mat_t src)¶
-
void
arb_mat_set_round_fmpz_mat(arb_mat_t dest, const fmpz_mat_t src, slong prec)¶
-
void
arb_mat_set_fmpq_mat(arb_mat_t dest, const fmpq_mat_t src, slong prec)¶ Sets dest to src. The operands must have identical dimensions.
Random generation¶
Input and output¶
Comparisons¶
-
int
arb_mat_equal(const arb_mat_t mat1, const arb_mat_t mat2)¶ Returns nonzero iff the matrices have the same dimensions and identical entries.
-
int
arb_mat_overlaps(const arb_mat_t mat1, const arb_mat_t mat2)¶ Returns nonzero iff the matrices have the same dimensions and each entry in mat1 overlaps with the corresponding entry in mat2.
-
int
arb_mat_contains_fmpz_mat(const arb_mat_t mat1, const fmpz_mat_t mat2)¶
-
int
arb_mat_contains_fmpq_mat(const arb_mat_t mat1, const fmpq_mat_t mat2)¶ Returns nonzero iff the matrices have the same dimensions and each entry in mat2 is contained in the corresponding entry in mat1.
-
int
arb_mat_eq(const arb_mat_t mat1, const arb_mat_t mat2)¶ Returns nonzero iff mat1 and mat2 certainly represent the same matrix.
-
int
arb_mat_ne(const arb_mat_t mat1, const arb_mat_t mat2)¶ Returns nonzero iff mat1 and mat2 certainly do not represent the same matrix.
Special matrices¶
Transpose¶
Norms¶
-
void
arb_mat_bound_inf_norm(mag_t b, const arb_mat_t A)¶ Sets b to an upper bound for the infinity norm (i.e. the largest absolute value row sum) of A.
Arithmetic¶
-
void
arb_mat_neg(arb_mat_t dest, const arb_mat_t src)¶ Sets dest to the exact negation of src. The operands must have the same dimensions.
-
void
arb_mat_add(arb_mat_t res, const arb_mat_t mat1, const arb_mat_t mat2, slong prec)¶ Sets res to the sum of mat1 and mat2. The operands must have the same dimensions.
-
void
arb_mat_sub(arb_mat_t res, const arb_mat_t mat1, const arb_mat_t mat2, slong prec)¶ Sets res to the difference of mat1 and mat2. The operands must have the same dimensions.
-
void
arb_mat_mul(arb_mat_t res, const arb_mat_t mat1, const arb_mat_t mat2, slong prec)¶ Sets res to the matrix product of mat1 and mat2. The operands must have compatible dimensions for matrix multiplication.
The threaded version splits the computation over the number of threads returned by flint_get_num_threads(). The default version automatically calls the threaded version if the matrices are sufficiently large and more than one thread can be used.
-
void
arb_mat_mul_entrywise(arb_mat_t C, const arb_mat_t A, const arb_mat_t B, slong prec)¶ Sets C to the entrywise product of A and B. The operands must have the same dimensions.
Scalar arithmetic¶
-
void
arb_mat_scalar_mul_2exp_si(arb_mat_t B, const arb_mat_t A, slong c)¶ Sets B to A multiplied by \(2^c\).
-
void
arb_mat_scalar_addmul_arb(arb_mat_t B, const arb_mat_t A, const arb_t c, slong prec)¶ Sets B to \(B + A \times c\).
Gaussian elimination and solving¶
-
int
arb_mat_lu(slong * perm, arb_mat_t LU, const arb_mat_t A, slong prec)¶ Given an \(n \times n\) matrix \(A\), computes an LU decomposition \(PLU = A\) using Gaussian elimination with partial pivoting. The input and output matrices can be the same, performing the decomposition in-place.
Entry \(i\) in the permutation vector perm is set to the row index in the input matrix corresponding to row \(i\) in the output matrix.
The algorithm succeeds and returns nonzero if it can find \(n\) invertible (i.e. not containing zero) pivot entries. This guarantees that the matrix is invertible.
The algorithm fails and returns zero, leaving the entries in \(P\) and \(LU\) undefined, if it cannot find \(n\) invertible pivot elements. In this case, either the matrix is singular, the input matrix was computed to insufficient precision, or the LU decomposition was attempted at insufficient precision.
-
void
arb_mat_solve_lu_precomp(arb_mat_t X, const slong * perm, const arb_mat_t LU, const arb_mat_t B, slong prec)¶ Solves \(AX = B\) given the precomputed nonsingular LU decomposition \(A = PLU\). The matrices \(X\) and \(B\) are allowed to be aliased with each other, but \(X\) is not allowed to be aliased with \(LU\).
-
int
arb_mat_solve(arb_mat_t X, const arb_mat_t A, const arb_mat_t B, slong prec)¶ Solves \(AX = B\) where \(A\) is a nonsingular \(n \times n\) matrix and \(X\) and \(B\) are \(n \times m\) matrices, using LU decomposition.
If \(m > 0\) and \(A\) cannot be inverted numerically (indicating either that \(A\) is singular or that the precision is insufficient), the values in the output matrix are left undefined and zero is returned. A nonzero return value guarantees that \(A\) is invertible and that the exact solution matrix is contained in the output.
-
int
arb_mat_inv(arb_mat_t X, const arb_mat_t A, slong prec)¶ Sets \(X = A^{-1}\) where \(A\) is a square matrix, computed by solving the system \(AX = I\).
If \(A\) cannot be inverted numerically (indicating either that \(A\) is singular or that the precision is insufficient), the values in the output matrix are left undefined and zero is returned. A nonzero return value guarantees that the matrix is invertible and that the exact inverse is contained in the output.
-
void
arb_mat_det(arb_t det, const arb_mat_t A, slong prec)¶ Computes the determinant of the matrix, using Gaussian elimination with partial pivoting. If at some point an invertible pivot element cannot be found, the elimination is stopped and the magnitude of the determinant of the remaining submatrix is bounded using Hadamard’s inequality.
Cholesky decomposition and solving¶
-
int
arb_mat_cho(arb_mat_t L, const arb_mat_t A, slong prec)¶ Computes the Cholesky decomposition of A, returning nonzero iff the symmetric matrix defined by the lower triangular part of A is certainly positive definite.
If a nonzero value is returned, then L is set to the lower triangular matrix such that \(A = L * L^T\).
If zero is returned, then either the matrix is not symmetric positive definite, the input matrix was computed to insufficient precision, or the decomposition was attempted at insufficient precision.
The underscore method computes L from A in-place, leaving the strict upper triangular region undefined.
-
void
arb_mat_solve_cho_precomp(arb_mat_t X, const arb_mat_t L, const arb_mat_t B, slong prec)¶ Solves \(AX = B\) given the precomputed Cholesky decomposition \(A = L L^T\). The matrices X and B are allowed to be aliased with each other, but X is not allowed to be aliased with L.
-
int
arb_mat_spd_solve(arb_mat_t X, const arb_mat_t A, const arb_mat_t B, slong prec)¶ Solves \(AX = B\) where A is a symmetric positive definite matrix and X and B are \(n \times m\) matrices, using Cholesky decomposition.
If \(m > 0\) and A cannot be factored using Cholesky decomposition (indicating either that A is not symmetric positive definite or that the precision is insufficient), the values in the output matrix are left undefined and zero is returned. A nonzero return value guarantees that the symmetric matrix defined through the lower triangular part of A is invertible and that the exact solution matrix is contained in the output.
-
void
arb_mat_inv_cho_precomp(arb_mat_t X, const arb_mat_t L, slong prec)¶ Sets \(X = A^{-1}\) where \(A\) is a symmetric positive definite matrix whose Cholesky decomposition L has been computed with
arb_mat_cho(). The inverse is calculated using the method of [Kri2013] which is more efficient than solving \(AX = I\) witharb_mat_solve_cho_precomp().
-
int
arb_mat_spd_inv(arb_mat_t X, const arb_mat_t A, slong prec)¶ Sets \(X = A^{-1}\) where A is a symmetric positive definite matrix. It is calculated using the method of [Kri2013] which computes fewer intermediate results than solving \(AX = I\) with
arb_mat_spd_solve().If A cannot be factored using Cholesky decomposition (indicating either that A is not symmetric positive definite or that the precision is insufficient), the values in the output matrix are left undefined and zero is returned. A nonzero return value guarantees that the symmetric matrix defined through the lower triangular part of A is invertible and that the exact inverse is contained in the output.
-
int
arb_mat_ldl(arb_mat_t res, const arb_mat_t A, slong prec)¶ Computes the \(LDL^T\) decomposition of A, returning nonzero iff the symmetric matrix defined by the lower triangular part of A is certainly positive definite.
If a nonzero value is returned, then res is set to a lower triangular matrix that encodes the \(L * D * L^T\) decomposition of A. In particular, \(L\) is a lower triangular matrix with ones on its diagonal and whose strictly lower triangular region is the same as that of res. \(D\) is a diagonal matrix with the same diagonal as that of res.
If zero is returned, then either the matrix is not symmetric positive definite, the input matrix was computed to insufficient precision, or the decomposition was attempted at insufficient precision.
The underscore methods compute res from A in-place, leaving the strict upper triangular region undefined. The default method uses algorithm 4.1.2 from [GVL1996].
-
void
arb_mat_solve_ldl_precomp(arb_mat_t X, const arb_mat_t L, const arb_mat_t B, slong prec)¶ Solves \(AX = B\) given the precomputed \(A = LDL^T\) decomposition encoded by L. The matrices X and B are allowed to be aliased with each other, but X is not allowed to be aliased with L.
-
void
arb_mat_inv_ldl_precomp(arb_mat_t X, const arb_mat_t L, slong prec)¶ Sets \(X = A^{-1}\) where \(A\) is a symmetric positive definite matrix whose \(LDL^T\) decomposition encoded by L has been computed with
arb_mat_ldl(). The inverse is calculated using the method of [Kri2013] which is more efficient than solving \(AX = I\) witharb_mat_solve_ldl_precomp().
Characteristic polynomial¶
-
void
arb_mat_charpoly(arb_poly_t cp, const arb_mat_t mat, slong prec)¶ Sets cp to the characteristic polynomial of mat which must be a square matrix. If the matrix has n rows, the underscore method requires space for \(n + 1\) output coefficients. Employs a division-free algorithm using \(O(n^4)\) operations.
Special functions¶
-
void
arb_mat_exp_taylor_sum(arb_mat_t S, const arb_mat_t A, slong N, slong prec)¶ Sets S to the truncated exponential Taylor series \(S = \sum_{k=0}^{N-1} A^k / k!\). Uses rectangular splitting to compute the sum using \(O(\sqrt{N})\) matrix multiplications. The recurrence relation for factorials is used to get scalars that are small integers instead of full factorials. As in [Joh2014b], all divisions are postponed to the end by computing partial factorials of length \(O(\sqrt{N})\). The scalars could be reduced by doing more divisions, but this appears to be slower in most cases.
-
void
arb_mat_exp(arb_mat_t B, const arb_mat_t A, slong prec)¶ Sets B to the exponential of the matrix A, defined by the Taylor series
\[\exp(A) = \sum_{k=0}^{\infty} \frac{A^k}{k!}.\]The function is evaluated as \(\exp(A/2^r)^{2^r}\), where \(r\) is chosen to give rapid convergence.
The elementwise error when truncating the Taylor series after N terms is bounded by the error in the infinity norm, for which we have
\[\left\|\exp(2^{-r}A) - \sum_{k=0}^{N-1} \frac{\left(2^{-r} A\right)^k}{k!} \right\|_{\infty} = \left\|\sum_{k=N}^{\infty} \frac{\left(2^{-r} A\right)^k}{k!}\right\|_{\infty} \le \sum_{k=N}^{\infty} \frac{(2^{-r} \|A\|_{\infty})^k}{k!}.\]We bound the sum on the right using
mag_exp_tail(). Truncation error is not added to entries whose values are determined by the sparsity structure of \(A\).
Sparsity structure¶
-
void
arb_mat_entrywise_is_zero(fmpz_mat_t dest, const arb_mat_t src)¶ Sets each entry of dest to indicate whether the corresponding entry of src is certainly zero. If the entry of src at row \(i\) and column \(j\) is zero according to
arb_is_zero()then the entry of dest at that row and column is set to one, otherwise that entry of dest is set to zero.
-
void
arb_mat_entrywise_not_is_zero(fmpz_mat_t dest, const arb_mat_t src)¶ Sets each entry of dest to indicate whether the corresponding entry of src is not certainly zero. This the complement of
arb_mat_entrywise_is_zero().
-
slong
arb_mat_count_is_zero(const arb_mat_t mat)¶ Returns the number of entries of mat that are certainly zero according to
arb_is_zero().
