Added FLUTE library as a top level standalone tool in Coriolis.

flute/CMakeLists.txt

@@ -0,0 +1,20 @@
+# -*- mode: CMAKE explicit-buffer-name: "CMakeLists.txt<flute>" -*-
+
+ set( CMAKE_LEGACY_CYGWIN_WIN32 0 )
+ project( FLUTE )
+
+ cmake_minimum_required( VERSION 2.8.9 )
+
+ set( ignoreVariables "${BUILD_DOC}" )
+
+ list( INSERT CMAKE_MODULE_PATH 0 "${DESTDIR}$ENV{CORIOLIS_TOP}/share/cmake/Modules/" )
+ find_package( Bootstrap REQUIRED )
+ setup_project_paths( CORIOLIS )
+ print_cmake_module_path()
+
+ set_cmake_policies()
+ check_distribution()
+ setup_sysconfdir( "${CMAKE_INSTALL_PREFIX}" )
+
+ add_subdirectory( src )
+ add_subdirectory( cmake_modules )

flute/cmake_modules/CMakeLists.txt

@@ -0,0 +1,2 @@
+
+ install( FILES FindFLUTE.cmake DESTINATION share/cmake/Modules ) 

flute/cmake_modules/FindFLUTE.cmake

@@ -0,0 +1,31 @@
+# - Find the FLUTE includes and libraries.
+# The following variables are set if Coriolis is found.  If FLUTE is not
+# found, FLUTE_FOUND is set to false.
+#  FLUTE_FOUND        - True when the Coriolis include directory is found.
+#  FLUTE_INCLUDE_DIR - the path to where the Coriolis include files are.
+#  FLUTE_LIBRARIES - The path to where the Coriolis library files are.
+
+
+set( FLUTE_INCLUDE_PATH_DESCRIPTION "The directory containing the FLUTE include files. E.g /usr/local/include/coriolis/flute/3.1 or /soc/coriolis2/include/coriolis/flute/3.1" )
+set( FLUTE_DIR_MESSAGE              "Set the FLUTE_INCLUDE_DIR cmake cache entry to the ${FLUTE_INCLUDE_PATH_DESCRIPTION}" )
+
+if( UNIX )
+  #
+  # Look for an installation.
+  #
+  find_path( FLUTE_INCLUDE_PATH NAMES         flute.h
+                                PATHS         ${CORIOLIS_DIR_SEARCH}
+                                PATH_SUFFIXES include/flute/3.1
+                                DOC "${FLUTE_INCLUDE_PATH_DESCRIPTION}"
+  )
+
+  find_library( FLUTE_LIBRARY_PATH NAMES         flute
+                                   PATHS         ${CORIOLIS_DIR_SEARCH}
+                                   PATH_SUFFIXES lib${LIB_SUFFIX}
+                                   DOC           "${FLUTE_INCLUDE_PATH_DESCRIPTION}"
+  )
+
+  set_libraries_path       (FLUTE FLUTE)
+  hurricane_check_libraries(FLUTE)
+
+endif( UNIX )

flute/src/3.1/CMakeLists.txt

@@ -0,0 +1,31 @@
+# -*- mode: CMAKE explicit-buffer-name: "CMakeLists.txt<flute/src/3.1>" -*-
+
+   include_directories ( ${FLUTE_SOURCE_DIR}/src/3.1
+                       )
+
+                   set ( includes      flute.h
+                                       dl.h
+                                       mst2.h
+                                       err.h
+                                       heap.h
+                                       dist.h
+                                       global.h
+                                       neighbors.h
+                       )
+                   set ( cpps          flute.cpp
+                                       flute_mst.cpp
+                                       dist.cpp
+                                       dl.cpp
+                                       err.cpp
+                                       mst2.cpp
+                                       heap.cpp
+                                       neighbors.cpp
+                       )
+
+           add_library ( flute         ${cpps} )
+ set_target_properties ( flute         PROPERTIES VERSION 3.1 SOVERSION 3 )
+
+               install ( TARGETS       flute       DESTINATION lib${LIB_SUFFIX} )
+               install ( FILES         ${includes} DESTINATION include/flute/3.1 ) 
+               install ( FILES         POST9.dat
+                                       POWV9.dat   DESTINATION share/flute/3.1 ) 

flute/src/3.1/ChangeLog.txt

@@ -0,0 +1,40 @@
+Version 1.0 [Apr. 2004]
+ * Wirelength estimation based on the idea of ICCAD 04 paper.
+
+Version 2.0 [Feb. 16, 2005]
+ * RSMT construction and Wirelength estimation based on the idea of ISPD
+   05 paper.
+
+Version 2.1 [Mar. 18, 2005]
+ * The program flute-ckt and associated files are added to compute the
+   FLUTE wirelength for a circuit in bookshelf format.
+
+Version 2.2 [Oct. 16, 2005]  
+ * The net breaking technique is improved by incorporating the HPWL based
+   estimation in ICCAD 04 paper.
+ * The option to remove duplicated pins is added.
+ * The code is reorganized.
+
+Version 2.3 [Sep. 23, 2006]  
+ * PORT is renamed to POST (Potentially Optimal Steiner Tree)
+ * The function readLUT() is rewritten to significantly cut the POWV 
+   and POST reading time.
+
+Version 2.4 [Jan. 16, 2007]  
+ * Internal parameters of flutes_wl_MD() and flutes_MD() are tuned to
+   improve runtime-quality tradeoff.
+ * The file size of POST9.dat is reduced by improved encoding.
+ * The function readLUT() is improved to speed up the reading of POST9.dat.
+ * Bookshelf files for ibm01 is included to demonstrate program flute-ckt.
+
+Version 2.5 [Jun. 20, 2007] 
+ * Local refinement technique is added to improve runtime-accuracy
+   tradeoff for high accuracy.
+
+Version 3.0 [Apr. 10, 2008]
+ * Implemented the net breaking and merging ideas of VLSIDAT 08 paper to
+   achieve much better scalability and accuracy for high-degree nets.
+
+Version 3.1 [Feb. 28, 2011]
+ * Functions in bookshelf_IO.c and memAlloc.c used by flute-ckt are
+   improved to speed up the reading of GSRC Bookshelf format circuits.

flute/src/3.1/Makefile

@@ -0,0 +1,26 @@
+CC = gcc
+CFLAGS = -O3 -I.
+
+SRC     = dist.c dl.c err.c heap.c mst2.c neighbors.c \
+	bookshelf_IO.c memAlloc.c flute.c flute_mst.c
+OBJ     = $(SRC:.c=.o)
+
+all: flute-net flute-ckt rand-pts
+
+flute-net: flute-net.c ${OBJ}
+	$(CC) $(CFLAGS) -o flute-net flute-net.c ${OBJ} -lm
+
+flute-ckt: flute-ckt.c ${OBJ}
+	$(CC) $(CFLAGS) -o flute-ckt flute-ckt.c ${OBJ} -lm
+
+rand-pts: rand-pts.c
+	$(CC) $(CFLAGS) -o rand-pts rand-pts.c 
+
+flute.o: flute.c flute.h
+	$(CC) $(CFLAGS) -c -o flute.o flute.c
+
+flute_mst.o: flute_mst.c flute.h
+	$(CC) $(CFLAGS) -c -o flute_mst.o flute_mst.c
+
+clean:
+	rm -f *.o core*

