Improve method choice in gr_mat_rank; wrap gr_mat_rank_lu in finite f…

…ield modules

src/fmpz_mod_mat/rank.c

@@ -11,26 +11,15 @@
     (at your option) any later version.  See <https://www.gnu.org/licenses/>.
 */
 
+#include "gr.h"
+#include "gr_mat.h"
 #include "fmpz_mod_mat.h"
 
 slong fmpz_mod_mat_rank(const fmpz_mod_mat_t A, const fmpz_mod_ctx_t ctx)
 {
-    slong m, n, rank;
-    slong *perm;
-    fmpz_mod_mat_t tmp;
-
-    m = A->r;
-    n = A->c;
-
-    if (m == 0 || n == 0)
-        return 0;
-
-    fmpz_mod_mat_init_set(tmp, A, ctx);
-    perm = flint_malloc(sizeof(slong) * m);
-
-    rank = fmpz_mod_mat_lu(perm, tmp, 0, ctx);
-
-    flint_free(perm);
-    fmpz_mod_mat_clear(tmp, ctx);
+    gr_ctx_t gr_ctx;
+    slong rank;
+    _gr_ctx_init_fmpz_mod_from_ref(gr_ctx, ctx);
+    GR_MUST_SUCCEED(gr_mat_rank_lu(&rank, (const gr_mat_struct *) A, gr_ctx));
     return rank;
 }

src/fq_mat_templates/rank.c

@@ -12,29 +12,18 @@
 
 #ifdef T
 
+#include "gr.h"
+#include "gr_mat.h"
 #include "templates.h"
 
 slong
 TEMPLATE(T, mat_rank) (const TEMPLATE(T, mat_t) A,
                        const TEMPLATE(T, ctx_t) ctx)
 {
-    slong m, n, rank;
-    slong *perm;
-    TEMPLATE(T, mat_t) tmp;
-
-    m = A->r;
-    n = A->c;
-
-    if (m == 0 || n == 0)
-        return 0;
-
-    TEMPLATE(T, mat_init_set) (tmp, A, ctx);
-    perm = flint_malloc(sizeof(slong) * m);
-
-    rank = TEMPLATE(T, mat_lu) (perm, tmp, 0, ctx);
-
-    flint_free(perm);
-    TEMPLATE(T, mat_clear) (tmp, ctx);
+    gr_ctx_t gr_ctx;
+    slong rank;
+    TEMPLATE3(_gr_ctx_init, T, from_ref)(gr_ctx, ctx);
+    GR_MUST_SUCCEED(gr_mat_rank_lu(&rank, (const gr_mat_struct *) A, gr_ctx));
     return rank;
 }
 

src/gr_mat/diagonalization.c

@@ -106,6 +106,8 @@ gr_mat_diagonalization_generic(gr_vec_t D, gr_mat_t L, gr_mat_t R, const gr_mat_
     if (gr_mat_is_square(A, ctx) != T_TRUE)
         return GR_DOMAIN;
 
+    gr_ctx_println(ctx);
+
     gr_ctx_init_fmpz(ZZ);
     gr_vec_init(eigenvalues, 0, ctx);
     gr_vec_init(mult, 0, ZZ);

src/gr_mat/rank.c

@@ -1,5 +1,5 @@
 /*
-    Copyright (C) 2022 Fredrik Johansson
+    Copyright (C) 2022, 2025 Fredrik Johansson
 
     This file is part of FLINT.
 
@@ -13,17 +13,108 @@
 #include "gr.h"
 #include "gr_mat.h"
 
-/* todo: different algorithms */
+int
+gr_mat_rank_fflu(slong * rank, const gr_mat_t A, gr_ctx_t ctx)
+{
+    slong n, m;
+    slong * P;
+    int status;
+    gr_mat_t T;
+    gr_ptr den;
+
+    n = gr_mat_nrows(A, ctx);
+    m = gr_mat_ncols(A, ctx);
+
+    if (n == 0 || m == 0)
+    {
+        *rank = 0;
+        return GR_SUCCESS;
+    }
+    else
+    {
+        GR_TMP_INIT(den, ctx);
+
+        gr_mat_init(T, n, m, ctx);
+        P = _perm_init(n);
+
+        status = gr_mat_fflu(rank, P, T, den, A, 0, ctx);
+
+        gr_mat_clear(T, ctx);
+        _perm_clear(P);
+
+        GR_TMP_CLEAR(den, ctx);
+
+        if (status != GR_SUCCESS)
+            status |= GR_UNABLE;
+
+        return status;
+    }
+}
+
+int
+gr_mat_rank_lu(slong * rank, const gr_mat_t A, gr_ctx_t ctx)
+{
+    slong n, m;
+    slong * P;
+    int status;
+    gr_mat_t T;
+
+    n = gr_mat_nrows(A, ctx);
+    m = gr_mat_ncols(A, ctx);
+
+    if (n == 0 || m == 0)
+    {
+        *rank = 0;
+        return GR_SUCCESS;
+    }
+    else
+    {
+        gr_mat_init(T, n, m, ctx);
+        P = _perm_init(n);
+
+        status = gr_mat_lu(rank, P, T, A, 0, ctx);
+
+        gr_mat_clear(T, ctx);
+        _perm_clear(P);
+
+        if (status != GR_SUCCESS)
+            status |= GR_UNABLE;
+
+        return status;
+    }
+}
+
 int
 gr_mat_rank(slong * rank, const gr_mat_t A, gr_ctx_t ctx)
 {
     truth_t dom;
 
-    dom = gr_ctx_is_integral_domain(ctx);
+    /* Sensible definition over any ring. */
+    if (gr_mat_nrows(A, ctx) == 0 || gr_mat_ncols(A, ctx) == 0)
+    {
+        *rank = 0;
+        return GR_SUCCESS;
+    }
+
+    dom = gr_ctx_is_field(ctx);
+
+    /* Prefer standard LU only over finite fields */
+    /* Prefer FFLU over non-finite fields as it often results in
+       smaller coefficients. TODO: this choice surely wants
+       tuning, e.g. when we have fast matrix multiplication.
+       For example, ca_mat currently uses LU for number fields. */
+    if (dom == T_TRUE && gr_ctx_is_finite(ctx) == T_TRUE)
+        return gr_mat_rank_lu(rank, A, ctx);
+
+    if (dom != T_TRUE)
+        dom = gr_ctx_is_integral_domain(ctx);
 
     if (dom == T_TRUE)
         return gr_mat_rank_fflu(rank, A, ctx);
 
+    /* There are ways to generalize the notation of rank to
+       non-integral domains, but we do not currently implement them;
+       for now we define the rank as undefined. */
     if (dom == T_FALSE)
         return GR_DOMAIN;