gr_mat.h

Initialises mat as a zero matrix with the given numbers of rows and columns, which may be zero.

Initialises res as a copy of mat.

Clears the entries and frees the memory of mat. A window must be released with gr_mat_window_clear instead.

Swaps mat1 and mat2 in constant time by exchanging the structures.

Swaps the entries of two matrices of the same shape one by one, leaving the structures in place. This is the way to exchange the contents of windows.

Return a pointer to the entry in row i and column j. The indices must be in bounds.

Return whether mat has no entries (zero rows or zero columns), respectively as many rows as columns. These predicates depend only on the shape and always return T_TRUE or T_FALSE.

gr_mat_window_init makes window the submatrix of mat formed by rows r_1 \le i < r_2 and columns c_1 \le j < c_2, sharing the entries of mat (writing to the window modifies mat). The indices must be within bounds. gr_mat_window_clear releases a window without touching the entries.

Write mat to a stream or to standard output, as a bracketed list of rows with one row per line.

Sets mat to the matrix described by s, of the form [[e11, ..., e1n], ..., [em1, ..., emn]], where each row is parsed as by gr_vec_set_str and each entry by gr_set_str. Whitespace around brackets and separators, including the line breaks written by gr_mat_write, is ignored.

Set res to the zero matrix, respectively to the matrix with ones on the main diagonal and zeros elsewhere (the identity if square).

Set res to a copy of mat, converting the entries from integers, rationals or the ring mat_ctx into R. A conversion that is not possible gives GR_DOMAIN or GR_UNABLE. The shapes must agree.

Set res to the matrix with c on the main diagonal and zeros elsewhere.

Returns whether the matrices are equal: T_FALSE if the shapes differ, otherwise the conjunction of the entrywise comparisons.

Return whether mat is zero, has 1 or -1 on the main diagonal and zeros elsewhere, or is a scalar matrix (a diagonal matrix with equal diagonal entries).

Return whether all entries off the main diagonal, below it, above it, or below the first subdiagonal (upper Hessenberg form) are zero. The matrix need not be square.

Set res to (A \mid B), respectively to A stacked above B. Return GR_DOMAIN if the numbers of rows (columns) of the inputs, or the shape of res, do not match.

Sets res to the transpose of mat. Aliasing is allowed for square matrices, in which case the entries are swapped in place.

Swap rows (columns) r and s of mat, and the entries r and s of the permutation perm if it is not NULL. The non-underscore versions return GR_DOMAIN for indices out of bounds; the underscore versions do not check them.

Reverse the order of the rows (columns): swap row i and row r - 1 - i for 0 \le i < r/2, and apply the same permutation to perm if it is not NULL.

Permute the rows (columns) by perm_act, or by its inverse for the inv versions, and apply the same permutation to perm_store if it is not NULL. gr_mat_permute_rows moves row perm_act[i] to row i (and perm_store[perm_act[i]] to perm_store[i]); the inv version moves row i to row perm_act[i], without computing the inverse permutation. perm_store and perm_act may be the same array.

Moves row i to position \mathit{new\_i}, shifting the rows in between by one: with i = 4 and \mathit{new\_i} = 7 the rows A_4, A_5, A_6, A_7 become A_5, A_6, A_7, A_4, and with i = 7 and \mathit{new\_i} = 4 they become A_7, A_4, A_5, A_6. Returns GR_DOMAIN if an index is out of bounds and otherwise succeeds.

Set each entry of res to f applied to the corresponding entry of mat, to the entries of mat1 and mat2, or to the entry of mat and the scalar c. Return GR_DOMAIN if the shapes do not match. These are declared only when gr.h is included before gr_mat.h.

Return whether the predicate f holds for all entries, for some entry, or for all pairs of corresponding entries (T_FALSE if the shapes differ).

Set res to -A, A + B or A - B.

Set res to A + xI or xI + A, adding x to the main diagonal.

Set res to A - xI or xI - A.

Set res to Ax or xA (the order matters over noncommutative rings).

Set res to A / x, dividing each entry with gr_div (GR_DOMAIN if a quotient does not exist in R).

Set res to \mathit{res} + Ac, respectively \mathit{res} - Ac.

Set res to AD or DA, where D is the diagonal matrix with the entries of the vector D on its diagonal: column j (row i) of A is multiplied by D_j (D_i).

Set res to the product AB. gr_mat_mul calls the ring method GR_METHOD_MAT_MUL; the default gr_mat_mul_generic is classical multiplication. Most rings with a dedicated matrix module override the method. Returns GR_DOMAIN if the shapes are incompatible: gr_mat_mul itself does not check them, so every ring method must (all methods in FLINT do). See Multiplying matrices over generic rings.

