interpolate-%: interpolate+%.o mgraph+%.o common.o
$(CC) $(CFLAGS) -o $@ $^ $(LIBGSL)
-minimise-%: energy+%.o graph+%.o mgraph+%.o minimise+%.o half+%.o common.o
- $(CXX) $(CXXFLAGS) -o $@ $^ $(LIBGSL)
+minimise-%: energy+%.o graph+%.o mgraph+%.o minimise+%.o \
+ half+%.o parallel.o common.o
+ $(CXX) $(CXXFLAGS) -o $@ $^ $(LIBGSL) -lpthread
# this ridiculous repetition is due to make being too lame
#include "common.h"
#include "minimise.h"
#include "mgraph.h"
+#include "parallel.h"
double vertex_areas[N], vertex_mean_edge_lengths[N], edge_lengths[N][V6];
static double best_energy= DBL_MAX;
static void addcost(double *energy, double tweight, double tcost, int pr);
-#define COST(weight, compute) addcost(&energy, (weight), (compute), printing)
+
+/*---------- main energy computation, weights, etc. ----------*/
+
+typedef double CostComputation(const Vertices vertices, int section);
+
+typedef struct {
+ double weight;
+ CostComputation *fn;
+} CostContribution;
+
+static const CostContribution costs[]= {
+#define COST(weight, compute) { (weight),(compute) },
+
+#if XBITS==3
+#define STOP_EPSILON 1e-6;
+ COST( 3e2, line_bending_cost)
+ COST( 1e3, edge_length_variation_cost)
+ COST( 0.2e3, rim_proximity_cost)
+ COST( 1e8, noncircular_rim_cost)
+#endif
+
+#if XBITS==4
+#define STOP_EPSILON 1e-5;
+ COST( 3e2, line_bending_cost)
+ COST( 3e3, edge_length_variation_cost)
+ COST( 3.8e1, rim_proximity_cost) // 5e1 is too much
+ // 2.5e1 is too little
+ // 0.2e1 grows compared to previous ?
+ // 0.6e0 shrinks compared to previous ?
+ COST( 1e12, noncircular_rim_cost)
+#endif
+};
+
+#define NCOSTS ((sizeof(costs)/sizeof(costs[0])))
void energy_init(void) {
+ stop_epsilon= STOP_EPSILON;
+}
+
+void compute_energy_separately(const struct Vertices *vs,
+ int section, void *energies_v, void *totals_v) {
+ double *energies= energies_v;
+ int ci;
+
+ compute_edge_lengths(vs->a, section);
+ compute_vertex_areas(vs->a, section);
+
+ for (ci=0; ci<NCOSTS; ci++)
+ energies[ci]= costs[ci].fn(vs->a, section);
}
-/*---------- main energy computation and subroutines ----------*/
+/*---------- energy computation machinery ----------*/
+
+void compute_energy_combine(const struct Vertices *vertices,
+ int section, void *energies_v, void *totals_v) {
+ int ci;
+
+ double *energies= energies_v;
+ double *totals= totals_v;
+
+ for (ci=0; ci<NCOSTS; ci++)
+ totals[ci] += energies[ci];
+}
double compute_energy(const struct Vertices *vs) {
static int bests_unprinted;
-
- double energy;
- int printing;
- compute_edge_lengths(vs->a);
- compute_vertex_areas(vs->a);
- energy= 0;
+ double totals[NCOSTS], energy;
+ int ci, printing;
printing= printing_check(pr_cost,0);
if (printing) printf("%15lld c>e |", evaluations);
- if (XBITS==3) {
- COST( 3e2, line_bending_cost(vs->a));
- COST( 1e3, edge_length_variation_cost(vs->a));
- COST( 0.2e3, rim_proximity_cost(vs->a));
- COST( 1e8, noncircular_rim_cost(vs->a));
- stop_epsilon= 1e-6;
- } else if (XBITS==4) {
- COST( 3e2, line_bending_cost(vs->a));
- COST( 3e3, edge_length_variation_cost(vs->a));
- COST( 3.8e1, rim_proximity_cost(vs->a)); // 5e1 is too much
- // 2.5e1 is too little
- // 0.2e1 grows compared to previous ?
- // 0.6e0 shrinks compared to previous ?
- COST( 1e12, noncircular_rim_cost(vs->a));
- stop_epsilon= 1e-5;
- } else {
- abort();
- }
+ inparallel(vs,
+ compute_energy_separately,
+ compute_energy_combine,
+ sizeof(totals) /* really, size of energies */,
+ totals);
+
+ energy= 0;
+
+ for (ci=0; ci<NCOSTS; ci++)
+ addcost(&energy, costs[ci].weight, totals[ci], printing);
if (printing) printf("| total %# e |", energy);
/*---------- Precomputations ----------*/
-void compute_edge_lengths(const Vertices vertices) {
+void compute_edge_lengths(const Vertices vertices, int section) {
int v1,e,v2;
- FOR_EDGE(v1,e,v2)
+ FOR_EDGE(v1,e,v2,OUTER)
edge_lengths[v1][e]= hypotD(vertices[v1],vertices[v2]);
}
-void compute_vertex_areas(const Vertices vertices) {
+void compute_vertex_areas(const Vertices vertices, int section) {
int v0,v1,v2, e1,e2;
// int k;
- FOR_VERTEX(v0) {
+ FOR_VERTEX(v0,OUTER) {
double total= 0.0, edges_total=0;
int count= 0;
* Q,e
*/
-double line_bending_cost(const Vertices vertices) {
+double line_bending_cost(const Vertices vertices, int section) {
static const double axb_epsilon= 1e-6;
static const double exponent_r= 3;
double a[D3], b[D3], axb[D3];
double total_cost= 0;
- FOR_EDGE(qi,e,ri) {
+ FOR_EDGE(qi,e,ri, OUTER) {
pi= EDGE_END2(qi,(e+3)%V6); if (pi<0) continue;
K a[k]= -vertices[pi][k] + vertices[qi][k];
* Q,e
*/
-double edge_length_variation_cost(const Vertices vertices) {
+double edge_length_variation_cost(const Vertices vertices, int section) {
double diff, cost= 0, exponent_r= 2;
int q, e,r, eback;
- FOR_EDGE(q,e,r) {
+ FOR_EDGE(q,e,r, OUTER) {
eback= edge_reverse(q,e);
diff= edge_lengths[q][e] - edge_lengths[q][eback];
cost += pow(diff,exponent_r);
oncircle[1] *= mult;
}
-double rim_proximity_cost(const Vertices vertices) {
+double rim_proximity_cost(const Vertices vertices, int section) {
double oncircle[3], cost=0;
int v;
- FOR_VERTEX(v) {
+ FOR_VERTEX(v, OUTER) {
int y= v >> YSHIFT;
int nominal_edge_distance= y <= Y/2 ? y : Y-1-y;
if (nominal_edge_distance==0) continue;
/*---------- noncircular rim cost ----------*/
-double noncircular_rim_cost(const Vertices vertices) {
+double noncircular_rim_cost(const Vertices vertices, int section) {
int vy,vx,v;
double cost= 0.0;
double oncircle[3];
- FOR_RIM_VERTEX(vy,vx,v) {
+ FOR_RIM_VERTEX(vy,vx,v, OUTER) {
find_nearest_oncircle(oncircle, vertices[v]);
double d2= hypotD2(vertices[v], oncircle);
#include "mgraph.h"
#include "minimise.h"
+#include "parallel.h"
static int sqdistances[N][N];
int v,e, current, future, dfuture;
buf_push= buf_pop= buffer;
- FOR_VERTEX(v) d[v]= -1;
+ FOR_VERTEX(v,INNER) d[v]= -1;
d[start]= 0;
*buf_push++= start;
assert(buf_pop==buf_push);
assert(buf_push <= buffer+sizeof(buffer)/sizeof(buffer[0]));
- FOR_VERTEX(v) {
+ FOR_VERTEX(v,INNER) {
assert(d[v] >= 0);
sqdistances_r[v]= d[v] * d[v];
}
//
}
-void graph_layout_prepare() {
+void graph_layout_prepare(void) {
int v1;
- FOR_VERTEX(v1)
+ FOR_VERTEX(v1,INNER)
breadth_first_search(v1, sqdistances[v1]);
alpha= 2;
}
-double graph_layout_cost(const Vertices v) {
+double graph_layout_cost(const Vertices v, int section) {
/* For each (vi,vj) computes shortest path s_ij = |vi..vj|
* along edges, and actual distance d_ij = |vi-vj|.
*
int v1,v2,e, nedges=0;
double totaledgelength=0, meanedgelength, meanedgelength2;
- FOR_EDGE(v1,e,v2) {
+ FOR_EDGE(v1,e,v2,INNER) {
totaledgelength += hypotD(v[v1], v[v2]);
nedges++;
}
meanedgelength2= meanedgelength * meanedgelength;
// printf("mean=%g mean^2=%g\n", meanedgelength, meanedgelength2);
- FOR_VERTEX(v1) {
- FOR_VERTEX(v2) {
+ FOR_VERTEX(v1,OUTER) {
+ FOR_VERTEX(v2,INNER) {
if (v1 == v2) continue;
double d2= hypotD2(v[v1],v[v2]);
bad= 0;
- FOR_VERTEX(vc)
+ FOR_VERTEX(vc,INNER)
FOR_COORD(kc) {
count[vc][kc]= 0;
ixfrom[vc][kc]= -1;
count[v][k]++;
);
- FOR_VERTEX(vc)
+ FOR_VERTEX(vc,INNER)
FOR_COORD(kc) {
if (count[vc][kc]==1) continue;
printf("BAD %03x %d count=%d ixfrom=%d\n",
#define V6 6
#define V3 3
-#define FOR_VERTEX(v) \
- for ((v)=0; (v)<N; (v)++)
+#define INNER(v,zero,n, precomp) \
+ for ((v)=(zero); precomp, (v)<(n); (v)++)
+
+#define FOR_VERTEX(v,loop) \
+ loop ((v), 0, N, NOTHING)
#define FOR_VPEDGE(e) \
for ((e)=0; (e)<V6; (e)++)
#define FOR_VEDGE(v1,e,v2) \
FOR_VEDGE_X(v1,e,v2,NOTHING,NOTHING)
-#define FOR_EDGE(v1,e,v2) \
- FOR_VERTEX((v1)) \
+#define FOR_EDGE(v1,e,v2, loop) \
+ FOR_VERTEX((v1), loop) \
FOR_VEDGE((v1),(e),(v2))
-#define FOR_RIM_VERTEX(vy,vx,v) \
- for ((vy)=0; (vy)<Y; (vy)+=Y-1) \
- for ((vx)=0; (v)= (vy)<<YSHIFT | (vx), (vx)<X; (vx)++)
+#define FOR_RIM_VERTEX(vy,vx,v, loop) \
+ for ((vy)=0; (vy)<Y; (vy)+=Y-1) \
+ loop ((vx), 0, X, (v)= (vy)<<YSHIFT | (vx))
int vertices_span_join_p(int v0, int v1);
double compute_energy(const struct Vertices *vs);
void energy_init(void);
-double graph_layout_cost(const Vertices v);
+double graph_layout_cost(const Vertices v, int section);
void graph_layout_prepare();
-void compute_vertex_areas(const Vertices vertices);
-void compute_edge_lengths(const Vertices vertices);
+void compute_vertex_areas(const Vertices vertices, int section);
+void compute_edge_lengths(const Vertices vertices, int section);
+
extern double vertex_areas[N], vertex_mean_edge_lengths[N], edge_lengths[N][V6];
-double line_bending_cost(const Vertices vertices);
-double noncircular_rim_cost(const Vertices vertices);
-double edge_length_variation_cost(const Vertices vertices);
-double rim_proximity_cost(const Vertices vertices);
+double line_bending_cost(const Vertices vertices, int section);
+double noncircular_rim_cost(const Vertices vertices, int section);
+double edge_length_variation_cost(const Vertices vertices, int section);
+double rim_proximity_cost(const Vertices vertices, int section);
extern const char *input_file, *best_file;
extern char *best_file_tmp;
static void compute_outvertices(void) {
int v0,k,side,ab,x,y;
- FOR_VERTEX(v0) {
+ FOR_VERTEX(v0, INNER) {
for (ab=0; ab<2; ab++) {
int v1= EDGE_END2(v0, ab?5:0);
int v2= EDGE_END2(v0, ab?0:1);
K Ok(ovAB[v0][ab][1], centroid[k] - normal[k]);
}
}
- FOR_VERTEX(v0) {
+ FOR_VERTEX(v0, INNER) {
int vw= EDGE_END2(v0,3);
int vnw= EDGE_END2(v0,2);
int vsw= EDGE_END2(v0,4);
K Ok(ovC[v0][side], in[v0][k] + adjust[k]);
}
}
- FOR_RIM_VERTEX(y,x,v0) {
+ FOR_RIM_VERTEX(y,x,v0, INNER) {
double rim[D3], inner[D3], radius_cos[D3], radius_sin[D3];
int vback, vfwd, around;
sin(angle) * radius_sin[k]);
}
}
- FOR_RIM_VERTEX(y,x,v0) {
+ FOR_RIM_VERTEX(y,x,v0, INNER) {
int vfwd= EDGE_END2(v0,0);
assert(vfwd >= 0);
int aroung;
static void outfacets(void) {
int v0,e,side,aroung;
- FOR_VERTEX(v0) {
+ FOR_VERTEX(v0, INNER) {
OutVertex *defs=0, *defs1=0;
int rimy=-1;
int_map *defs1aroundmap= 0;
--- /dev/null
+/*
+ * Parallel processing
+ */
+
+#include <pthread.h>
+
+#include "mgraph.h"
+#include "parallel.h"
+
+typedef struct {
+ const struct Vertices *vertices;
+ Computation *separately;
+ void *gendata;
+} ForAllThreads;
+
+typedef struct {
+ ForAllThreads *allthreads;
+ int section;
+ void *secdata;
+ pthread_t thread;
+} PerThread;
+
+static void *routine(void *thread_v) {
+ PerThread *t= thread_v;
+ ForAllThreads *a= t->allthreads;
+
+ a->separately(a->vertices, t->section, t->secdata, a->gendata);
+
+ return 0;
+}
+
+void inparallel(const struct Vertices *vertices,
+ Computation *separately,
+ Computation *combine,
+ size_t secdatasz, void *gendata) {
+ typedef struct { unsigned char secdata[secdatasz]; } SecData;
+
+ ForAllThreads allthreads;
+ SecData secdatas[nsections];
+ PerThread threads[nsections];
+
+ int s, r;
+
+ allthreads.vertices= vertices;
+ allthreads.separately= separately;
+ allthreads.gendata= gendata;
+
+ for (s=0; s<nsections; s++) {
+ threads[s].allthreads= &allthreads;
+ threads[s].section= s;
+ threads[s].secdata= secdatas[s].secdata;
+ r= pthread_create(&threads[s].thread,0,routine,&threads[s]);
+ if (r) diee("pthread_create");
+ }
+
+ for (s=0; s<nsections; s++) {
+ r= pthread_join(threads[s].thread, 0);
+ if (r) diee("pthread_join");
+ combine(vertices, s, threads[s].secdata, gendata);
+ }
+}
--- /dev/null
+/*
+ * Parallel processing
+ */
+
+#ifndef PARALLEL_H
+#define PARALLEL_H
+
+#ifdef NPROCESSORS
+# define NSECTIONS NPROCESSORS
+#else
+# define NSECTIONS 1
+#endif
+
+/* used anamorphically: section, nsections */
+
+#define OUTER_PERSECTION_BASE(zero,n, sect) \
+ ((zero) + sect * (((n)-(zero) + NSECTIONS-1) / NSECTIONS))
+#define OUTER(v,zero,n, precomp) \
+ for ((v)= OUTER_PERSECTION_BASE((zero),(n), section); \
+ precomp, \
+ (v) < OUTER_PERSECTION_BASE((zero),(n), section + 1) && (v) < (n); \
+ (v)++)
+
+#define nsections NSECTIONS
+
+typedef void Computation(const struct Vertices *vertices,
+ int section, void *secdata, void *gendata);
+
+void inparallel(const struct Vertices *vertices,
+ Computation *separately,
+ Computation *combine,
+ size_t secdatasz, void *gendata);
+
+#endif /*PARALLEL_H*/
int v, k;
- FOR_VERTEX(v) {
+ FOR_VERTEX(v,INNER) {
input_gsl.data= &conformation[v][0];
GA( gsl_blas_dgemv(CblasNoTrans, 1.0,&transform_gsl, &input_gsl,
0.0, &result_gsl) );
int vb, ve[V6], e;
ntris= 0;
- FOR_VERTEX(vb) {
+ FOR_VERTEX(vb,INNER) {
/* We use the two triangles in the parallelogram vb, vb+e5, vb+e0, vb+e1.
* We go round each triangle clockwise (although our surface is non-
* orientable so it shouldn't matter). Picking the parallelogram
if (!vertex_in_triangles_checked) {
int v, expd;
- FOR_VERTEX(v) {
+ FOR_VERTEX(v,INNER) {
expd= RIM_VERTEX_P(v) ? 3 : 6;
if (vertex_in_triangles[v] != expd) {
fprintf(stderr,"vertex %02x used for %d triangles, expected %d\n",
static void topocheck(void) {
int v1,e,v2,eprime,v1prime, count;
- FOR_EDGE(v1,e,v2) {
+ FOR_EDGE(v1,e,v2, INNER) {
count= 0;
FOR_VEDGE(v2,eprime,v1prime)
if (v1prime==v1) count++;
~ian / moebius2.git / commitdiff