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.
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,
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));
* (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];
* Generally useful stuff.
*/
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_blas.h>
+
#include "common.h"
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]) {
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];
}
~ian / moebius2.git / commitdiff