Computes the m \times n by n \times p product with the classical algorithm, as mp dot products of length n (_gr_vec_dot), with the right factor transposed shallowly so that its columns are contiguous.

One level of Strassen's algorithm [Str1969] with the 7-multiplication, 15-addition sequence of Winograd [Win1971] in the order given by Bodrato [Bodrato2010], which saves further operations when squaring. The seven products of half size are computed with gr_mat_mul, so the recursion continues only if the ring method dispatches back to this function above a cutoff; then the complexity is O(n^{\log_2 7}) = O(n^{2.81}). See Strassen–Winograd multiplication.

Compute AB with about half the multiplications of the classical algorithm over a commutative ring, at the cost of more additions, which pays off when multiplications are much more expensive than additions (large integers, polynomials). Waksman's algorithm [Waksman1970] (refining Winograd's inner product trick [Win1968]) needs exact division by 2, implemented as gr_mul_2exp_si by -1. Rosowski's algorithm [Rosowski2023] needs no division and uses fewer multiplications than Waksman's when the inner dimension is odd.

Sets res to A^2 (currently gr_mat_mul of A by itself).

Set res to A^e for a square matrix, by repeated squaring in the matrix ring (gr_generic_pow_ui, gr_generic_pow_fmpz); a negative exponent requires an invertible matrix. A^0 is the identity. Returns GR_DOMAIN if the matrices are not square of the same size.

Set res to f(A) for a scalar polynomial f and a square matrix A, by rectangular splitting [PS1973]: about 2\sqrt{\mathit{len}} matrix products and O(\mathit{len}) scalar multiplications.

Set res to the max norm \max_{i,j} |a_{ij}|, the 1-norm (largest absolute column sum) \max_j \sum_i |a_{ij}|, the \infty-norm (largest absolute row sum) \max_i \sum_j |a_{ij}|, or the Frobenius norm (\sum_{i,j} |a_{ij}|^2)^{1/2}. The ring must provide absolute values, comparisons and (for the Frobenius norm) square roots; the result is an element of R.

Solve LX = B for a lower triangular L, or UX = B for an upper triangular U. Only the relevant triangle is read. If unit is set, the diagonal is taken to be all ones and is not read (it may hold other data, such as the U factor of a packed LU decomposition). Otherwise the diagonal entries are inverted if possible and divided by otherwise; a division that fails gives GR_DOMAIN, even if the system has a solution.

Look for an entry that is provably nonzero in column column, in rows \mathit{start\_row} \le i < \mathit{end\_row}. On success, set *pivot_row and return GR_SUCCESS. Return GR_DOMAIN if all these entries are provably zero, and GR_UNABLE if none is provably nonzero but some cannot be decided.

Compute a generalised LU decomposition PA = LU of an m \times n matrix A with row permutation P, writing the rank to rank, the permutation to the array P of length m, and L and U packed into LU. See LU decompositions and linear systems.

Computes a fraction-free LU decomposition of A by Bareiss' algorithm, over an integral domain with exact division, writing the final pivot (a minor of A) to den. See LU decompositions and linear systems.

Solve AX = B for a square matrix A. Return GR_DOMAIN if A is singular, even if the system has solutions, and if X does not have entries in R. See LU decompositions and linear systems.

Solve AX = B given an LU or fraction-free LU decomposition of A with full rank, as computed by gr_mat_lu or gr_mat_fflu. The data are assumed to be valid, which is not checked. The fflu version computes \mathit{den} \cdot X, where den is the denominator returned by gr_mat_fflu.

Solve AX = B over the fraction field of an integral domain R: compute X and \mathit{den} \in R such that A X = \mathit{den} \cdot B, with \mathit{den} = \pm \det A from fraction-free elimination (\mathit{den} = 1 for 0 \times 0 matrices). Return GR_DOMAIN if A is singular. See LU decompositions and linear systems.

Solves AX = B for an arbitrary (possibly rectangular or singular) A over a field. If the system has several solutions, an arbitrary one is returned; GR_DOMAIN means that it has none. See LU decompositions and linear systems.

Set res to the determinant of a square matrix (GR_DOMAIN otherwise), with the division-free Berkowitz algorithm, cofactor expansion, LU decomposition or fraction-free LU decomposition. See Determinants and permanents.

