chiark / gitweb /
compiles
authorIan Jackson <ian@davenant.relativity.greenend.org.uk>
Sun, 30 Dec 2007 13:30:17 +0000 (13:30 +0000)
committerIan Jackson <ian@davenant.relativity.greenend.org.uk>
Sun, 30 Dec 2007 13:30:17 +0000 (13:30 +0000)
.bzrignore [new file with mode: 0644]
Makefile
bgl.cpp
common.c
common.h
energy.c
mgraph.h

diff --git a/.bzrignore b/.bzrignore
new file mode 100644 (file)
index 0000000..22710ca
--- /dev/null
@@ -0,0 +1 @@
+minimise
index 9e3aa83..d0e7f2e 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -1,11 +1,16 @@
 
 TARGETS=minimise
 
-CFLAGS=        -Wall -Wwrite-strings -Wpointer-arith -Werror
+CWARNS=        -Wall -Wwrite-strings -Wpointer-arith -Werror
+
+OPTIMISE=      -O2
+CFLAGS=                $(CWARNS) $(OPTIMISE) -g
+CXXFLAGS=      $(CWARNS) $(OPTIMISE) -g
 
 all:           $(TARGETS)
 
 minimise:      energy.o bgl.o common.o mgraph.o
+               $(CXX) $(CXXFLAGS) -o $@ $^ -lgsl -lgslcblas
 
 clean:
                rm -f *.o $(TARGETS) *.new
