cleaned up; ready for release
authorLawrence Cayton <lcayton@tuebingen.mpg.de>
Fri, 19 Nov 2010 19:47:29 +0000 (20:47 +0100)
committerLawrence Cayton <lcayton@tuebingen.mpg.de>
Fri, 19 Nov 2010 19:47:29 +0000 (20:47 +0100)
Makefile
defs.h
driver.cu
kernels.cu
rbc.cu
readme.txt
versions.txt [new file with mode: 0644]

index 61ccaf8..750a214 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -1,7 +1,7 @@
 CC=gcc
 NVCC=nvcc
 CCFLAGS=
-NVCCFLAGS= --ptxas-options=-v -arch=sm_20 
+NVCCFLAGS= --ptxas-options=-v
 #other flags: -deviceemu -arch=sm_20 --compiler-bindir=/usr/bin/gcc-4.3
 LINKFLAGS=-lcuda -lgsl -lgslcblas -lm
 #other linkflags: 
diff --git a/defs.h b/defs.h
index c2a6be0..79d99b5 100644 (file)
--- a/defs.h
+++ b/defs.h
@@ -5,7 +5,7 @@
 
 #define FLOAT_TOL 1e-7
 
-#define K 32 //for k-NN
+#define K 32 //for k-NN.  Do not change!
 #define BLOCK_SIZE 16 //must be a power of 2 (current 
 // implementation of findRange requires a power of 4, in fact)
 