Set res to the permanent \operatorname{perm}(A) = \sum_{\sigma} \prod_i a_{i,\sigma(i)} of a square matrix, the determinant without signs. Computing it is #P-hard [Val1979], and all algorithms take exponential time. See Permanents.

Sets res to the sum of the diagonal entries of a square matrix. Returns GR_DOMAIN if the matrix is not square.

Set rank to the rank of A, by LU decomposition or by fraction-free LU decomposition. gr_mat_rank uses LU over finite fields (and finite rings pretending to be fields), fraction-free LU over other integral domains, which keeps the entries small over \mathbb{Z} and \mathbb{Q}, and returns GR_DOMAIN for rings that are known not to be integral domains, where the rank is not defined here. The lu and fflu versions do not check the ring; any failure of the decomposition (an impossible inverse, an undecidable pivot) is reported with GR_UNABLE.

Set R to the reduced row echelon form of A (pivots equal to 1, zeros above and below them) and rank to the rank, using LU decomposition or fraction-free LU decomposition followed by back substitution. This requires a field, or at least that all pivots are invertible. gr_mat_rref uses the LU version.

Compute R and \mathit{den} such that R / \mathit{den} is the reduced row echelon form of A over the fraction field of an integral domain, with entries of R in the ring. The denominator is a common, not necessarily minimal, denominator: the pivots of R all equal \mathit{den}.

Set the columns of X to a basis of the right nullspace \{x : Ax = 0\} of the m \times n matrix A, read off from the reduced row echelon form: one basis vector per non-pivot column. gr_mat_nullspace resizes X to n rows and as many columns as the nullity; the no_resize version requires X with n rows and at least that many columns, zero-pads the rest and writes the nullity to nullity.

Sets X to a nullspace basis given a reduced row echelon form A of rank rank, as computed by gr_mat_rref, or by gr_mat_rref_den with the denominator Aden (NULL for an ordinary rref). X must have n rows and n - \mathit{rank} columns. Returns GR_UNABLE if a pivot position cannot be identified because a zero test is undecidable.

Sets res to the inverse of a square matrix in the ring, computed with gr_mat_nonsingular_solve and the identity matrix. Returns GR_DOMAIN if the matrix is shown to be not invertible over R (this is not reliable over rings with zero divisors) and GR_UNABLE if invertibility cannot be decided. For the inverse over the fraction field, use gr_mat_nonsingular_solve_den or gr_mat_adjugate.

Set adj to the adjugate of A and det to its determinant, so that \operatorname{adj}(A) A = A \operatorname{adj}(A) = \det(A) I. The cofactor version computes n^2 minors; the charpoly version evaluates (-1)^{n-1}(A^{n-1} + c_{n-1} A^{n-2} + \cdots + c_1 I) from the characteristic polynomial x^n + c_{n-1} x^{n-1} + \cdots + c_0 by rectangular splitting, with O(\sqrt{n}) matrix products. The default uses cofactors for n \le 5 and the characteristic polynomial otherwise. All versions are division-free when the characteristic polynomial is.

Set res to the characteristic polynomial \det(xI - A) of a square matrix (GR_DOMAIN otherwise). _gr_mat_charpoly calls the ring method GR_METHOD_MAT_CHARPOLY; the default _gr_mat_charpoly_generic uses Berkowitz's algorithm, which works over any commutative ring. Rings with fast arithmetic override it: nmod and fmpz_mod use Danilevsky's algorithm above dimension 8 (16 for fmpz_mod), falling back to Berkowitz when it meets a non-invertible pivot; fq and its variants use Berkowitz up to dimension 12 and Danilevsky above; fmpz and fmpq call fmpz_mat_charpoly and fmpq_mat_charpoly; mpn_mod has its own. The underscore versions write n + 1 coefficients. See Characteristic and minimal polynomials.

The characteristic polynomial by Berkowitz's division-free algorithm [Ber1984], with O(n^4) ring operations (as O(n) matrix-vector products for each of n leading principal submatrices). Since it needs no divisions and no zero tests, it succeeds whenever the ring arithmetic does.

The characteristic polynomial with O(n^3) operations, by Danilevsky's reduction to a block companion (Frobenius) form [Dan1937], or by reduction to Hessenberg form with Gaussian elimination or with Householder reflections followed by gr_mat_charpoly_from_hessenberg. The inplace version overwrites its input.