flute/src/3.1/Readme

@@ -0,0 +1,40 @@
+-------------- FLUTE - Version 3.1 -----------------
+                       by
+                Chris C.-N. Chu
+       Dept. of ECpE, Iowa State University
+             Copyright (c) - 2005
+  Iowa State University Research Foundation, Inc.
+----------------------------------------------------
+
+This package contains the following files:
+
+ flute.c -- The rectilinear Steiner minimal tree and wirelength estimation
+	    algorithm described in the ICCAD 04 and ISPD 05 papers with
+	    some improvements described in TCAD 07 paper.
+ flute.h -- The interface to use flute.
+ flute_mst.c -- The net breaking and merging techniques described in the
+	    VLSIDAT 08 paper.
+ dist.[ch], dl.[ch], err.[ch], heap.[ch], mst2.[ch], neighbors.[ch],
+	    global.h -- Utility functions used by flute_mst.c
+ POWV9.dat -- The lookup-table of optimal POWVs up to degree 9.
+ POST9.dat -- The lookup-table for optimal Steiner trees up to degree 9.
+ flute-net.c -- A program to evaluate the wirelength of a net. It takes
+	    input from stdin as a list of points.
+ rand-pts.c -- A program to generate a list of random points.
+ flute-ckt.c -- A program to find FLUTE and half-perimeter wirelength
+	    of a circuit in bookshelf format.
+ bookshelf_IO.[ch] -- Functions for flute-ckt.c to read bookshelf files.
+ memAlloc.[ch] -- Functions for flute-ckt.c to allocate memory.
+ ibm01/ibm01.* -- ibm01 bookshelf files that can be read by flute-ckt.c
+ license.txt -- License agreement.
+ ChangeLog.txt
+ Makefile
+ Readme
+
+To run the programs, first do a 'make'. POWV9.dat and POST9.dat is assume
+to be in the current directory.  Some example commands:
+
+ rand-pts | flute-net
+ rand-pts 20 | flute-net	    // 20-pin nets
+ rand-pts -r 20 | flute-net	    // randomized
+ flute-ckt ibm01 ibm01.aux ibm01/ibm01.pl

flute/src/3.1/bookshelf_IO.h

@@ -0,0 +1,166 @@
+/* -----------  FastPlace - Version 3.1 ----------------
+                       by 
+   Natarajan Viswanathan and Chris C.-N. Chu
+     Dept. of ECpE, Iowa State University
+          Copyright (c) - 2004 
+Iowa State University Research Foundation, Inc.
+-------------------------------------------------------- */
+/* --------------------------------------------------------------------------
+   Header file used in bookshelf_IO.c 
+----------------------------------------------------------------------------*/
+
+#ifndef FLUTE_BOOKSHELF_IO_H
+#define FLUTE_BOOKSHELF_IO_H
+
+namespace Flute {
+
+#define BUFFERSIZE 1000
+#define LINESIZE 2000
+
+    /* -----------------------------------------------------------------------------
+        Reads the .nodes file and creates a hash linking cell name to cell index for 
+        all the nodes in the circuit (movable nodes + fixed nodes + I/O pads)
+
+        creates extern vars:
+            cellName[i]     (i = 1..movableNodes + numTerminals)
+    -------------------------------------------------------------------------------- */
+    extern void createHash(char benchmarkPath[], char nodesFile[]);
+    extern void freeHash();
+
+    /* -----------------------------------------------------------------------------
+        Reads the .aux file to get the other file names
+  
+        creates extern vars:
+            nodesFile[], netsFile[], wtsFile[], sclFile[], plFile[], benchmarkName[]
+    -------------------------------------------------------------------------------- */
+    extern void readAuxFile(char benchmarkPath[], char auxFile[]);
+
+    /* -----------------------------------------------------------------------------
+        Reads the .nodes file to get cell widths and heights
+        
+        creates extern vars: 
+            movableNodes, numTerminals, averageCellWidth, cellWidth[], cellHeight[]
+            cellArea[]
+    -------------------------------------------------------------------------------- */
+    extern void readNodesFile(char benchmarkPath[], char nodesFile[]);
+
+    /* -----------------------------------------------------------------------------
+        Reads the .nets file to get the netlist information
+   
+        creates extern vars: 
+            numNets, numPins, netlist[], netlistIndex[], xPinOffset[], yPinOffset[]
+    -------------------------------------------------------------------------------- */
+    extern void readNetsFile(char benchmarkPath[], char netsFile[]);
+
+    /* -----------------------------------------------------------------------------
+        Reads the .pl file to get coordinates of all the nodes and the placement 
+        boundary based on the position of the I/O pads
+  
+        creates extern vars:
+            xCellCoord[], yCellCoord[]
+    -------------------------------------------------------------------------------- */
+    extern void readPlFile(char benchmarkPath[], char plFile[]);
+    
+    /* -----------------------------------------------------------------------------
+        Reads the .scl file to get placement (core) region information
+  
+        creates extern vars:
+            siteOriginX, siteEndX, siteOriginY, siteEndY
+            minX, maxX, minY, maxY, numRows
+            coreRowHeight, coreWidth, coreHeight 
+            core_row data structure
+    -------------------------------------------------------------------------------- */
+    extern void readSclFile(char benchmarkPath[], char sclFile[]);
+    
+    /* -----------------------------------------------------------
+        Reads the .wts file to get net-weights
+        
+        creates extern vars:
+            netWeight[]
+    ----------------------------------------------------------- */
+    extern void readWtsFile(char benchmarkPath[], char wtsFile[], int net_weights);
+    
+    /* -----------------------------------------------------------------------------
+        writes out a bookshelf format .pl file
+    -------------------------------------------------------------------------------- */
+    extern void writePlFile(char outputDir[], char fileName[], float xCoord[], float yCoord[], 
+                            int totalNodes);
+    
+    
+    /*--------------  Extern Variables  ------------------*/
+
+    extern char **cellName;
+    extern char *nodesFile, *netsFile, *wtsFile, *sclFile, *plFile, *benchmarkName;
+    
+    extern int movableNodes, numTerminals;
+    extern float averageCellWidth, *cellWidth, *cellHeight, *cellArea; 
+    
+    extern int numNets, numPins, *netlist, *netlistIndex;
+    extern float *xPinOffset, *yPinOffset;
+    extern char **netName;
+       
+    extern float *xCellCoord, *yCellCoord;
+    extern char **cellOrientation;
+    
+    extern int numRows, numRowBlockages;
+    extern float siteOriginY, siteEndY, coreHeight;
+    extern float siteOriginX, siteEndX, coreWidth;
+    extern float coreRowHeight, minX, maxX, minY, maxY;
+    extern float *xRowBlockage, *yRowBlockage, *widthRowBlockage;
+    typedef struct row_data ROW;
+    struct row_data {
+        float x_low;
+        float x_high;
+        float y_low;
+        float height;
+        float siteWidth;
+        float siteSpacing;
+        float subrow_origin[100], subrow_end[100];
+        int subrow_count;
+        char siteOrient[3], siteSymmetry[3];
+    };
+    extern ROW *core_row;
+    
+    extern float *netWeight;
+
+}  // Flute namespace.
+
+#endif /* FLUTE_BOOKSHELF_IO_H */ 
+
+/* -----------------------------------------------------------------------------------------------
+                                Description of Extern Variables
+
+    cellName[i]         =   cell name corresponding to cell index "i" 
+                            (i = 1..movableNodes+numTerminals)
+
+    movableNodes        =   number of movable nodes,
+    numTerminals        =   number of fixed nodes (I/O Pads + Fixed Macros)
+    averageCellWidth    =   avg width of movable nodes,
+    cellWidth[i]        =   width of cell "i"   (i = 1..movableNodes+numTerminals)
+    cellHeight[i]       =   height of cell "i"  (i = 1..movableNodes+numTerminals)
+    cellArea[i]         =   area of cell "i"    (i = 1..movableNodes+numTerminals)
+
+    numNets             =   number of nets
+    numPins             =   number of pins
+    netlist[]           =   netlist of the circuit
+    xPinOffset[]        =   x-offset of the pins from the center of the cell
+    yPinOffset[]        =   y-offset (      "       )
+    netlistIndex[]      =   index to the netlist and offset vectors
+
+    xCellCoord[i]       =   x-coordinate of cell "i"  (i = 1..movableNodes+numTerminals)
+    yCellCoord[i]       =   y-coordinate of cell "i", 
+    minX, maxX          =   left and right boundaries of the chip (Note: not the placement region) 
+    minY, maxY          =   bottom and top boundaries of the chip
+    
+    siteOriginX         =   left boundary of the placement region
+    siteEndX            =   right boundary of the placement region   
+    siteOriginY         =   bottom boundary of the placement region    
+    siteEndY            =   top boundary of the placement region   
+    siteWidth           =   width of a placement site within a row 
+    siteSpacing         =   the space b/w adjacent placement sites within a row 
+    numRows             =   number of placement rows 
+    coreRowHeight       =   row Height 
+    coreWidth           =   siteEndX - siteOriginX 
+    coreHeight          =   siteEndY - siteOriginY 
+
+---------------------------------------------------------------------------------------------------*/

