diff --git a/knik/src/CMakeLists.txt b/knik/src/CMakeLists.txt index 0fa21064..e4019571 100644 --- a/knik/src/CMakeLists.txt +++ b/knik/src/CMakeLists.txt @@ -2,7 +2,7 @@ # include ( ${QT_USE_FILE} ) include_directories ( ${KNIK_SOURCE_DIR}/src - ${KNIK_SOURCE_DIR}/src/flute-2.4/src + ${KNIK_SOURCE_DIR}/src/flute-3.1/src ${HURRICANE_INCLUDE_DIR} ${CORIOLIS_INCLUDE_DIR} ${UTILITIES_INCLUDE_DIR} @@ -39,13 +39,29 @@ KnikEngine.cpp GraphicKnikEngine.cpp ) - set ( fluteIncludes flute-2.4/src/knik/flute.h ) - set ( fluteCpps flute-2.4/src/flute.cpp ) + set ( fluteIncludes flute-3.1/src/knik/flute.h + flute-3.1/src/knik/dl.h + flute-3.1/src/knik/mst2.h + flute-3.1/src/knik/err.h + flute-3.1/src/knik/heap.h + flute-3.1/src/knik/dist.h + flute-3.1/src/knik/global.h + flute-3.1/src/knik/neighbors.h + ) + set ( fluteCpps flute-3.1/src/flute.cpp + flute-3.1/src/flute_mst.cpp + flute-3.1/src/dist.cpp + flute-3.1/src/dl.cpp + flute-3.1/src/err.cpp + flute-3.1/src/mst2.cpp + flute-3.1/src/heap.cpp + flute-3.1/src/neighbors.cpp + ) qtX_wrap_cpp ( mocCpps ${mocIncludes} ) add_library ( flute ${fluteCpps} ) - set_target_properties ( flute PROPERTIES VERSION 2.4 SOVERSION 2 ) + set_target_properties ( flute PROPERTIES VERSION 3.1 SOVERSION 3 ) target_link_libraries ( flute ${HURRICANE_LIBRARIES} ${CORIOLIS_LIBRARIES} ) @@ -71,5 +87,5 @@ install ( FILES ${includes} ${mocIncludes} ${fluteIncludes} DESTINATION include/coriolis2/knik ) - install ( FILES flute-2.4/etc/POST9.dat - flute-2.4/etc/POWV9.dat DESTINATION share/coriolis2/flute-2.4 ) + install ( FILES flute-3.1/etc/POST9.dat + flute-3.1/etc/POWV9.dat DESTINATION share/coriolis2/flute-3.1 ) diff --git a/knik/src/Graph.cpp b/knik/src/Graph.cpp index ac0b65ff..433e008e 100644 --- a/knik/src/Graph.cpp +++ b/knik/src/Graph.cpp @@ -39,7 +39,6 @@ #define __USE_MATRIXVERTEX__ #define EPSILON 10e-4 -#define FLUTE_LIMIT 150 #define HISTORIC_INC 1.5 // define the increment of historic cost for ripup & reroute namespace Knik { @@ -1547,13 +1546,6 @@ void Graph::UpdateEstimateCongestion ( bool create ) { if ( _vertexes_to_route.size() < 2 ) return; - if ( _vertexes_to_route.size() >= FLUTE_LIMIT ) { - if ( create ) - cerr << Warning( "Graph::UpdateEstimateCongestion(): %s\n" - " has more than %d vertex/terminals and cannot be handled by FLUTE." - , getString(_working_net).c_str(), FLUTE_LIMIT ) << endl; - return; - } //cerr << "Running FLUTE for net : " << _working_net << endl; auto_ptr flutetree ( createFluteTree() ); diff --git a/knik/src/flute-2.4/LICENSE.txt b/knik/src/flute-3.1/LICENSE.txt similarity index 100% rename from knik/src/flute-2.4/LICENSE.txt rename to knik/src/flute-3.1/LICENSE.txt diff --git a/knik/src/flute-2.4/etc/POST9.dat b/knik/src/flute-3.1/etc/POST9.dat similarity index 100% rename from knik/src/flute-2.4/etc/POST9.dat rename to knik/src/flute-3.1/etc/POST9.dat diff --git a/knik/src/flute-2.4/etc/POWV9.dat b/knik/src/flute-3.1/etc/POWV9.dat similarity index 100% rename from knik/src/flute-2.4/etc/POWV9.dat rename to knik/src/flute-3.1/etc/POWV9.dat diff --git a/knik/src/flute-3.1/src/dist.cpp b/knik/src/flute-3.1/src/dist.cpp new file mode 100644 index 00000000..06f51b75 --- /dev/null +++ b/knik/src/flute-3.1/src/dist.cpp @@ -0,0 +1,44 @@ +#include "knik/global.h" + +/*********************************************************************/ +/* + Return the Manhattan distance between two points +*/ + +long dist( + Point p, + Point q +) +{ + long dx, dy; + + dx = (p.x) - (q.x); + if( dx < 0 ) dx = -dx; + dy = (p.y) - (q.y); + if( dy < 0 ) dy = -dy; + + return dx + dy; +} + +/*********************************************************************/ +/* + Return the Manhattan distance between two points +*/ + +long dist2( + Point* p, + Point* q +) +{ + long dx, dy; + + dx = (p->x) - (q->x); + if( dx < 0 ) dx = -dx; + dy = (p->y) - (q->y); + if( dy < 0 ) dy = -dy; + + return dx + dy; +} + +/*********************************************************************/ +/*********************************************************************/ diff --git a/knik/src/flute-3.1/src/dl.cpp b/knik/src/flute-3.1/src/dl.cpp new file mode 100644 index 00000000..350dc953 --- /dev/null +++ b/knik/src/flute-3.1/src/dl.cpp @@ -0,0 +1,161 @@ +#include "knik/dl.h" +#include +#include + +dl_t dl_alloc() +{ + dl_t dl = (dl_t)malloc(sizeof(dl_s)); + if (!dl) { + printf("Out of memory!!\n"); + } else { + dl->first = dl->last = 0; dl->count = 0; + } + return dl; +} + +void dl_delete(dl_t dl, dl_el *el) +{ + if (dl->first == el) { + dl->first = el->next; + } + if (dl->last == el) { + dl->last = el->prev; + } + if (el->next) { + el->next->prev = el->prev; + } + if (el->prev) { + el->prev->next = el->next; + } + free(el); dl->count--; +} + +void dl_clear(dl_t dl) +{ + dl_el *el, *next; + if (dl->count > 0) { + for (el=dl->first; el; el=next) { + next = el->next; + free(el); + } + } + dl->first = dl->last = 0; + dl->count = 0; +} + +void dl_concat(dl_t first_list, dl_t second_list) +{ + if (first_list->count <= 0) { + *first_list = *second_list; + } else if (second_list->count > 0) { + first_list->last->next = second_list->first; + second_list->first->prev = first_list->last; + first_list->last = second_list->last; + first_list->count += second_list->count; + } + + free(second_list); +} + +static void dl_insertion_sort(dl_t dl, size_t el_size, + int(*compar)(void *, void *)) +{ + char *buf; + void *curr_d, *srch_d; + dl_el *curr, *srch; + + if (dl_length(dl) <= 1) { + return; + } + + buf = (char*)malloc(el_size); + + for (curr=dl->first; curr!=dl->last; curr=curr->next) { + curr_d = (void*)(((dl_el*)curr)+1); + + for (srch=dl->last; srch!=curr; srch=srch->prev) { + srch_d = (void*)(((dl_el*)srch)+1); + if (compar(curr_d, srch_d) > 0) { + memcpy((void*)buf, curr_d, el_size); + memcpy(curr_d, srch_d, el_size); + memcpy(srch_d, (void*)buf, el_size); + } + } + } + + + free(buf); +} + +void dl_sort(dl_t dl, size_t el_size, int(*compar)(void *, void *)) +{ + dl_el *el, *first_head, *second_head; + dl_s first_list, second_list; + void *first_item, *second_item; + int i, len; + + if (dl_length(dl) <= 25) { + dl_insertion_sort(dl, el_size, compar); + return; + } + + len = dl_length(dl)/2; + for (i=0, el=dl->first; inext; + } + + first_list.first = dl->first; + first_list.last = el->prev; + first_list.count = len; + first_list.last->next = 0; + + second_list.first = el; + second_list.last = dl->last; + second_list.count = dl_length(dl)-len; + second_list.first->prev = 0; + + dl_sort(&first_list, el_size, compar); + dl_sort(&second_list, el_size, compar); + + /* in-place merging */ + first_head = first_list.first; + second_head = second_list.first; + + first_item = (void*)(((dl_el*)first_head)+1); + second_item = (void*)(((dl_el*)second_head)+1); + if (compar(first_item, second_item) <= 0) { + dl->first = el = first_head; + first_head = first_head->next; + } else { + dl->first = el = second_head; + second_head = second_head->next; + } + + while (1) { + first_item = (void*)(((dl_el*)first_head)+1); + second_item = (void*)(((dl_el*)second_head)+1); + if (compar(first_item, second_item) <= 0) { + el->next = first_head; + first_head->prev = el; + el = first_head; + first_head = first_head->next; + if (!first_head) { + el->next = second_head; + second_head->prev = el; + dl->last = second_list.last; + break; + } + } else { + el->next = second_head; + second_head->prev = el; + el = second_head; + second_head = second_head->next; + if (!second_head) { + el->next = first_head; + first_head->prev = el; + dl->last = first_list.last; + break; + } + } + } +} diff --git a/knik/src/flute-3.1/src/err.cpp b/knik/src/flute-3.1/src/err.cpp new file mode 100644 index 00000000..59cb4ff1 --- /dev/null +++ b/knik/src/flute-3.1/src/err.cpp @@ -0,0 +1,28 @@ +#include +#include + +/**************************************************************************/ +/* + print error message and continue +*/ + +void err_msg( +char* msg +) +{ + fprintf(stderr, "%s\n", msg); +} + +/**************************************************************************/ +/* + print error message and exit +*/ + +void err_exit( +char* msg +) +{ + fprintf(stderr, "%s\n", msg); + exit(1); +} + diff --git a/knik/src/flute-2.4/src/flute.cpp b/knik/src/flute-3.1/src/flute.cpp similarity index 86% rename from knik/src/flute-2.4/src/flute.cpp rename to knik/src/flute-3.1/src/flute.cpp index b9dddf3b..ebb734ef 100644 --- a/knik/src/flute-2.4/src/flute.cpp +++ b/knik/src/flute-3.1/src/flute.cpp @@ -17,11 +17,6 @@ using CRL::AllianceFramework; #include #include "knik/flute.h" -#define max(x,y) ((x)>(y)?(x):(y)) -#define min(x,y) ((x)<(y)?(x):(y)) -#define abs(x) ((x)<0?(-x):(x)) -#define ADIFF(x,y) ((x)>(y)?(x-y):(y-x)) // Absolute difference - #if D<=7 #define MGROUP 5040/4 // Max. # of groups, 7! = 5040 #define MPOWV 15 // Max. # of POWVs per group @@ -37,13 +32,19 @@ int numgrp[10]={0,0,0,0,6,30,180,1260,10080,90720}; struct csoln { unsigned char parent; - unsigned char seg[11]; // add: 0..i, Sub: j..10; seg[i+1]=seg[j-1]=0 + unsigned char seg[11]; // Add: 0..i, Sub: j..10; seg[i+1]=seg[j-1]=0 unsigned char rowcol[D-2]; // row = rowcol[]/16, col = rowcol[]%16, unsigned char neighbor[2*D-2]; }; struct csoln *LUT[D+1][MGROUP]; // storing 4 .. D int numsoln[D+1][MGROUP]; +struct point +{ + DTYPE x, y; + int o; +}; + void readLUT(); DTYPE flute_wl(int d, DTYPE x[], DTYPE y[], int acc); DTYPE flutes_wl_LD(int d, DTYPE xs[], DTYPE ys[], int s[]); @@ -56,10 +57,12 @@ FTree flutes_RDP(int d, DTYPE xs[], DTYPE ys[], int s[], int acc); FTree dmergetree(FTree t1, FTree t2); FTree hmergetree(FTree t1, FTree t2, int s[]); FTree vmergetree(FTree t1, FTree t2); +void local_refinement(FTree *tp, int p); DTYPE wirelength(FTree t); void printtree(FTree t); void plottree(FTree t); + void readLUT() { unsigned char charnum[256], line[32], *linep, c; @@ -151,15 +154,13 @@ void readLUT() fclose ( fprt ); } + DTYPE flute_wl(int d, DTYPE x[], DTYPE y[], int acc) { DTYPE xs[MAXD], ys[MAXD], minval, l, xu, xl, yu, yl; int s[MAXD]; - int i, j, /*k,*/ minidx; - struct point { - DTYPE x, y; - int o; - } pt[MAXD], *ptp[MAXD], *tmpp; + int i, j, k, minidx; + struct point pt[MAXD], *ptp[MAXD], *tmpp; if (d==2) l = ADIFF(x[0], x[1]) + ADIFF(y[0], y[1]); @@ -248,7 +249,7 @@ DTYPE flute_wl(int d, DTYPE x[], DTYPE y[], int acc) // xs[] and ys[] are coords in x and y in sorted order // s[] is a list of nodes in increasing y direction // if nodes are indexed in the order of increasing x coord -// i.e., s[i] = s_i in defined in paper +// i.e., s[i] = s_i as defined in paper // The points are (xs[s[i]], ys[i]) for i=0..d-1 // or (xs[i], ys[si[i]]) for i=0..d-1 @@ -405,13 +406,13 @@ DTYPE flutes_wl_MD(int d, DTYPE xs[], DTYPE ys[], int s[], int acc) if (lb < 2) lb = 2; ub=d-1-lb; -// compute scores +// Compute scores #define AAWL 0.6 #define BBWL 0.3 float CCWL = 7.4/((d+10.)*(d-3.)); float DDWL = 4.8/(d-1); - // compute penalty[] + // Compute penalty[] dx = CCWL*(xs[d-2]-xs[1]); dy = CCWL*(ys[d-2]-ys[1]); for (r = d/2, pnlty = 0; r>=0; r--, pnlty += dx) @@ -422,7 +423,7 @@ DTYPE flutes_wl_MD(int d, DTYPE xs[], DTYPE ys[], int s[], int acc) // for (r=0; rx < pb->x) return -1; + if (pa->x > pb->x) return 1; + return 0; +} + +static int ordery(const void *a, const void *b) +{ + struct point *pa, *pb; + + pa = *(struct point**)a; + pb = *(struct point**)b; + + if (pa->y < pb->y) return -1; + if (pa->y > pb->y) return 1; + return 0; +} + FTree flute(int d, DTYPE x[], DTYPE y[], int acc) { - DTYPE xs[MAXD], ys[MAXD], minval; - int s[MAXD]; - int i, j, /*k,*/ minidx; - struct point { - DTYPE x, y; - int o; - } pt[MAXD], *ptp[MAXD], *tmpp; + DTYPE *xs, *ys, minval; + int *s; + int i, j, k, minidx; + struct point *pt, **ptp, *tmpp; FTree t; if (d==2) { @@ -591,27 +613,37 @@ FTree flute(int d, DTYPE x[], DTYPE y[], int acc) t.branch[1].n = 1; } else { + xs = (DTYPE *)malloc(sizeof(DTYPE)*(d)); + ys = (DTYPE *)malloc(sizeof(DTYPE)*(d)); + s = (int *)malloc(sizeof(int)*(d)); + pt = (struct point *)malloc(sizeof(struct point)*(d+1)); + ptp = (struct point **)malloc(sizeof(struct point*)*(d+1)); + for (i=0; ix; - minidx = i; - for (j=i+1; j ptp[j]->x) { - minval = ptp[j]->x; - minidx = j; + if (d<200) { + for (i=0; ix; + minidx = i; + for (j=i+1; j ptp[j]->x) { + minval = ptp[j]->x; + minidx = j; + } } + tmpp = ptp[i]; + ptp[i] = ptp[minidx]; + ptp[minidx] = tmpp; } - tmpp = ptp[i]; - ptp[i] = ptp[minidx]; - ptp[minidx] = tmpp; + } else { + qsort(ptp, d, sizeof(struct point *), orderx); } - + #if REMOVE_DUPLICATE_PIN==1 ptp[d] = &pt[d]; ptp[d]->x = ptp[d]->y = -999999; @@ -632,31 +664,46 @@ FTree flute(int d, DTYPE x[], DTYPE y[], int acc) } // sort y to find s[] - for (i=0; iy; - minidx = i; - for (j=i+1; j ptp[j]->y) { - minval = ptp[j]->y; - minidx = j; + if (d<200) { + for (i=0; iy; + minidx = i; + for (j=i+1; j ptp[j]->y) { + minval = ptp[j]->y; + minidx = j; + } } + ys[i] = ptp[minidx]->y; + s[i] = ptp[minidx]->o; + ptp[minidx] = ptp[i]; + } + ys[d-1] = ptp[d-1]->y; + s[d-1] = ptp[d-1]->o; + } else { + qsort(ptp, d, sizeof(struct point *), ordery); + for (i=0; iy; + s[i] = ptp[i]->o; } - ys[i] = ptp[minidx]->y; - s[i] = ptp[minidx]->o; - ptp[minidx] = ptp[i]; } - ys[d-1] = ptp[d-1]->y; - s[d-1] = ptp[d-1]->o; t = flutes(d, xs, ys, s, acc); + + free(xs); + free(ys); + free(s); + free(pt); + free(ptp); } + return t; } // xs[] and ys[] are coords in x and y in sorted order // s[] is a list of nodes in increasing y direction // if nodes are indexed in the order of increasing x coord -// i.e., s[i] = s_i in defined in paper +// i.e., s[i] = s_i as defined in paper // The points are (xs[s[i]], ys[i]) for i=0..d-1 // or (xs[i], ys[si[i]]) for i=0..d-1 @@ -917,13 +964,13 @@ FTree flutes_MD(int d, DTYPE xs[], DTYPE ys[], int s[], int acc) if (lb < 2) lb = 2; ub=d-1-lb; -// compute scores +// Compute scores #define AA 0.6 // 2.0*BB #define BB 0.3 float CC = 7.4/((d+10.)*(d-3.)); float DD = 4.8/(d-1); - // compute penalty[] + // Compute penalty[] dx = CC*(xs[d-2]-xs[1]); dy = CC*(ys[d-2]-ys[1]); for (r = d/2, pnlty = 0; r>=2; r--, pnlty += dx) @@ -939,7 +986,7 @@ FTree flutes_MD(int d, DTYPE xs[], DTYPE ys[], int s[], int acc) // for (r=0; rdeg; + root = tp->branch[p].n; + +// Reverse edges to point to root + prev = root; + curr = tp->branch[prev].n; + next = tp->branch[curr].n; + while (curr != next) { + tp->branch[curr].n = prev; + prev = curr; + curr = next; + next = tp->branch[curr].n; + } + tp->branch[curr].n = prev; + tp->branch[root].n = root; + +// Find Steiner nodes that are at pins + for (i=d; i<=2*d-3; i++) + SteinerPin[i] = -1; + for (i=0; ibranch[i].n; + if (tp->branch[i].x == tp->branch[next].x && + tp->branch[i].y == tp->branch[next].y) + SteinerPin[next] = i; // Steiner 'next' at Pin 'i' + } + SteinerPin[root] = p; + +// Find pins that are directly connected to root + dd = 0; + for (i=0; ibranch[i].n; + if (SteinerPin[curr] == i) + curr = tp->branch[curr].n; + while (SteinerPin[curr] < 0) + curr = tp->branch[curr].n; + if (curr == root) { + x[dd] = tp->branch[i].x; + if (SteinerPin[tp->branch[i].n] == i && tp->branch[i].n != root) + index[dd++] = tp->branch[i].n; // Steiner node + else index[dd++] = i; // Pin + } + } + + if (4 <= dd && dd <= D) { +// Find Steiner nodes that are directly connected to root + ii=dd; + for (i=0; ibranch[index[i]].n; + while (SteinerPin[curr] < 0) { + index[ii++] = curr; + SteinerPin[curr] = INT_MAX; + curr = tp->branch[curr].n; + } + } + index[ii] = root; + + for (ii=0; iibranch[index[ii]].y; + } + + tt = flutes_LD(dd, xs, ys, ss); + +// Find new wirelength + tp->length += tt.length; + for (ii=0; ii<2*dd-3; ii++) { + i = index[ii]; + j = tp->branch[i].n; + tp->length -= ADIFF(tp->branch[i].x, tp->branch[j].x) + + ADIFF(tp->branch[i].y, tp->branch[j].y); + } + +// Copy tt into t + for (ii=0; iibranch[index[ii]].n = index[tt.branch[ii].n]; + } + for (; ii<=2*dd-3; ii++) { + tp->branch[index[ii]].x = tt.branch[ii].x; + tp->branch[index[ii]].y = tt.branch[ii].y; + tp->branch[index[ii]].n = index[tt.branch[ii].n]; + } + free(tt.branch); + } + + return; +} + DTYPE wirelength(FTree t) { int i, j; diff --git a/knik/src/flute-3.1/src/flute_mst.cpp b/knik/src/flute-3.1/src/flute_mst.cpp new file mode 100644 index 00000000..c585ac77 --- /dev/null +++ b/knik/src/flute-3.1/src/flute_mst.cpp @@ -0,0 +1,2780 @@ +#include +#include +#include +#include +#include +#include +#include "knik/dl.h" +#include "knik/flute.h" +#include "knik/mst2.h" + +#define INFNTY INT_MAX + +#define D2M D2(1) // Max net degree that flute_mr will handle +#define MR_FOR_SMALL_CASES_ONLY 1 +#if MR_FOR_SMALL_CASES_ONLY +#define MAXPART D2M // max partition of an MST +#define MAXT (D2M/D*2) +#else +#define MAXPART (d/9*2) //(MAXD/THD*2) // max partition of an MST +#define MAXPART2 ((t1.deg+t2.deg)/9*2) +#define MAXT (d/5) +#endif + +int D3=INFNTY; + +int FIRST_ROUND=2; // note that num of total rounds = 1+FIRST_ROUND +int EARLY_QUIT_CRITERIA=1; + +#define DEFAULT_QSIZE (3+min(d,1000)) + +#define USE_HASHING 1 +#if USE_HASHING +#define new_ht 1 +//int new_ht=1; +dl_t ht[D2M+1]; // hash table of subtrees indexed by degree +#endif + +unsigned int curr_mark=0; + +FTree wmergetree(FTree t1, FTree t2, int *order1, int *order2, DTYPE cx, DTYPE cy, int acc); +FTree xmergetree(FTree t1, FTree t2, int *order1, int *order2, DTYPE cx, DTYPE cy); +void color_tree(FTree t, int *color); +int longest_path_edge(int i, int j, int *e, int *p, int *es); +void preprocess_edges(int num_edges, int *edges, DTYPE *len, + int *e, int *p, int *es); + +#define init_queue(q) { q[1] = 2; } +inline void enqueue(int **q, int e) +{ + int _qsize; + if ((*q)[0]==(*q)[1]) { + _qsize=2*((*q)[0]+1); + (*q)=(int*)realloc((*q), _qsize*sizeof(int)); + (*q)[0]=_qsize; + } + (*q)[(*q)[1]++]=e; +} + +#define MAX_HEAP_SIZE (MAXD*2) +DTYPE **hdist; +typedef struct node_pair_s { // pair of nodes representing an edge + int node1, node2; +} node_pair; +node_pair *heap; //heap[MAXD*MAXD]; +int heap_size=0; +int max_heap_size = MAX_HEAP_SIZE; + +int in_heap_order(int e1, int e2) +{ + if (hdist[heap[e1].node1][heap[e1].node2] < + hdist[heap[e2].node1][heap[e2].node2]) { + return 1; + } else { + return 0; + } +} + +void sift_up(int i) +{ + node_pair tmp; + int j; + + for (j=i/2; j>=1 && in_heap_order(i, j); i=j, j/=2) { + tmp = heap[j]; + heap[j] = heap[i]; + heap[i] = tmp; + } +} + +void sift_down(int i) +{ + int left, right, j; + node_pair tmp; + + left = i*2; + right = left+1; + + while (left <= heap_size) { + if (left == heap_size || in_heap_order(left, right)) { + j = left; + } else { + j = right; + } + if (in_heap_order(j, i)) { + tmp = heap[j]; + heap[j] = heap[i]; + heap[i] = tmp; + i = j; + left = i*2; + right = left+1; + } else { + break; + } + } +} + +void insert_heap(node_pair *np) +{ + if (heap_size >= max_heap_size) { + max_heap_size *= 2; + heap = (node_pair*)realloc(heap, sizeof(node_pair)*(max_heap_size+1)); + } + heap[++heap_size] = *np; + sift_up(heap_size); +} + +void extract_heap(node_pair *np) +{ + // caller has to make sure heap is not empty + *np = heap[1]; + heap[1] = heap[heap_size--]; + sift_down(1); +} + +void init_param() +{ + int i; + + heap = (node_pair*)malloc(sizeof(node_pair)*(max_heap_size+1)); +} + +FTree reftree; // reference for qsort +int cmp_branch(const void *a, const void *b) { + int n; + DTYPE x1, x2, x3; + + x1 = reftree.branch[*(int*)a].x; + n = reftree.branch[*(int*)a].n; + x3 = reftree.branch[n].x; + if (x3 < x1) x1=x3; + + x2 = reftree.branch[*(int*)b].x; + n = reftree.branch[*(int*)b].n; + x3 = reftree.branch[n].x; + if (x3 < x2) x2=x3; + + return (x1 <= x2) ? -1 : 1; +} + +void update_dist2(FTree t, DTYPE **dist, DTYPE longest, + int *host, int *min_node1, int *min_node2, + int **nb) +{ + int i, j, m, n, dd, node1, node2, node3, node4, p1, p2, pi, pn; + DTYPE min_dist, smallest; + DTYPE x1, x2, x3, x4, y1, y2, y3, y4; + DTYPE threshold_x, threshold_y; + DTYPE md = dist[*min_node1][*min_node2]; + +#if MR_FOR_SMALL_CASES_ONLY + int isPin_base[D2M], *isPin, id[2*D2M]; + int u, v, b[D2M*2]; +#else + int *isPin_base, *isPin, *id; + int u, v, *b; + + isPin_base = (int*)malloc(sizeof(int)*t.deg); + id = (int*)malloc(sizeof(int)*t.deg*2); + b = (int*)malloc(sizeof(int)*t.deg*2); +#endif + + isPin = &(isPin_base[0]) - t.deg; + dd = t.deg*2-2; + + for (i=0; i=0 || n==i) { + continue; + } + + if (id[i] < id[n]) { + id[n] = id[i]; + } else { + id[i] = id[n]; + } + + } + + for (i=0; ithreshold_x || + ADIFF(t.branch[i].y, t.branch[j].y)>threshold_y) + continue; + m = t.branch[j].n; + node3 = host[j]; + node4 = host[m]; + if (node3 < 0 && node4 < 0) { + continue; + } + + if (t.branch[j].x <= t.branch[m].x) { + x1 = t.branch[j].x; + x2 = t.branch[m].x; + } else { + x1 = t.branch[m].x; + x2 = t.branch[j].x; + } + if (t.branch[j].y <= t.branch[m].y) { + y1 = t.branch[j].y; + y2 = t.branch[m].y; + } else { + y1 = t.branch[m].y; + y2 = t.branch[j].y; + } + + if (x2 < x3) { + min_dist = x3 - x2; + } else if (x4 < x1) { + min_dist = x1 - x4; + } else { + min_dist = 0; + } + + if (min_dist >= threshold_x) { + break; + } + + if (y2 < y3) { + min_dist += y3 - y2; + } else if (y4 < y1) { + min_dist += y1 - y4; + } + + if (min_dist >= longest) { + continue; + } + + p1 = (node1 < 0) ? node2 : ((node2 < 0) ? node1 : -1); + p2 = (node3 < 0) ? node4 : ((node4 < 0) ? node3 : -1); + + if (p1 >= 0 && p2 < 0) { + dist[p1][node3] = + ADIFF(t.branch[p1].x, t.branch[node3].x) + + ADIFF(t.branch[p1].y, t.branch[node3].y); + dist[p1][node4] = + ADIFF(t.branch[p1].x, t.branch[node4].x) + + ADIFF(t.branch[p1].y, t.branch[node4].y); + p2 = (dist[p1][node3] <= dist[p1][node4]) ? node3 : node4; + } else if (p1 < 0 && p2 >= 0) { + dist[node1][p2] = + ADIFF(t.branch[node1].x, t.branch[p2].x) + + ADIFF(t.branch[node1].y, t.branch[p2].y); + dist[node2][p2] = + ADIFF(t.branch[node2].x, t.branch[p2].x) + + ADIFF(t.branch[node2].y, t.branch[p2].y); + p1 = (dist[node1][p2] <= dist[node2][p2]) ? node1 : node2; + } else if (p1 < 0 && p2 < 0) { + // all 4 nodes are real, pick the closest pair + + dist[node1][node3] = + ADIFF(t.branch[node1].x, t.branch[node3].x) + + ADIFF(t.branch[node1].y, t.branch[node3].y); + dist[node1][node4] = + ADIFF(t.branch[node1].x, t.branch[node4].x) + + ADIFF(t.branch[node1].y, t.branch[node4].y); + dist[node2][node3] = + ADIFF(t.branch[node2].x, t.branch[node3].x) + + ADIFF(t.branch[node2].y, t.branch[node3].y); + dist[node2][node4] = + ADIFF(t.branch[node2].x, t.branch[node4].x) + + ADIFF(t.branch[node2].y, t.branch[node4].y); + + p1 = node1; + p2 = node3; + smallest = dist[p1][p2]; + + if (dist[node2][node3] < smallest) { + p1 = node2; + smallest = dist[p1][p2]; + } + if (dist[node1][node4] < smallest) { + p1 = node1; + p2 = node4; + smallest = dist[p1][p2]; + } + if (dist[node2][node4] < smallest) { + p1 = node2; + p2 = node4; + smallest = dist[p1][p2]; + } + } else { + dist[p1][p2] = + ADIFF(t.branch[p1].x, t.branch[p2].x) + + ADIFF(t.branch[p1].y, t.branch[p2].y); + } + + if (min_dist < dist[p1][p2]) { + dist[p1][p2] = dist[p2][p1] = min_dist; + enqueue(&nb[p1], p2); + enqueue(&nb[p2], p1); + + if (min_dist < md) { + md = min_dist; + *min_node1 = p1; + *min_node2 = p2; + } + } + } + } + +#if !(MR_FOR_SMALL_CASES_ONLY) + free(isPin_base); free(id); free(b); +#endif +} + +void mst_from_heap(int d, DTYPE **dist, int node1, int node2, int **nb, + int *edges, int *tree_id) +{ + int i, j, itr, idx; + node_pair e; + + hdist = dist; + heap_size=0; + + for (i=0; i1; j--) { + i=nb[node1][j]; + if (tree_id[i]) continue; + { + e.node2 = i; + insert_heap(&e); + } + } + for (itr=d-2; itr>=1; itr--) { + e.node1 = node2; + + for (j=nb[node2][1]-1; j>1; j--) { + i=nb[node2][j]; + if (tree_id[i]) continue; + { + e.node2 = i; + insert_heap(&e); + } + } + + do { + //assert(heap_size>0); + //extract_heap(&e); + e = heap[1]; + while (tree_id[heap[heap_size].node2]) { + heap_size--; + } + heap[1] = heap[heap_size--]; + sift_down(1); + node2 = e.node2; + } while (tree_id[node2]); + node1 = e.node1; + tree_id[node2] = tree_id[node1]; + edges[idx++] = node1; + edges[idx++] = node2; + } + //assert(idx==2*d-2); +} + +void build_rmst(long d, DTYPE *x, DTYPE *y, int *edges, int *inMST) +{ + int i, j, idx, n; + int *map = (int*)calloc(d, sizeof(int)); + Point *pt = (Point*)calloc( 4*d, sizeof(Point) ); + long *parent = (long*)calloc( 4*d, sizeof(long) ); + long *par = (long*)calloc( 4*d, sizeof(long) ); + int *size = (int*)calloc(d, sizeof(int)); + + for (i=0; i=0; j--) { + if (pt[j].x==pt[i].x && pt[j].y==pt[i].y && par[j]>=0) { + par[i]=j; + break; + } + } + } + } + + for (i=0; i=0); + size[parent[i]]++; + } + + idx = 2*d-3; + for (i=0; i0) { + //assert(!inMST[j]); + inMST[j] = 1; + edges[idx--] = j; + edges[idx--] = j = parent[j]; + size[j]--; + } + } + //assert(idx==-1); + + inMST[edges[0]]=1; + + free(pt); + free(map); + free(parent); + free(par); + free(size); +} + + +/* cached version */ +FTree flutes_c(int d, DTYPE *xs, DTYPE *ys, int *s, int acc) +{ + int i; + //int orig_ht_flag; + FTree t, tdup; + +#if USE_HASHING + dl_forall(FTree, ht[d], t) { + for (i=0; i= d) { // found a match + tdup = t; + tdup.branch = (Branch*)malloc(sizeof(Branch)*(2*d-2)); + for (i=2*d-3; i>=0; i--) { + tdup.branch[i] = t.branch[i]; + } + return tdup; + } + } dl_endfor; + + //orig_ht_flag = new_ht; + //new_ht = 0; +#endif + + t = flutes_LMD(d, xs, ys, s, acc); + +#if USE_HASHING + //new_ht = orig_ht_flag; + + tdup = t; + tdup.branch = (Branch*)malloc(sizeof(Branch)*(2*d-2)); + for (i=2*d-3; i>=0; i--) { + tdup.branch[i] = t.branch[i]; + } + dl_prepend(FTree, ht[d], tdup); +#endif + + return t; +} + +FTree flute_mr(int d, DTYPE *xs, DTYPE *ys, int *s, + int acc, int rounds, DTYPE **dist, + DTYPE *threshold_x, DTYPE *threshold_y, + DTYPE *threshold, + int *best_round, + int *min_node1, int *min_node2, + int **nb) +{ + int i, j, k, m, n, itr, node1, node2; + DTYPE min_dist, longest; + DTYPE dist1, dist2; + FTree t, best_t, *subtree, ttmp; + DTYPE min_x, max_x; + +#if MR_FOR_SMALL_CASES_ONLY + int num_subtree, subroot[MAXPART], suproot[MAXPART], isSuperRoot[D2M]; + int tree_id[D2M], tid, tree_size[D2M], edges[2*D2M]; + int idx[MAXPART], offset[MAXPART], *order[MAXT], order_base[MAXT*D2M]; + DTYPE x[D2M+MAXPART], y[D2M+MAXPART]; + int new_s[D2M+MAXPART], si[D2M], xmap[D2M+MAXPART]; +#else + int num_subtree, *subroot, *suproot, *isSuperRoot; + int *tree_id, tid, *tree_size, *edges; + int *idx, *offset, **order, *order_base; + DTYPE *x, *y; + int *new_s, *si, *xmap; + + subroot = (int*)malloc(sizeof(int)*MAXPART); + suproot = (int*)malloc(sizeof(int)*MAXPART); + idx = (int*)malloc(sizeof(int)*MAXPART); + offset = (int*)malloc(sizeof(int)*MAXPART); + order = (int**)malloc(sizeof(int*)*MAXT); + + isSuperRoot = (int*)malloc(sizeof(int)*d); + tree_id = (int*)malloc(sizeof(int)*d); + tree_size = (int*)malloc(sizeof(int)*d); + edges = (int*)malloc(sizeof(int)*d*2); + order_base = (int*)malloc(sizeof(int)*d*MAXT); + new_s = (int*)malloc(sizeof(int)*(d+MAXPART)); + si = (int*)malloc(sizeof(int)*d); + xmap = (int*)malloc(sizeof(int)*(d+MAXPART)); + x = (DTYPE*)malloc(sizeof(DTYPE)*(d+MAXPART)); + y = (DTYPE*)malloc(sizeof(DTYPE)*(d+MAXPART)); +#endif + +#if USE_HASHING + if (new_ht) { + for (i=0; i<=d; i++) { + ht[i] = dl_alloc(); + } + } +#endif + + best_t.branch = NULL; + best_t.length = INFNTY; + + for (i=0; i=0) { + for (i=0; i=0; ) { + node2 = edges[i--]; + node1 = edges[i--]; + j = tree_size[node1]+tree_size[node2]; + //Chris + if (j >= TAU(acc)) { + isSuperRoot[node1] = 1; + suproot[num_subtree] = node1; + subroot[num_subtree++] = node2; + } else { + tree_size[node1] = j; + } + } + + //assert(num_subtree<=MAXT); + + for (i=1; i= EARLY_QUIT_CRITERIA) { + if (t.branch) { + free(t.branch); + } + break; + } + } else if (best_t.length == t.length) { + *best_round = rounds; + } else if (best_t.length > t.length) { + if (best_t.branch) { + free(best_t.branch); + } + best_t = t; + *best_round = rounds; + } + + if (rounds > 0) { + for (i=0; i=0; ) { + node2 = edges[i--]; + node1 = edges[i--]; + j = tree_size[node1]+tree_size[node2]; + //if (j >= d/2) { + if (j >= d/2 && num_subtree<2) { + isSuperRoot[node1] = 1; + suproot[num_subtree] = node1; + subroot[num_subtree++] = node2; + } else { + tree_size[node1] = j; + } + } + */ + + for (i=2*d-3; i>=0; ) { + node2 = edges[i--]; + node1 = edges[i--]; + tree_size[node1] += tree_size[node2]; + } + + j = 0; + smallest_gap = ADIFF(tree_size[j], d/2); + for (i=1; i=0; ) { + node2 = edges[i--]; + node1 = edges[i--]; + if (node2==j) { + isSuperRoot[node1] = 1; + suproot[num_subtree] = node1; + subroot[num_subtree++] = node2; + tree_size[subroot[0]] -= tree_size[j]; + break; + } + } + + //assert(num_subtree<=MAXT); + + for (i=1; i= INFNTY && d > 1000) { + D3 = (d <= 10000) ? 1000 : 10000; + } + t = flute_am(d, xs, ys, s, A, + &threshold_x, &threshold_y, &threshold); + if (orig_D3 >= INFNTY) { + D3 = orig_D3; + } + } + + return t; +} + +int pickWin(FTree t, DTYPE cx, DTYPE cy, int inWin[]) +{ +#if MR_FOR_SMALL_CASES_ONLY + int i, j, n, dd, cnt, stack[D2M*2], top=0, prev, curr, next; + int isPin_base[D2M], *isPin, nghbr_base[D2M], *nghbr, q[2*D2M]; +#else + int i, j, n, dd, cnt, *stack, top=0, prev, curr, next; + int *isPin_base, *isPin, *nghbr_base, *nghbr, *q; + + stack = (int*)malloc(sizeof(int)*t.deg*2); + isPin_base = (int*)malloc(sizeof(int)*t.deg); + nghbr_base = (int*)malloc(sizeof(int)*t.deg); + q = (int*)malloc(sizeof(int)*t.deg*2); +#endif + + if (t.deg <= 3) { + for (i=0; i 0); + + while (top > 0) { + i = stack[--top]; + if (inWin[i]) { + continue; + } + inWin[i] = 1; + if (i < t.deg) { + cnt++; + continue; + } + n = isPin[i]; + if (n >= 0) { + if (!inWin[n]) { + inWin[n] = 1; + cnt++; + } + } else { + stack[top++] = t.branch[i].n; + for (j=nghbr[i]; j0); +#if !(MR_FOR_SMALL_CASES_ONLY) + free(stack); free(isPin_base); free(nghbr_base); free(q); +#endif + return cnt; +} + +/* merge tree t2 into tree t1, which have shared common nodes. The intention + is to add the non-common tree nodes of t2 into t1, as well as the + associated steiner nodes */ +FTree merge_into(FTree t1, FTree t2, int common[], int nc, int *o1, int *o2) +{ + FTree t; + DTYPE cx, cy; +#if MR_FOR_SMALL_CASES_ONLY + int i, j, k, d, n, offset, map[2*D2M], reachable[2*D2M]; + int o[D2M+MAXPART]; +#else + int i, j, k, d, n, offset, *map, *reachable; + int *o; + + map = (int*)malloc(sizeof(int)*(t1.deg+t2.deg)*2); + reachable = (int*)malloc(sizeof(int)*(t1.deg+t2.deg)*2); + o = (int*)malloc(sizeof(int)*(t1.deg+t2.deg+MAXPART2)); +#endif + + if (t2.deg <= nc) { + free(t2.branch); +#if !(MR_FOR_SMALL_CASES_ONLY) + free(map); free(reachable); free(o); +#endif + return t1; + } + + t.deg = t1.deg + t2.deg - nc; + t.branch = (Branch *) malloc((2*t.deg-2)*sizeof(Branch)); + + offset = t2.deg - nc; + + for (i=t1.deg; i=0; i--) { + if (!common[i] && t2.branch[i].n!=i) { + reachable[t2.branch[i].n] = 1; + } + } + + for (i=2*t2.deg-3; i>=0; i--) { + map[i] = -1; + } + + d = t1.deg*2 - 2; + + /* a slow method; could be sped up here */ + for (i=0; i=t2.deg) { + for (; i=t1.deg) { + for (; k=t1.deg); + t.branch[j].n = map[n] + offset; + o[j] = o2[k]; + j++; + } + } else if (o1[i] < o2[k]) { + t.branch[j] = t1.branch[i]; + t.branch[j].n = t1.branch[i].n + offset; + o[j] = o1[i]; + j++; + i++; + } else { + t.branch[j] = t2.branch[k]; + n = t2.branch[k].n; + //assert(map[n]>=t1.deg); + t.branch[j].n = map[n] + offset; + o[j] = o2[k]; + j++; + do { + k++; + } while (k=t1.deg); + t.branch[map[i]+offset].n = map[n] + offset; + j++; + } + } + + j = d+offset; + //assert(j <= t.deg*2-2); + while (j < t.deg*2-2) { + /* redundant steiner nodes */ + t.branch[j] = t2.branch[0]; + t.branch[j].n = j; + j++; + } + + /* + for (i=0; i=t.deg); + } + */ + + t.length = wirelength(t); + + free(t1.branch); + free(t2.branch); + +#if !(MR_FOR_SMALL_CASES_ONLY) + free(map); free(reachable); free(o); +#endif + return t; +} + +/* simply merge two trees at their common node */ +FTree smergetree(FTree t1, FTree t2, int *o1, int *o2, + DTYPE cx, DTYPE cy) +{ + FTree t; + int d, i, j, k, n, ci, cn, mapped_cn, prev, curr, next, offset; +#if MR_FOR_SMALL_CASES_ONLY + int o[D2M+MAXPART], map[2*D2M]; +#else + int *o, *map; + + map = (int*)malloc(sizeof(int)*(t1.deg+t2.deg)*2); + o = (int*)malloc(sizeof(int)*(t1.deg+t2.deg+MAXPART2)); +#endif + + t.deg = t1.deg + t2.deg - 1; + t.branch = (Branch *) malloc((2*t.deg-2)*sizeof(Branch)); + + offset = t2.deg - 1; + d = t1.deg*2-2; + + for (i=0; i=t2.deg) { + for (; i=t1.deg) { + for (; k2 && t2.deg>2) { + if (d= t2.deg) { + for (; i1= t1.deg) { + for (; i2 2) { + for (i=0; i0) { + dl_pop_first(FTreeNode*, subtree, p); + p->e = p; + grandp = p->parent; + if (grandp) { + p->len = p->blen = ADIFF(p->x, grandp->x) + ADIFF(p->y, grandp->y); + if (p->len < grandp->len) { + p->len = grandp->len; + p->e = grandp->e; + } + } else { + p->len = 0; + } + + if (id) { + p->id = id; + } + + dl_forall(FTreeNode*, p->children, child) { + dl_prepend(FTreeNode*, subtree, child); + } dl_endfor; + } + + dl_free(subtree); +} + +FTreeNode *createRootedFTree(FTree t, int *order, int id, dl_t list_of_nodes) +{ + int i, dd, n; + FTreeNode *root=0, **nodes, *p; + + dd = t.deg*2-2; + nodes = (FTreeNode**)malloc(sizeof(FTreeNode*)*dd); + for (i=0; imark = curr_mark; + nodes[i]->children = dl_alloc(); + } + + curr_mark++; + for (i=0; imark = curr_mark; + n = t.branch[i].n; + if (i==n) { + if (i < t.deg) { + //assert(root==0); + nodes[i]->parent = 0; + root = nodes[i]; + } else { /* must be redundant */ + dl_free(nodes[i]->children); + free(nodes[i]); + nodes[i] = 0; + continue; + } + } else { + p = nodes[n]; + nodes[i]->parent = p; + dl_append(FTreeNode*, p->children, nodes[i]); + } + nodes[i]->order = (i < t.deg) ? order[i] : -1; + nodes[i]->id = id; + nodes[i]->x = t.branch[i].x; + nodes[i]->y = t.branch[i].y; + + /* len will be computed in update_subtree + nodes[i]->blen = + ADIFF(t.branch[i].x, t.branch[n].x)+ADIFF(t.branch[i].y, t.branch[n].y); + + nodes[i]->e = nodes[i]; + nodes[i]->len = + ADIFF(t.branch[i].x, t.branch[n].x)+ADIFF(t.branch[i].y, t.branch[n].y); + */ + + dl_append(FTreeNode*, list_of_nodes, nodes[i]); + } + + //assert(root); + + update_subtree(root, 0); + + for (i=0; imark!=curr_mark) { + dl_free(nodes[i]->children); + free(nodes[i]); + } + } + + free(nodes); + return root; +} + +void freeFTree(FTreeNode *t) +{ + FTreeNode *child; + dl_forall(FTreeNode *, t->children, child) { + freeFTree(child); + } dl_endfor; + dl_free(t->children); + free(t); +} + +int cmpNodeByYX(const void* a, const void* b) +{ + DTYPE ay = (*(FTreeNode**)a)->y; + DTYPE by = (*(FTreeNode**)b)->y; + DTYPE ax, bx; + + if (ay < by) return -1; + if (ay > by) return 1; + + ax = (*(FTreeNode**)a)->x; + bx = (*(FTreeNode**)b)->x; + + if (ax < bx) return -1; + if (ax > bx) return 1; + return 0; +} + +int cmpNodeByXY(const void* a, const void* b) +{ + DTYPE ax = (*(FTreeNode**)a)->x; + DTYPE bx = (*(FTreeNode**)b)->x; + DTYPE ay, by; + + if (ax < bx) return -1; + if (ax > bx) return 1; + + ay = (*(FTreeNode**)a)->y; + by = (*(FTreeNode**)b)->y; + + if (ay < by) return -1; + if (ay > by) return 1; + return 0; +} + +void remove_child(dl_t children_list, FTreeNode* c) +{ + FTreeNode *child; + dl_forall(FTreeNode*, children_list, child) { + if (child==c) { + dl_delete_current(); + break; + } + } dl_endfor; +} + +void cleanFTree(FTreeNode* tn) +{ + /* + FTreeNode *c, *p; + + dl_forall(FTreeNode*, tn->children, c) { + cleanFTree(c); + } dl_endfor; + + p = tn->parent; + if (!p) return; + + if (tn->order >= 0) return; // don't clean pin nodes + + if (dl_length(tn->children)<=0) { + remove_child(p->children, tn); + dl_free(tn->children); + free(tn); + } else if (dl_length(tn->children)<=1) { + c = dl_first(FTreeNode*, tn->children); + c->parent = p; + dl_append(FTreeNode*, p->children, c); + remove_child(p->children, tn); + dl_free(tn->children); + free(tn); + } + */ + + // non-recursive version + FTreeNode *c, *p; + dl_t nlist=dl_alloc(); + + dl_append(FTreeNode*, nlist, tn); + + while (dl_length(nlist)>0) { + dl_pop_first(FTreeNode*, nlist, tn); + dl_forall(FTreeNode*, tn->children, c) { + dl_append(FTreeNode*, nlist, c); + } dl_endfor; + + p = tn->parent; + if (p && tn->order < 0) { + if (dl_length(tn->children)<=0) { + remove_child(p->children, tn); + dl_free(tn->children); + free(tn); + } else if (dl_length(tn->children)<=1) { + c = dl_first(FTreeNode*, tn->children); + c->parent = p; + dl_append(FTreeNode*, p->children, c); + remove_child(p->children, tn); + dl_free(tn->children); + free(tn); + } + } + } + + dl_free(nlist); +} + +int cmpNodeByOrder(void* a, void* b) +{ + int ax = (*(FTreeNode**)a)->order; + int bx = (*(FTreeNode**)b)->order; + + if (ax < bx) return -1; + if (ax > bx) return 1; + return 0; +} + +FTree mergeRootedFTrees(FTreeNode *tn1, FTreeNode *tn2, int *order1) +{ + int i, n, redundant; + FTree t; + FTreeNode *child, *p; + dl_t list_of_nodes=dl_alloc(); + dl_t pin_nodes=dl_alloc(), steiner_nodes=dl_alloc(); + + //assert(tn1->x==tn2->x && tn1->y==tn2->y); + + /* merge tn2 to tn1 */ + while (dl_length(tn2->children)>0) { + dl_pop_first(FTreeNode*, tn2->children, child); + child->parent = tn1; + dl_append(FTreeNode*, tn1->children, child); + } + dl_free(tn2->children); + free(tn2); + + cleanFTree(tn1); + + /* convert tn1 back to a FTree */ + + dl_append(FTreeNode*, list_of_nodes, tn1); + do { + dl_pop_first(FTreeNode*, list_of_nodes, child); + if (child->order < 0) { + if (dl_length(child->children)==1) { /* redundant steiner node */ + p = dl_first(FTreeNode*, child->children); + p->parent = child->parent; + /* note that p->parent's children list is already gone */ + dl_append(FTreeNode*, list_of_nodes, p); + dl_free(child->children); + free(child); + continue; + } else if (dl_length(child->children)==0) { + dl_free(child->children); + free(child); + continue; + } + dl_append(FTreeNode*, steiner_nodes, child); + } else { + dl_append(FTreeNode*, pin_nodes, child); + } + dl_concat(list_of_nodes, child->children); + } while (dl_length(list_of_nodes)>0); + dl_free(list_of_nodes); + + dl_sort(pin_nodes, sizeof(FTreeNode*), cmpNodeByOrder); + + i=0; + dl_forall(FTreeNode*, pin_nodes, child) { + child->id = i++; + } dl_endfor; + + t.deg = i; + + dl_forall(FTreeNode*, steiner_nodes, child) { + child->id = i++; + } dl_endfor; + + //assert(i<=2*t.deg-2); + + t.branch = (Branch*)malloc(sizeof(Branch)*(t.deg*2-2)); + + redundant = i; + for (; i<2*t.deg-2; i++) { + t.branch[i].n = i; + t.branch[i].x = tn1->x; + t.branch[i].y = tn1->y; + } + + t.branch[tn1->id].n = -1; + + dl_forall(FTreeNode*, pin_nodes, child) { + i = child->id; + if (child->order >= 0) { + order1[i] = child->order; + } + t.branch[i].x = child->x; + t.branch[i].y = child->y; + p = child->parent; + if (p) { + if (p->id >= t.deg) { + t.branch[i].n = p->id; + } else { + //assert(p==tn1); + //assert(redundantid].n = redundant; + t.branch[redundant].x = p->x; + t.branch[redundant].y = p->y; + redundant++; + } + } + } dl_endfor; + dl_forall(FTreeNode*, steiner_nodes, child) { + i = child->id; + if (child->order >= 0) { + order1[i] = child->order; + } + t.branch[i].x = child->x; + t.branch[i].y = child->y; + p = child->parent; + if (p->id < t.deg) { // must be the root + if (t.branch[p->id].n < 0) { + t.branch[p->id].n = i; + t.branch[i].n = i; + } else { + n = t.branch[p->id].n; + if (t.branch[p->id].x==t.branch[n].x && + t.branch[p->id].y==t.branch[n].y) { + t.branch[i].n = n; + } else { + //assert(redundantid].x; + t.branch[redundant].y = t.branch[p->id].y; + t.branch[redundant].n = t.branch[p->id].n; + t.branch[p->id].n = redundant; + t.branch[i].n = redundant; + redundant++; + } + } + } else { + t.branch[i].n = p->id; + } + } dl_endfor; + + dl_forall(FTreeNode*, pin_nodes, child) { + free(child); + } dl_endfor; + dl_free(pin_nodes); + + dl_forall(FTreeNode*, steiner_nodes, child) { + free(child); + } dl_endfor; + dl_free(steiner_nodes); + + t.length = wirelength(t); + return t; +} + +void collect_nodes(FTreeNode* tn, dl_t nlist) +{ + /* + FTreeNode* c; + + dl_append(FTreeNode*, nlist, tn); + dl_forall(FTreeNode*, tn->children, c) { + collect_nodes(c, nlist); + }dl_endfor; + */ + // non-recursive version + FTreeNode* c; + dl_el *curr; + + dl_append(FTreeNode*, nlist, tn); + + for (curr=nlist->last; curr; curr=curr->next) { + tn = dl_data(FTreeNode*, curr); + dl_forall(FTreeNode*, tn->children, c) { + dl_append(FTreeNode*, nlist, c); + }dl_endfor; + } +} + +typedef struct { + FTreeNode *n1, *n2; + DTYPE new_x, new_y, gain; +} xdata; + +int cmpXdata(void *a, void *b) +{ + DTYPE ga = (*(xdata*)a).gain; + DTYPE gb = (*(xdata*)b).gain; + if (ga > gb) return -1; + if (ga < gb) return 1; + return 0; +} + +inline FTreeNode *cedge_lca(FTreeNode* n1, FTreeNode* n2, DTYPE *len, int *n2ton1) +{ + int i; + FTreeNode *curr, *lca, *e; + + curr_mark++; + + curr = n1; + while (curr) { + curr->mark = curr_mark; + curr = curr->parent; + } + + lca = n2; + while (lca && lca->mark!=curr_mark) { + lca->mark = curr_mark; + lca = lca->parent; + } + + if (!lca) { + n1 = n1->parent; + if (n1 && n1!=lca && (n1->len > n2->len)) { + *n2ton1 = 0; + *len = n1->len; + return n1->e; + } else { + *n2ton1 = 1; + *len = n2->len; + return n2->e; + } + } + + if (lca==n1 || lca==n1->parent || lca==n2) { + if (lca!=n2) { + *n2ton1 = 1; + *len = n2->blen; + e = n2; + curr = n2->parent; + } else { + *n2ton1 = 0; + *len = n1->blen; + e = n1; + curr = n1->parent; + } + while (curr != lca) { + if (*len < curr->blen) { + *len = curr->blen; + e = curr; + } + curr = curr->parent; + } + return e; + } + + /* lca is above both n1 and n2 */ + *n2ton1 = 0; + n1 = n1->parent; + *len = n1->blen; + e = n1; + curr = n1; + for (i=0; i<2; i++, curr=n2) { + while (curr != lca) { + if (*len < curr->blen) { + if (i>0) { + *n2ton1 = 1; + } + *len = curr->blen; + e = curr; + } + curr = curr->parent; + } + } + + return e; + +} + +FTreeNode *critical_edge(FTreeNode* n1, FTreeNode* n2, DTYPE *len, int *n2ton1) +{ + if (n1->id != n2->id) { + n1 = n1->parent; + if (n1 && (n1->len > n2->len)) { + *n2ton1 = 0; + *len = n1->len; + return n1->e; + } else { + *n2ton1 = 1; + *len = n2->len; + return n2->e; + } + } + + return cedge_lca(n1, n2, len, n2ton1); +} + +void splice2(FTreeNode *n1, FTreeNode *n2, FTreeNode *e) +{ + FTreeNode *curr, *prev, *next, *s; + + //assert(n2->parent); + //assert(e->id==n2->id); + + prev = n2; + curr = n2->parent; + next = curr->parent; + while (prev != e) { + remove_child(curr->children, prev); + curr->parent = prev; + dl_append(FTreeNode*, prev->children, curr); + prev = curr; + curr = next; + next = curr->parent; + } + remove_child(curr->children, prev); + + n2->parent = n1; + dl_append(FTreeNode*, n1->children, n2); + + update_subtree(n1, n1->parent->id); +} + +void cut_and_splice(FTreeNode *n1, FTreeNode *n2, + DTYPE new_x, DTYPE new_y, + DTYPE *x1, DTYPE *y1, DTYPE *x2, DTYPE *y2, + FTreeNode *e, int n2ton1) +{ + FTreeNode *p1, *node, *s; + + /* new steiner node */ + p1 = n1->parent; + remove_child(p1->children, n1); + node = (FTreeNode*)malloc(sizeof(FTreeNode)); + node->x = new_x; + node->y = new_y; + node->mark = curr_mark; + + node->parent = p1; + dl_append(FTreeNode*, p1->children, node); + n1->parent = node; + node->children = dl_alloc(); + dl_append(FTreeNode*, node->children, n1); + node->order = -1; + + node->e = n1->e; + node->id = n1->id; + + if (*x1==n1->x) { + *x2 = new_x; + } else { + *x1 = new_x; + } + if (*y1==n1->y) { + *y2 = new_y; + } else { + *y1 = new_y; + } + + if (n2->order >= 0) { + /* n2 is a pin, need to replicate a steiner node */ + s = n2->parent; + if (s->x!=n2->x || s->y!=n2->y) { + s = (FTreeNode*)malloc(sizeof(FTreeNode)); + s->mark = curr_mark; + s->order = -1; + s->id = n2->id; + s->x = n2->x; + s->y = n2->y; + s->e = n2->e; + if (s->e == n2) { + s->e = s; + } + if (e == n2) { + e = s; + } + s->len = n2->len; + s->blen = n2->blen; + n2->blen = 0; + + remove_child(n2->parent->children, n2); + dl_append(FTreeNode*, n2->parent->children, s); + s->parent = n2->parent; + n2->parent = s; + s->children = dl_alloc(); + dl_append(FTreeNode*, s->children, n2); + } + n2 = s; + } + + if (n2ton1) { + splice2(node, n2, e); + } else { + splice2(n2, node, e); + } +} + +typedef struct { + FTreeNode *n1, *n2; + DTYPE min_dist, new_x, new_y; + int n2ton1; +} splice_info; + +DTYPE exchange_branches_order_x(int num_nodes, FTreeNode **nodes, + DTYPE threshold_x, DTYPE threshold_y, + DTYPE max_len) +{ + int n2ton1; + FTreeNode *n1, *p1, *n2, *p2, *node, *e, *s; + DTYPE x1, x2, y1, y2, min_dist, new_x, new_y, len; + DTYPE gain=0; + int i, j, curr_row, next_header, num_rows, start, end, mid; + int *header=(int*)malloc(sizeof(int)*(num_nodes+1)); + dl_t batch_list=dl_alloc(); + splice_info sinfo; + + int batch_mode = (num_nodes >= D3); + + header[0]=0; + + y1 = nodes[0]->y; + for (i=num_rows=1; iy == y1) { + continue; + } + header[num_rows++] = i; + y1 = nodes[i]->y; + } + header[num_rows] = i; + + curr_row = 0; + next_header = header[1]; + for (i=0; i= next_header) { + curr_row++; + next_header = header[curr_row+1]; + } + n1 = nodes[i]; + p1 = n1->parent; + if (!p1) { + continue; + } + if (p1->x == n1->x && p1->y == n1->y) { + continue; + } + if (n1->x <= p1->x) { + x1 = n1->x; x2 = p1->x; + } else { + x1 = p1->x; x2 = n1->x; + } + if (n1->y <= p1->y) { + y1 = n1->y; y2 = p1->y; + } else { + y1 = p1->y; y2 = n1->y; + } + + if (curr_row > 0) { + for (j=curr_row-1; j>0; j--) { + if (y1 - threshold_y > nodes[header[j]]->y) { + j++; + break; + } + } + } else { + j = 0; + } + for (;j < num_rows && nodes[header[j]]->y <= y2 + threshold_y; j++) { + /* find the closest node on row j */ + start = header[j]; + end = header[j+1]; + while (start < end) { + mid = (start+end)/2; + if (nodes[mid]->x <= x1) { + start = mid + 1; + } else { + end = mid; + } + } + //assert(start==end); + + if (start >= header[j+1]) { + continue; + } + n2 = nodes[start]; + + if (batch_mode && n1->id==n2->id) continue; + + if (!n2->parent) { + continue; + } + + min_dist = n2->x - x2; + + if (abs(min_dist) > threshold_x) { + continue; + } else if (min_dist < 0) { + min_dist = 0; + new_x = n2->x; + } else { + new_x = x2; + } + + if (n2->y < y1) { + min_dist += y1 - n2->y; + new_y = y1; + } else if (n2->y > y2) { + min_dist += n2->y - y2; + new_y = y2; + } else { + new_y = n2->y; + } + + if (min_dist ==0 || min_dist > max_len) { + continue; + } + + e = critical_edge(n1, n2, &len, &n2ton1); + if (min_dist < len && e!=n1) { + if (batch_mode) { + sinfo.n1 = n1; + sinfo.n2 = n2; + sinfo.min_dist = min_dist; + sinfo.new_x = new_x; + sinfo.new_y = new_y; + sinfo.n2ton1 = n2ton1; + dl_append(splice_info, batch_list, sinfo); + } else { + gain += len - min_dist; + cut_and_splice(n1, n2, new_x, new_y, &x1, &y1, &x2, &y2, e, n2ton1); + } + } + } + } + + dl_forall(splice_info, batch_list, sinfo) { + n1 = sinfo.n1; + n2 = sinfo.n2; + n2ton1 = sinfo.n2ton1; + min_dist = sinfo.min_dist; + + e = critical_edge(n1, n2, &len, &n2ton1); + if (min_dist < len && e!=n1) { + gain += len - min_dist; + cut_and_splice(n1, n2, sinfo.new_x, sinfo.new_y, + &x1, &y1, &x2, &y2, e, n2ton1); + } + } dl_endfor; + + dl_free(batch_list); + + free(header); + + return gain; +} + +DTYPE exchange_branches_order_y(int num_nodes, FTreeNode **nodes, + DTYPE threshold_x, DTYPE threshold_y, + DTYPE max_len) +{ + int n2ton1; + FTreeNode *n1, *p1, *n2, *p2, *node, *e, *s; + DTYPE x1, x2, y1, y2, min_dist, new_x, new_y, len; + DTYPE gain=0; + int i, j, curr_row, next_header, num_rows, start, end, mid; + int *header=(int*)malloc(sizeof(int)*(num_nodes+1)); + dl_t batch_list=dl_alloc(); + splice_info sinfo; + + int batch_mode = (num_nodes >= D3); + + header[0]=0; + + x1 = nodes[0]->x; + for (i=num_rows=1; ix == x1) { + continue; + } + header[num_rows++] = i; + x1 = nodes[i]->x; + } + header[num_rows] = i; + + curr_row = 0; + next_header = header[1]; + for (i=0; i= next_header) { + curr_row++; + next_header = header[curr_row+1]; + } + n1 = nodes[i]; + p1 = n1->parent; + if (!p1) { + continue; + } + if (p1->x == n1->x && p1->y == n1->y) { + continue; + } + if (n1->x <= p1->x) { + x1 = n1->x; x2 = p1->x; + } else { + x1 = p1->x; x2 = n1->x; + } + if (n1->y <= p1->y) { + y1 = n1->y; y2 = p1->y; + } else { + y1 = p1->y; y2 = n1->y; + } + + if (curr_row > 0) { + for (j=curr_row-1; j>0; j--) { + if (x1 - threshold_x > nodes[header[j]]->x) { + j++; + break; + } + } + } else { + j = 0; + } + for (;j < num_rows && nodes[header[j]]->x <= x2 + threshold_x; j++) { + /* find the closest node on row j */ + start = header[j]; + end = header[j+1]; + while (start < end) { + mid = (start+end)/2; + if (nodes[mid]->y <= y1) { + start = mid + 1; + } else { + end = mid; + } + } + //assert(start==end); + if (start >= header[j+1]) { + continue; + } + n2 = nodes[start]; + + if (batch_mode && n1->id==n2->id) continue; + + if (!n2->parent) { + continue; + } + + min_dist = n2->y - y2; + + if (abs(min_dist) > threshold_y) { + continue; + } else if (min_dist < 0) { + min_dist = 0; + new_y = n2->y; + } else { + new_y = y2; + } + + if (n2->x < x1) { + min_dist += x1 - n2->x; + new_x = x1; + } else if (n2->x > x2) { + min_dist += n2->x - x2; + new_x = x2; + } else { + new_x = n2->x; + } + + if (min_dist ==0 || min_dist > max_len) { + continue; + } + + e = critical_edge(n1, n2, &len, &n2ton1); + if (min_dist < len && e!=n1) { + if (batch_mode) { + sinfo.n1 = n1; + sinfo.n2 = n2; + sinfo.min_dist = min_dist; + sinfo.new_x = new_x; + sinfo.new_y = new_y; + sinfo.n2ton1 = n2ton1; + dl_append(splice_info, batch_list, sinfo); + } else { + gain += len - min_dist; + cut_and_splice(n1, n2, new_x, new_y, &x1, &y1, &x2, &y2, e, n2ton1); + } + } + } + } + + dl_forall(splice_info, batch_list, sinfo) { + n1 = sinfo.n1; + n2 = sinfo.n2; + n2ton1 = sinfo.n2ton1; + min_dist = sinfo.min_dist; + + e = critical_edge(n1, n2, &len, &n2ton1); + if (min_dist < len && e!=n1) { + gain += len - min_dist; + cut_and_splice(n1, n2, sinfo.new_x, sinfo.new_y, + &x1, &y1, &x2, &y2, e, n2ton1); + } + } dl_endfor; + + dl_free(batch_list); + + free(header); + + return gain; +} + +/* cross exchange branches after merging */ +FTree xmergetree(FTree t1, FTree t2, int *order1, int *order2, + DTYPE cx, DTYPE cy) +{ + int i, num, cnt, order_by_x=1; + FTree t; + FTreeNode *tn1, *tn2, *n1, *p1, **nodes; + dl_t list_of_nodes=dl_alloc(); + DTYPE threshold_x, threshold_y; + DTYPE min_x, max_x, max_len, len, gain; + + if (t1.deg <= 0) { + for (i=0; ix; + for (i=0; iparent; + if (p1) { + len = ADIFF(n1->x, p1->x) + ADIFF(n1->y, p1->y); + if (len > max_len) { + max_len = len; + } + } + if (n1->x < min_x) { + min_x = n1->x; + } else if (n1->x > max_x) { + max_x = n1->x; + } + } + + threshold_x = (max_x - min_x)/4; + threshold_y = (nodes[num-1]->y - nodes[0]->y)/4; + + threshold_x = min(threshold_x, max_len); + threshold_y = min(threshold_y, max_len); + + for (cnt=(t1.deg+t2.deg)/2; cnt>0; cnt--) { + gain = (order_by_x) ? + exchange_branches_order_x(num, nodes, threshold_x, threshold_y, max_len): + exchange_branches_order_y(num, nodes, threshold_x, threshold_y, max_len); + + //assert(gain>=0); + + if (gain <= 0 && !order_by_x) { + break; + } + if (cnt>1) { + collect_nodes(tn1, list_of_nodes); + num = dl_length(list_of_nodes); + if (num <= 1) { + break; + } + + collect_nodes(tn2, list_of_nodes); + if (dl_length(list_of_nodes)-num <= 1) { + break; + } + + free(nodes); + num = dl_length(list_of_nodes); + nodes = (FTreeNode**)malloc(sizeof(FTreeNode*)*num); + i = 0; + dl_forall(FTreeNode*, list_of_nodes, n1) { + nodes[i++] = n1; + } dl_endfor; + dl_clear(list_of_nodes); + + if (order_by_x) { + order_by_x = 0; + qsort(nodes, num, sizeof(FTreeNode*), cmpNodeByXY); + } else { + order_by_x = 1; + qsort(nodes, num, sizeof(FTreeNode*), cmpNodeByYX); + } + } + } + + dl_free(list_of_nodes); + free(nodes); + + t = mergeRootedFTrees(tn1, tn2, order1); + + free(t1.branch); + free(t2.branch); + + return t; +} + diff --git a/knik/src/flute-3.1/src/heap.cpp b/knik/src/flute-3.1/src/heap.cpp new file mode 100644 index 00000000..f075490e --- /dev/null +++ b/knik/src/flute-3.1/src/heap.cpp @@ -0,0 +1,177 @@ +/****************************************************************************/ +/* + Binary heap routines for use in Prim's algorithm, + with points are numbered from 0 to n-1 +*/ + +#include +#include "knik/heap.h" +#include "knik/err.h" + + +Heap* _heap = (Heap*)NULL; +long _max_heap_size = 0; +long _heap_size = 0; + +/****************************************************************************/ +/* +*/ + +void allocate_heap( long n ) +{ + if( _max_heap_size < n ) + { + _heap = (Heap*)realloc( (void*)_heap, (size_t)(n+1)*sizeof(Heap) ); + if( ! _heap ) + { + err_exit( "Cannot reallocate memory in allocate_heap!" ); + } + _max_heap_size = n; + } +} +/****************************************************************************/ +/* +*/ + +void deallocate_heap() +{ + _max_heap_size = 0; + if( _heap ) + { + free( (void*)_heap ); + _heap = (Heap*)NULL; + } +} + +/****************************************************************************/ + +void heap_init( long n ) +{ + register long p; + + allocate_heap( n ); + _heap_size = 0; + for( p = 0; p < n; p++ ) + { + heap_idx( p ) = 0; + } + +} /* END heap_init() */ + +/****************************************************************************/ + +void heap_insert( + long p, + long key +) +{ + register long k; /* hole in the heap */ + register long j; /* parent of the hole */ + register long q; /* heap_elt(j) */ + + heap_key( p ) = key; + + if( _heap_size == 0 ) + { + _heap_size = 1; + heap_elt( 1 ) = p; + heap_idx( p ) = 1; + return; + } + + k = ++ _heap_size; + j = k >> 1; /* k/2 */ + + while( (j > 0) && (heap_key(q=heap_elt(j)) > key) ) { + + heap_elt( k ) = q; + heap_idx( q ) = k; + k = j; + j = k>>1; /* k/2 */ + + } + + /* store p in the position of the hole */ + heap_elt( k ) = p; + heap_idx( p ) = k; + +} /* END heap_insert() */ + + +/****************************************************************************/ + +void heap_decrease_key +( + long p, + long new_key +) +{ + register long k; /* hole in the heap */ + register long j; /* parent of the hole */ + register long q; /* heap_elt(j) */ + + heap_key( p ) = new_key; + k = heap_idx( p ); + j = k >> 1; /* k/2 */ + + if( (j > 0) && (heap_key(q=heap_elt(j)) > new_key) ) { /* change is needed */ + do { + + heap_elt( k ) = q; + heap_idx( q ) = k; + k = j; + j = k>>1; /* k/2 */ + + } while( (j > 0) && (heap_key(q=heap_elt(j)) > new_key) ); + + /* store p in the position of the hole */ + heap_elt( k ) = p; + heap_idx( p ) = k; + } + +} /* END heap_decrease_key() */ + + +/****************************************************************************/ + +long heap_delete_min() +{ + long min, last; + register long k; /* hole in the heap */ + register long j; /* child of the hole */ + register long l_key; /* key of last point */ + + if( _heap_size == 0 ) /* heap is empty */ + return( -1 ); + + min = heap_elt( 1 ); + last = heap_elt( _heap_size -- ); + l_key = heap_key( last ); + + k = 1; j = 2; + while( j <= _heap_size ) { + + if( heap_key(heap_elt(j)) > heap_key(heap_elt(j+1)) ) + j++; + + if( heap_key(heap_elt(j)) >= l_key) + break; /* found a position to insert 'last' */ + + /* else, sift hole down */ + heap_elt(k) = heap_elt(j); /* Note that j <= _heap_size */ + heap_idx( heap_elt(k) ) = k; + k = j; + j = k << 1; + } + + heap_elt( k ) = last; + heap_idx( last ) = k; + + heap_idx( min ) = -1; /* mark the point visited */ + return( min ); + +} /* END heap_delete_min() */ + + +/****************************************************************************/ + diff --git a/knik/src/flute-3.1/src/knik/dist.h b/knik/src/flute-3.1/src/knik/dist.h new file mode 100644 index 00000000..69eb97e7 --- /dev/null +++ b/knik/src/flute-3.1/src/knik/dist.h @@ -0,0 +1,16 @@ +#ifndef _DIST_H_ +#define _DIST_H_ + +#include "global.h" + +long dist( + Point p, + Point q +); + +long dist2( + Point* p, + Point* q +); + +#endif diff --git a/knik/src/flute-3.1/src/knik/dl.h b/knik/src/flute-3.1/src/knik/dl.h new file mode 100644 index 00000000..afa71614 --- /dev/null +++ b/knik/src/flute-3.1/src/knik/dl.h @@ -0,0 +1,180 @@ +#ifndef DL_H +#define DL_H + +#include +#include + +typedef struct dl_el_s { + struct dl_el_s *prev, *next; +} dl_el; + +typedef struct { + dl_el *first, *last; + unsigned int count; +} dl_s, *dl_t; + +dl_t dl_alloc(void); +void dl_delete(dl_t dl, dl_el *el); +void dl_clear(dl_t dl); +void dl_concat(dl_t list1, dl_t list2); +void dl_sort(dl_t dl, size_t el_size, int(*compar)(void *, void *)); + +#define dl_length(dl) (dl)->count + +#define dl_empty(dl) ((dl)->count <= 0) + +#define dl_data(type, el) \ + *(type*)(((dl_el*)(el))+1) + +#define dl_data_p(type, el) \ + ((type*)(((dl_el*)(el))+1)) + +#define dl_forall(type, dl, data) \ +{ \ + dl_el *_el, *_next; \ + dl_t _curr_dl = (dl); \ + for (_el=_curr_dl->first; _el; _el=_next) { \ + _next = _el->next; \ + (data) = dl_data(type, _el); + +#define dl_forall_p(type, dl, data_p) \ +{ \ + dl_el *_el, *_next; \ + dl_t _curr_dl = (dl); \ + for (_el=_curr_dl->first; _el; _el=_next) { \ + _next = _el->next; \ + (data_p) = dl_data_p(type, _el); + +#define dl_current() _el +#define dl_delete_current() dl_delete(_curr_dl, _el) + +#define dl_endfor \ + } \ +} + +#define dl_forall_reverse(type, dl, data) \ +{ \ + dl_el *_el, *_next; \ + dl_t _curr_dl = (dl); \ + for (_el=_curr_dl->last; _el; _el=_next) { \ + _next = _el->prev; \ + (data) = dl_data(type, _el); + +#define dl_forall_reverse_p(type, dl, data_p) \ +{ \ + dl_el *_el, *_next; \ + dl_t _curr_dl = (dl); \ + for (_el=_curr_dl->last; _el; _el=_next) { \ + _next = _el->prev; \ + (data_p) = dl_data_p(type, _el); + +#define dl_first(type, dl) \ + dl_data(type, (dl)->first) + + +#define dl_first_element(dl) (dl)->first + + +#define dl_last(type, dl) \ + dl_data(type, (dl)->last) + +#define dl_pop_first(type, dl, data) \ +{ \ + (data) = dl_first(type, dl); \ + dl_delete((dl), (dl)->first); \ +} + +#define dl_pop_last(type, dl, data) \ +{ (data) = dl_last(type, dl); dl_delete((dl), (dl)->last); } + +#define dl_insert_before(type, dl, element, data) \ +{ \ + if ((element) == (dl)->first) { \ + dl_prepend(type, dl, data); \ + } else { \ + dl_el *_el = (dl_el*) malloc(sizeof(dl_el)+sizeof(type)); \ + if (!_el) { \ + printf("Out of memory!!\n"); \ + } else { \ + memcpy(_el+1, &(data), sizeof(type)); \ + _el->prev = (element)->prev; _el->next = (element); \ + (element)->prev->next = _el; (element)->prev = _el; \ + (dl)->count++; \ + } \ + } \ +} + +#define dl_insert_after(type, dl, element, data) \ +{ \ + if ((element) == (dl)->last) { \ + dl_append(type, dl, data); \ + } else { \ + dl_el *_el = (dl_el*) malloc(sizeof(dl_el)+sizeof(type)); \ + if (!_el) { \ + printf("Out of memory!!\n"); \ + } else { \ + memcpy(_el+1, &(data), sizeof(type)); \ + _el->next = (element)->next; _el->prev = (element); \ + (element)->next->prev = _el; (element)->next = _el; \ + (dl)->count++; \ + } \ + } \ +} + +#define dl_append(type, dl, data) \ +{ \ + dl_el *_el = (dl_el*) malloc(sizeof(dl_el)+sizeof(type)); \ + if (!_el) { \ + printf("Out of memory!!\n"); \ + } else { \ + memcpy(_el+1, &(data), sizeof(type)); \ + _el->next = 0; \ + if ((dl)->count <= 0) { \ + _el->prev = 0; \ + (dl)->first = (dl)->last = _el; \ + (dl)->count = 1; \ + } else { \ + _el->prev = (dl)->last; \ + (dl)->last->next = _el; \ + (dl)->last = _el; \ + (dl)->count++; \ + } \ + } \ +} + +#define dl_prepend(type, dl, data) \ +{ \ + dl_el *_el = (dl_el*) malloc(sizeof(dl_el)+sizeof(type)); \ + if (!_el) { \ + printf("Out of memory!!\n"); \ + } else { \ + memcpy(_el+1, &(data), sizeof(type)); \ + _el->prev = 0; \ + if ((dl)->count <= 0) { \ + _el->next = 0; \ + (dl)->first = (dl)->last = _el; \ + (dl)->count = 1; \ + } else { \ + _el->next = (dl)->first; \ + (dl)->first->prev = _el; \ + (dl)->first = _el; \ + (dl)->count++; \ + } \ + } \ +} + +#define dl_free(dl) \ +{ \ + dl_clear(dl); free(dl); dl = 0; \ +} + +#define dl_duplicate(dest, src, type) \ +{ \ + dest = dl_alloc(); \ + type _data_el; \ + dl_forall(type, src, _data_el) { \ + dl_append(type, dest, _data_el); \ + } dl_endfor; \ +} + +#endif diff --git a/knik/src/flute-3.1/src/knik/err.h b/knik/src/flute-3.1/src/knik/err.h new file mode 100644 index 00000000..a35724e6 --- /dev/null +++ b/knik/src/flute-3.1/src/knik/err.h @@ -0,0 +1,12 @@ +#ifndef _ERR_H_ +#define _ERR_H_ + +void err_msg( +char* msg +); + +void err_exit( +char* msg +); + +#endif diff --git a/knik/src/flute-2.4/src/knik/flute.h b/knik/src/flute-3.1/src/knik/flute.h similarity index 53% rename from knik/src/flute-2.4/src/knik/flute.h rename to knik/src/flute-3.1/src/knik/flute.h index 5de45c5f..95ec6ca7 100644 --- a/knik/src/flute-2.4/src/knik/flute.h +++ b/knik/src/flute-3.1/src/knik/flute.h @@ -1,19 +1,43 @@ #ifndef _KNIK_FLUTE_H #define _KNIK_FLUTE_H -#define POWVFILE "POWV9.dat" // LUT for POWV (Wirelength Vector) -#define POSTFILE "POST9.dat" // LUT for POST (Steiner FTree) -#define MAXD 150 // max. degree of a net that can be handled - // setting MAXD to more than 150 is not recommended -#define D 9 // LUT is used for d <= D, D <= 9 -#define ROUTING 1 // 1 to construct routing, 0 to estimate WL only -#define REMOVE_DUPLICATE_PIN 0 // remove dup. pin for flute_wl() & flute() +/*****************************/ +/* User-Defined Parameters */ +/*****************************/ +#define MAXD 150 // max. degree that can be handled #define ACCURACY 3 // Default accuracy +#define ROUTING 1 // 1 to construct routing, 0 to estimate WL only +#define LOCAL_REFINEMENT 0 // Suggestion: Set to 1 if ACCURACY >= 5 +#define REMOVE_DUPLICATE_PIN 0 // Remove dup. pin for flute_wl() & flute() #ifndef DTYPE // Data type for distance #define DTYPE int #endif + +/*****************************/ +/* User-Callable Functions */ +/*****************************/ +// void readLUT(); +// DTYPE flute_wl(int d, DTYPE x[], DTYPE y[], int acc); +// DTYPE flutes_wl(int d, DTYPE xs[], DTYPE ys[], int s[], int acc); +// FTree flute(int d, DTYPE x[], DTYPE y[], int acc); +// FTree flutes(int d, DTYPE xs[], DTYPE ys[], int s[], int acc); +// DTYPE wirelength(FTree t); +// void printtree(FTree t); +// void plottree(FTree t); + + +/*************************************/ +/* Internal Parameters and Functions */ +/*************************************/ +#define POWVFILE "POWV9.dat" // LUT for POWV (Wirelength Vector) +#define POSTFILE "POST9.dat" // LUT for POST (Steiner FTree) +#define D 9 // LUT is used for d <= D, D <= 9 +#define TAU(A) (8+1.3*(A)) +#define D1(A) (25+120/((A)*(A))) // flute_mr is used for D1 < d <= D2 +#define D2(A) ((A)<=6 ? 500 : 75+5*(A)) + typedef struct { DTYPE x, y; // starting point of the branch @@ -27,7 +51,7 @@ struct FTree Branch *branch; // array of tree branches }; -// Major functions +// User-Callable Functions extern void readLUT(); extern DTYPE flute_wl(int d, DTYPE x[], DTYPE y[], int acc); //Macro: DTYPE flutes_wl(int d, DTYPE xs[], DTYPE ys[], int s[], int acc); @@ -37,13 +61,14 @@ extern DTYPE wirelength(FTree t); extern void printtree(FTree t); extern void plottree(FTree t); - // Other useful functions +extern void init_param(); extern DTYPE flutes_wl_LD(int d, DTYPE xs[], DTYPE ys[], int s[]); extern DTYPE flutes_wl_MD(int d, DTYPE xs[], DTYPE ys[], int s[], int acc); extern DTYPE flutes_wl_RDP(int d, DTYPE xs[], DTYPE ys[], int s[], int acc); extern FTree flutes_LD(int d, DTYPE xs[], DTYPE ys[], int s[]); extern FTree flutes_MD(int d, DTYPE xs[], DTYPE ys[], int s[], int acc); +extern FTree flutes_HD(int d, DTYPE xs[], DTYPE ys[], int s[], int acc); extern FTree flutes_RDP(int d, DTYPE xs[], DTYPE ys[], int s[], int acc); #if REMOVE_DUPLICATE_PIN==1 @@ -55,12 +80,20 @@ extern FTree flutes_RDP(int d, DTYPE xs[], DTYPE ys[], int s[], int acc); #endif #define flutes_wl_ALLD(d, xs, ys, s, acc) flutes_wl_LMD(d, xs, ys, s, acc) -#define flutes_ALLD(d, xs, ys, s, acc) flutes_LMD(d, xs, ys, s, acc) -//#define flutes_ALLD(d, xs, ys, s, acc) (d<=D ? flutes_LD(d, xs, ys, s) : (d<=D2 ? flutes_MD(d, xs, ys, s, acc) : flutes_HD(d, xs, ys, s, acc))) +#define flutes_ALLD(d, xs, ys, s, acc) \ + (d<=D ? flutes_LD(d, xs, ys, s) \ + : (d<=D1(acc) ? flutes_MD(d, xs, ys, s, acc) \ + : flutes_HD(d, xs, ys, s, acc))) #define flutes_wl_LMD(d, xs, ys, s, acc) \ (d<=D ? flutes_wl_LD(d, xs, ys, s) : flutes_wl_MD(d, xs, ys, s, acc)) #define flutes_LMD(d, xs, ys, s, acc) \ (d<=D ? flutes_LD(d, xs, ys, s) : flutes_MD(d, xs, ys, s, acc)) +#define max(x,y) ((x)>(y)?(x):(y)) +#define min(x,y) ((x)<(y)?(x):(y)) +#define abs(x) ((x)<0?(-x):(x)) +#define ADIFF(x,y) ((x)>(y)?(x-y):(y-x)) // Absolute difference + #endif + diff --git a/knik/src/flute-3.1/src/knik/global.h b/knik/src/flute-3.1/src/knik/global.h new file mode 100644 index 00000000..96f8f89b --- /dev/null +++ b/knik/src/flute-3.1/src/knik/global.h @@ -0,0 +1,19 @@ +#ifndef _GLOBAL_H_ +#define _GLOBAL_H_ + +#include + +#define TRUE 1 +#define FALSE 0 +#define MAXLONG 0x7fffffffL + +struct point +{ + long x, y; +}; + +typedef struct point Point; + +typedef long nn_array[8]; + +#endif /* _GLOBAL_H_ */ diff --git a/knik/src/flute-3.1/src/knik/heap.h b/knik/src/flute-3.1/src/knik/heap.h new file mode 100644 index 00000000..eaa51294 --- /dev/null +++ b/knik/src/flute-3.1/src/knik/heap.h @@ -0,0 +1,31 @@ +#ifndef _HEAP_H_ +#define _HEAP_H_ + +#include "global.h" + +struct heap_info +{ + long key; + long idx; + long elt; +}; + +typedef struct heap_info Heap; + +extern Heap* _heap; + +#define heap_key( p ) ( _heap[p].key ) +#define heap_idx( p ) ( _heap[p].idx ) +#define heap_elt( k ) ( _heap[k].elt ) + +#define in_heap( p ) ( heap_idx(p) > 0 ) +#define never_seen( p ) ( heap_idx(p) == 0 ) + +void allocate_heap( long n ); +void deallocate_heap(); +void heap_init( long n ); +void heap_insert( long p, long key ); +void heap_decrease_key( long p, long new_key ); +long heap_delete_min(); + +#endif /* _HEAP_H_ */ diff --git a/knik/src/flute-3.1/src/knik/mst2.h b/knik/src/flute-3.1/src/knik/mst2.h new file mode 100644 index 00000000..e4650515 --- /dev/null +++ b/knik/src/flute-3.1/src/knik/mst2.h @@ -0,0 +1,11 @@ +#ifndef _MST2_H_ +#define _MST2_H_ + +#include "global.h" + +void mst2_package_init( long n ); +void mst2_package_done(); +void mst2( long n, Point* pt, long* parent ); + +#endif + diff --git a/knik/src/flute-3.1/src/knik/neighbors.h b/knik/src/flute-3.1/src/knik/neighbors.h new file mode 100644 index 00000000..491a0a94 --- /dev/null +++ b/knik/src/flute-3.1/src/knik/neighbors.h @@ -0,0 +1,19 @@ +#include "global.h" + +void allocate_nn_arrays( long n ); +void deallocate_nn_arrays(); + +void brute_force_nearest_neighbors +( + long n, + Point* pt, + nn_array* nn +); + +void dq_nearest_neighbors +( + long n, + Point* pt, + nn_array* nn +); + diff --git a/knik/src/flute-3.1/src/mst2.cpp b/knik/src/flute-3.1/src/mst2.cpp new file mode 100644 index 00000000..a3fa0e9b --- /dev/null +++ b/knik/src/flute-3.1/src/mst2.cpp @@ -0,0 +1,92 @@ +#include +#include +#include +#include "knik/global.h" +#include "knik/neighbors.h" +#include "knik/dist.h" +#include "knik/heap.h" +#include "knik/err.h" + + + +void mst2_package_init( long n ) +{ + allocate_heap( n ); + allocate_nn_arrays( n ); +} + +/****************************************************************************/ +/* +*/ + +void mst2_package_done() +{ + deallocate_heap(); + deallocate_nn_arrays(); +} + +/****************************************************************************/ +/* +*/ + +void mst2 +( + long n, + Point* pt, + long* parent +) +{ + long i, k, nn1; + long d; + long oct; + long root = 0; + extern nn_array* nn; + +// brute_force_nearest_neighbors( n, pt, nn ); + dq_nearest_neighbors( n, pt, nn ); + + /* + Binary heap implementation of Prim's algorithm. + Runs in O(n*log(n)) time since at most 8n edges are considered + */ + + heap_init( n ); + heap_insert( root, 0 ); + parent[root] = root; + + for( k = 0; k < n; k++ ) /* n points to be extracted from heap */ + { + i = heap_delete_min(); + + if (i<0) break; +#ifdef DEBUG + assert( i >= 0 ); +#endif + + /* + pt[i] entered the tree, update heap keys for its neighbors + */ + for( oct = 0; oct < 8; oct++ ) + { + nn1 = nn[i][oct]; + if( nn1 >= 0 ) + { + d = dist( pt[i], pt[nn1] ); + if( in_heap(nn1) && (d < heap_key(nn1)) ) + { + heap_decrease_key( nn1, d ); + parent[nn1] = i; + } + else if( never_seen(nn1) ) + { + heap_insert( nn1, d ); + parent[nn1] = i; + } + } + } + } +} + +/****************************************************************************/ +/****************************************************************************/ + diff --git a/knik/src/flute-3.1/src/neighbors.cpp b/knik/src/flute-3.1/src/neighbors.cpp new file mode 100644 index 00000000..d17fec1d --- /dev/null +++ b/knik/src/flute-3.1/src/neighbors.cpp @@ -0,0 +1,527 @@ +#include +#include +#include +#include "knik/global.h" +#include "knik/err.h" +#include "knik/dist.h" + +long octant +( + Point from, + Point to +); + +static Point* _pt; + +/***************************************************************************/ +/* + For efficiency purposes auxiliary arrays are allocated as globals +*/ + +long max_arrays_size = 0; +nn_array* nn = (nn_array*)NULL; +Point* sheared = (Point*)NULL; +long* sorted = (long*)NULL; +long* aux = (long*)NULL; + +/***************************************************************************/ +/* + resize the auxiliary arrays to fit the specified number of points +*/ + +void allocate_nn_arrays( long n ) +{ + if( max_arrays_size < n ) + { + nn = (nn_array*)realloc( (void*)nn, (size_t)n*sizeof(nn_array) ); + sheared = (Point*)realloc( (void*)sheared, (size_t)n*sizeof(Point) ); + sorted = (long*)realloc( (void*)sorted, (size_t)n*sizeof(long) ); + aux = (long*)realloc( (void*)aux, (size_t)n*sizeof(long) ); + if( !nn || !sheared || !sorted || !aux ) + { + err_exit( "Cannot allocate memory in allocate_nn_arrays!" ); + } + max_arrays_size = n; + } +} + +/***************************************************************************/ +/* + free memory used by auxiliary arrays +*/ + +void deallocate_nn_arrays() +{ + max_arrays_size = 0; + if( nn ) + { + free( (void*)nn ); + nn = (nn_array*)NULL; + } + if( sheared ) + { + free( (void*)sheared ); + sheared = (Point*)NULL; + } + if( sorted ) + { + free( (void*)sorted ); + sorted = (long*)NULL; + } + if( aux ) + { + free( (void*)aux ); + aux = (long*)NULL; + } + +} + +/***************************************************************************/ +/* + comparison function for use in quicksort +*/ + +static int compare_x +( + const void* i, + const void* j +) +{ + /* + points with the same x must appear in increasing order of y + */ + if( sheared[*((long*)i)].x == sheared[*((long*)j)].x) + { + return sheared[*((long*)i)].y - sheared[*((long*)j)].y; + } + else + { + return sheared[*((long*)i)].x - sheared[*((long*)j)].x; + } +} + + +/***************************************************************************/ +/* + Combine step of the Guibas-Stolfi divide-and-conquer NE nearest neighbor + algorithm. For efficiency purposes SW nearest neighbors are computed + at the same time. +*/ + +void ne_sw_combine +( + long left, + long mid, + long right, + Point* pt, + long* sorted, + long* aux, + long oct, + nn_array* nn +) +{ + long i, j, k, y2; + long i1; + long i2; + long best_i2; /* index of current best nearest-neighbor */ + long best_dist; /* distance to best nearest-neighbor */ + long d; + +#ifdef DEBUG + assert( right > mid ); + assert( mid > left ); +#endif + + /* + update north-east nearest neighbors accross the mid-line + */ + + i1 = left; + i2 = mid; y2 = pt[ sorted[i2] ].y; + + while( (i1 < mid) && (pt[ sorted[i1] ].y >= y2) ) + { + i1++; + } + + if( i1 < mid ) + { + best_i2 = i2; + best_dist = dist2( pt + sorted[i1], pt + sorted[best_i2] ); + i2++; + + while( (i1 < mid) && (i2 < right) ) + { + if( pt[ sorted[i1] ].y < pt[ sorted[i2] ].y ) + { + d = dist2( pt + sorted[i1], pt + sorted[i2] ); + if( d < best_dist ) + { + best_i2 = i2; + best_dist = d; + } + i2++; + } + else + { + if( (nn[ sorted[i1] ][oct] == -1) || + ( best_dist < dist2( pt + sorted[i1], pt + nn[ sorted[i1] ][oct]) ) + ) + { + nn[ sorted[i1] ][oct] = sorted[best_i2]; + } + i1++; + if( i1 < mid ) + { + best_dist = dist2( pt + sorted[i1], pt + sorted[best_i2] ); + } + } + } + + while( i1 < mid ) + { + if( (nn[ sorted[i1] ][oct] == -1) || + ( dist2( pt + sorted[i1], pt + sorted[best_i2] ) < + dist2( pt + sorted[i1], pt + nn[ sorted[i1] ][oct]) ) + ) + { + nn[ sorted[i1] ][oct] = sorted[best_i2]; + } + i1++; + } + } + /* + repeat for south-west nearest neighbors + */ + + oct = (oct + 4) % 8; + + i1 = right - 1; + i2 = mid - 1; y2 = pt[ sorted[i2] ].y; + + while( (i1 >= mid) && (pt[ sorted[i1] ].y <= y2) ) + { + i1--; + } + + if( i1 >= mid ) + { + best_i2 = i2; + best_dist = dist2( pt + sorted[i1], pt + sorted[best_i2] ); + i2--; + + while( (i1 >= mid) && (i2 >= left) ) + { + if( pt[ sorted[i1] ].y > pt[ sorted[i2] ].y ) + { + d = dist2( pt + sorted[i1], pt + sorted[i2] ); + if( d < best_dist ) + { + best_i2 = i2; + best_dist = d; + } + i2--; + } + else + { + if( (nn[ sorted[i1] ][oct] == -1) || + ( best_dist < dist2( pt + sorted[i1], pt + nn[ sorted[i1] ][oct]) ) + ) + { + nn[ sorted[i1] ][oct] = sorted[best_i2]; + } + i1--; + if( i1 >= mid ) + { + best_dist = dist2( pt + sorted[i1], pt + sorted[best_i2] ); + } + } + } + + while( i1 >= mid ) + { + if( (nn[ sorted[i1] ][oct] == -1) || + ( dist2( pt + sorted[i1], pt + sorted[best_i2] ) < + dist2( pt + sorted[i1], pt + nn[ sorted[i1] ][oct]) ) + ) + { + nn[ sorted[i1] ][oct] = sorted[best_i2]; + } + i1--; + } + } + + /* + merge sorted[left..mid-1] with sorted[mid..right-1] by y-coordinate + */ + + i = left; /* first unprocessed element in left list */ + j = mid; /* first unprocessed element in right list */ + k = left; /* first free available slot in output list */ + + while( (i < mid) && (j < right) ) + { + if( pt[ sorted[i] ].y >= pt[ sorted[j] ].y ) + { + aux[k++] = sorted[i++]; + } + else + { + aux[k++] = sorted[j++]; + } + } + + /* + copy leftovers + */ + while( i < mid ) { aux[k++] = sorted[i++]; } + while( j < right ) { aux[k++] = sorted[j++]; } + + /* + now copy sorted points from 'aux' to 'sorted' + */ + + for( i = left; i < right; i++ ) { sorted[i] = aux[i]; } + +#if 0 + memcpy( (void*)(sorted+left), /* destination */ + (void*)(aux+left), /* source */ + (size_t)(right-left)*sizeof(long) /* number of bytes */ + ); +#endif + +} + +/***************************************************************************/ +/* + compute north-east and south-west nearest neighbors for points indexed + by {sorted[left],...,sorted[right-1]} +*/ + +void ne_sw_nearest_neighbors +( + long left, + long right, + Point* pt, + long* sorted, + long* aux, + long oct, + nn_array* nn +) +{ + long mid; + +#ifdef DEBUG + assert( right > left ); +#endif + + if( right == left + 1 ) + { + nn[ sorted[left] ][oct] = nn[ sorted[left]][(oct+4) % 8] = -1; + } + else + { + mid = (left + right) / 2; + ne_sw_nearest_neighbors( left, mid, pt, sorted, aux, oct, nn ); + ne_sw_nearest_neighbors( mid, right, pt, sorted, aux, oct, nn ); + ne_sw_combine( left, mid, right, pt, sorted, aux, oct, nn ); + } +} + +/***************************************************************************/ +/* + Guibas-Stolfi algorithm for computing nearest NE neighbors +*/ + +void dq_nearest_neighbors +( + long n, + Point* pt, + nn_array* nn +) +{ + long i, oct; + void check_nn( long, Point*, nn_array* ); + + long shear[4][4] = { + {1, -1, 0, 2}, + {2, 0, -1, 1}, + {1, 1, -2, 0}, + {0, 2, -1, -1} + }; + + + +_pt = pt; + + for( oct = 0; oct < 4; oct++ ) + { + for( i = 0; i < n; i++ ) + { + sheared[i].x = shear[oct][0]*pt[i].x + shear[oct][1]*pt[i].y; + sheared[i].y = shear[oct][2]*pt[i].x + shear[oct][3]*pt[i].y; + sorted[i] = i; + } + + qsort( sorted, n, sizeof(long), compare_x ); + ne_sw_nearest_neighbors( 0, n, sheared, sorted, aux, oct, nn ); + } + +#ifdef DEBUG + check_nn( n, pt, nn ); +#endif + +} + +/***************************************************************************/ +/***************************************************************************/ +/* + Brute-force nearest-neighbor computation for debugging purposes +*/ + +/***************************************************************************/ +/* + Half-open octants are numbered from 0 to 7 in anti-clockwise order + starting from ( dx >= dy > 0 ). +*/ + +#define sgn(x) ( x>0 ? 1 : (x < 0 ? -1 : 0) ) + +long octant +( + Point from, + Point to +) +{ + long dx = to.x - from.x; + long dy = to.y - from.y; + long sgn1 = sgn(dx)*sgn(dy); + long sgn2 = sgn(dx+dy)*sgn(dx-dy); + long oct = 0x0; + + + if( (dy < 0) || ((dy==0) && (dx>0)) ) oct += 4; + if( (sgn1 < 0) || (dy==0) ) oct += 2; + if( (sgn1*sgn2 < 0) || (dy==0) || (dx==0) ) oct += 1; + + return oct; +} + +/***************************************************************************/ +/* + O(n^2) algorithm for computing all nearest neighbors +*/ + +void brute_force_nearest_neighbors +( + long n, + Point* pt, + nn_array* nn +) +{ + long i, j, oct; + long d; + + /* + compute nearest neighbors by inspecting all pairs of points + */ + for( i = 0; i < n; i++ ) + { + for( oct = 0; oct < 8; oct++ ) + { + nn[i][oct] = -1; + } + } + + for( i = 0; i < n; i++ ) + { + for( j = i+1; j < n; j++ ) + { + d = dist(pt[i], pt[j]); + + oct = octant( pt[i], pt[j] ); + if( ( nn[i][oct] == -1 ) || + ( d < dist(pt[i], pt[ nn[i][oct] ]) ) + ) + { + nn[i][oct] = j; + } + + oct = (oct + 4) % 8; + if( ( nn[j][oct] == -1 ) || + ( d < dist(pt[j], pt[ nn[j][oct] ]) ) + ) + { + nn[j][oct] = i; + } + } + } +} + + +/***************************************************************************/ +/* + compare nearest neighbors against those computed by brute force +*/ + +void check_nn +( + long n, + Point* pt, + nn_array* nn +) +{ + long i, j, oct; + nn_array* nn1; + + nn1 = (nn_array*)calloc( (size_t)n, (size_t)sizeof(nn_array) ); + brute_force_nearest_neighbors( n, pt, nn1 ); + + for( i = 0; i < n; i++ ) + { + for( oct = 0; oct < 8; oct++ ) + { + if( nn[i][oct] == -1 ) + { + assert( nn1[i][oct] == -1 ); + } + else + { + assert( nn1[i][oct] != -1 ); + + if( octant(pt[i], pt[ nn[i][oct] ]) != oct ) + { + printf( "WRONG OCTANT!\noct=%ld\n", oct ); + printf( "i=%ld, x=%ld, y=%ld\n", i, pt[i].x, pt[i].y ); + j = nn[i][oct]; + printf( "nn=%ld, x=%ld, y=%ld, dist = %ld\n", j, pt[j].x, pt[j].y, + dist(pt[i], pt[j ]) ); + } +// assert( octant(pt[i], pt[ nn[i][oct] ]) == oct ); + + assert( octant(pt[i], pt[ nn1[i][oct] ]) == oct ); + + if( dist(pt[i], pt[ nn[i][oct] ]) != + dist(pt[i], pt[ nn1[i][oct] ]) ) + { + printf( "NNs DON'T MATCH!\noct=%ld\n", oct ); + printf( "i=%ld, x=%ld, y=%ld\n", i, pt[i].x, pt[i].y ); + j = nn[i][oct]; + printf( "nn=%ld, x=%ld, y=%ld, dist = %ld\n", j, pt[j].x, pt[j].y, + dist(pt[i], pt[j ]) ); + j = nn1[i][oct]; + printf( "nn1=%ld, x=%ld, y=%ld, dist = %ld\n", j, pt[j].x, pt[j].y, + dist(pt[i], pt[ j ]) ); + } +// assert( dist(pt[i], pt[ nn[i][oct] ]) == +// dist(pt[i], pt[ nn1[i][oct] ]) ); + } + } + } + + free( nn1 ); +} + +/***************************************************************************/ +/***************************************************************************/ +