The characteristic polynomial by the Faddeev-LeVerrier algorithm [FadSom1972], which computes M_k = A M_{k-1} + c_{n-k+1} I and c_{n-k} = -\operatorname{tr}(A M_k) / k, and, if adj is not NULL, the adjugate matrix at no extra cost. The bsgs version is the baby-step giant-step variant of Preparata and Sarwate [PreSar1978], with O(n^{3.5}) operations instead of O(n^4) and \sqrt{n} temporary matrices. Both divide by 1, \ldots, n and fail (with GR_DOMAIN or GR_UNABLE) in small positive characteristic or without division.

The characteristic polynomial of a matrix in upper Hessenberg form (not checked), by the recurrence on its leading principal submatrices, with O(n^3) operations and no divisions.

Set res to an upper Hessenberg matrix similar to A, computed with Gaussian elimination (similarity transformations by elementary matrices) or with Householder reflections; the default uses Gaussian elimination. They need divisions and zero tests, and the Householder version also square roots and complex conjugation, and return GR_DOMAIN or GR_UNABLE if these fail.

Sets res to the minimal polynomial of a square matrix over a field, by computing the minimal polynomials of Krylov sequences v, Av, A^2 v, \ldots of unit vectors with incremental Gaussian elimination (gr_mat_reduce_row) and taking their least common multiple, with O(n^3) field operations. See Minimal polynomial.

Sets the n \times n matrix res to the companion matrix of a polynomial of degree n, whose characteristic polynomial is the monic associate of the polynomial. The leading coefficient must be invertible (not checked). The underscore version reads n + 1 coefficients.

Set res_num and the scalar res_den such that \mathit{res\_num} / \mathit{res\_den} is the companion matrix of a polynomial of degree n, without dividing: res_den is the leading coefficient. The underscore version reads n + 1 coefficients.

Replaces the n \times n matrix M by P^{-1} M P in place, where P is the identity matrix with the zero entries of row r replaced by d (an elementary similarity transformation). Requires 0 \le r < n; returns GR_DOMAIN if M is not square or r is out of range. Such transformations preserve the determinant and the characteristic and minimal polynomials.

Reduces row m of A by the previous rows, which are assumed to be in Gauss form though not necessarily in order, as needed for incremental elimination (minimal polynomials). Entry i of P is the row of A with a pivot in column i, or -1; entry i of L is the number of possibly nonzero columns of row i (all equal to the number of columns unless A is zero on the right), with nondecreasing entries. On return, column is the column of the new pivot of row m (the first available column from the left, also recorded in P), or -1 if the row reduces to zero.

Find the eigenvalues of a square matrix in the ring ctx, without repetition, with their multiplicities in mult, as the roots of the characteristic polynomial (gr_poly_roots, respectively gr_poly_roots_other for a matrix over another ring mat_ctx). As for gr_poly_roots, the multiplicities may add up to less than n if the ring is not algebraically closed. Over exact real and complex numbers, see Linear algebra over exact numbers.

Compute a diagonalisation L A R = D of a square matrix: D is the vector of eigenvalues repeated according to multiplicity, the columns of R are right eigenvectors and the rows of L = R^{-1} left eigenvectors, so that A = R D R^{-1}. L or R may be NULL to skip computing it.

Computes the Jordan block structure of a square matrix over a field containing its eigenvalues: the distinct eigenvalues in lambda, the number of blocks in num_blocks, and for each block the index of its eigenvalue in block_lambda and its size in block_size (arrays of length at least n). The block sizes are determined from the ranks of the powers of A - \lambda I.

Sets mat to the Jordan matrix with the given blocks (each block has its eigenvalue on the diagonal and ones on the superdiagonal).

Given the block structure computed by gr_mat_jordan_blocks, sets mat to an invertible matrix P with A = P J P^{-1}, where J is the Jordan matrix.

Computes the Jordan form J of a square matrix and, if P is not NULL, a transformation with A = P J P^{-1}. The eigenvalues must lie in the ring, which must be a field with decidable equality for the result to be meaningful (exact fields such as qqbar and ca). The algorithm and its use over exact numbers are described in linear algebra over exact numbers.

Compute the matrix function f(A) = P f(J) P^{-1} from the Jordan decomposition A = P J P^{-1}. The function is given by jet_func, with signature int jet_func(gr_ptr r, gr_srcptr x, slong n, gr_ctx_t ctx), which writes the Taylor coefficients f(x), f'(x), \ldots, f^{(n-1)}(x) / (n-1)! to r. The param version takes a function int jet_func(gr_ptr r, gr_srcptr x, slong n, gr_srcptr c, gr_ctx_t ctx) of x and an extra parameter c, which may in fact point to arbitrary data understood by jet_func.