flute/src/3.1/dist.cpp

@@ -0,0 +1,48 @@
+#include "global.h"
+
+namespace Flute {
+
+/*********************************************************************/
+/*
+   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; 
+}
+
+/*********************************************************************/
+/*********************************************************************/
+
+}  // Flute namespace.

flute/src/3.1/dist.h

@@ -0,0 +1,20 @@
+#ifndef FLUTE_DIST_H
+#define FLUTE_DIST_H
+
+#include "global.h"
+
+namespace Flute {
+
+long  dist(
+  Point  p,
+  Point  q
+);
+
+long  dist2(   
+  Point*  p,
+  Point*  q
+);
+
+}  // Flute namespace.
+
+#endif

flute/src/3.1/dl.cpp

@@ -0,0 +1,165 @@
+#include <dl.h>
+#include <assert.h>
+#include <stdio.h>
+
+namespace Flute {
+
+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; i<len; i++) {
+    el = el->next;
+  }
+
+  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;
+      }
+    }
+  }
+}
+
+}  // Flute namespace.

flute/src/3.1/dl.h

@@ -0,0 +1,184 @@
+#ifndef FLUTE_DL_H
+#define FLUTE_DL_H
+
+#include <string.h>
+#include <stdlib.h>
+
+namespace Flute {
+
+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; \
+}
+
+}  // Flute namespace.
+
+#endif

flute/src/3.1/err.cpp

