LU factorization with full pivoting

[lampretl]

Hey, I'm just learning how to use FLINT through Julia's FLINT_jll. So far, I have 3 questions:

(1) Does the library allow computing the LU factorization of a not-necessarily-square not-necessarily-full-rank matrix, say, over the rationals (fmpz_mat) or over a finite field? If the rank is not full, it can happen during the algorithm that a zero row or column appears, which would mean we need row and column permutation (or transformation matrices). I'm not sure if fmpz_mat_fflu or fmpq_mat_rref offer this.

(2) Are small finite fields implemented to avoid unnecessary RAM usage. E.g. do elements from fields Z/2, Z/3, Z/5, Z/7, Z/11, Z/13 require 1 byte (which is enough) or more? More generally, do elements from prime fields Z/p, where p<sqrt(typemax(UIntN)) for N=8,16,32,64,128, require N/8 bytes (which suffices), or more?

(3) Is the LU factorization of a matrix with entries from a finite field optimized to use multithreading, cache-blocking, vectorization?

[albinahlback]

As for (2), Fredrik has implemented nmod8 type stuff (and nmod16 etc.), which is present in the gr module. However, we do not treat powers of two in a special manner (not that you asked, just that this can be implemented faster).

As for (3), I'm not an expert, but it does not look like it uses multithreading. I'm not so sure about cache-blocking (you'd have to look more carefully about this), but vectorization should be implemented, at least in parts of the code (but perhaps not optimally).

[albinahlback]

@fredrik-johansson is the real expert here though.

[lampretl]

@albinahlback thanks for your reply! I still hope to get info on (1) and (3). I'd appreciate if @fredrik-johansson could chip in. : )

[fredrik-johansson]

(1) Yes, functions like fmpz_mat_fflu accept matrices of any shape and rank.

(2) It is possible to use compact representations, e.g. gr_ctx_init_nmod8 initializes a prime finite field using one byte per element. However, this representation is not currently as optimized for speed as nmod which use 8 bytes per element.

(3) Yes, but our implementations of these techniques are far from optimal. Note that FLINT can be configured to use BLAS for matrix multiplication. This gives the best performance if you want to compute LU factorizations of large matrices over finite fields.

[lampretl]

Wow, I was not aware that BLAS can be used on matrices over finite fields! Could you give me a code snippet that does an in-place matrix update C += αA*B?

[fredrik-johansson]

BLAS itself doesn't multiply matrices over finite fields. Internally we convert finite field elements to floating-point numbers, multiply with BLAS, convert back to integers and reduce mod p.

[lampretl]

Aaa, conversion to floats. Fascinating, how much of modern engineering is reliant on a library, written in 1987 in Fortran. :) This is probably the source you refer to: https://github.com/flintlib/flint/blob/0b4dfbc580679a68df9778edb829efcda5e3a798/src/fmpz_mat/mul_blas.c

In the approach you mentioned, are there ever any problems with float weirdness (will the decimal part of entries in the product matrix always be 0 or must we round/floor/ceil)?

If we are computing the product AB of matrices of size (l,m),(m,n) over T∈(Float16,Float32,Float64) modulo p (so entries in A and B are 0 ≤ _ < p), is it sufficient that m(p-1)^2 ≤ floatmax(T)? For instance, if m≤655, can I compute modulo 11 over Float16, since 655*10^2 ≤ 6.55e4 = floatmax(Float16)? Or must the condition be more restrictive to avoid rounding and overflow? I'm not very knowledgeable on how floats work.

[fredrik-johansson]

It is sufficient that $m(p-1)^2 &lt; 2^p$ where $p$ is the mantissa precision of the type, though one can extend the range with intermediate reductions. There are several papers where you can read about this, e.g. https://membres-ljk.imag.fr/Jean-Guillaume.Dumas/Publications/Field_blas.pdf