In ArithProg bootstrap: we need to LU decompose the whole Krylov matr…

…ix, not just the first N rows. Fixes #381. Add a regression test
linbox-team · Sep 6, 2023 · e08aea9 · e08aea9
1 parent ee1bbde
commit e08aea9
Show file tree

Hide file tree

Showing 2 changed files with 31 additions and 8 deletions.
diff --git 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
@@ -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;
 }