@@ -0,0 +1,31 @@
+#include <stdio.h>
+#include <stdlib.h>
+
+namespace Flute {
+
+/**************************************************************************/
+/*
+  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);
+}
+
+}  // Flute namespace.

flute/src/3.1/err.h

@@ -0,0 +1,16 @@
+#ifndef FLUTE_ERR_H
+#define FLUTE_ERR_H
+
+namespace Flute {
+
+void  err_msg(
+char* msg
+);
+
+void  err_exit(
+char* msg
+);
+
+}  // Flute namespace.
+
+#endif

flute/src/3.1/flute-ckt.cpp

@@ -0,0 +1,98 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "bookshelf_IO.h"
+#include "memAlloc.h"
+#include "flute.h"
+
+using namespace Flute;
+
+float HPwl();
+float flutewl();
+
+int main (int argc, char *argv[])
+{
+    char benchmarkPath[BUFFERSIZE], auxFile[BUFFERSIZE], placefile[BUFFERSIZE];
+
+    if(argc != 4) {
+        printf("Usage: %s <benchmark_dir> <aux_file> <placement_file>\n",
+               argv[0]);
+        printf("    <benchmark_dir> is the benchmark file directory.\n");
+        printf("    <aux_file> is the bookshelf format auxiliary file");
+        printf(" (assume in <benchmark_dir>).\n");
+        printf("    <placement_file> is the placement file");
+        printf(" (assume in current directory).\n");
+        exit(1);
+    }    
+
+    strcpy(benchmarkPath, argv[1]);
+    strcpy(auxFile, argv[2]);
+    strcpy(placefile, argv[3]);
+
+    readAuxFile(benchmarkPath, auxFile);
+    createHash(benchmarkPath, nodesFile);
+    readNodesFile(benchmarkPath, nodesFile);
+    readNetsFile(benchmarkPath, netsFile);
+    readPlFile(".", placefile);
+    freeHash();
+
+    readLUT();
+
+    printf("Half-perimeter wirelength: %.2f\n", HPwl());
+    printf("FLUTE wirelength         : %.2f\n", flutewl());
+}
+
+float HPwl()
+{
+    float totx, toty, xu, xl, yu, yl, xOffset, yOffset;
+    int i, j, k;
+
+    totx = 0.0; toty = 0.0;
+    for (j=1; j<=numNets; j++) {
+        xl = yl = 1e9;
+        xu = yu = -1e9;
+        for (k=netlistIndex[j]; k<netlistIndex[j+1]; k++) {
+            i = netlist[k];
+            xOffset = xPinOffset[k];
+            yOffset = yPinOffset[k];
+
+            if (xCellCoord[i]+xOffset < xl) xl = xCellCoord[i]+xOffset;
+            if (xu < xCellCoord[i]+xOffset) xu = xCellCoord[i]+xOffset;
+            if (yCellCoord[i]+yOffset < yl) yl = yCellCoord[i]+yOffset;
+            if (yu < yCellCoord[i]+yOffset) yu = yCellCoord[i]+yOffset;
+        }
+        totx += xu - xl;
+        toty += yu - yl;
+    }
+    
+    return totx + toty;
+}
+
+float flutewl()
+{
+    Tree t;
+    DTYPE totwl;
+    DTYPE x[MAXD], y[MAXD];
+    float xOffset, yOffset;
+    int i, j, k, r, d;
+
+    totwl = 0;
+    for (j=1; j<=numNets; j++) {
+        d = netlistIndex[j+1] - netlistIndex[j];
+        k = netlistIndex[j]; 
+        for (r=0; r<d; r++) {
+            i = netlist[k+r];
+            xOffset = xPinOffset[k+r];
+            yOffset = yPinOffset[k+r];
+            x[r] = (DTYPE) xCellCoord[i]+xOffset;
+            y[r] = (DTYPE) yCellCoord[i]+yOffset;
+        }
+#if ROUTING==1
+        t = flute(d, x, y, ACCURACY); totwl += t.length;
+#else        
+        totwl += flute_wl(d, x, y, ACCURACY);
+#endif        
+    }
+    
+    return (float) totwl;
+}    

flute/src/3.1/flute-net.cpp

@@ -0,0 +1,25 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include "flute.h"
+
+using namespace Flute;
+
+int main()
+{
+    int d=0;
+    int x[MAXD], y[MAXD];
+    Tree flutetree;
+    int flutewl;
+    
+    while (!feof(stdin)) {
+        scanf("%d %d\n", &x[d], &y[d]);
+        d++;
+    }
+    readLUT();
+
+    flutetree = flute(d, x, y, ACCURACY);
+    printf("FLUTE wirelength = %d\n", flutetree.length);
+
+    flutewl = flute_wl(d, x, y, ACCURACY);
+    printf("FLUTE wirelength (without RSMT construction) = %d\n", flutewl);
+}

flute/src/3.1/flute.h

@@ -0,0 +1,108 @@
+
+#ifndef FLUTE_FLUTE_H
+#define FLUTE_FLUTE_H
+
+#include <string>
+
+namespace Flute {
+
+  using std::string;
+  
+
+/*****************************/
+/*  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(string);
+// DTYPE flute_wl(int d, DTYPE x[], DTYPE y[], int acc);
+// DTYPE flutes_wl(int d, DTYPE xs[], DTYPE ys[], int s[], int acc);
+// Tree flute(int d, DTYPE x[], DTYPE y[], int acc);
+// Tree flutes(int d, DTYPE xs[], DTYPE ys[], int s[], int acc);
+// DTYPE wirelength(Tree t);
+// void printtree(Tree t);
+// void plottree(Tree t);
+
+
+/*************************************/
+/* Internal Parameters and Functions */
+/*************************************/
+#define POWVFILE "POWV9.dat"        // LUT for POWV (Wirelength Vector)
+#define POSTFILE "POST9.dat"        // LUT for POST (Steiner Tree)
+#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
+    int n;   // index of neighbor
+} Branch;
+
+struct Tree
+{
+    int deg;   // degree
+    DTYPE length;   // total wirelength
+    Branch *branch;   // array of tree branches
+};
+
+// User-Callable Functions
+extern void readLUT(string directory);
+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);
+extern Tree flute(int d, DTYPE x[], DTYPE y[], int acc);
+//Macro: Tree flutes(int d, DTYPE xs[], DTYPE ys[], int s[], int acc);
+extern DTYPE wirelength(Tree t);
+extern void printtree(Tree t);
+extern void plottree(Tree 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 Tree flutes_LD(int d, DTYPE xs[], DTYPE ys[], int s[]);
+extern Tree flutes_MD(int d, DTYPE xs[], DTYPE ys[], int s[], int acc);
+extern Tree flutes_HD(int d, DTYPE xs[], DTYPE ys[], int s[], int acc);
+extern Tree flutes_RDP(int d, DTYPE xs[], DTYPE ys[], int s[], int acc);
+
+#if REMOVE_DUPLICATE_PIN==1
+  #define flutes_wl(d, xs, ys, s, acc) flutes_wl_RDP(d, xs, ys, s, acc) 
+  #define flutes(d, xs, ys, s, acc) flutes_RDP(d, xs, ys, s, acc) 
+#else
+  #define flutes_wl(d, xs, ys, s, acc) flutes_wl_ALLD(d, xs, ys, s, acc) 
+  #define flutes(d, xs, ys, s, acc) flutes_ALLD(d, xs, ys, s, 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) \
+    (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
+
+}  // Flute namespace.
+
+#endif /* FLUTE_FLUTE_H */

flute/src/3.1/global.h

@@ -0,0 +1,23 @@
+#ifndef  FLUTE_GLOBAL_H
+#define  FLUTE_GLOBAL_H
+
+#include <stdio.h>
+
+namespace Flute {
+
+#define  TRUE   1
+#define  FALSE  0
+#define  MAXLONG  0x7fffffffL
+
+struct point
+{
+  long    x, y;
+};
+
+typedef  struct point Point;
+
+typedef  long nn_array[8];
+
+}  // Flute namespace.
+
+#endif  /* FLUTE_GLOBAL_H */

flute/src/3.1/heap.cpp

@@ -0,0 +1,179 @@
+/****************************************************************************/
+/*
+  Binary heap routines for use in Prim's algorithm, 
+  with points are numbered from 0 to n-1
+*/
+
+#include <stdlib.h>
+#include "heap.h"
+#include "err.h"
+
+namespace Flute {
+
+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() */
+
+
+/****************************************************************************/
+
+}  // Flute namespace.

flute/src/3.1/heap.h

@@ -0,0 +1,35 @@
+#ifndef  FLUTE_HEAP_H
+#define  FLUTE_HEAP_H
+
+#include "global.h"
+
+namespace Flute {
+
+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();
+
+}  // Flute namespace.
+
+#endif /* FLUTE_HEAP_H */

flute/src/3.1/license.txt

@@ -0,0 +1,50 @@
+READ THIS LICENSE AGREEMENT CAREFULLY BEFORE USING THIS PRODUCT. BY USING
+THIS PRODUCT YOU INDICATE YOUR ACCEPTANCE OF THE TERMS OF THE FOLLOWING
+AGREEMENT. THESE TERMS APPLY TO YOU AND ANY SUBSEQUENT LICENSEE OF THIS
+PRODUCT.
+
+License Agreement for FLUTE
+
+Copyright (c) 2004 by Dr. Chris C. N. Chu
+All rights reserved
+
+ATTRIBUTION ASSURANCE LICENSE (adapted from the original BSD license)
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the conditions below are
+met. These conditions require a modest attribution to Dr. Chris C. N. Chu
+(the "Author").
+
+1. Redistributions of the source code, with or without modification (the
+   "Code"), must be accompanied by any documentation and, each time
+   the resulting executable program or a program dependent thereon is
+   launched, a prominent display (e.g., splash screen or banner text) of
+   the Author's attribution information, which includes:
+   (a) Dr. Chris C. N. Chu ("AUTHOR"),
+   (b) Iowa State University ("PROFESSIONAL IDENTIFICATION"), and
+   (c) http://home.engineering.iastate.edu/~cnchu/ ("URL").
+
+2. Users who intend to use the Code for commercial purposes will notify
+   Author prior to such commercial use.
+
+3. Neither the name nor any trademark of the Author may be used to
+   endorse or promote products derived from this software without
+   specific prior written permission.
+
+4. Users are entirely responsible, to the exclusion of the Author and any
+   other persons, for compliance with (1) regulations set by owners or
+   administrators of employed equipment, (2) licensing terms of any other
+   software, and (3) local, national, and international regulations
+   regarding use, including those regarding import, export, and use of
+   encryption software.
+
+THIS FREE SOFTWARE IS PROVIDED BY THE AUTHOR "AS IS" AND ANY EXPRESS OR
+IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
+OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
+IN NO EVENT SHALL THE AUTHOR OR ANY CONTRIBUTOR BE LIABLE FOR ANY DIRECT,
+INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+(INCLUDING, BUT NOT LIMITED TO, EFFECTS OF UNAUTHORIZED OR MALICIOUS
+NETWORK ACCESS; PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
+THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

flute/src/3.1/memAlloc.cpp

@@ -0,0 +1,468 @@
+/* --------------------------------------------------------------------------
+   Public domain memory allocation and de-allocation routines.
+   Taken from Appendix B of: 
+   Numerical Recipes in C: The Art of Scientific Computing, Second Edition,
+   Cambridge University Press, 1992
+----------------------------------------------------------------------------*/
+#include <stdio.h>
+#include <stddef.h>
+#include <stdlib.h>
+
+#include "memAlloc.h"
+
+namespace Flute {
+
+#define MEM_END 1
+#define FREE_ARG char*
+
+void runtimeError(char error_text[])
+/* error handler */
+{
+  fprintf(stderr, "ERROR: %s \n", error_text);
+  fprintf(stderr, "Aborting !! \n");
+  fflush(stdout);
+  fflush(stderr);
+  exit(1);
+}
+
+float *vector(long nl, long nh)
+/* allocate a float vector with subscript range v[nl..nh] */
+{
+  float *v;
+  v=(float *)malloc((size_t) ((nh-nl+1+MEM_END)*sizeof(float)));
+  if (!v) runtimeError("allocation failure in vector()");
+  return v-nl+MEM_END;
+}
+
+int *ivector(long nl, long nh)
+/* allocate an int vector with subscript range v[nl..nh] */
+{
+  int *v;
+  v=(int *)malloc((size_t) ((nh-nl+1+MEM_END)*sizeof(int)));
+  if (!v) runtimeError("allocation failure in ivector()");
+  return v-nl+MEM_END;
+}
+
+unsigned char *cvector(long nl, long nh)
+/* allocate an unsigned char vector with subscript range v[nl..nh] */
+{
+  unsigned char *v;
+  v=(unsigned char *)malloc((size_t) ((nh-nl+1+MEM_END)*sizeof(unsigned char)));
+  if (!v) runtimeError("allocation failure in cvector()");
+  return v-nl+MEM_END;
+}
+
+unsigned long *lvector(long nl, long nh)
+/* allocate an unsigned long vector with subscript range v[nl..nh] */
+{
+  unsigned long *v;
+  v=(unsigned long *)malloc((size_t) ((nh-nl+1+MEM_END)*sizeof(long)));
+  if (!v) runtimeError("allocation failure in lvector()");
+  return v-nl+MEM_END;
+}
+
+double *dvector(long nl, long nh)
+/* allocate a double vector with subscript range v[nl..nh] */
+{
+  double *v;
+  v=(double *)malloc((size_t) ((nh-nl+1+MEM_END)*sizeof(double)));
+  if (!v) runtimeError("allocation failure in dvector()");
+  return v-nl+MEM_END;
+}
+
+float **matrix(long nrl, long nrh, long ncl, long nch)
+/* allocate a float matrix with subscript range m[nrl..nrh][ncl..nch] */
+{
+  long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
+  float **m;
+  
+  /* allocate pointers to rows */
+  m=(float **) malloc((size_t)((nrow+MEM_END)*sizeof(float*)));
+  if (!m) runtimeError("allocation failure 1 in matrix()");
+  m += MEM_END;
+  m -= nrl;
+  
+  /* allocate rows and set pointers to them */
+  m[nrl]=(float *) malloc((size_t)((nrow*ncol+MEM_END)*sizeof(float)));
+  if (!m[nrl]) runtimeError("allocation failure 2 in matrix()");
+  m[nrl] += MEM_END;
+  m[nrl] -= ncl;
+  for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
+  /* return pointer to array of pointers to rows */
+  return m;
+}
+
+double **dmatrix(long nrl, long nrh, long ncl, long nch)
+/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
+{
+  long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
+  double **m;
+
+  /* allocate pointers to rows */
+  m=(double **) malloc((size_t)((nrow+MEM_END)*sizeof(double*)));
+  if (!m) runtimeError("allocation failure 1 in dmatrix()");
+  m += MEM_END;
+  m -= nrl;
+  
+  /* allocate rows and set pointers to them */
+  m[nrl]=(double *) malloc((size_t)((nrow*ncol+MEM_END)*sizeof(double)));
+  if (!m[nrl]) runtimeError("allocation failure 2 in dmatrix()");
+  m[nrl] += MEM_END;
+  m[nrl] -= ncl;
+  
+  for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
+  
+  /* return pointer to array of pointers to rows */
+  return m;
+}
+
+int **imatrix(long nrl, long nrh, long ncl, long nch)
+/* allocate a int matrix with subscript range m[nrl..nrh][ncl..nch] */
+{
+  long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
+  int **m;
+
+  /* allocate pointers to rows */
+  m=(int **) malloc((size_t)((nrow+MEM_END)*sizeof(int*)));
+  if (!m) runtimeError("allocation failure 1 in imatrix()");
+  m += MEM_END;
+  m -= nrl;
+  
+  /* allocate rows and set pointers to them */
+  m[nrl]=(int *) malloc((size_t)((nrow*ncol+MEM_END)*sizeof(int)));
+  if (!m[nrl]) runtimeError("allocation failure 2 in imatrix()");
+  m[nrl] += MEM_END;
+  m[nrl] -= ncl;
+  
+  for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
+  
+  /* return pointer to array of pointers to rows */
+  return m;
+}
+
+
+char **cmatrix(long nrl, long nrh, long ncl, long nch)
+/* allocate a char matrix with subscript range m[nrl..nrh][ncl..nch] */
+{
+  long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
+  char **m;
+
+  /* allocate pointers to rows */
+  m=(char **) malloc((size_t)((nrow+MEM_END)*sizeof(char*)));
+  if (!m) runtimeError("allocation failure 1 in cmatrix()");
+  m += MEM_END;
+  m -= nrl;
+  
+  /* allocate rows and set pointers to them */
+  m[nrl]=(char *) malloc((size_t)((nrow*ncol+MEM_END)*sizeof(char)));
+  if (!m[nrl]) runtimeError("allocation failure 2 in cmatrix()");
+  m[nrl] += MEM_END;
+  m[nrl] -= ncl;
+  
+  for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
+  
+  /* return pointer to array of pointers to rows */
+  return m;
+}
+
+
+unsigned long **lmatrix(long nrl, long nrh, long ncl, long nch)
+/* allocate a int matrix with subscript range m[nrl..nrh][ncl..nch] */
+{
+  long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
+  unsigned long **m;
+
+  /* allocate pointers to rows */
+  m=(unsigned long **) malloc((size_t)((nrow+MEM_END)*sizeof(long*)));
+  if (!m) runtimeError("allocation failure 1 in lmatrix()");
+  m += MEM_END;
+  m -= nrl;
+  
+  /* allocate rows and set pointers to them */
+  m[nrl]=(unsigned long *) malloc((size_t)((nrow*ncol+MEM_END)*sizeof(long)));
+  if (!m[nrl]) runtimeError("allocation failure 2 in lmatrix()");
+  m[nrl] += MEM_END;
+  m[nrl] -= ncl;
+  
+  for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
+  
+  /* return pointer to array of pointers to rows */
+  return m;
+}
+
+
+float **submatrix(float **a, long oldrl, long oldrh, long oldcl, long oldch, 
+  long newrl, long newcl)
+/* point a submatrix [newrl..][newcl..] to a[oldrl..oldrh][oldcl..oldch] */
+{
+  long i,j,nrow=oldrh-oldrl+1,ncol=oldcl-newcl;
+  float **m;
+  
+  /* allocate array of pointers to rows */
+  m=(float **) malloc((size_t) ((nrow+MEM_END)*sizeof(float*)));
+  if (!m) runtimeError("allocation failure in submatrix()");
+  m += MEM_END;
+  m -= newrl;
+  
+  /* set pointers to rows */
+  for(i=oldrl,j=newrl;i<=oldrh;i++,j++) m[j]=a[i]+ncol;
+  
+  /* return pointer to array of pointers to rows */
+  return m;
+}
+
+float **convert_matrix(float *a, long nrl, long nrh, long ncl, long nch)
+/* allocate a float matrix m[nrl..nrh][ncl..nch] that points to the matrix
+declared in the standard C manner as a[nrow][ncol], where nrow=nrh-nrl+1
+and ncol=nch-ncl+1. The routine should be called with the address
+&a[0][0] as the first argument. */
+{
+  long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1;
+  float **m;
+
+  /* allocate pointers to rows */
+  m=(float **) malloc((size_t) ((nrow+MEM_END)*sizeof(float*)));
+  if (!m) runtimeError("allocation failure in convert_matrix()");
+  m += MEM_END;
+  m -= nrl;
+  
+  /* set pointers to rows */
+  m[nrl]=a-ncl;
+  for(i=1,j=nrl+1;i<nrow;i++,j++) m[j]=m[j-1]+ncol;
+  
+  /* return pointer to array of pointers to rows */
+  return m;
+}
+
+float ***f3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh)
+/* allocate a float 3tensor with range t[nrl..nrh][ncl..nch][ndl..ndh] */
+{
+  long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1,ndep=ndh-ndl+1;
+  float ***t;
+
+  /* allocate pointers to pointers to rows */
+  t=(float ***) malloc((size_t)((nrow+MEM_END)*sizeof(float**)));
+  if (!t) runtimeError("allocation failure 1 in f3tensor()");
+  t += MEM_END;
+  t -= nrl;
+  
+  /* allocate pointers to rows and set pointers to them */
+  t[nrl]=(float **) malloc((size_t)((nrow*ncol+MEM_END)*sizeof(float*)));
+  if (!t[nrl]) runtimeError("allocation failure 2 in f3tensor()");
+  t[nrl] += MEM_END;
+  t[nrl] -= ncl;
+  
+  /* allocate rows and set pointers to them */
+  t[nrl][ncl]=(float *) malloc((size_t)((nrow*ncol*ndep+MEM_END)*sizeof(float)));
+  if (!t[nrl][ncl]) runtimeError("allocation failure 3 in f3tensor()");
+  t[nrl][ncl] += MEM_END;
+  t[nrl][ncl] -= ndl;
+  
+  for(j=ncl+1;j<=nch;j++) t[nrl][j]=t[nrl][j-1]+ndep;
+  for(i=nrl+1;i<=nrh;i++) {
+    t[i]=t[i-1]+ncol;
+    t[i][ncl]=t[i-1][ncl]+ncol*ndep;
+    for(j=ncl+1;j<=nch;j++) t[i][j]=t[i][j-1]+ndep;
+  }
+  
+  /* return pointer to array of pointers to rows */
+  return t;
+}
+
+
+unsigned long ***lmatrix3D(long nrl, long nrh, long ncl, long nch, long ndl, long ndh)
+/* allocate an unsigned long 3D matrix with range t[nrl..nrh][ncl..nch][ndl..ndh] */
+{
+  long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1,ndep=ndh-ndl+1;
+  unsigned long ***t;
+  long d1, d2, d3;
+  
+  /* allocate pointers to pointers to rows */
+  t=(unsigned long ***) malloc((size_t)((nrow+MEM_END)*sizeof(long**)));
+  if (!t) runtimeError("allocation failure 1 in lmatrix3D()");
+  t += MEM_END;
+  t -= nrl;
+  
+  /* allocate pointers to rows and set pointers to them */
+  t[nrl]=(unsigned long **) malloc((size_t)((nrow*ncol+MEM_END)*sizeof(long*)));
+  if (!t[nrl]) runtimeError("allocation failure 2 in lmatrix3D()");
+  t[nrl] += MEM_END;
+  t[nrl] -= ncl;
+  
+  /* allocate rows and set pointers to them */
+  t[nrl][ncl]=(unsigned long *) malloc((size_t)((nrow*ncol*ndep+MEM_END)*sizeof(long)));
+  if (!t[nrl][ncl]) runtimeError("allocation failure 3 in lmatrix3D()");
+  t[nrl][ncl] += MEM_END;
+  t[nrl][ncl] -= ndl;
+  
+  for(j=ncl+1;j<=nch;j++) t[nrl][j]=t[nrl][j-1]+ndep;
+  for(i=nrl+1;i<=nrh;i++) {
+    t[i]=t[i-1]+ncol;
+    t[i][ncl]=t[i-1][ncl]+ncol*ndep;
+    for(j=ncl+1;j<=nch;j++) t[i][j]=t[i][j-1]+ndep;
+  }  
+  
+  /* return pointer to array of pointers to rows */
+  return t;
+}
+
+
+int ***imatrix3D(int nrl, int nrh, int ncl, int nch, int ndl, int ndh)
+/* allocate an int 3D matrix with range t[nrl..nrh][ncl..nch][ndl..ndh] */
+{
+  int i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1,ndep=ndh-ndl+1;
+  int ***t;
+  int d1, d2, d3;
+  
+  /* allocate pointers to pointers to rows */
+  t=(int ***) malloc((size_t)((nrow+MEM_END)*sizeof(int**)));
+  if (!t) runtimeError("allocation failure 1 in imatrix3D()");
+  t += MEM_END;
+  t -= nrl;
+  
+  /* allocate pointers to rows and set pointers to them */
+  t[nrl]=(int **) malloc((size_t)((nrow*ncol+MEM_END)*sizeof(int*)));
+  if (!t[nrl]) runtimeError("allocation failure 2 in imatrix3D()");
+  t[nrl] += MEM_END;
+  t[nrl] -= ncl;
+  
+  /* allocate rows and set pointers to them */
+  t[nrl][ncl]=(int *) malloc((size_t)((nrow*ncol*ndep+MEM_END)*sizeof(int)));
+  if (!t[nrl][ncl]) runtimeError("allocation failure 3 in imatrix3D()");
+  t[nrl][ncl] += MEM_END;
+  t[nrl][ncl] -= ndl;
+  
+  for(j=ncl+1;j<=nch;j++) t[nrl][j]=t[nrl][j-1]+ndep;
+  for(i=nrl+1;i<=nrh;i++) {
+    t[i]=t[i-1]+ncol;
+    t[i][ncl]=t[i-1][ncl]+ncol*ndep;
+    for(j=ncl+1;j<=nch;j++) t[i][j]=t[i][j-1]+ndep;
+  }  
+  
+  /* return pointer to array of pointers to rows */
+  return t;
+}
+
+
+void free_vector(float *v, long nl, long nh)
+/* free a float vector allocated with vector() */
+{
+  free((FREE_ARG) (v+nl-MEM_END));
+}
+
+void free_ivector(int *v, long nl, long nh)
+/* free an int vector allocated with ivector() */
+{
+  free((FREE_ARG) (v+nl-MEM_END));
+}
+
+void free_cvector(unsigned char *v, long nl, long nh)
+/* free an unsigned char vector allocated with cvector() */
+{
+  free((FREE_ARG) (v+nl-MEM_END));
+}
+
+void free_lvector(unsigned long *v, long nl, long nh)
+/* free an unsigned long vector allocated with lvector() */
+{
+  free((FREE_ARG) (v+nl-MEM_END));
+}
+
+void free_dvector(double *v, long nl, long nh)
+/* free a double vector allocated with dvector() */
+{
+  free((FREE_ARG) (v+nl-MEM_END));
+}
+
+void free_matrix(float **m, long nrl, long nrh, long ncl, long nch)
+/* free a float matrix allocated by matrix() */
+{
+  free((FREE_ARG) (m[nrl]+ncl-MEM_END));
+  free((FREE_ARG) (m+nrl-MEM_END));
+}
+
+void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch)
+/* free a double matrix allocated by dmatrix() */
+{
+  free((FREE_ARG) (m[nrl]+ncl-MEM_END));
+  free((FREE_ARG) (m+nrl-MEM_END));
+}
+
+void free_imatrix(int **m, long nrl, long nrh, long ncl, long nch)
+/* free an int matrix allocated by imatrix() */
+{
+  free((FREE_ARG) (m[nrl]+ncl-MEM_END));
+  free((FREE_ARG) (m+nrl-MEM_END));
+}
+
+void free_cmatrix(char **m, long nrl, long nrh, long ncl, long nch)
+/* free a char matrix allocated by imatrix() */
+{
+  free((FREE_ARG) (m[nrl]+ncl-MEM_END));
+  free((FREE_ARG) (m+nrl-MEM_END));
+}
+
+void free_lmatrix(unsigned long **m, long nrl, long nrh, long ncl, long nch)
+/* free an unsigned long matrix allocated by lmatrix() */
+{
+  free((FREE_ARG) (m[nrl]+ncl-MEM_END));
+  free((FREE_ARG) (m+nrl-MEM_END));
+}
+
+void free_submatrix(float **b, long nrl, long nrh, long ncl, long nch)
+/* free a submatrix allocated by submatrix() */
+{
+  free((FREE_ARG) (b+nrl-MEM_END));
+}
+
+void free_convert_matrix(float **b, long nrl, long nrh, long ncl, long nch)
+/* free a matrix allocated by convert_matrix() */
+{
+  free((FREE_ARG) (b+nrl-MEM_END));
+}
+
+void free_f3tensor(float ***t, long nrl, long nrh, long ncl, long nch,
+  long ndl, long ndh)
+/* free a float f3tensor allocated by f3tensor() */
+{
+  free((FREE_ARG) (t[nrl][ncl]+ndl-MEM_END));
+  free((FREE_ARG) (t[nrl]+ncl-MEM_END));
+  free((FREE_ARG) (t+nrl-MEM_END));
+}
+
+
+void free_lmatrix3D(unsigned long ***t, long nrl, long nrh, long ncl, long nch,
+  long ndl, long ndh)
+/* free an unsigned long 3D matrix allocated by lmatrix3D() */
+{
+  free((FREE_ARG) (t[nrl][ncl]+ndl-MEM_END));
+  free((FREE_ARG) (t[nrl]+ncl-MEM_END));
+  free((FREE_ARG) (t+nrl-MEM_END));
+}
+
+
+void free_imatrix3D(int ***t, int nrl, int nrh, int ncl, int nch,
+  int ndl, int ndh)
+/* free an unsigned int 3D matrix allocated by imatrix3D() */
+{
+  free((FREE_ARG) (t[nrl][ncl]+ndl-MEM_END));
+  free((FREE_ARG) (t[nrl]+ncl-MEM_END));
+  free((FREE_ARG) (t+nrl-MEM_END));
+}
+
+
+int comp_float(const void *i, const void *j) {
+    if(*(float *)i < *(float *)j)
+        return -1;
+    else if(*(float *)i > *(float *)j)
+        return 1;
+    else
+        return 0;
+}
+
+
+int comp_int(const void *i, const void *j) {
+    return *(int *)i - *(int *)j;
+}
+
+}  // Flute namespace.