diff --git a/bgl.cpp b/bgl.cpp
index 1beae8c..ae36ea4 100644 (file)
--- a/bgl.cpp
+++ b/bgl.cpp
@@ -48,12 +48,25 @@ extern "C" {
 
 class Graph { }; // this is a dummy as our graph has no actual representation
 
-struct OutEdgeIncrable {
+using namespace boost;
+
+struct OutEdgeIterator :
+  public iterator_facade<
+    OutEdgeIterator,
+    int const,
+    forward_traversal_tag
+> {
   int f;
-  OutEdgeIncrable& operator++() { f += 1<<ESHIFT; return *this; }
-  OutEdgeIncrable(int v, int e) : f(v | (e << ESHIFT)) { }
+  void increment() { f += 1<<ESHIFT; }
+  bool equal(OutEdgeIterator const& other) const { return f == other.f; }
+  int const& dereference() const { return f; }
+  OutEdgeIterator() { }
+  OutEdgeIterator(int _f) : f(_f) { }
+  OutEdgeIterator(int v, int e) : f(e << ESHIFT | v) { }
 };
  
+typedef counting_iterator<int> VertexIterator;
+
 namespace boost {
   // We make Graph a model of various BGL Graph concepts.
   // This mainly means that graph_traits<Graph> has lots of stuff.
@@ -64,44 +77,44 @@ namespace boost {
     public virtual incidence_graph_tag,
     public virtual vertex_list_graph_tag,
     public virtual edge_list_graph_tag { };
-  
-  struct graph_traits<Graph> {
 
+  struct graph_traits<Graph> {
     // Concept Graph:
-  
     typedef int vertex_descriptor; /* vertex number, -1 => none */
     typedef int edge_descriptor; /* see above */
     typedef undirected_tag directed_category;
     typedef disallow_parallel_edge_tag edge_parallel_category;
     typedef layout_graph_traversal_category traversal_category;
-    inline int null_vertex() { return -1; }
 
     // Concept IncidenceGraph:
-    
-    typedef counting_iterator<OutEdgeIncrable,
-      std::forward_iterator_tag> out_edge_iterator;
+    typedef OutEdgeIterator out_edge_iterator;
     typedef unsigned degree_size_type;
-    
-    inline int source(int f, const Graph&) { return f&VMASK; }
-    inline int target(int f, const Graph&) { return EDGE_END2(f&VMASK, f>>ESHIFT); }
-    inline std::pair<out_edge_iterator,out_edge_iterator>
-    out_edges(int v, const Graph&) {
-      return std::make_pair(out_edge_iterator(OutEdgeIncrable(v, VE_MIN(v))),
-                           out_edge_iterator(OutEdgeIncrable(v, VE_MAX(v))));
-    }
-    inline unsigned out_degree(int v, const Graph&) {
-      return VE_MAX(v) - VE_MIN(v);
-    }
 
     // Concept VertexListGraph:
-    typedef counting_iterator<int> vertex_iterator;
+    typedef VertexIterator vertex_iterator;
     typedef unsigned vertices_size_type;
-    inline std::pair<vertex_iterator,vertex_iterator>
-    vertices(const Graph&) {
-      return std::make_pair(vertex_iterator(0), vertex_iterator(N));
-    }
-    inline unsigned num_vertices(const Graph&) { return N; }
   };
+    
+  // Concept Graph:
+  inline int null_vertex() { return -1; }
+
+  // Concept IncidenceGraph:
+  inline int source(int f, const Graph&) { return f&VMASK; }
+  inline int target(int f, const Graph&) { return EDGE_END2(f&VMASK, f>>ESHIFT); }
+  inline std::pair<OutEdgeIterator,OutEdgeIterator>
+  out_edges(int v, const Graph&) {
+    return std::make_pair(OutEdgeIterator(v, VE_MIN(v)),
+                         OutEdgeIterator(v, VE_MAX(v)));
+  }
+  inline unsigned out_degree(int v, const Graph&) {
+    return VE_MAX(v) - VE_MIN(v);
+  }
+
+  // Concept VertexListGraph:
+  inline std::pair<VertexIterator,VertexIterator> vertices(const Graph&) {
+    return std::make_pair(VertexIterator(0), VertexIterator(N));
+  }
+  inline unsigned num_vertices(const Graph&) { return N; }
 };
 
 static void single_source_shortest_paths(int v1,
@@ -109,7 +122,7 @@ static void single_source_shortest_paths(int v1,
                                         double vertex_distances[/*v*/]) {
   Graph g;
 
-  boost::dijkstra_shortest_paths(g, v1,
+  dijkstra_shortest_paths(g, v1,
      weight_map(edge_weights).
      vertex_index_map(identity_property_map()).
      distance_map(vertex_distances));
@@ -132,15 +145,16 @@ double graph_layout_cost(const Vertices v, const double vertex_areas[N]) {
    * (In practice we compute d^2+epsilon and use it for the
    *  divisions, to avoid division by zero.)
    */
-  static const d2_epsilon= 1e-6;
+  static const double d2_epsilon= 1e-6;
   
-  double edge_weights[N*V6], vertex_distances[N], total_cost;
+  double edge_weights[N*V6], vertex_distances[N], total_cost=0;
   int v1,v2,e,f;
 
-  FOR_VEDGE_X(v1,e,v2,
-             f= v1 | e << ESHIFT,
-             edge_weights[f]= NaN)
-    edge_weights[f]= hypotD(v[v1], v[v2]);
+  FOR_VERTEX(v1)
+    FOR_VEDGE_X(v1,e,v2,
+               f= v1 | e << ESHIFT,
+               edge_weights[f]= NAN)
+      edge_weights[f]= hypotD(v[v1], v[v2]);
 
   FOR_VERTEX(v1) {
     double a1= vertex_areas[v1];
index 47edc4d..9001025 100644 (file)
--- a/common.c
+++ b/common.c
@@ -2,6 +2,9 @@
  * Generally useful stuff.
  */
 
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_blas.h>
+
 #include "common.h"
 
 double magnD(const double pq[D3]) {
@@ -9,10 +12,10 @@ double magnD(const double pq[D3]) {
 
   v.size= D3;
   v.stride= 1;
-  v.vector.data= pq;
+  v.data= (double*)pq;
   /* owner and block ought not to be used */
 
-  return gsl_blas_snrm2(&v);
+  return gsl_blas_dnrm2(&v);
 }
 
 double hypotD(const double p[D3], const double q[D3]) {
@@ -20,22 +23,25 @@ double hypotD(const double p[D3], const double q[D3]) {
   double pq[D3];
 
   K pq[k]= p[k] - q[k];
-  return hypotD1(pq);
+  return magnD(pq);
 }
 
 double hypotD2(const double p[D3], const double q[D3]) {
+  int k;
   double d2= 0;
+  
   K d2= ffsqa(p[k] - q[k], d2);
   return d2;
 }
 
 double hypotD2plus(const double p[D3], const double q[D3], double d2) {
+  int k;
   K d2= ffsqa(p[k] - q[k], d2);
   return d2;
 }
 
 void xprod(double r[D3], const double a[D3], const double b[D3]) {
-  r[0]= a[1]*b[2] - a[2]*b[1]);
-  r[1]= a[2]*b[0] - a[0]*b[2]);
-  r[2]= a[0]*b[1] - a[1]*b[0]);
+  r[0]= a[1]*b[2] - a[2]*b[1];
+  r[1]= a[2]*b[0] - a[0]*b[2];
+  r[2]= a[0]*b[1] - a[1]*b[0];
 }
index 3978b5b..f890f84 100644 (file)
--- a/common.h
+++ b/common.h
@@ -25,6 +25,11 @@ double hypotD2plus(const double p[D3], const double q[D3], double add);
 double magnD(const double pq[D3]);
 void xprod(double r[D3], const double a[D3], const double b[D3]);
 
+#define FOR_COORD(k) \
+  for ((k)=0; (k)<D3; (k)++)
+
+#define K FOR_COORD(k)
+
 #ifdef FP_FAST_FMA
 # define fma_fast fma
 #else
index f28135a..e0ebbe2 100644 (file)
--- a/energy.c
+++ b/energy.c
@@ -257,7 +257,7 @@ static double edgewise_vertex_displacement_cost(const Vertices vertices) {
   static const double l3_epsilon= 1e-6;
 
   int pi,e,qi,ri,si, k;
-  double m[D3], mprime[D3], b, d2, l, sigma_bd2_l3;
+  double m[D3], mprime[D3], b, d2, l, sigma_bd2_l3=0;
 
   FOR_EDGE(pi,e,qi) {
     ri= EDGE_END2(pi,(e+1)%V6); if (ri<0) continue;
index 24b02dc..85577e9 100644 (file)
--- a/mgraph.h
+++ b/mgraph.h
@@ -98,11 +98,6 @@ extern int edge_end2(unsigned v1, int e);
   FOR_VERTEX((v1))                             \
     FOR_VEDGE((v1),(e),(v2))
 
-#define FOR_COORD(k) \
-  for ((k)=0; (k)<D3; (k)++)
-
-#define K FOR_COORD(k)
-
 #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)++)