Commit 70d9e5e6 authored by Charles Bouillaguet's avatar Charles Bouillaguet
Browse files

column renumbering to avoid empty columns

parent cfffd0f2
......@@ -167,15 +167,19 @@ void parser_store_polynomial(void *opaque, int line)
* Degree-d monomials are in monomials[Nd[d]:Nd[d + 1]].
* The number of monomials of degree <= d is Nd[d + 1].
* This is NOT a monomial order : a<b but a*b > b*b
*
* "bad" monomial contain more than one special variable.
* bad_renumbering[j] == -1 if j is a good monomial
* bad_renumbering[j] == i (>= 0) if j is the i-th bad monomial
*/
typedef u128 monomial_t;
monomial_t *monomials; /* array to rank/unrank monomials */
int *mul_table; /* multiplication table */
bool *bad_monomial; /* contains more than one special variable */
monomial_t vmask; /* monomial that contains all special variables */
int nbad; /* total number of bad monomials */
int *bad_renumbering; /* renumber the bad monomials consecutively */
static inline monomial_t mkvar(int i)
{
......@@ -327,9 +331,9 @@ static inline bool monomial_is_bad(monomial_t m)
void test_bad_monomials()
{
assert(!bad_monomial[0]);
assert(bad_renumbering[0] < 0);
for (int i = 0; i < n; i++) {
assert(!bad_monomial[i + 1]);
assert(bad_renumbering[i + 1] < 0);
}
for (int i = 0; i < v; i++)
for (int j = 0; j < i; j++)
......@@ -365,13 +369,12 @@ void monomial_setup_list()
test_ranking();
log(1, " - Singling out bad monomials\n");
bad_monomial = malloc(Nd[D + 1] * sizeof(*bad_monomial));
if (bad_monomial == NULL)
bad_renumbering = malloc(Nd[D + 1] * sizeof(*bad_renumbering));
if (bad_renumbering == NULL)
err(1, "cannot allocate list of bad monomials");
#pragma omp parallel for reduction(+:nbad)
for (int j = 0; j < Nd[D + 1]; j++) {
bool bad = monomial_is_bad(monomials[j]);
bad_monomial[j] = bad;
bad_renumbering[j] = bad ? nbad : -1;
nbad += bad;
}
log(2, " - %d bad monomials\n", nbad);
......@@ -423,7 +426,7 @@ mzd_t *M; /* dense matrix for degree D-2 multiples */
void matrix_setup()
{
ncols = Nd[D + 1];
ncols = nbad;
row_capacity = 10;
entry_capacity = 100;
Up = malloc((row_capacity + 1) * sizeof(*Up));
......@@ -438,7 +441,7 @@ void matrix_setup()
scratch[j] = 0;
}
/* Eliminate duplicates and good monomials */
/* Eliminate duplicates and good monomials, renumbers columns */
void matrix_store_row(int *row, int size)
{
if (row_capacity < nrows + 1) {
......@@ -453,17 +456,19 @@ void matrix_store_row(int *row, int size)
if (Uj == NULL)
errx(1, "Reallocating entries failed");
}
for (int k = 0; k < size; k++) {
for (int k = 0; k < size; k++) { // scan one, update scratch
int j = row[k];
if (!bad_monomial[j])
continue;
scratch[j] ^= 1;
int jj = bad_renumbering[j];
if (jj < 0)
continue; // skip good monomial
scratch[jj] ^= 1;
}
for (int k = 0; k < size; k++) {
for (int k = 0; k < size; k++) { // scan again, read scratch
int j = row[k];
if (scratch[j]) {
scratch[j] = 0;
Uj[nnz] = j;
int jj = bad_renumbering[j];
if (jj >= 0 && scratch[jj]) {
scratch[jj] = 0;
Uj[nnz] = jj;
nnz += 1;
}
}
......@@ -683,14 +688,26 @@ int main(int argc, char **argv)
hnrows, hncols, hnnz, (double) nnz / nrows, skipped);
}
/* check empty cols */
for (int i = 0; i < nrows; i++)
for (u64 p = Up[i]; p < Up[i+1]; p++)
/* check empty cols & rows */
int nempty_rows = 0;
for (int i = 0; i < nrows; i++) {
if (Up[i] == Up[i+1])
nempty_rows++;
for (u64 p = Up[i]; p < Up[i+1]; p++) {
assert(0 <= Uj[p]);
assert(Uj[p] < ncols);
scratch[Uj[p]] = 1;
int ntouched = 0;
}
}
int ntouched_cols = 0;
for (int j = 0; j < ncols; j++)
ntouched += scratch[j];
log(1, " - %d / %d columns non-empty\n", ntouched, ncols);
ntouched_cols += scratch[j];
log(1, " - %d / %d empty rows / columns\n", ncols - ntouched_cols, nempty_rows);
if (nrows > ncols)
log(0, "Sanity check: %d extra rows (OK)\n", nrows - ncols);
else
errx(1, "ERROR : more columns than rows (%d extra columns)", ncols - nrows);
if (out_filename) {
log(1, "saving...\n");
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment