From: Ian Jackson <ian@davenant.relativity.greenend.org.uk>
Date: Sat, 27 Sep 2008 23:12:47 +0000 (+0100)
Subject: Merge multiprocessor support; fix up Makefile runes
X-Git-Url: http://www.chiark.greenend.org.uk/ucgi/~ian/git?p=moebius2.git;a=commitdiff_plain;h=c9f085069bc4e555689dffedbfa08af250772f2a

Merge multiprocessor support; fix up Makefile runes
---

c9f085069bc4e555689dffedbfa08af250772f2a
diff --cc Makefile
index b41c543,9443d70..927f45d
--- a/Makefile
+++ b/Makefile
@@@ -11,9 -11,12 +11,10 @@@ SGTATHAM=sgtatha
  CWARNS=	-Wall -Wwrite-strings -Wpointer-arith -Werror -Wshadow
  CXXWARNS= $(CWARNS) -Wno-shadow -Wno-error
  
 -NPROCCFLAGS := -DNPROCESSORS=$(shell ./nprocessors)
 -
  OPTIMISE=	-O2
- CFLAGS=		-MMD $(OPTIMISE) -g $(CWARNS) $(DIMCFLAGS)
+ CFLAGS_UNIPROC=	-MMD $(OPTIMISE) -g $(CWARNS) $(DIMCFLAGS)
  CXXFLAGS=	-MMD $(OPTIMISE) -g $(CXXWARNS)
+ CFLAGS=		$(CFLAGS_UNIPROC) $(NPROCCFLAGS)
  
  LIBGSL= -lgsl -lgslcblas
  
@@@ -31,6 -34,12 +32,19 @@@ minimise:	energy.o graph.o common.o mgr
  primer:		primer.o common.o
  		$(CC) $(CFLAGS) -o $@ $^ $(LIBGSL)
  
+ nprocessors:	nprocessors.o common.o
+ 		$(CC) $(CFLAGS_UNIPROC) -o $@ $^ $(LIBGSL)
+ 
+ common.o nprocessors.o: %.o: %.c
+ 		$(CC) -c $(CPPFLAGS) $(CFLAGS_UNIPROC) $< -o $@
+ 
++.nprocessors.make: ./nprocessors
++	set -e; n=`./nprocessors`; \
++	echo "NPROCCFLAGS := -DNPROCESSORS=$$n" $o
++
++include .nprocessors.make
++
++
  prime.data:	primer
  		./$^ $o
  
diff --cc energy.c
index 540f6a9,5aa1c60..e43eced
--- a/energy.c
+++ b/energy.c
@@@ -11,12 -12,69 +12,77 @@@ double vertex_areas[N], vertex_mean_edg
  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)
++#define STOP_EPSILON 1e-6
++    COST(  3e3,   line_bending_cost)
++    COST(  3e3,   edge_length_variation_cost)
++    COST( 0.4e3,  rim_proximity_cost)
++    COST(  1e6,   edge_angle_cost)
++                  #define EDGE_ANGLE_COST_CIRCCIRCRAT (0.5/1.7)
++//    COST(  1e1,   small_triangles_cost)
++    COST(  1e12,   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
++#define STOP_EPSILON 1e-6
++    COST(  3e5,   line_bending_cost)
++    COST( 10e2,   edge_length_variation_cost)
++    COST( 9.0e1,  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)
++
++    COST(  1e12,   edge_angle_cost)
++                  #define EDGE_ANGLE_COST_CIRCCIRCRAT (0.5/1.3)
++    COST(  1e18,   noncircular_rim_cost)
+ #endif
+ };
+ 
+ #define NCOSTS ((sizeof(costs)/sizeof(costs[0])))
+ 
++const double edge_angle_cost_circcircrat= EDGE_ANGLE_COST_CIRCCIRCRAT;
 +
  void energy_init(void) {
+   stop_epsilon= STOP_EPSILON;
  }
  
- /*---------- main energy computation and subroutines ----------*/
+ 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);
+ }
+ 
+ /*---------- 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;
@@@ -92,10 -135,10 +143,10 @@@ static void addcost(double *energy, dou
  
  /*---------- 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)
++  FOR_EDGE(v1,e,v2, OUTER)
      edge_lengths[v1][e]= hypotD(vertices[v1],vertices[v2]);
  }
  
@@@ -103,7 -146,7 +154,7 @@@ void compute_vertex_areas(const Vertice
    int v0,v1,v2, e1,e2;
  //  int k;
  
-   FOR_VERTEX(v0) {
 -  FOR_VERTEX(v0,OUTER) {
++  FOR_VERTEX(v0, OUTER) {
      double total= 0.0, edges_total=0;
      int count= 0;
  
@@@ -165,19 -208,17 +216,19 @@@
     *          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;
 +  static const double exponent_r= 4;
  
    int pi,e,qi,ri, k;
    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;
  
 +//if (!(qi&XMASK)) fprintf(stderr,"%02x-%02x-%02x (%d)\n",pi,qi,ri,e);
 +
      K a[k]= -vertices[pi][k] + vertices[qi][k];
      K b[k]= -vertices[qi][k] + vertices[ri][k];
  
@@@ -261,128 -305,3 +312,128 @@@ double noncircular_rim_cost(const Verti
    }
    return cost;
  }
 +
 +/*---------- overly sharp edge cost ----------*/
 +
 +  /*
 +   *
 +   *   	            Q `-_
 +   *              / |    `-_		               P'Q' ------ S'
 +   *             /  |       `-.		             _,' `.       .
 +   *            /   |           S	          _,'     :      .
 +   *           /    |      _,-'                _,'        :r    .r
 +   *          /     |  _,-'		   R' '           `.   .
 +   *         /    , P '			       ` .   r     :  .
 +   *        /  ,-'			            `  .   : 
 +   *       /,-'				                 ` C'
 +   *      /'
 +   *   	 R
 +   *
 +   *
 +   *
 +   *  Let delta =  angle between two triangles' normals
 +   *
 +   *  Giving energy contribution:
 +   *
 +   *                                   2
 +   *    E             =  F   .  delta
 +   *     vd, edge PQ      vd
 +   */
 +
- double edge_angle_cost(const Vertices vertices, double circcircrat) {
++double edge_angle_cost(const Vertices vertices, int section) {
 +  double pq1[D3], rp[D3], ps[D3], rp_2d[D3], ps_2d[D3], rs_2d[D3];
 +  double a,b,c,s,r;
 +  const double minradius_base= 0.2;
 +
 +  int pi,e,qi,ri,si, k;
 +//  double our_epsilon=1e-6;
 +  double total_cost= 0;
 +
-   FOR_EDGE(pi,e,qi) {
++  FOR_EDGE(pi,e,qi, OUTER) {
 +//    if (!(RIM_VERTEX_P(pi) || RIM_VERTEX_P(qi))) continue;
 +
 +    si= EDGE_END2(pi,(e+V6-1)%V6);  if (si<0) continue;
 +    ri= EDGE_END2(pi,(e   +1)%V6);  if (ri<0) continue;
 +
 +    K {
 +      pq1[k]= -vertices[pi][k] + vertices[qi][k];
 +      rp[k]=  -vertices[ri][k] + vertices[pi][k];
 +      ps[k]=  -vertices[pi][k] + vertices[si][k];
 +    }
 +
 +    normalise(pq1,1,1e-6);
 +    xprod(rp_2d, rp,pq1); /* projects RP into plane normal to PQ */
 +    xprod(ps_2d, ps,pq1); /* likewise PS */
 +    K rs_2d[k]= rp_2d[k] + ps_2d[k];
 +    /* radius of circumcircle of R'P'S' from Wikipedia
 +     * `Circumscribed circle' */
 +    a= magnD(rp_2d);
 +    b= magnD(ps_2d);
 +    c= magnD(rs_2d);
 +    s= 0.5*(a+b+c);
 +    r= a*b*c / sqrt((a+b+c)*(a-b+c)*(b-c+a)*(c-a+b) + 1e-6);
 +
-     double minradius= minradius_base + circcircrat*(a+b);
++    double minradius= minradius_base + edge_angle_cost_circcircrat*(a+b);
 +    double deficit= minradius - r;
 +    if (deficit < 0) continue;
 +    double cost= deficit*deficit;
 +
 +    total_cost += cost;
 +  }
 +
 +  return total_cost;
 +}
 +
 +/*---------- small triangles cost ----------*/
 +
 +  /*
 +   *
 +   *   	            Q `-_
 +   *              / |    `-_
 +   *             /  |       `-.
 +   *            /   |           S
 +   *           /    |      _,-'
 +   *          /     |  _,-'
 +   *         /    , P '
 +   *        /  ,-'
 +   *       /,-'
 +   *      /'
 +   *   	 R
 +   *
 +   *  Let delta =  angle between two triangles' normals
 +   *
 +   *  Giving energy contribution:
 +   *
 +   *                                   2
 +   *    E             =  F   .  delta
 +   *     vd, edge PQ      vd
 +   */
 +
- double small_triangles_cost(const Vertices vertices) {
++double small_triangles_cost(const Vertices vertices, int section) {
 +  double pq[D3], ps[D3];
 +  double x[D3];
 +  int pi,e,qi,si, k;
 +//  double our_epsilon=1e-6;
 +  double total_cost= 0;
 +
-   FOR_EDGE(pi,e,qi) {
++  FOR_EDGE(pi,e,qi, OUTER) {
 +//    if (!(RIM_VERTEX_P(pi) || RIM_VERTEX_P(qi))) continue;
 +
 +    si= EDGE_END2(pi,(e+V6-1)%V6);  if (si<0) continue;
 +
 +    K {
 +      pq[k]= vertices[qi][k] - vertices[pi][k];
 +      ps[k]= vertices[si][k] - vertices[pi][k];
 +    }
 +    xprod(x, pq,ps);
 +
 +    double cost= 1/(magnD2(x) + 0.01);
 +
 +//double cost= pow(magnD(spqxpqr), 3);
 +//assert(dot>=-1 && dot <=1);
 +//double cost= 1-dot;
 +    total_cost += cost;
 +  }
 +
 +  return total_cost;
 +}
diff --cc minimise.h
index a8a2288,097b8ac..c2e13e7
--- a/minimise.h
+++ b/minimise.h
@@@ -10,19 -10,18 +10,22 @@@
  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 edge_angle_cost(const Vertices vertices, double circcircrat);
- double small_triangles_cost(const Vertices vertices);
++extern const double edge_angle_cost_circcircrat;
++
+ 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);
++double edge_angle_cost(const Vertices vertices, int section);
++double small_triangles_cost(const Vertices vertices, int section);
  
  extern const char *input_file, *best_file;
  extern char *best_file_tmp;
diff --cc primer.c
index 70cbcfe,70cbcfe..76939a6
--- a/primer.c
+++ b/primer.c
@@@ -17,7 -17,7 +17,7 @@@ int main(int argc, const char **argv) 
  	 X*Y, N, X, Y, D3,
  	 prec+5,prec);
  
--  FOR_VERTEX(vi) {
++  FOR_VERTEX(vi, INNER) {
      int x= vi & XMASK; /* distance along strip */
      int y= vi >> YSHIFT; /* distance across strip */
      double u= y * 1.0 / (Y-1); /* SGT's u runs 0..1 across the strip */