flute/src/3.1/memAlloc.h

@@ -0,0 +1,94 @@
+/* --------------------------------------------------------------------------
+   Public domain memory allocation and de-allocation routine header file.
+   Taken from Appendix B of: 
+   Numerical Recipes in C: The Art of Scientific Computing, Second Edition,
+   Cambridge University Press, 1992
+----------------------------------------------------------------------------*/
+#ifndef FLUTE_MEMALLOC_H
+#define FLUTE_MEMALLOC_H
+
+namespace Flute {
+
+static float sqrarg;
+#define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 : sqrarg*sqrarg)
+
+static double dsqrarg;
+#define DSQR(a) ((dsqrarg=(a)) == 0.0 ? 0.0 : dsqrarg*dsqrarg)
+
+static double dmaxarg1,dmaxarg2;
+#define DMAX(a,b) (dmaxarg1=(a),dmaxarg2=(b),(dmaxarg1) > (dmaxarg2) ?\
+(dmaxarg1) : (dmaxarg2))
+
+static double dminarg1,dminarg2;
+#define DMIN(a,b) (dminarg1=(a),dminarg2=(b),(dminarg1) < (dminarg2) ?\
+(dminarg1) : (dminarg2))
+
+static float maxarg1,maxarg2;
+#define FMAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1) > (maxarg2) ?\
+(maxarg1) : (maxarg2))
+
+static float minarg1,minarg2;
+#define FMIN(a,b) (minarg1=(a),minarg2=(b),(minarg1) < (minarg2) ?\
+(minarg1) : (minarg2))
+
+static long lmaxarg1,lmaxarg2;
+#define LMAX(a,b) (lmaxarg1=(a),lmaxarg2=(b),(lmaxarg1) > (lmaxarg2) ?\
+(lmaxarg1) : (lmaxarg2))
+
+static long lminarg1,lminarg2;
+#define LMIN(a,b) (lminarg1=(a),lminarg2=(b),(lminarg1) < (lminarg2) ?\
+(lminarg1) : (lminarg2))
+
+static int imaxarg1,imaxarg2;
+#define IMAX(a,b) (imaxarg1=(a),imaxarg2=(b),(imaxarg1) > (imaxarg2) ?\
+(imaxarg1) : (imaxarg2))
+
+static int iminarg1,iminarg2;
+#define IMIN(a,b) (iminarg1=(a),iminarg2=(b),(iminarg1) < (iminarg2) ?\
+(iminarg1) : (iminarg2))
+
+#define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
+#define MAX(a,b) ((a)>(b) ? (a) : (b))
+#define MIN(a,b) ((a)<(b) ? (a) : (b))
+
+void runtimeError(char error_text[]);
+float *vector(long nl, long nh);
+int *ivector(long nl, long nh);
+unsigned char *cvector(long nl, long nh);
+unsigned long *lvector(long nl, long nh);
+double *dvector(long nl, long nh);
+float **matrix(long nrl, long nrh, long ncl, long nch);
+double **dmatrix(long nrl, long nrh, long ncl, long nch);
+int **imatrix(long nrl, long nrh, long ncl, long nch);
+char **cmatrix(long nrl, long nrh, long ncl, long nch);
+unsigned long **lmatrix(long nrl, long nrh, long ncl, long nch);
+float **submatrix(float **a, long oldrl, long oldrh, long oldcl, long oldch,
+                  long newrl, long newcl);
+float **convert_matrix(float *a, long nrl, long nrh, long ncl, long nch);
+float ***f3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh);
+unsigned long ***lmatrix3D(long nrl, long nrh, long ncl, long nch, long ndl, long ndh);
+int ***imatrix3D(int nrl, int nrh, int ncl, int nch, int ndl, int ndh);
+
+void free_vector(float *v, long nl, long nh);
+void free_ivector(int *v, long nl, long nh);
+void free_cvector(unsigned char *v, long nl, long nh);
+void free_lvector(unsigned long *v, long nl, long nh);
+void free_dvector(double *v, long nl, long nh);
+void free_matrix(float **m, long nrl, long nrh, long ncl, long nch);
+void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch);
+void free_imatrix(int **m, long nrl, long nrh, long ncl, long nch);
+void free_cmatrix(char **m, long nrl, long nrh, long ncl, long nch);
+void free_lmatrix(unsigned long **m, long nrl, long nrh, long ncl, long nch);
+void free_submatrix(float **b, long nrl, long nrh, long ncl, long nch);
+void free_convert_matrix(float **b, long nrl, long nrh, long ncl, long nch);
+void free_f3tensor(float ***t, long nrl, long nrh, long ncl, long nch,
+                   long ndl, long ndh);
+void free_lmatrix3D(unsigned long ***t, long nrl, long nrh, long ncl, long nch,
+                   long ndl, long ndh);
+void free_imatrix3D(int ***t, int nrl, int nrh, int ncl, int nch, int ndl, int ndh);
+int comp_float(const void *i, const void *j);
+int comp_int(const void *i, const void *j);
+
+}  // Flute namespace.
+
+#endif /* FLUTE_MEMALLOC_H */

