From e08aea90f4f85d6b81018358d75bab235f9972ad Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Cl=C3=A9ment=20Pernet?= <clement.pernet@gmail.com>
Date: Wed, 6 Sep 2023 18:41:01 +0200
Subject: [PATCH 1/4] In ArithProg bootstrap: we need to LU decompose the whole
 Krylov matrix, not just the first N rows. Fixes #381. Add a regression test

---
 fflas-ffpack/ffpack/ffpack_frobenius.inl | 16 ++++++++--------
 tests/regression-check.C                 | 23 +++++++++++++++++++++++
 2 files changed, 31 insertions(+), 8 deletions(-)

diff --git a/fflas-ffpack/ffpack/ffpack_frobenius.inl b/fflas-ffpack/ffpack/ffpack_frobenius.inl
index a42fb5497..268a65cc8 100644
--- a/fflas-ffpack/ffpack/ffpack_frobenius.inl
+++ b/fflas-ffpack/ffpack/ffpack_frobenius.inl
@@ -93,7 +93,6 @@ namespace FFPACK { namespace Protected {
         size_t *dK = FFLAS::fflas_new<size_t>(noc*degree);
         for (size_t i=0; i<noc; ++i)
             dK[i]=0;
-
 #ifdef __FFLASFFPACK_ARITHPROG_PROFILING
         Givaro::Timer timkrylov, timelim, timsimilar, timrest;
         timkrylov.start();
@@ -128,14 +127,15 @@ namespace FFPACK { namespace Protected {
         timelim.start();
 #endif
         size_t * Pk = FFLAS::fflas_new<size_t>(N);
-        size_t * Qk = FFLAS::fflas_new<size_t>(N);
-        for (size_t i=0; i<N; ++i)
+        size_t nrows = noc*degree;
+        size_t * Qk = FFLAS::fflas_new<size_t>(nrows);
+        for (size_t i=0; i<nrows; ++i)
             Qk[i] = 0;
         for (size_t i=0; i<N; ++i)
             Pk[i] = 0;
 
             // @todo: replace by PLUQ
-        size_t R = LUdivine(F, FFLAS::FflasNonUnit, FFLAS::FflasNoTrans, N, N, K, ldk, Pk, Qk);
+        size_t R = LUdivine(F, FFLAS::FflasNonUnit, FFLAS::FflasNoTrans, noc*degree, N, K, ldk, Pk, Qk);
         size_t row_idx = 0;
         size_t ii=0;
         size_t dold = degree;
@@ -198,8 +198,9 @@ namespace FFPACK { namespace Protected {
         size_t Ma = Mk;
         size_t Ncurr = R;
         size_t offset = Ncurr-1;
-        for (size_t i=Mk-1; i>=nb_full_blocks+1;  --i){
-            if (dK[i] >= 1){
+        for (size_t ip1=Mk; ip1 > nb_full_blocks;  --ip1){
+                size_t i = ip1-1;
+                if (dK[i] >= 1){
                 for (size_t j = offset+1; j<R; ++j)
                     if (!F.isZero(*(K4 + i*ldk + j))){
                         //std::cerr<<"FAIL C != 0 in preconditionning"<<std::endl;
@@ -286,7 +287,6 @@ namespace FFPACK { namespace Protected {
         ldb = Ma;
         Nb = Ncurr;
         B = FFLAS::fflas_new (F, Ncurr, ldb);
-
         for (size_t j=0; j<Ma; ++j)
             FFLAS::fassign(F, Ncurr, K4+j*ldk, 1, B+j, ldb);
         FFLAS::fflas_delete (dA, dK, K4);
@@ -299,7 +299,7 @@ namespace FFPACK { namespace Protected {
                const size_t N, typename PolRing::Domain_t::Element_ptr A, const size_t lda,
                const size_t degree)
     {
-
+        if (!N) return frobeniusForm;
         typedef typename PolRing::Domain_t Field;
         typedef typename PolRing::Element Polynomial;
         const Field& F = PR.getdomain();
diff --git a/tests/regression-check.C b/tests/regression-check.C
index 3747e295f..7632a85aa 100644
--- a/tests/regression-check.C
+++ b/tests/regression-check.C
@@ -120,6 +120,28 @@ bool gf2ModularBalanced(){
     return true;
 }
 
+bool charpolyArithProg(){
+
+    Modular<double> F(37);
+    Poly1Dom<Modular<double> >PolRing(F,'X');
+
+    // Reading the matrix from a file
+    double* A;
+    size_t m, n;
+    std::string file("data/regression_charpoly.sms");
+    ReadMatrix(file.c_str(), F, m, n, A);
+    // here m=n=35
+    Poly1Dom<Modular<double> >::Element charp(n+1);
+
+    CharPoly(PolRing, charp, n, A, n);
+
+    fflas_delete(A);
+    bool pass = F.isOne(charp[n]);
+    for (size_t i = 0; i<n; i++)
+            pass = pass && (F.isZero(charp[i]));
+    if (!pass) std::cerr<<"charpolyArithProg FAILED"<<std::endl;
+    return pass;
+}
 
 int main() {
     bool pass = true ;
@@ -129,6 +151,7 @@ int main() {
     pass = pass && checkZeroDimCharpoly();
     pass = pass && checkZeroDimMinPoly();
     pass = pass && gf2ModularBalanced();
+    pass = pass && charpolyArithProg();
     return !pass;
 }
 

From 452c2ba9e645970c997306057f095ef89f4e1b36 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Cl=C3=A9ment=20Pernet?= <clement.pernet@gmail.com>
Date: Wed, 6 Sep 2023 18:41:17 +0200
Subject: [PATCH 2/4] matrix for the regression test

---
 tests/data/regression_charpoly.sms | 97 ++++++++++++++++++++++++++++++
 1 file changed, 97 insertions(+)
 create mode 100644 tests/data/regression_charpoly.sms

diff --git a/tests/data/regression_charpoly.sms b/tests/data/regression_charpoly.sms
new file mode 100644
index 000000000..477e99ecb
--- /dev/null
+++ b/tests/data/regression_charpoly.sms
@@ -0,0 +1,97 @@
+35 35 M
+1 2 113
+1 4 229
+1 33 173
+2 19 113
+2 21 113
+2 27 229
+2 29 173
+3 9 113
+3 24 229
+4 25 229
+4 27 113
+5 8 113
+5 10 113
+5 14 113
+5 26 229
+6 8 229
+6 12 113
+6 18 113
+6 20 173
+6 29 113
+7 6 113
+7 16 229
+7 28 173
+7 33 113
+8 20 113
+8 22 113
+9 15 113
+10 9 229
+10 22 113
+10 23 113
+11 32 113
+12 5 113
+12 16 113
+12 22 229
+13 5 229
+13 6 113
+13 8 173
+13 27 113
+13 34 113
+14 20 113
+14 23 113
+14 32 229
+15 24 113
+16 8 113
+17 8 113
+17 28 113
+17 32 173
+18 5 113
+18 17 113
+18 20 229
+18 26 173
+18 31 113
+19 4 113
+19 13 113
+19 16 173
+19 30 229
+20 11 113
+20 26 113
+21 7 113
+21 12 229
+21 13 113
+21 17 173
+22 26 113
+23 3 113
+23 15 229
+23 26 113
+25 3 229
+25 10 113
+26 9 113
+26 32 113
+27 10 229
+27 12 113
+27 30 113
+28 15 173
+28 20 113
+28 35 113
+29 11 173
+29 16 113
+29 17 113
+30 5 113
+30 23 229
+30 25 113
+31 9 173
+31 11 229
+31 14 113
+31 28 113
+32 15 113
+33 29 113
+33 35 173
+34 14 229
+34 18 113
+34 22 173
+34 30 113
+35 11 113
+35 24 173
+0 0 0

From af017f372e00d5a1a8d76ee0b4eba7ba5176b896 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Cl=C3=A9ment=20Pernet?= <clement.pernet@gmail.com>
Date: Wed, 6 Sep 2023 19:07:16 +0200
Subject: [PATCH 3/4] restoring the proper switch in Auto mode depending on the
 threshold macors

---
 fflas-ffpack/ffpack/ffpack_charpoly.inl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/fflas-ffpack/ffpack/ffpack_charpoly.inl b/fflas-ffpack/ffpack/ffpack_charpoly.inl
index a8210280d..3098dbdaf 100644
--- a/fflas-ffpack/ffpack/ffpack_charpoly.inl
+++ b/fflas-ffpack/ffpack/ffpack_charpoly.inl
@@ -65,9 +65,9 @@ namespace FFPACK {
 
         FFPACK_CHARPOLY_TAG tag = CharpTag;
         if (tag == FfpackAuto){
-            if (N < degree)
+            if (N < __FFLASFFPACK_CHARPOLY_Danilevskii_LUKrylov_THRESHOLD)
                 tag = FfpackDanilevski;
-            else if (N < degree)
+            else if (N < __FFLASFFPACK_CHARPOLY_LUKrylov_ArithProg_THRESHOLD)
                 tag = FfpackLUK;
             else
                 tag = FfpackArithProgKrylovPrecond;
@@ -110,7 +110,7 @@ namespace FFPACK {
                         else
                             return CharPoly (R, charp, N, A, lda, G, FfpackLUK);
                     }
-                                   } while (cont);
+                } while (cont);
                 return charp;
             }
             case FfpackArithProg:

From 175ec1e1ddfd707e73c90c3650ecbb7d5704142c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Cl=C3=A9ment=20Pernet?= <clement.pernet@gmail.com>
Date: Thu, 7 Sep 2023 16:04:22 +0200
Subject: [PATCH 4/4] fixed number of completed blocks

---
 fflas-ffpack/ffpack/ffpack_frobenius.inl | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/fflas-ffpack/ffpack/ffpack_frobenius.inl b/fflas-ffpack/ffpack/ffpack_frobenius.inl
index 268a65cc8..532d41df6 100644
--- a/fflas-ffpack/ffpack/ffpack_frobenius.inl
+++ b/fflas-ffpack/ffpack/ffpack_frobenius.inl
@@ -164,7 +164,6 @@ namespace FFPACK { namespace Protected {
         std::cerr <<"  LU (Krylov)              : "<<timelim.usertime()<<std::endl;
         timsimilar.start();
 #endif
-
         // Selection of the last iterate of each block
 
         typename Field::Element_ptr K3 = FFLAS::fflas_new (F, Mk, N);
@@ -198,7 +197,11 @@ namespace FFPACK { namespace Protected {
         size_t Ma = Mk;
         size_t Ncurr = R;
         size_t offset = Ncurr-1;
-        for (size_t ip1=Mk; ip1 > nb_full_blocks;  --ip1){
+        // How many non-full blocks are completed :
+        //  - if no full_blocks: all of them
+        //  - if there is at least one full block, then the first non-full will be further eliminated later on and should not be pushed as completed
+        size_t last_completed = nb_full_blocks + nb_full_blocks?1:0;
+        for (size_t ip1=Mk; ip1 > last_completed;  --ip1){
                 size_t i = ip1-1;
                 if (dK[i] >= 1){
                 for (size_t j = offset+1; j<R; ++j)