X-Git-Url: http://www.chiark.greenend.org.uk/ucgi/~ianmdlvl/git?p=matchsticks-search.git;a=blobdiff_plain;f=main.c;h=57a979654f705bd01912fe773cd01668a92ddb07;hp=5eb9e19d6b0b7caa3324e9c515561b895639a14e;hb=7f7db5029057d9d6cfd74de0f464e1913cb80cff;hpb=1b6dda6020dc3991bf1e642d9bee77cd0b6dff91 diff --git a/main.c b/main.c index 5eb9e19..57a9796 100644 --- a/main.c +++ b/main.c @@ -34,6 +34,7 @@ #include #include #include +#include #include #include #include @@ -41,6 +42,10 @@ #include +#ifndef VERSION +#define VERSION "(unknown-version)" +#endif + /* * Algorithm. * @@ -53,7 +58,8 @@ * * We search all possible adjacency matrices, and for each one we run * GLPK's simplex solver. We represent the adjacency matrix as an - * array of bitmaps. + * array of bitmaps: one word per input stick, with one bit per output + * stick. * * However, there are a couple of wrinkles: * @@ -71,11 +77,15 @@ * nondecreasing in array order. * * Once we have a solution, we also avoid considering any candidate - * which involves dividing one of the output sticks into so many + * which involves dividing one of the input sticks into so many * fragment that the smallest fragment would necessarily be no bigger * than our best solution. That is, we reject candidates where any of * the hamming weights of the adjacency bitmap words are too large. * + * We further winnow the set of possible adjacency matrices, by + * ensuring the same bit is not set in too many entries of adjmatrix + * (ie, as above, only considering output sticks). + * * And, we want to do the search in order of increasing maximum * hamming weight. This is because in practice optimal solutions tend * to have low hamming weight, and having found a reasonable solution @@ -86,6 +96,8 @@ typedef uint32_t AdjWord; #define PRADJ "08"PRIx32 +#define FOR_BITS(j,m) for (j=0, j##bit=1; j < (m); j++, j##bit<<=1) + static int n, m, maxhamweight; static AdjWord *adjmatrix; static AdjWord adjall; @@ -94,6 +106,9 @@ static double best; static glp_prob *best_prob; static AdjWord *best_adjmatrix; +static int n_max_frags, m_max_frags; +static int *weight; + static unsigned printcounter; static void iterate(void); @@ -107,6 +122,22 @@ static void progress_eol(void) { fflush(stderr); } +static void set_best(double new_best) { + best = new_best; + /* + * When computing n_max_frags, we want to set a value that will skip + * anything that won't provide strictly better solutions. So we + * want + * frags < n / best + * _ _ + * <=> frags < | n / best | + * _ _ + * <=> frags <= | n / best | - 1 + */ + n_max_frags = ceil(n / best) - 1; + m_max_frags = ceil(m / best) - 1; +} + /*----- multicore support -----*/ /* @@ -146,10 +177,11 @@ typedef struct { pid_t pid; } Worker; static Worker *mc_us; +static bool mc_am_generator; static void multicore_check_for_new_best(void); -#define MAX_NIOVS 3 +#define MAX_NIOVS 4 static AdjWord mc_iter_min; static int mc_niovs; static size_t mc_iovlen; @@ -170,6 +202,7 @@ static void mc_rwvsetup_outer(void) { IOV(maxhamweight, 1); IOV(mc_iter_min, 1); IOV(*adjmatrix, multicore_iteration_boundary); + IOV(*weight, m); } static void mc_rwvsetup_full(void) { @@ -202,6 +235,8 @@ static void mc_awaitpid(int wnum, pid_t pid) { } static void multicore_outer_iteration(int i, AdjWord min) { + static unsigned check_counter; + assert(i == multicore_iteration_boundary); mc_iter_min = min; mc_rwvsetup_outer(); @@ -209,6 +244,9 @@ static void multicore_outer_iteration(int i, AdjWord min) { assert(r == mc_iovlen); /* effectively, this writev arranges to transfers control * to some worker's instance of iterate_recurse via mc_iterate_worker */ + + if (!(check_counter++ & 0xff)) + multicore_check_for_new_best(); } static void mc_iterate_worker(void) { @@ -231,8 +269,8 @@ static void mc_iterate_worker(void) { iterate_recurse(mc_org_it_bound, mc_iter_min); multicore_iteration_boundary = mc_org_it_bound; } - LPRINTF("worker %2d reporting",mc_us->w); if (best_adjmatrix) { + LPRINTF("worker %2d reporting",mc_us->w); adjmatrix = best_adjmatrix; mc_rwvsetup_full(); ssize_t r = writev(fileno(mc_us->results), mc_iov, mc_niovs); @@ -275,6 +313,7 @@ static void multicore(void) { genpid = fork(); assert(genpid >= 0); if (!genpid) { + mc_am_generator = 1; LPRINTF("generator running"); iterate(); exit(0); @@ -290,13 +329,15 @@ static void multicore(void) { LPRINTF("reading report from %2d",w); ssize_t sr = preadv(fileno(mc_workers[w].results), mc_iov, mc_niovs, 0); if (!sr) continue; + LPRINTF("got report from %2d",w); maxhamweight = 0; optimise(1); } } static void multicore_check_for_new_best(void) { - if (!ncpus) return; + if (!(mc_us || mc_am_generator)) + return; for (;;) { double msg; @@ -304,13 +345,14 @@ static void multicore_check_for_new_best(void) { if (!got) break; assert(got == sizeof(msg)); if (msg > best) - best = msg; + set_best(msg); mc_bus_read += sizeof(msg); } } static void multicore_found_new_best(void) { - if (!ncpus) return; + if (!mc_us) + return; if (mc_us /* might be master */) fprintf(stderr," w%-2d ",mc_us->w); ssize_t wrote = write(mc_bus, &best, sizeof(best)); @@ -328,16 +370,22 @@ static void prep(void) { adjmatrix = xalloc_adjmatrix(); glp_term_out(GLP_OFF); setlinebuf(stderr); + weight = calloc(sizeof(*weight), m); assert(weight); + n_max_frags = INT_MAX; + m_max_frags = INT_MAX; } +#if 0 static AdjWord one_adj_bit(int bitnum) { return (AdjWord)1 << bitnum; } +#endif static int count_set_adj_bits(AdjWord w) { - int j, total; - for (j=0, total=0; j best; + return maxhamweight <= m_max_frags; } static bool preconsider_ok(int nwords, bool doprint) { @@ -358,14 +405,13 @@ static bool preconsider_ok(int nwords, bool doprint) { bool had_max = 0; for (i=0, totalfrags=0; i= maxhamweight); - totalfrags += frags; PRINTF("%"PRADJ" ", adjmatrix[i]); - double maxminsize = (double)m / frags; - if (maxminsize <= best) { + if (frags > m_max_frags) { PRINTF(" too fine"); goto out; } + had_max += (frags >= maxhamweight); + totalfrags += frags; } if (!had_max) { /* Skip this candidate as its max hamming weight is lower than @@ -385,6 +431,7 @@ static void optimise(bool doprint) { /* Consider the best answer (if any) for a given adjacency matrix */ glp_prob *prob = 0; int i, j; + AdjWord jbit; /* * Up to a certain point, optimise() can be restarted. We use this @@ -470,8 +517,8 @@ static void optimise(bool doprint) { glp_set_obj_coef(prob, X_minimum, 1); for (i=0; i= n) { printcounter++; optimise(!(printcounter & 0xfff)); @@ -589,8 +639,20 @@ static void iterate_recurse(int i, AdjWord min) { if (i == 0 && (adjmatrix[i] & (1+adjmatrix[i]))) goto again; + FOR_BITS(j,m) + if (adjmatrix[i] & jbit) + weight[j]++; + for (int j = 0; j < m; j++) + if (weight[j] >= n_max_frags) + goto takeout; + iterate_recurse(i+1, adjmatrix[i]); + takeout: + FOR_BITS(j,m) + if (adjmatrix[i] & jbit) + weight[j]--; + again: if (adjmatrix[i] == adjall) return; @@ -622,7 +684,7 @@ static void report(void) { continue; a[x][y] = min + glp_get_col_prim(best_prob, i); } - printf("%d into %d: min fragment %g\n", n, m, min); + printf("%d into %d: min fragment %g [%s]\n", n, m, min, VERSION); for (i = 0; i < n; i++) { for (j = 0; j < m; j++) { if (a[i][j])