flute/src/3.1/mst2.cpp

@@ -0,0 +1,93 @@
+#include  <stdlib.h>
+#include  <stdio.h>
+#include   <assert.h>
+#include  "global.h"
+#include  "neighbors.h"
+#include  "dist.h"
+#include  "heap.h"
+#include  "err.h"
+
+namespace Flute {
+
+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;
+        }
+      }
+    }
+  }
+}
+
+/****************************************************************************/
+/****************************************************************************/
+
+}  // Flute namespace.

flute/src/3.1/mst2.h

@@ -0,0 +1,15 @@
+#ifndef FLUTE_MST2_H
+#define FLUTE_MST2_H
+
+#include "global.h"
+
+namespace Flute {
+
+void  mst2_package_init( long  n );
+void  mst2_package_done();
+void  mst2( long n, Point* pt, long* parent );
+
+}  // Flute namespace.
+
+#endif 
+

flute/src/3.1/neighbors.cpp

@@ -0,0 +1,530 @@
+#include  <assert.h>
+#include  <string.h>
+#include  <stdlib.h>
+#include  "global.h"
+#include  "err.h"
+#include  "dist.h"
+
+namespace Flute {
+
+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 );
+}
+
+/***************************************************************************/
+/***************************************************************************/
+
+}  // Flute namespace.

