Commit 91261ba5 authored by Quentin HAMMERER's avatar Quentin HAMMERER
Browse files

le script qui utilise cado-nfs utilise son chemin, donc il faut l'ajuster si...

le script qui utilise cado-nfs utilise son chemin, donc il faut l'ajuster si tu fais tourner sur ta machine, toujours dans ce script, il prend en compte que l'output de macaulay_gen sera 'mat' (oui je n'ai jamais touché au shell'
parent 1b712504
CC = gcc
CFLAGS = -g -Wall -Wextra -O3 -Wno-unused-parameter -Wno-maybe-uninitialized -Werror -march=native
# OpenMP
CFLAGS += -fopenmp
LDFLAGS += -fopenmp
LDLIBS += -lm
all: macaulay_gen reconstruct_monomials
macaulay_gen: parser.o tools.o macaulay_gen.o
macaulay_gen: LDLIBS += `pkg-config --libs m4ri`
reconstruct_monomials: parser.o tools.o reconstruct_monomials.o
reconstruct_monomials: LDLIBS += `pkg-config --libs m4ri`
parser.o: parser.h
tools.o: tools.h
macaulay_gen.o: parser.h tools.h
macaulay_gen.o: CFLAGS += `pkg-config --cflags m4ri`
reconstruct_monomials.o: parser.h tools.h
reconstruct_monomials.o: CFLAGS += `pkg-config --cflags m4ri`
clean:
rm -rf *.o my_mat/ *log.bin *w.bin
rm -f parser_demo moebius_demo parser_demo_alt monica monica_vector macaulay_gen lanczos reconstruct_monomials
import os
D = 3
v = 7
input_file = "large.in"
output = "mat"
new_system_name = "new_system.in"
os.system("./macaulay_gen --in " + input_file + " --out " + output + " --degree " + str(D) + " --inner-hybridation " + str(v))
os.system("sh get_kernel.sh")
os.system("./reconstruct_monomials --in " + input_file + " --out " + new_system_name + " --degree " + str(D) + " --inner-hybridation " + str(v) + " --ilog " + output + ".ilog.bin " + " --jlog " + output + ".jlog.bin")
\ No newline at end of file
#!/bin/bash
echo "Start script"
rm -rf my_mat/ kernel.bin
mkdir my_mat
mv mat.bin my_mat/
mv mat.cw.bin my_mat/
mv mat.rw.bin my_mat/
~/stage/cado-nfs-master/build/poste-almasty-02/linalg/bwc/./bwc.pl :complete balancing_options="reorder=columns" matrix=`pwd`/my_mat/mat.bin mn=64
cp my_mat/mat.bin-1x1/W kernel.bin
echo "End script"
\ No newline at end of file
File added
This diff is collapsed.
// TODO : il faudrait tester ceci
#define _POSIX_C_SOURCE 200809L // for strdup
#include <stdio.h>
#include <stdbool.h>
#include <stdlib.h>
#include <err.h>
#include <getopt.h>
#include <inttypes.h>
#include <assert.h>
#include <string.h> // strcmp, strdup
#include <sys/time.h> // gettimeofday
#include <m4ri/m4ri.h>
#include "parser.h"
#include "tools.h"
/*
* This program generates a degree-D macaulay matrix from a quadratic system,
* where columns corresponding to monomials having degree at most one in certain
* variables have been removed. Finding vectors in the left-kernel of the matrix
* is the basis of the crossbred algorithm.
*/
// #define MAXDEG 8
int n; /* number of variables */
int v; /* number of non-enumerated variables (the first ones) */
int m; /* number of polynomials */
int D;
int Nd[MAXDEG]; /* number of degree-d monomials */
monomial_t *monomials; /* array to rank/unrank monomials */
// monomial_t vmask; /* monomial that contains all special variables */
int nbad; /* total number of bad monomials */
int *bad_renumbering; /* renumber the bad monomials consecutively */
// int *inv_bad_renumbering; /* renumber the bad monomials consecutively */
struct input_poly_t * input_system;
/********************************** output buffering ********************************/
int buffer_capacity;
int buffer_size;
int * buffer;
u64 words_written;
int rows_written;
FILE * buffer_stream;
void output_setup()
{
rows_written = 0;
words_written = 0;
buffer_capacity = 1024*256;
buffer_size = 0;
buffer = alloc_array(buffer_capacity, sizeof(*buffer),"output buffer");
}
void flush_buffer()
{
if (buffer_stream != NULL) {
/* flush the old buffer to disk */
size_t check = fwrite(buffer, sizeof(*buffer), buffer_size, buffer_stream);
if (check != (size_t) buffer_size)
errx(1, "error when writing to disk");
}
buffer_size = 0;
}
/* write a row to the matrix file; return the row number */
int buffer_write_row(int row_size, int *row)
{
int i = -1;
#pragma omp critical(buffer)
{
i = rows_written;
rows_written += 1;
/* getting this out of the critical section requires the
use of an RCU-like mechanism */
if (buffer_size + row_size + 1 < buffer_capacity)
flush_buffer();
buffer[buffer_size] = row_size;
buffer_size += 1;
for (int k = 0; k < row_size; k++)
buffer[buffer_size + k] = row[k];
buffer_size += row_size;
words_written += row_size + 1;
}
return i;
}
/********************************** sparse matrix ***********************************/
mzd_t *M; /* dense matrix for degree D-2 multiples */
/* sparse matrix */
int nrows;
int ncols;
bool *LM; /* leading monomoials of smaller row-reduced macaulay matrix (one per row of M) */
int skipped; /* multiples skipped to avoid trivial linear dependency */
int valid_multipliers; /* number of monomial multipliers used of the final matrix */
int *rw; /* row weights */
int *cw; /* columns weights */
int *ilog; /* k-th row in matrix file is m_j * f_i, with i = ilog[k] and j = jlog(k) */
int *jlog;
/* nrows and ncols have been setup */
void matrix_setup()
{
assert(nrows > 0);
assert(ncols > 0);
skipped = 0;
/* warning, scratch should be per-thread */
rw = alloc_array(nrows, sizeof(*rw), "row weights");
cw = alloc_array(ncols, sizeof(*cw), "col weights");
ilog = alloc_array(nrows, sizeof(*ilog), "ilog");
jlog = alloc_array(nrows, sizeof(*jlog), "jlog");
for (int j = 0; j < ncols; j++)
cw[j] = 0;
for (int i = 0; i < nrows; i++)
rw[i] = -1;
}
void matrix_clean()
{
free(rw);
free(cw);
free(ilog);
free(jlog);
}
/* Store m_j * f_i in the sparse matrix, and update metadata.
Eliminate duplicates and good monomials, renumbers columns.
Scratch should be a zero-filled array. This destroys row. */
void matrix_store_row(int i, int j, int *row, int size, bool *scratch)
{
for (int k = 0; k < size; k++) { // scan one, update scratch
int p = row[k];
int jj = bad_renumbering[p];
row[k] = jj;
if (jj < 0)
continue; // skip good monomial
scratch[jj] ^= 1;
}
int newrow[9000];
int newsize = 0;
for (int k = 0; k < size; k++) { // scan again, read scratch
int jj = row[k];
if (jj >= 0 && scratch[jj]) {
scratch[jj] = 0;
newrow[newsize] = jj;
newsize += 1;
#pragma omp atomic update
cw[jj] += 1;
}
}
int u = buffer_write_row(newsize, newrow);
rw[u] = newsize;
ilog[u] = i;
jlog[u] = j;
}
static inline void mzd_flip_bit (mzd_t *M, rci_t row, rci_t col)
{
mzd_write_bit(M, row, col, mzd_read_bit(M, row, col) ^ 1);
}
static inline int mzd_ffs(mzd_t *M, rci_t row)
{
word const *truerow = mzd_row(M, row);
for (int i = 0; i < M->width; i++) {
if (truerow[i] == 0)
continue;
int j = m4ri_radix * i;
while (!mzd_read_bit(M, row, j))
j += 1;
return j;
}
return M->ncols;
}
void macaulay_dense(int d)
{
valid_multipliers = Nd[d + 1];
if (d < 2)
return;
log(0, "Generating degree-%d (dense) Macaulay matrix \n", d);
log(1, " - matrix size: %d x %d\n", m * Nd[d - 1], Nd[d + 1]);
M = mzd_init(m * Nd[d - 1], Nd[d + 1]);
LM = alloc_array(Nd[d + 1], sizeof(*LM), "leading monomials");
int row = 0;
/* compute m_j * f_i, for all monomial m_j of degree <= d */
for (int j = 0; j < Nd[d - 1]; j++) {
/* generate multiplication table for m_j */
int mmul[9000];
for (int k = 0; k < Nd[3]; k++) {
monomial_t product = monomials[j] | monomials[k];
mmul[k] = monomial_rank(product);
}
for (int i = 0; i < m; i++) {
for (int k = 0; k < input_system[i].nterms; k++) {
int l = mmul[input_system[i].terms[k]];
mzd_flip_bit(M, row, M->ncols - 1 - l);
}
row += 1;
}
}
assert(row == M->nrows);
log(1, " - echelonization\n");
int rank = mzd_echelonize(M, 0);
if (rank < M->nrows)
errx(1, "Bizarre, rank defect!");
for (int i = 0; i < Nd[d + 1]; i++)
LM[i] = 0;
for (int i = 0; i < rank; i++) {
int j = M->ncols - 1 - mzd_ffs(M, i);
assert(j >= 0);
if (j == 0)
errx(1, "Polynomial system has been proven inconsistent!");
if (j < n + 1)
errx(1, "New linear equation found! (decrease D)");
LM[j] = 1;
valid_multipliers -= 1;
}
}
/* multiply by degree-d monomials. May run inside a parallel region */
void macaulay_sparse(int d)
{
bool * scratch = alloc_array(ncols, sizeof(*scratch), "scratch (per thread)");
for (int j = 0; j < ncols; j++)
scratch[j] = 0;
#pragma omp for
for (int j = Nd[d]; j < Nd[d + 1]; j++) {
if (LM && LM[j]) {
#pragma omp atomic update
skipped += 1;
continue;
}
/* generate multiplication table for m_j */
int mmul[9000];
for (int k = 0; k < Nd[3]; k++) {
monomial_t product = monomials[j] | monomials[k];
mmul[k] = monomial_rank(product);
}
for (int i = 0; i < m; i++) {
/* compute and store m_j * f_i */
int row[9000];
int size = input_system[i].nterms;
for (int k = 0; k < size; k++)
row[k] = mmul[input_system[i].terms[k]];
matrix_store_row(i, j, row, size, scratch);
}
}
free(scratch);
}
/************************************* main ***********************************/
int main(int argc, char **argv)
{
// test code, not useful for the actual computation
test_shift_128();
test_HW();
/* process command-line options */
struct option longopts[5] = {
{"in", required_argument, NULL, 'i'},
{"out", required_argument, NULL, 'o'},
{"degree", required_argument, NULL, 'D'},
{"inner-hybridation", required_argument, NULL, 'v'},
{NULL, 0, NULL, 0}
};
char *in_filename = NULL;
char *out_basename = NULL;
v = -1;
signed char ch;
while ((ch = getopt_long(argc, argv, "", longopts, NULL)) != -1) {
switch (ch) {
case 'i':
in_filename = optarg;
break;
case 'o':
out_basename = optarg;
break;
case 'D':
D = atoi(optarg);
break;
case 'v':
v = atoi(optarg);
break;
default:
errx(1, "Unknown option\n");
}
}
/* validate arguments */
if (D == 0)
errx(1, "the --degree argument is mandatory");
if (D < 2)
errx(1, "the --degree must be at least 2");
if (D > MAXDEG)
errx(1, "increase MAXDEG");
if (v < 0)
errx(1, "the --inner-hybridation argument is mandatory");
FILE *input_file = stdin;
if (in_filename) {
input_file = fopen(in_filename, "r");
if (input_file == NULL)
err(1, "Cannot open %s", in_filename);
}
log(0, "Reading input system\n");
parser(input_file, NULL, &parser_setup, &parser_store_monomial, &parser_store_polynomial, NULL);
log(1, " - Read %d polynomials in %d variables\n", m, n);
log(1, " - Generating degree-%d Macaulay matrix\n", D);
log(1, " - Truncating columns linear in the first %d variables\n", v);
monomial_setup_list();
convert_input_system();
log(1, " - %d total monomials of degree <= %d\n", Nd[D+1], D);
/* prepare monomial multiples and setup sparse matrix */
macaulay_dense(D - 2);
log(1, " - %d monomial multipliers\n", valid_multipliers);
nrows = m * valid_multipliers;
ncols = nbad;
if (ncols < nrows)
log(1, " - %d extra rows (OK)\n", nrows - ncols);
else
errx(1, "ERROR : more columns than rows (%d extra columns)", ncols - nrows);
log(0, "Matrix production (pass 1/2)\n");
matrix_setup();
output_setup();
char hnrows[10];
char hncols[10];
human_format(nrows, hnrows);
human_format(ncols, hncols);
log(1, " - Final size: %s rows, %s columns\n", hnrows, hncols);
/* actual stuff -- two passes because of potential empty columns */
for (int d = 2; d <= D; d++) {
log(1, " - Expanding to degree %d\n", d);
#pragma omp parallel
macaulay_sparse(d - 2);
char hnnz[10];
human_format(rows_written, hnrows);
human_format(words_written, hnnz);
log(2, " - %s rows, %s cols, %s nnz (%.0f per row), %d poly skipped\n",
hnrows, hncols, hnnz, (double) words_written / rows_written, skipped);
}
log(0, "Detecting empty columns\n");
/* check empty cols & rows */
int nempty_rows = 0;
u64 nnz = 0;
for (int i = 0; i < nrows; i++) {
if (rw[i] <= 0)
nempty_rows++;
nnz += rw[i];
}
if (nempty_rows == 0)
log(1, " - no empty rows (OK)\n");
else
errx(1, " - %d empty rows (KO)!\n", nempty_rows);
int ncols_before = ncols;
ncols = 0;
u64 nnz_check = 0;
/* check empty columns ; enumerate all monomials again */
for (int j = 0; j < Nd[D + 1]; j++) {
int jj = bad_renumbering[j];
assert(jj < 0 || ncols <= jj);
if (jj < 0)
continue;
bad_renumbering[j] = ncols;
nnz_check += cw[jj];
if (cw[jj] != 0)
ncols += 1;
}
assert(nnz == nnz_check);
assert(ncols <= ncols_before);
log(1, " - %d empty columns\n", ncols_before - ncols);
if (out_basename == NULL) {
log(0, "Not saving matrix because no --out option given\n");
exit(EXIT_SUCCESS);
}
log(1, " - exact dimensions %d x %d (%ld nnz)\n", nrows, ncols, nnz);
char matrix_filename[256];
char rw_filename[256];
char cw_filename[256];
char ilog_filename[256];
char jlog_filename[256];
if (strlen(out_basename) > 240)
errx(1, "output base name too long");
snprintf(matrix_filename, 256, "%s.bin", out_basename);
snprintf(rw_filename, 256, "%s.rw.bin", out_basename);
snprintf(cw_filename, 256, "%s.cw.bin", out_basename);
snprintf(ilog_filename, 256, "%s.ilog.bin", out_basename);
snprintf(jlog_filename, 256, "%s.jlog.bin", out_basename);
log(1, " - saving to %s / %s / %s / %s / %s\n", matrix_filename,
rw_filename, cw_filename, ilog_filename, jlog_filename);
buffer_stream = fopen(matrix_filename, "w");
if (buffer_stream == NULL)
errx(1, "cannot open %s", matrix_filename);
matrix_clean();
log(0, "Matrix production (pass 2/2)\n");
/* output pass 2/2 */
matrix_setup();
output_setup();
for (int d = 2; d <= D; d++) {
log(1, " - Expanding to degree %d\n", d);
#pragma omp parallel
macaulay_sparse(d - 2);
char hnnz[10];
human_format(rows_written, hnrows);
human_format(words_written, hnnz);
log(2, " - %s rows, %s cols, %s nnz (%.0f per row), %d poly skipped\n",
hnrows, hncols, hnnz, (double) words_written / rows_written, skipped);
}
log(0, "Finalization\n");
flush_buffer();
fclose(buffer_stream);
write_array(rw_filename, rw, nrows);
write_array(cw_filename, cw, ncols);
write_array(ilog_filename, ilog, nrows);
write_array(jlog_filename, jlog, nrows);
exit(EXIT_SUCCESS);
}
\ No newline at end of file
0000000000000000
80cf271899cba635
80cf271899cba635
0000000000000000
0000000000000000
80cf271899cba635
0000000000000000
0000000000000000
0000000000000000
80cf271899cba635
80cf271899cba635
80cf271899cba635
80cf271899cba635
80cf271899cba635
80cf271899cba635
80cf271899cba635
0000000000000000
80cf271899cba635
80cf271899cba635
80cf271899cba635
80cf271899cba635
80cf271899cba635
80cf271899cba635
0000000000000000
80cf271899cba635
80cf271899cba635
0000000000000000
0000000000000000
80cf271899cba635
0000000000000000
80cf271899cba635
80cf271899cba635
80cf271899cba635
0000000000000000
0000000000000000
0000000000000000
80cf271899cba635
0000000000000000
80cf271899cba635
0000000000000000
80cf271899cba635
0000000000000000
0000000000000000
0000000000000000
80cf271899cba635
80cf271899cba635
0000000000000000
0000000000000000
0000000000000000
0000000000000000
80cf271899cba635
0000000000000000
0000000000000000
80cf271899cba635
80cf271899cba635
0000000000000000
80cf271899cba635
0000000000000000
80cf271899cba635
80cf271899cba635
80cf271899cba635
0000000000000000
80cf271899cba635
0000000000000000
0000000000000000
80cf271899cba635
80cf271899cba635
80cf271899cba635
0000000000000000
0000000000000000
0000000000000000