output of fmpz_poly_mat_fflu

[Hjalte Frellesvig]

Dear Flint developers.

I am struggling to understand the output of fmpz_poly_mat_fflu
slong fmpz_poly_mat_fflu(fmpz_poly_mat_t B, fmpz_poly_t den, slong *perm, const fmpz_poly_mat_t A, int rank_check)
which is supposed to generate a LU decomposition of the input matrix A.

How do I extract "L" and "U" from (I assume) B and den?
Or am I misunderstanding what the function does?

Furthermore please let me know if you have a preferred way for users to ask such questions, since it might not be my last.

Best wishes, Hjalte Frellesvig

[tho...@gmail.com]

If you don't mind a bit of julia code, you can have a look at how we do it in Nemo for fmpz_mat_fflu here: here: https://github.com/Nemocas/Nemo.jl/blob/97ee5f6da894b9d1a295f3f731dad88aac35bf40/src/flint/fmpz_mat.jl#L620-L657

I think it should be similar for fmpz_poly_mat_fflu.

Best wishes

Tommy

[Claus Fieker]

On Fri, Dec 01, 2023 at 01:43:12AM -0800, tho...@gmail.com wrote:
> If you don't mind a bit of julia code, you can have a look at how we do it
> in Nemo for fmpz_mat_fflu here:
> here: https://github.com/Nemocas/Nemo.jl/blob/97ee5f6da894b9d1a295f3f731dad88aac35bf40/src/flint/fmpz_mat.jl#L620-L657
>
> I think it should be similar for fmpz_poly_mat_fflu.
>

I had the problem last week, so

#takes the result of lu! and turns it into a proper lu
function real_lu(P::Perm, U, r::Int)
m = nrows(U)
n = ncols(U)
R = base_ring(U)
L = similar(U, m, m)

for i = 1:m
for j = 1:n
if i > j
L[i, j] = U[i, j]
U[i, j] = R(0)
elseif i == j
L[i, j] = R(1)
elseif j <= m
L[i, j] = R(0)
end
end
end
return r, P, L, U
end

> Best wishes
> Tommy
>
> hjalte.f...@gmail.com schrieb am Montag, 20. November 2023 um 19:47:55
> UTC+1:
>
> > Dear Flint developers.
> >
> > I am struggling to understand the output of fmpz_poly_mat_fflu
> > slong fmpz_poly_mat_fflu(fmpz_poly_mat_t B, fmpz_poly_t den, slong *perm,
> > const fmpz_poly_mat_t A, int rank_check)
> > which is supposed to generate a LU decomposition of the input matrix A.
> >
> > How do I extract "L" and "U" from (I assume) B and den?
> > Or am I misunderstanding what the function does?
> >
> > Furthermore please let me know if you have a preferred way for users to
> > ask such questions, since it might not be my last.
> >
> > Best wishes, Hjalte Frellesvig
> >
> >
> >
>

> --
>
> ---
> You received this message because you are subscribed to the Google Groups "flint-devel" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to flint-devel...@googlegroups.com.
> To view this discussion on the web, visit https://groups.google.com/d/msgid/flint-devel/66474ffc-4bfa-4f87-bca2-cac16d09d505n%40googlegroups.com.

[Hjalte Frellesvig]

Dear Tommy, Claus and All.
Thank you for the feedback.
I don't think the Julia algorithm you refer to is exactly what is used here.
Using that algorithm to compute P, L, and U, and then comparing PLU to A I get, for a two-by-two example, one row to agree, but the other is off by a multiplicative factor.
I add a short test program containing a simple example illustrating this.
Best wishes, Hjalte

You received this message because you are subscribed to a topic in the Google Groups "flint-devel" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/flint-devel/VniPUD9mB18/unsubscribe.
To unsubscribe from this group and all its topics, send an email to flint-devel...@googlegroups.com.
To view this discussion on the web, visit https://groups.google.com/d/msgid/flint-devel/ZWm0QxbkMWXUZb1R%40mathematik.uni-kl.de.

[Claus Fieker]

On Thu, Dec 07, 2023 at 11:23:20AM +0100, Hjalte Frellesvig wrote:
> Dear Tommy, Claus and All.
> Thank you for the feedback.
> I don't think the Julia algorithm you refer to is exactly what is used here.
> Using that algorithm to compute P, L, and U, and then comparing PLU to A I
> get, for a two-by-two example, one row to agree, but the other is off by a
> multiplicative factor.
> I add a short test program containing a simple example illustrating this.
> Best wishes, Hjalte

From the Nemo Docu:

Return a tuple r, d, p, L, U consisting of the rank of A, a denominator d, a
permutation p of A belonging to P, a lower triangular matrix L and an upper
triangular matrix U such that p(A) = LDU, where p(A) stands for the matrix
whose rows are the given permutation p of the rows of A and such that D is
the diagonal matrix diag(p_1, p_1p_2, \ldots, p_{n-2}p_{n-1}, p_{n-1}p_n)
where the p_i are the inverses of the diagonal entries of L. The denominator
d is set to \pm \mathrm{det}(S) where S is an appropriate submatrix of A (S
= A if A is square and nonsingular) and the sign is decided by the parity of
the permutation.


In this case:

julia> matrix(Qx, 2, 2, [2*x^2, -3x, 2x-1, 5x+3])
[ 2*x^2 -3*x]
[2*x - 1 5*x + 3]

julia> fflu(ans)
(2, 10*x^3 + 12*x^2 - 3*x, (), [2*x^2 0; 2*x-1 10*x^3+12*x^2-3*x], [2*x^2 -3*x; 0 10*x^3+12*x^2-3*x])

julia> Qs,s = rational_function_field(QQ, "s")
(Rational function field over QQ, s)

julia> matrix(Qs, 2, 2, [2*s^2, -3s, 2s-1, 5s+3])
[ 2*s^2 -3*s]
[2*s - 1 5*s + 3]

julia> fflu(ans)
(2, 10*s^3 + 12*s^2 - 3*s, (), [2*s^2 0; 2*s-1 10*s^3+12*s^2-3*s], [2*s^2 -3*s; 0 10*s^3+12*s^2-3*s])

julia> L, U = ans[4], ans[5]
([2*s^2 0; 2*s-1 10*s^3+12*s^2-3*s], [2*s^2 -3*s; 0 10*s^3+12*s^2-3*s])