flute/src/3.1/neighbors.h

@@ -0,0 +1,28 @@
+
+#ifndef FLUTE_NEIGHBORS_H 
+#define FLUTE_NEIGHBORS_H 
+
+#include "global.h"
+
+namespace Flute {
+
+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
+);
+
+}  // Flute namespace.
+
+#endif /* FLUTE_NEIGHBORS_H */

flute/src/3.1/rand-pts.cpp

@@ -0,0 +1,35 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <unistd.h>
+
+int main(int ac, char *av[])
+{
+    int d=10, tmp, i;
+    int PNUM = 0;
+
+    for (i=1; i<ac; i++) {
+        if (strcmp(av[i], "-r")==0)  // random
+            srandom((int) getpid());
+        else if (strncmp(av[i], "-s", 2)==0)  // set random seed
+            srandom(atoi(av[i]+2));
+        else if (strcmp(av[i], "-n")==0)  // print # of points first
+            PNUM=1;
+        else if (sscanf(av[i], "%d", &tmp))  // set # of points
+            d = tmp;
+        else {
+            printf("Usage: %s [-r] [-s<S>] [-n] [<D>]\n", av[0]);
+            printf("  Output <D> random points ");
+            printf("as <D> lines of coordinate pairs.\n");
+            printf("  Default <D> is 10.\n");
+            printf("  -r\t Randomize. Use getpid() as seed.\n");
+            printf("  -s<S>\t Set random seed to <S>.\n");
+            printf("  -n\t Write <D> first before the random points.\n");
+            exit(-1);
+        }
+    }
+    
+    if (PNUM)
+        printf("%d\n", d);
+    for (i=1; i<=d; i++)
+        printf("%4d %4d\n", (int) random()%10000, (int) random()%10000);
+}

flute/src/CMakeLists.txt

@@ -0,0 +1,2 @@
+
+ add_subdirectory(3.1)