Skip to content

Fixes a bug the bootstrap of Charpoly ArithProgPrecond method #382

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 4 commits into
base: master
Choose a base branch
from
Open

Fixes a bug the bootstrap of Charpoly ArithProgPrecond method #382

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
e08aea9
In ArithProg bootstrap: we need to LU decompose the whole Krylov matr…
ClementPernet Sep 6, 2023
452c2ba
matrix for the regression test
ClementPernet Sep 6, 2023
af017f3
restoring the proper switch in Auto mode depending on the threshold m…
ClementPernet Sep 6, 2023
175ec1e
fixed number of completed blocks
ClementPernet Sep 7, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 3 additions & 3 deletions fflas-ffpack/ffpack/ffpack_charpoly.inl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -110,7 +110,7 @@ namespace FFPACK {
else
return CharPoly (R, charp, N, A, lda, G, FfpackLUK);
}
} while (cont);
} while (cont);
return charp;
}
case FfpackArithProg:
Expand Down
21 changes: 12 additions & 9 deletions fflas-ffpack/ffpack/ffpack_frobenius.inl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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();
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -198,8 +197,13 @@ 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){
// 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)
if (!F.isZero(*(K4 + i*ldk + j))){
//std::cerr<<"FAIL C != 0 in preconditionning"<<std::endl;
Expand Down Expand Up @@ -286,7 +290,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);
Expand All @@ -299,7 +302,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();
Expand Down
97 changes: 97 additions & 0 deletions tests/data/regression_charpoly.sms
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
23 changes: 23 additions & 0 deletions tests/regression-check.C
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 ;
Expand All @@ -129,6 +151,7 @@ int main() {
pass = pass && checkZeroDimCharpoly();
pass = pass && checkZeroDimMinPoly();
pass = pass && gf2ModularBalanced();
pass = pass && charpolyArithProg();
return !pass;
}

Expand Down