julia> d = diagonal_matrix([inv(L[1,1]), inv(L[1,1]*L[2,2])])
[1//2//s^2 0]
[ 0 1//20//(s^5 + 6//5*s^4 - 3//10*s^3)]

julia> L*d*U
[ 2*s^2 -3*s]
[2*s - 1 5*s + 3]

the fflu is fraction free so there is some "displacement"

> To view this discussion on the web, visit https://groups.google.com/d/msgid/flint-devel/CAFTk_RsuG0NtcZy7qDxY3BGFCAz3bsdnkgYrwPaHxLgEwP1W_g%40mail.gmail.com.

>
> #include <stdlib.h>
> #include <stdio.h>
>
> #include "flint/fmpz_poly_mat.h"
>
>
>
> void testfflu()
> {
> char *output;
>
> output = (char*) malloc(500*sizeof(char));
>
> fmpz_poly_mat_t A, B;
> fmpz_poly_t det;
>
>
> slong rank, rr=2;
>
> fmpz_poly_mat_init(A, rr, rr);
> fmpz_poly_mat_init(B, rr, rr);
> fmpz_poly_init(det);
>
>
> // 2*x^2
> fmpz_poly_set_coeff_si(fmpz_poly_mat_entry(A,0,0), 2, 2);
> // -3*x
> fmpz_poly_set_coeff_si(fmpz_poly_mat_entry(A,0,1), 1, -3);
> // 2*x-1
> fmpz_poly_set_coeff_si(fmpz_poly_mat_entry(A,1,0), 1, 2);
> fmpz_poly_set_coeff_si(fmpz_poly_mat_entry(A,1,0), 0, -1);
> // 5*x+3
> fmpz_poly_set_coeff_si(fmpz_poly_mat_entry(A,1,1), 1, 5);
> fmpz_poly_set_coeff_si(fmpz_poly_mat_entry(A,1,1), 0, 3);
>
> printf("A:\n");
> fmpz_poly_mat_print(A, "x");
>
> slong perm[2] = {0, 1};
>
>
>
>
> rank = fmpz_poly_mat_fflu(B, det, perm, A, 0);
>
> printf("B:\n");
> fmpz_poly_mat_print(B, "x");
>
> slong i, j;
>
> printf("det = ");
> fmpz_poly_print_pretty(det, "x");
> printf("\n");
>
> printf("rank = %ld\n", rank);
>
> printf("perm = ");
> for(i=0; i<rr; ++i){
> printf("%ld", perm[i]);
> }
> printf("\n\n");
>
>
> fmpz_poly_mat_t P, L, U, PLU;
>
> fmpz_poly_mat_init(P, rr, rr);
> fmpz_poly_mat_init(L, rr, rr);
> fmpz_poly_mat_init(U, rr, rr);
> fmpz_poly_mat_init(PLU, rr, rr);
>
> fmpz_poly_mat_zero(P);
> for(i=0; i<rr; ++i){
> fmpz_poly_one(fmpz_poly_mat_entry(P,i,perm[i]));
> }
>
> fmpz_poly_mat_zero(L);
> fmpz_poly_mat_set(U, B);
>
> for(i=0; i<rr; ++i){
> for(j=0; j<rr; ++j){
> if(i>j){
> fmpz_poly_set(fmpz_poly_mat_entry(L,i,j), fmpz_poly_mat_entry(U,i,j));
> fmpz_poly_zero(fmpz_poly_mat_entry(U,i,j));
> }else if(i==j){
> fmpz_poly_one(fmpz_poly_mat_entry(L,i,j));
> }else if(j<=rr){
> fmpz_poly_zero(fmpz_poly_mat_entry(L,i,j));
> }
> }
> }
>
> fmpz_poly_mat_mul(PLU, P, L);
> fmpz_poly_mat_mul(PLU, PLU, U);
>
> printf("P:\n");
> fmpz_poly_mat_print(P, "x");
> printf("L:\n");
> fmpz_poly_mat_print(L, "x");
> printf("U:\n");
> fmpz_poly_mat_print(U, "x");
>
> printf("PLU:\n");
> fmpz_poly_mat_print(PLU, "x");
>
> fmpz_poly_mat_sub(PLU, PLU, A);
> printf("PLU - A:\n");
> fmpz_poly_mat_print(PLU, "x");
> // For this difference I get the second row to be zeros, but the first is not.
> // PLU and A differ by a multiplicative factor.
>
>
> fmpz_poly_mat_clear(A);
> fmpz_poly_mat_clear(B);
> fmpz_poly_mat_clear(P);
> fmpz_poly_mat_clear(L);
> fmpz_poly_mat_clear(U);
> fmpz_poly_mat_clear(PLU);
> fmpz_poly_clear(det);
> }
>
>
>
> int main(int argc, char* argv[])
> {
> testfflu();
>
> return 0;
> }
>

[Claus Fieker]

On Thu, Dec 07, 2023 at 11:23:20AM +0100, Hjalte Frellesvig wrote:

> Dear Tommy, Claus and All.
> Thank you for the feedback.
> I don't think the Julia algorithm you refer to is exactly what is used here.
> Using that algorithm to compute P, L, and U, and then comparing PLU to A I
> get, for a two-by-two example, one row to agree, but the other is off by a
> multiplicative factor.
> I add a short test program containing a simple example illustrating this.
> Best wishes, Hjalte
>

I'd say the flint-doc is slightly wrong...or at least incomplete

[Hjalte Frellesvig]

Dear Claus and all.

I agree in the end that the flint output should be interpreted as done by the nemo code linked above.

At least I get agreement in all the cases I have tried.

I also agree that this probably should be explained in the flint documentation.
Thanks a lot for the help!

To view this discussion on the web, visit https://groups.google.com/d/msgid/flint-devel/ZXHlkY1tlFduWhke%40mathematik.uni-kl.de.