index 72b1495..358020b 100644 (file)
--- a/driver.cu
+++ b/driver.cu
@@ -28,7 +28,7 @@ int main(int argc, char**argv){
   real *data;
   matrix x, q;
   unint *NNs;
-  intMatrix NNsK, NNsKCPU, kNNsRBC;
+  intMatrix NNsK, kNNsRBC;
   unint i;
   struct timeval tvB,tvE;
   cudaError_t cE;
@@ -72,47 +72,24 @@ int main(int argc, char**argv){
   orgData(data, (n+m), d, x, q);
   free(data);
 
-  for(i=0;i<m;i++)
-    NNs[i]=DUMMY_IDX;
-  
   NNsK.r=q.r; NNsK.pr=q.pr; NNsK.pc=NNsK.c=K; NNsK.ld=NNsK.pc;
-  NNsKCPU.r=q.r; NNsKCPU.pr=q.pr; NNsKCPU.pc=NNsKCPU.c=K; NNsKCPU.ld=NNsKCPU.pc;
   kNNsRBC.r=q.r; kNNsRBC.pr=q.pr; kNNsRBC.pc=kNNsRBC.c=K; kNNsRBC.ld=kNNsRBC.pc;
   kNNsRBC.mat = (unint*)calloc(kNNsRBC.pr*kNNsRBC.pc, sizeof(*kNNsRBC.mat));
   NNsK.mat = (unint*)calloc(NNsK.pr*NNsK.pc, sizeof(*NNsK.mat));
-  NNsKCPU.mat = (unint*)calloc(NNsKCPU.pr*NNsKCPU.pc, sizeof(*NNsKCPU.mat));
-
+  
   /* printf("running k-brute force..\n"); */
   /* gettimeofday(&tvB,NULL); */
   /* bruteK(x,q,NNsK); */
   /* gettimeofday(&tvE,NULL); */
   /* printf("\t.. time elapsed = %6.4f \n",timeDiff(tvB,tvE)); */
 
-  /* printf("running regular brute force..\n"); */
-  /* gettimeofday(&tvB,NULL); */
-  /* bruteSearch(x,q,NNs); */
-  /* gettimeofday(&tvE,NULL); */
-  /* printf("\t.. time elapsed = %6.4f \n",timeDiff(tvB,tvE)); */
-
-
-  /* printf("running cpu version..\n"); */
-  /* gettimeofday(&tvB,NULL); */
-  /* bruteKCPU(x,q,NNsKCPU); */
-  /* gettimeofday(&tvE,NULL); */
-  /* printf("\t.. time elapsed = %6.4f \n",timeDiff(tvB,tvE)); */
-
-
   printf("\nrunning rbc..\n");
   gettimeofday(&tvB,NULL);
   buildRBC(x, &rbcS, numReps, s);
   gettimeofday(&tvE,NULL);
   printf("\t.. build time for rbc = %6.4f \n",timeDiff(tvB,tvE));
 
-  /* gettimeofday(&tvB,NULL); */
-  /* queryRBC(q, rbcS, NNs); */
-  /* gettimeofday(&tvE,NULL); */
-  /* printf("\t.. query time for rbc = %6.4f \n",timeDiff(tvB,tvE)); */
-
+  //This finds the 32-NN; if you are only interested in the 1-NN, use queryRBC(..) instead
   gettimeofday(&tvB,NULL);
   kqueryRBC(q, rbcS, kNNsRBC);
   gettimeofday(&tvE,NULL);
@@ -128,26 +105,9 @@ int main(int argc, char**argv){
   
   evalKNNerror(x,q,kNNsRBC);
 
-/* for(i=0;i<q.r;i++){ */
-    /* int j; */
- /*    for(j=0;j<K;j++){ */
- /*      if(NNsK.mat[IDX(i,j,NNsK.ld)] != kNNsRBC.mat[IDX(i,j,kNNsRBC.ld)]) */
- /*    printf("%d %d: (%d: %6.2f %d: %6.2f) ",i,j,NNsK.mat[IDX(i,j,NNsK.ld)],distVec(q,x,i,NNsK.mat[IDX(i,j,NNsK.ld)]),kNNsRBC.mat[IDX(i,j,kNNsRBC.ld)] ,distVec(q,x,i,kNNsRBC.mat[IDX(i,j,kNNsRBC.ld)])); */
- /*    } */
- /* /\*    printf("\n"); *\/ */
- /*  } */
-  /* if(outFile){ */
-  /*   FILE* fp = fopen(outFile, "a"); */
-  /*   fprintf( fp, "%d %d %6.5f %6.5f \n", numReps, s, mean, sqrt(var) ); */
-  /*   fclose(fp); */
-  /* } */
-
-  /* free(ranges); */
-  /* free(cnts); */
   free(NNs);
   free(NNsK.mat);
-  free(NNsKCPU.mat);
+  free(kNNsRBC.mat);
   free(x.mat);
   free(q.mat);
 }
@@ -302,6 +262,15 @@ void evalNNerror(matrix x, matrix q, unint *NNs){
   }
   printf("\tavg rank = %6.4f; std dev = %6.4f \n\n", mean, sqrt(var));
   printf("(range count took %6.4f) \n", timeDiff(tvB, tvE));
+  
+  if(outFile){
+    FILE* fp = fopen(outFile, "a");
+    fprintf( fp, "%d %d %6.5f %6.5f \n", numReps, s, mean, sqrt(var) );
+    fclose(fp);
+  }
+
+  free(ranges);
+  free(cnts);
 }
 
 
@@ -315,7 +284,6 @@ void evalKNNerror(matrix x, matrix q, intMatrix NNs){
   
   unint *ol = (unint*)calloc( q.r, sizeof(*ol) );
   
-
   intMatrix NNsB;
   NNsB.r=q.r; NNsB.pr=q.pr; NNsB.c=NNsB.pc=32; NNsB.ld=NNsB.pc;
   NNsB.mat = (unint*)calloc( NNsB.pr*NNsB.pc, sizeof(*NNsB.mat) );
@@ -344,8 +312,13 @@ void evalKNNerror(matrix x, matrix q, intMatrix NNs){
     var += (((double)ol[i])-mean)*(((double)ol[i])-mean)/((double)m);
   }
   printf("\tavg overlap = %6.4f/%d; std dev = %6.4f \n", mean, K, sqrt(var));
-  
-  
+
+  FILE* fp;
+  if(outFile){
+    fp = fopen(outFile, "a");
+    fprintf( fp, "%d %d %6.5f %6.5f ", numReps, s, mean, sqrt(var) );
+  }
+
   real *ranges = (real*)calloc(q.pr,sizeof(*ranges));
   for(i=0;i<q.r;i++){
     ranges[i] = distVec(q,x,i,NNs.mat[IDX(i, K-1, NNs.ld)]);
@@ -364,10 +337,13 @@ void evalKNNerror(matrix x, matrix q, intMatrix NNs){
     var += (((double)cnts[i])-mean)*(((double)cnts[i])-mean)/((double)m);
   }
   printf("\tavg actual rank of 32nd NN returned by the RBC = %6.4f; std dev = %6.4f \n\n", mean, sqrt(var));
-
-
   printf("(brute k-nn took %6.4f) \n", timeDiff(tvB, tvE));
 
+  if(outFile){
+    fprintf( fp, "%6.5f %6.5f \n", mean, sqrt(var) );
+    fclose(fp);
+  }
+
   free(cnts);
   free(ol);
   free(NNsB.mat);
index 3524b06..bdef21e 100644 (file)
@@ -87,6 +87,9 @@ __global__ void planNNKernel(const matrix Q, const unint *qMap, const matrix X,
 }
 
 
+//This is indentical to the planNNkernel, except that it maintains a list of 32-NNs.  At 
+//each iteration-chunk, the next 16 distances are computed, then sorted, then merged 
+//with the previously computed 32-NNs.
 __global__ void planKNNKernel(const matrix Q, const unint *qMap, const matrix X, const intMatrix xMap, matrix dMins, intMatrix dMinIDs, compPlan cP,  unint qStartPos ){
   unint qB = qStartPos + blockIdx.y * BLOCK_SIZE;  //indexes Q
   unint xB; //X (DB) Block;
@@ -158,9 +161,8 @@ __global__ void planKNNKernel(const matrix Q, const unint *qMap, const matrix X,
 }
 
 
-
+//The basic 1-NN search kernel.
 __global__ void nnKernel(const matrix Q, unint numDone, const matrix X, real *dMins, unint *dMinIDs){
-
   unint qB = blockIdx.y * BLOCK_SIZE + numDone;  //indexes Q
   unint xB; //indexes X;
   unint cB; //colBlock
@@ -220,8 +222,11 @@ __global__ void nnKernel(const matrix Q, unint numDone, const matrix X, real *dM
 }
 
 
+//Computes the 32-NNs for each query in Q.  It is similar to nnKernel above, but maintains a 
+//list of the 32 currently-closest points in the DB, instead of just the single NN.  After each 
+//batch of 16 points is processed, it sorts these 16 points according to the distance from the 
+//query, then merges this list with the other list.
 __global__ void knnKernel(const matrix Q, unint numDone, const matrix X, matrix dMins, intMatrix dMinIDs){
-  
   unint qB = blockIdx.y * BLOCK_SIZE + numDone;  //indexes Q
   unint xB; //indexes X;
   unint cB; //colBlock
@@ -276,7 +281,7 @@ __global__ void knnKernel(const matrix Q, unint numDone, const matrix X, matrix
   
 }
 
-
+//Computes all pairs of distances between Q and X.
 __global__ void dist1Kernel(const matrix Q, unint qStart, const matrix X, unint xStart, matrix D){
   unint c, i, j;
 
@@ -312,9 +317,9 @@ __global__ void dist1Kernel(const matrix Q, unint qStart, const matrix X, unint
 }
 
 
-
+//This function is used by the rbc building routine.  It find an appropriate range 
+//such that roughly cntWant points fall within this range.  D is a matrix of distances.
 __global__ void findRangeKernel(const matrix D, unint numDone, real *ranges, unint cntWant){
-  
   unint row = blockIdx.y*(BLOCK_SIZE/4)+threadIdx.y + numDone;
   unint ro = threadIdx.y;
   unint co = threadIdx.x;
@@ -466,53 +471,11 @@ __global__ void rangeCountKernel(const matrix Q, unint numDone, const matrix X,
 }
 
 
-__device__ void sort16off(real x[][48], unint xi[][48]){
-  int i = threadIdx.x;
-  int j = threadIdx.y;
-
-  if(i%2==0)
-    mmGateI( x[j]+K+i, x[j]+K+i+1, xi[j]+K+i, xi[j]+K+i+1 );
-  __syncthreads();
-
-  if(i%4<2)
-    mmGateI( x[j]+K+i, x[j]+K+i+2, xi[j]+K+i, xi[j]+K+i+2 );
-  __syncthreads();
-
-  if(i%4==1)
-    mmGateI( x[j]+K+i, x[j]+K+i+1, xi[j]+K+i, xi[j]+K+i+1 );
-  __syncthreads();
-  
-  if(i%8<4)
-    mmGateI( x[j]+K+i, x[j]+K+i+4, xi[j]+K+i, xi[j]+K+i+4 );
-  __syncthreads();
-  
-  if(i%8==2 || i%8==3)
-    mmGateI( x[j]+K+i, x[j]+K+i+2, xi[j]+K+i, xi[j]+K+i+2 );
-  __syncthreads();
-
-  if( i%2 && i%8 != 7 ) 
-    mmGateI( x[j]+K+i, x[j]+K+i+1, xi[j]+K+i, xi[j]+K+i+1 );
-  __syncthreads();
-  
-  //0-7; 8-15 now sorted.  merge time.
-  if( i<8)
-    mmGateI( x[j]+K+i, x[j]+K+i+8, xi[j]+K+i, xi[j]+K+i+8 );
-  __syncthreads();
-  
-  if( i>3 && i<8 )
-    mmGateI( x[j]+K+i, x[j]+K+i+4, xi[j]+K+i, xi[j]+K+i+4 );
-  __syncthreads();
-  
-  int os = (i/2)*4+2 + i%2;
-  if(i<6)
-    mmGateI( x[j]+K+os, x[j]+K+os+2, xi[j]+K+os, xi[j]+K+os+2 );
-  __syncthreads();
-  
-  if( i%2 && i<15)
-    mmGateI( x[j]+K+i, x[j]+K+i+1, xi[j]+K+i, xi[j]+K+i+1 );
-}
-
+//**************************************************************************
+// The following functions are an implementation of Batcher's sorting network.  
+// All computations take place in (on-chip) shared memory.
 
+// The function name is descriptive; it sorts each row of x, whose indices are xi.
 __device__ void sort16(real x[][16], unint xi[][16]){
   int i = threadIdx.x;
   int j = threadIdx.y;
@@ -561,6 +524,9 @@ __device__ void sort16(real x[][16], unint xi[][16]){
 }
 
 
+// This function takes an array of lists, each of length 48. It is assumed
+// that the first 32 numbers are sorted, and the last 16 numbers.  The 
+// routine then merges these lists into one sorted list of length 48.
 __device__ void merge32x16(real x[][48], unint xi[][48]){
   int i = threadIdx.x;
   int j = threadIdx.y;
@@ -598,7 +564,57 @@ __device__ void merge32x16(real x[][48], unint xi[][48]){
 
 }
 
+//This is the same as sort16, but takes as input lists of length 48
+//and sorts the last 16 entries.  This cleans up some of the NN code, 
+//though it is inelegant.
+__device__ void sort16off(real x[][48], unint xi[][48]){
+  int i = threadIdx.x;
+  int j = threadIdx.y;
+
+  if(i%2==0)
+    mmGateI( x[j]+K+i, x[j]+K+i+1, xi[j]+K+i, xi[j]+K+i+1 );
+  __syncthreads();
+
+  if(i%4<2)
+    mmGateI( x[j]+K+i, x[j]+K+i+2, xi[j]+K+i, xi[j]+K+i+2 );
+  __syncthreads();
+
+  if(i%4==1)
+    mmGateI( x[j]+K+i, x[j]+K+i+1, xi[j]+K+i, xi[j]+K+i+1 );
+  __syncthreads();
+  
+  if(i%8<4)
+    mmGateI( x[j]+K+i, x[j]+K+i+4, xi[j]+K+i, xi[j]+K+i+4 );
+  __syncthreads();
+  
+  if(i%8==2 || i%8==3)
+    mmGateI( x[j]+K+i, x[j]+K+i+2, xi[j]+K+i, xi[j]+K+i+2 );
+  __syncthreads();
+
+  if( i%2 && i%8 != 7 ) 
+    mmGateI( x[j]+K+i, x[j]+K+i+1, xi[j]+K+i, xi[j]+K+i+1 );
+  __syncthreads();
+  
+  //0-7; 8-15 now sorted.  merge time.
+  if( i<8)
+    mmGateI( x[j]+K+i, x[j]+K+i+8, xi[j]+K+i, xi[j]+K+i+8 );
+  __syncthreads();
+  
+  if( i>3 && i<8 )
+    mmGateI( x[j]+K+i, x[j]+K+i+4, xi[j]+K+i, xi[j]+K+i+4 );
+  __syncthreads();
+  
+  int os = (i/2)*4+2 + i%2;
+  if(i<6)
+    mmGateI( x[j]+K+os, x[j]+K+os+2, xi[j]+K+os, xi[j]+K+os+2 );
+  __syncthreads();
+  
+  if( i%2 && i<15)
+    mmGateI( x[j]+K+i, x[j]+K+i+1, xi[j]+K+i, xi[j]+K+i+1 );
+}
 
+//min-max gate: it sets the minimum of x and y into x, the maximum into y, and 
+//exchanges the indices (xi and yi) accordingly.
 __device__ void mmGateI(real *x, real *y, unint *xi, unint *yi){
   int ti = MINi( *x, *y, *xi, *yi );
   *yi = MAXi( *x, *y, *xi, *yi );
diff --git a/rbc.cu b/rbc.cu
index 4cfef0f..7629e03 100644 (file)
--- a/rbc.cu
+++ b/rbc.cu
@@ -66,7 +66,8 @@ void queryRBC(const matrix q, const rbcStruct rbcS, unint *NNs){
   free(groupCountQ);
 }
 
-
+//This function is very similar to queryRBC, with a couple of basic changes to handle
+//k-nn.  
 void kqueryRBC(const matrix q, const rbcStruct rbcS, intMatrix NNs){
   unint m = q.r;
   unint numReps = rbcS.dr.r;
@@ -99,6 +100,9 @@ void kqueryRBC(const matrix q, const rbcStruct rbcS, intMatrix NNs){
   // Setup the computation matrix.  Currently, the computation matrix is 
   // just the identity matrix: each query assigned to a particular 
   // representative is compared only to that representative's points.  
+
+  // NOTE: currently, idIntersection is the *only* computation matrix 
+  // that will work properly with k-nn search (this is not true for 1-nn above).
   idIntersection(cM);
 
   initCompPlan(&dcP, cM, groupCountQ, rbcS.groupCount, numReps);
@@ -119,7 +123,6 @@ void kqueryRBC(const matrix q, const rbcStruct rbcS, intMatrix NNs){
 
 
 void buildRBC(const matrix x, rbcStruct *rbcS, unint numReps, unint s){
-  //const matrix dr, intMatrix xmap, unint *counts, unint s){
   unint n = x.pr;
   intMatrix xmap;
 
@@ -335,7 +338,6 @@ void computeKNNs(matrix dx, intMatrix dxMap, matrix dq, unint *dqMap, compPlan d
   planKNNWrap(dq, dqMap, dx, dxMap, dMins, dMinIDs, dcP, compLength);
   cudaMemcpy( NNs.mat, dMinIDs.mat, dq.r*K*sizeof(*NNs.mat), cudaMemcpyDeviceToHost);
 
-  
   cudaFree(dMins.mat);
   cudaFree(dMinIDs.mat);
 }
index 0d0244f..7209623 100644 (file)
@@ -1,4 +1,4 @@
-***Random Ball Cover (RBC)***
+***Random Ball Cover (RBC) v0.2***
 Lawrence Cayton
 lcayton@tuebingen.mpg.de
 
@@ -70,13 +70,13 @@ may wish to call the readData function on separate files.  There is
 also a readDataText function in the driver.cu for your convenience.  
 
 The core of the implementation is in rbc.cu and in the kernel files.
-There is a buildRBC function and a queryRBC function which should
-suffice for basic use of the data structure.
+There is a buildRBC function, a queryRBC function, and a kqueryRBC
+function, which together should suffice for basic use of the data
+structure.   
 
-Currently, the kernel functions are not seriously optimized.  Indeed,
-the results appearing in the ADMS paper came from a slightly more
-optimized version than this one.  A tuned version will be available at
-some point.   
+Currently, the kernel functions are reasonably optimized, but can be
+improved.  Indeed, the results appearing in the ADMS paper came from a
+slightly more optimized version than this one. 
 
 
 ---------------------------------------------------------------------
@@ -89,6 +89,14 @@ MISC NOTES ON THE CODE
   but more complex metrics will require some aditional work.  The L_2
   metric (standard Euclidean distance) is already defined in defs.h.  
 
+* The k-NN code is currently hard-coded for K=32.  The k-nn code
+  contains a manual implementation of a sorting network, which is why
+  this is hard-coded.  This design allows all sorting to take place
+  in on-chip (shared) memory, and is highly efficient.  Note that
+  the NNs are returned in sorted order, so that if one wants only,
+  say, 5 NNs, one can simply ignore the last 27 returned indices.  For
+  k>32, contact the author.
+
 * The code requires that the entire DB and query set fit into the 
   device memory.  
 
diff --git a/versions.txt b/versions.txt
new file mode 100644 (file)
index 0000000..ba27e64
--- /dev/null
@@ -0,0 +1,5 @@
+0.2 nov 2010
+    Added support for k-NN queries, for k <= 32.  
+
+0.1 aug 2010 
+    Initial release.