Set res to the exponential or the (principal) logarithm of a square matrix. The default methods use the Jordan decomposition; arb and acb override the exponential with arb_mat_exp and acb_mat_exp (scaling and squaring of Taylor series with error bounds). The Jordan versions need the eigenvalues in the ring and are meant for exact fields such as qqbar and ca; elsewhere they typically return GR_UNABLE.

Set res to A^c for a scalar or rational exponent, using the Jordan decomposition (principal branch). The non-Jordan versions first check for small integer exponents and use gr_mat_pow_si for them.

Set res to a square root A^{1/2} or reciprocal square root A^{-1/2}, computed from the Jordan decomposition (the only algorithm implemented). This often fails with GR_UNABLE when the ring is not an algebraically closed field, even if a square root exists.

Returns whether A is square with A A^T = A^T A = I. The ring is assumed to be commutative.

Return whether A A^T (row versions) or A^T A (column versions) is diagonal (orthogonal) or the identity (orthonormal).

Sets the square matrix A to a random orthogonal matrix, using Cayley's construction (I + S)^{-1}(I - S) with S skew-symmetric, or a permutation matrix if that fails, and then negates random rows. Returns GR_DOMAIN if A is not square. The ring is assumed to be commutative.

Compute a reduced QR decomposition A = QR of an m \times n matrix of full column rank (m \ge n; Q is m \times n with orthonormal columns and R is n \times n upper triangular), or the transposed LQ decomposition A = LQ of a matrix of full row rank (m \le n; L is m \times m lower triangular and Q is m \times n with orthonormal rows).

Return whether the rows of A are LLL-reduced with parameters (\delta, \eta), by Gram-Schmidt orthogonalisation in the ring and checking the conditions row by row. The with_removal version checks only the first newd rows, as needed when rows whose Gram-Schmidt norm exceeds the bound gs_B are to be removed (as in fmpz_lll).

Sets all entries of res to one.

Sets res to a Pascal matrix of binomial coefficients: for triangular 0 the symmetric matrix with entries \binom{i+j}{i} (the rows of Pascal's triangle along the antidiagonals), for 1 the upper triangular matrix with entries \binom{j}{i}, and for -1 the lower triangular matrix with entries \binom{i}{j}.

Sets res to the matrix of Stirling numbers: unsigned of the first kind for kind 0, signed of the first kind for 1, and of the second kind for 2.

Sets res to the Hilbert matrix with entries 1 / (i + j + 1), which requires these inverses in R.

Sets the square matrix res to a Hadamard matrix of its size if possible. Returns GR_DOMAIN if none exists (the size must be 0, 1, 2 or a multiple of 4), and GR_UNABLE if the implementation cannot construct one.

Sets res to a random matrix with a nonuniform distribution meant for testing (sparse matrices and special values are favoured).

Randomises mat in place by at most count random additions or subtractions of one row (column) to another, which preserves the determinant and the rank.

Sets mat to a random row and column permutation of the diagonal matrix whose first n diagonal entries are given by the vector diag, writing the parity (0 or 1) of the permutation to parity. Returns GR_DOMAIN if the diagonal of mat is shorter than n.

Sets mat to a random sparse matrix with the given rank and exactly rank nonzero entries (meaningful over integral domains). A dense matrix of the same rank is obtained by calling gr_mat_randops afterwards.

Randomises the square matrix mat in place by at most opcount conjugations by random elementary matrices, which preserve the characteristic and minimal polynomials.

Solve the system of linear differential equations Y'(t) = A(t) Y(t) with Y(0) = Y_0, where A(t) is a matrix of rational functions given by a matrix of polynomials A_numerator (over the polynomial ring A_poly_ctx) and a common denominator A_denominator with nonzero constant term, and set Y to the power series solution truncated to length len. The entries of Y are polynomials over the ring sol_poly_ctx.

Helper for the contexts of matrix rings and spaces (here ctx is such a context, not the coefficient ring): resizes mat to r \times c if the context allows matrices of all sizes, and otherwise returns GR_DOMAIN unless r \times c is the shape of the context.

Test an implementation of gr_mat_mul, gr_mat_lu, gr_mat_det, _gr_mat_charpoly, gr_mat_nonsingular_solve_tril or gr_mat_nonsingular_solve_triu by comparing with the generic algorithms and checking algebraic identities.

Test the accuracy of a matrix multiplication over an approximate (floating-point) ring: the first checks \|C - AB\| \le \|A\| \|B\| \cdot \mathit{rel\_tol} in the max norm against classical multiplication; the second generates matrices with nonnegative entries and checks that every entry has relative error at most \mathit{rel\_tol}.