Cleaned up, debugged. Ready for 1st release
authorLawrence Cayton <lcayton@tuebingen.mpg.de>
Wed, 18 Aug 2010 17:02:48 +0000 (19:02 +0200)
committerLawrence Cayton <lcayton@tuebingen.mpg.de>
Wed, 18 Aug 2010 17:02:48 +0000 (19:02 +0200)
21 files changed:
Makefile
brute.cu
brute.h
defs.h
deprecated.cu
driver.cu
kernelWrap.cu
kernelWrap.h
kernels.cu
kernels.h
rbc.cu
rbc.h
readme.txt
sKernel.cu
sKernel.h
sKernelWrap.cu [new file with mode: 0644]
sKernelWrap.h [new file with mode: 0644]
utils.cu
utils.h
utilsGPU.cu
utilsGPU.h

index 30a7799..6c748aa 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -1,12 +1,12 @@
 CC=gcc
 NVCC=nvcc
 CCFLAGS=
-NVCCFLAGS= --ptxas-options=-v
+NVCCFLAGS= --ptxas-options=-v -arch=sm_20
 #other flags: -deviceemu -arch=sm_20 --compiler-bindir=/usr/bin/gcc-4.3
 LINKFLAGS=-lcuda 
 #other linkflags: 
 SOURCES=
-CUSOURCES=driver.cu utils.cu utilsGPU.cu rbc.cu kernels.cu brute.cu kernelWrap.cu sKernel.cu
+CUSOURCES= driver.cu utils.cu utilsGPU.cu rbc.cu brute.cu kernels.cu kernelWrap.cu sKernel.cu sKernelWrap.cu
 OBJECTS=$(SOURCES:.c=.o)
 CUOBJECTS=$(CUSOURCES:.cu=.o)
 EXECUTABLE=testRBC
index ca55c54..d2377c0 100644 (file)
--- a/brute.cu
+++ b/brute.cu
@@ -1,3 +1,7 @@
+/* This file is part of the Random Ball Cover (RBC) library.
+ * (C) Copyright 2010, Lawrence Cayton [lcayton@tuebingen.mpg.de]
+ */
+
 #ifndef BRUTE_CU
 #define BRUTE_CU
 
 #include<stdio.h>
 #include<cuda.h>
 
-void bruteRangeCount(matrix x, matrix q, real *ranges, int *cnts){
+void bruteRangeCount(matrix x, matrix q, real *ranges, unint *cnts){
   matrix dx, dq;
   real *dranges;
-  int *dcnts;
+  unint *dcnts;
   
   copyAndMove(&dx, &x);
   copyAndMove(&dq, &q);
 
-  cudaMalloc( (void**)&dranges, q.pr*sizeof(*dranges) );
+  checkErr( cudaMalloc( (void**)&dranges, q.pr*sizeof(*dranges) ) );
   cudaMemcpy( dranges, ranges, q.r*sizeof(*dranges), cudaMemcpyHostToDevice );
 
-  cudaMalloc( (void**)&dcnts, q.pr*sizeof(*dcnts) );
+  checkErr( cudaMalloc( (void**)&dcnts, q.pr*sizeof(*dcnts) ) );
   
   rangeCountWrap(dq, dx, dranges, dcnts);
   
@@ -35,19 +39,19 @@ void bruteRangeCount(matrix x, matrix q, real *ranges, int *cnts){
 }
 
 
-void bruteSearch(matrix x, matrix q, int *NNs){
+void bruteSearch(matrix x, matrix q, unint *NNs){
   real *dMins;
-  int *dMinIDs;
+  unint *dMinIDs;
   matrix dx, dq;
 
   
   dx.r=x.r; dx.pr=x.pr; dx.c=x.c; dx.pc=x.pc; dx.ld=x.ld;
   dq.r=q.r; dq.pr=q.pr; dq.c=q.c; dq.pc=q.pc; dq.ld=q.ld;
 
-  cudaMalloc((void**)&dMins, q.pr*sizeof(*dMins));
-  cudaMalloc((void**)&dMinIDs, q.pr*sizeof(*dMinIDs));
-  cudaMalloc((void**)&dx.mat, dx.pr*dx.pc*sizeof(*dx.mat));
-  cudaMalloc((void**)&dq.mat, dq.pr*dq.pc*sizeof(*dq.mat));
+  checkErr( cudaMalloc((void**)&dMins, q.pr*sizeof(*dMins)) );
+  checkErr( cudaMalloc((void**)&dMinIDs, q.pr*sizeof(*dMinIDs)) );
+  checkErr( cudaMalloc((void**)&dx.mat, dx.pr*dx.pc*sizeof(*dx.mat)) );
+  checkErr( cudaMalloc((void**)&dq.mat, dq.pr*dq.pc*sizeof(*dq.mat)) );
 
   cudaMemcpy(dx.mat,x.mat,x.pr*x.pc*sizeof(*dx.mat),cudaMemcpyHostToDevice);
   cudaMemcpy(dq.mat,q.mat,q.pr*q.pc*sizeof(*dq.mat),cudaMemcpyHostToDevice);
@@ -63,11 +67,11 @@ void bruteSearch(matrix x, matrix q, int *NNs){
 }
 
 
-void bruteCPU(matrix X, matrix Q, int *NNs){
+void bruteCPU(matrix X, matrix Q, unint *NNs){
   real *dtoNNs; 
   real temp;
 
-  int i, j;
+  unint i, j;
 
   dtoNNs = (real*)calloc(Q.r,sizeof(*dtoNNs));
   
@@ -75,7 +79,7 @@ void bruteCPU(matrix X, matrix Q, int *NNs){
     dtoNNs[i] = MAX_REAL;
     NNs[i] = 0;
     for(j=0; j<X.r; j++ ){
-      temp = distL1( Q, X, i, j );
+      temp = distVec( Q, X, i, j );
       if( temp < dtoNNs[i]){
        NNs[i] = j;
        dtoNNs[i] = temp;
diff --git a/brute.h b/brute.h
index 534cded..cb25e9f 100644 (file)
--- a/brute.h
+++ b/brute.h
@@ -1,10 +1,14 @@
+/* This file is part of the Random Ball Cover (RBC) library.
+ * (C) Copyright 2010, Lawrence Cayton [lcayton@tuebingen.mpg.de]
+ */
+
 #ifndef BRUTE_H
 #define BRUTE_H
 
 #include "defs.h"
 
-void bruteRangeCount(matrix,matrix,real*,int*);
-void bruteSearch(matrix,matrix,int*);
-void bruteCPU(matrix,matrix,int*);
+void bruteRangeCount(matrix,matrix,real*,unint*);
+void bruteSearch(matrix,matrix,unint*);
+void bruteCPU(matrix,matrix,unint*);
 
 #endif
diff --git a/defs.h b/defs.h
index ba308f0..2cbd4d3 100644 (file)
--- a/defs.h
+++ b/defs.h
@@ -4,11 +4,19 @@
 #include<float.h>
 
 #define FLOAT_TOL 1e-7
-#define BLOCK_SIZE 16 //must be a power of 2
+#define BLOCK_SIZE 16 //must be a power of 2 (current 
+// implementation of findRange requires a power of 4, in fact)
 
 #define MAX_BS 65535 //max block size (specified by CUDA)
 #define SCAN_WIDTH 1024
 
+#define MEM_USED_IN_SCAN(n) ( 2*( (n) + SCAN_WIDTH-1 )/SCAN_WIDTH*sizeof(unint))
+
+//The distance measure that is used.  This macro returns the 
+//distance for a single coordinate.
+#define DIST(i,j) ( fabs((i)-(j)) )  // L_1
+//#define DIST(i,j) ( ( (i)-(j) )*( (i)-(j) ) )  // L_2
+
 // Format that the data is manipulated in:
 typedef float real;
 #define MAX_REAL FLT_MAX
@@ -22,7 +30,7 @@ typedef float real;
 //Percentage of device mem to use
 #define MEM_USABLE .95
 
-#define DUMMY_IDX INT_MAX
+#define DUMMY_IDX UINT_MAX
 
 //Row major indexing
 #define IDX(i,j,ld) (((i)*(ld))+(j))
@@ -37,53 +45,53 @@ typedef float real;
 
 #define MIN(i,j) ((i) < (j) ? (i) : (j))
 
+typedef unsigned int unint;
+
+
 typedef struct {
   real *mat;
-  int r; //rows
-  int c; //cols
-  int pr; //padded rows
-  int pc; //padded cols
-  int ld; //the leading dimension (in this code, this is the same as pc)
+  unint r; //rows
+  unint c; //cols
+  unint pr; //padded rows
+  unint pc; //padded cols
+  unint ld; //the leading dimension (in this code, this is the same as pc)
 } matrix;
 
 
 typedef struct {
   char *mat;
-  int r;
-  int c;
-  int pr;
-  int pc;
-  int ld;
+  unint r;
+  unint c;
+  unint pr;
+  unint pc;
+  unint ld;
 } charMatrix;
 
-typedef struct {
-  int *mat;
-  int r;
-  int c;
-  int pr;
-  int pc;
-  int ld;
-} intMatrix;
-
 
 typedef struct {
-  int numComputeSegs;
-  int normSegSize;//The number of points handled in one computation,
-                     //though there will always be one leftover segment
-                     //with (possibly) a different number of points.
-  int lastSegSize;//.. and this is it.
-} memPlan;
+  unint *mat;
+  unint r;
+  unint c;
+  unint pr;
+  unint pc;
+  unint ld;
+} intMatrix;
 
 
 typedef struct{
-  int *numGroups; //The number of groups of DB points to be examined.
-  int sNumGroups; //size of the above array.
-  int *groupCountX; //The number of elements in each group.
-  int sGroupCountX;
-  int *groupOff;
-  int sGroupOff;
-  int *qToGroup; //map from query to group #.
-  int sQToGroup;
-  int ld; //the width of memPos and groupCount (= max over numGroups)
+  unint *numGroups; //The number of groups of DB points to be examined.
+  unint *groupCountX; //The number of elements in each DB group.
+  unint *qToQGroup; //map from query to query group #.
+  unint *qGroupToXGroup; //map from query group to DB gruop
+  unint ld; //the width of memPos and groupCount (= max over numGroups)
 } compPlan;
+
+
+typedef struct {
+  matrix dx;
+  intMatrix dxMap;
+  matrix dr;
+  unint *groupCount;
+} rbcStruct;
+
 #endif
index 3dfc315..f2c7d8e 100644 (file)
@@ -1,3 +1,10 @@
+/* This file contains functions that were in use at one point but 
+   are not currently used.  As far as I know, everything here is
+   debugged, and so can be plugged in reasonably safely.
+*/
+
+
+
 /* This is the more complex version of computeReps.  It can be used even if X doesn't
 fit into the device memory.  It does not work in emulation mode, because it checks to see
 how much mem is available on the device.  Thus for debugging purposes we currently 
@@ -314,3 +321,146 @@ void brute2Step(matrix x, matrix q, intMatrix NNs){
   cudaFree(dMins.mat);
 
 }
+
+
+memPlan createMemPlan(int,int);
+
+memPlan createMemPlan(unint nPts, unint memPerPt){
+  memPlan mp;
+  unsigned int memFree, memTot;
+  unint ptsAtOnce;
+
+  cuMemGetInfo(&memFree, &memTot);
+  memFree = (unint)(((float)memFree)*MEM_USABLE);
+  printf("memfree = %d \n",memFree);
+  ptsAtOnce = DPAD(memFree/memPerPt); //max number of pts that can be processed at once
+  printf("ptsAtOnce = %d \n",ptsAtOnce);
+  mp.numComputeSegs = nPts/ptsAtOnce + ((nPts%ptsAtOnce==0) ? 0 : 1);
+  mp.normSegSize=PAD(nPts/mp.numComputeSegs); 
+  mp.lastSegSize=PAD(nPts) - mp.normSegSize*(mp.numComputeSegs-1);
+  //Note that lastSegSize is automatically padded if nPts is.
+  return mp;
+}
+
+typedef struct {
+  unint numComputeSegs;
+  unint normSegSize;//The number of points handled in one computation,
+                     //though there will always be one leftover segment
+                     //with (possibly) a different number of points.
+  unint lastSegSize;//.. and this is it.
+} memPlan;
+
+
+void blockIntersection(charMatrix,matrix,real*,real*);
+
+void blockIntersection(charMatrix cM, matrix dr, real *radiiX, real *radiiQ){
+  matrix dD;
+  real *dradiiX, *dradiiQ;
+  unint pnR = dr.pr;
+  charMatrix dcM;
+  
+  dD.r=dD.c=dr.r; dD.pr=dD.pc=dD.ld=dr.pr;
+  dcM.r=cM.r; dcM.c=cM.c; dcM.pr=cM.pr; dcM.pc=cM.pc; dcM.ld=cM.ld;
+  
+  checkErr( cudaMalloc((void**)&dD.mat, pnR*pnR*sizeof(*dD.mat)) );
+  checkErr( cudaMalloc((void**)&dradiiX, pnR*sizeof(*dradiiX)) );
+  checkErr( cudaMalloc((void**)&dradiiQ, pnR*sizeof(*dradiiQ)) );
+  checkErr( cudaMalloc((void**)&dcM.mat, dcM.pr*dcM.pc*sizeof(*dcM.mat)) );
+  
+  // Copying over the radii. Note that everything after the first dr.r places 
+  // on the device variables is undefined.
+  cudaMemcpy(dradiiX,radiiX,dr.r*sizeof(*dradiiX),cudaMemcpyHostToDevice);
+  cudaMemcpy(dradiiQ,radiiQ,dr.r*sizeof(*dradiiQ),cudaMemcpyHostToDevice);
+  
+  dist1Wrap(dr, dr, dD);
+  pruneWrap(dcM, dD, dradiiX, dradiiQ);
+
+  cudaMemcpy(cM.mat,dcM.mat,pnR*pnR*sizeof(*dcM.mat),cudaMemcpyDeviceToHost);
+  
+  cudaFree(dcM.mat);
+  cudaFree(dradiiQ);
+  cudaFree(dradiiX);
+  cudaFree(dD.mat);
+}
+
+
+
+void groupPoints(matrix,unint*,unint*,unint);
+
+// This function sorts points by their repID.  It makes two passes through the 
+// matrix x; one to count the bucket sizes, the next to place points in the 
+// correct bucket.  Note that this function allocates a temporary
+// matrix the size of x, then copies the results over to x at the end.  The 
+// sort could be done in place, eg by doing numReps passes through x instead of 2.
+void groupPoints(matrix x, unint *xID, unint *repIDs, unint numReps){
+  matrix y;
+  unint n=x.r;
+  unint d=x.c;
+  unint i;
+  unint *gS; //groupSize
+  unint *yID;
+
+  yID = (unint*)calloc(n,sizeof(*yID));
+  y.mat = (real*)calloc(n*d,sizeof(*y.mat));
+  gS = (unint*)calloc(numReps+1,sizeof(*gS));
+
+  y.r=n; y.pr=n; y.c=d; y.pc=d; y.ld=d;
+
+  for(i=0;i<n;i++)
+    gS[repIDs[i]+1]++;
+  for(i=1;i<numReps;i++)
+    gS[i]=gS[i-1]+gS[i];
+  
+  for(i=0;i<n;i++){
+    copyVector(&y.mat[IDX(gS[repIDs[i]],0,y.ld)], &x.mat[IDX(i,0,x.ld)],d);
+    yID[gS[repIDs[i]]]=xID[i];
+    gS[repIDs[i]]++;
+  }
+  
+  for(i=0;i<n;i++){
+    copyVector(&x.mat[IDX(i,0,x.ld)],&y.mat[IDX(i,0,y.ld)],d);
+    xID[i]=yID[i];
+  }
+  
+  free(yID);
+  free(gS);
+  free(y.mat);
+}
+
+__global__ void pruneKernel(const matrix,const real*,const real*,charMatrix);
+
+
+__global__ void pruneKernel(const matrix D, const real *radiiX, const real *radiiQ, charMatrix cM){
+  unint offX = threadIdx.x;
+  unint offQ = threadIdx.y;
+
+  unint blockX = blockIdx.x * BLOCK_SIZE;
+  unint blockQ = blockIdx.y * BLOCK_SIZE;
+  
+  __shared__ real sD[BLOCK_SIZE][BLOCK_SIZE];
+  __shared__ real sRQ[BLOCK_SIZE];
+  __shared__ real sRX[BLOCK_SIZE];
+
+  sD[offQ][offX]=D.mat[IDX(blockQ+offQ,blockX+offX,D.ld)];
+  
+  if(offQ==0)
+    sRX[offX]=radiiX[blockX+offX];
+  if(offX==0)
+    sRQ[offQ]=radiiQ[blockQ+offQ];
+  
+  __syncthreads();
+  
+  if(blockQ+offQ < D.r && blockX+offX < D.c){
+    cM.mat[IDX(blockQ+offQ,blockX+offX,cM.ld)] = (sD[offQ][offX]-sRX[offX]-2*sRQ[offQ] <= 0) ? 1 : 0;
+  }
+}
+
+void pruneWrap(charMatrix dcM, matrix dD, real *dradiiX, real *dradiiQ){
+  dim3 block(BLOCK_SIZE,BLOCK_SIZE);
+  dim3 grid(dD.pr/BLOCK_SIZE,dD.pc/BLOCK_SIZE);
+  
+  pruneKernel<<<grid,block>>>(dD,dradiiX,dradiiQ,dcM);
+  cudaThreadSynchronize();
+}
+
+void pruneWrap(charMatrix,matrix,real*,real*);
index 6ec580f..7bd2711 100644 (file)
--- a/driver.cu
+++ b/driver.cu
@@ -1,3 +1,7 @@
+/* This file is part of the Random Ball Cover (RBC) library.
+ * (C) Copyright 2010, Lawrence Cayton [lcayton@tuebingen.mpg.de]
+ */
+
 #include<stdio.h>
 #include<stdlib.h>
 #include<cuda.h>
 #include "sKernel.h"
 
 void parseInput(int,char**);
-void readData(char*,int,int,real*);
-void orgData(real*,int,int,matrix,matrix);
+void readData(char*,unint,unint,real*);
+void readDataText(char*,unint,unint,real*);
+void orgData(real*,unint,unint,matrix,matrix);
 
 
 char *dataFile, *outFile;
-int n=0, m=0, d=0, numReps=0, s=0;
-int deviceNum=0;
+unint n=0, m=0, d=0, numReps=0, s=0;
+unint deviceNum=0;
 int main(int argc, char**argv){
   real *data;
   matrix x, q;
-  int *NNs, *NNsBrute;
-  int i;
+  unint *NNs, *NNsBrute;
+  unint i;
   struct timeval tvB,tvE;
   cudaError_t cE;
+  rbcStruct rbcS;
 
   printf("*****************\n");
   printf("RANDOM BALL COVER\n");
@@ -39,7 +45,6 @@ int main(int argc, char**argv){
     exit(1);
   }
   
-  
   unsigned int memFree, memTot;
   CUcontext pctx;
   unsigned int flags=0;
@@ -57,23 +62,15 @@ int main(int argc, char**argv){
   x.r = n; x.c = d; x.pr = PAD(n); x.pc = PAD(d); x.ld = x.pc;
   q.r = m; q.c = d; q.pr = PAD(m); q.pc = PAD(d); q.ld = q.pc;
 
-  NNs = (int*)calloc( m, sizeof(*NNs) );
+  NNs = (unint*)calloc( m, sizeof(*NNs) );
   for(i=0; i<m; i++)
-    NNs[i]=-1;
-  NNsBrute = (int*)calloc( m, sizeof(*NNsBrute) );
+    NNs[i]=DUMMY_IDX;
+  NNsBrute = (unint*)calloc( m, sizeof(*NNsBrute) );
 
   readData(dataFile, (n+m), d, data);
   orgData(data, (n+m), d, x, q);
   free(data);
 
-
-  
-  /* printf("db:\n"); */
-  /* printMat(x); */
-  /* printf("\nqueries: \n"); */
-  /* printMat(q); */
-  /* printf("\n\n"); */
-  
   for(i=0;i<m;i++)
     NNs[i]=NNsBrute[i]=DUMMY_IDX;
   
@@ -83,12 +80,18 @@ int main(int argc, char**argv){
   /* gettimeofday(&tvE,NULL); */
   /* printf("\t.. time elapsed = %6.4f \n",timeDiff(tvB,tvE)); */
   
-  
   printf("\nrunning rbc..\n");
   gettimeofday(&tvB,NULL);
-  rbc(x,q,numReps,s,NNs); 
+  buildRBC(x, &rbcS, numReps, s);
   gettimeofday(&tvE,NULL);
-  printf("\t.. total time elapsed for rbc = %6.4f \n",timeDiff(tvB,tvE));
+  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));
+
+  destroyRBC(&rbcS);
   printf("finished \n");
   
   cE = cudaGetLastError();
@@ -98,11 +101,12 @@ int main(int argc, char**argv){
 
   printf("\nComputing error rates (this might take a while)\n");
   real *ranges = (real*)calloc(q.pr,sizeof(*ranges));
-  for(i=0;i<q.r;i++)
-    ranges[i] = distL1(q,x,i,NNs[i]) - 10e-6;
-  
+  for(i=0;i<q.r;i++){
+    if(NNs[i]>n) printf("error");
+    ranges[i] = distVec(q,x,i,NNs[i]) - 10e-6;
+  }
 
-  int *cnts = (int*)calloc(q.pr,sizeof(*cnts));
+  unint *cnts = (unint*)calloc(q.pr,sizeof(*cnts));
   gettimeofday(&tvB,NULL);
   bruteRangeCount(x,q,ranges,cnts);
   gettimeofday(&tvE,NULL);
@@ -119,6 +123,7 @@ int main(int argc, char**argv){
   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) );
@@ -138,6 +143,15 @@ void parseInput(int argc, char **argv){
   int i=1;
   if(argc <= 1){
     printf("\nusage: \n  testRBC -f datafile (bin) -n numPts (DB) -m numQueries -d dim -r numReps -s numPtsPerRep [-o outFile] [-g GPU num]\n\n");
+    printf("\tdatafile     = binary file containing the data\n");
+    printf("\tnumPts       = size of database\n");
+    printf("\tnumQueries   = number of queries\n");
+    printf("\tdim          = dimensionailty\n");
+    printf("\tnumReps      = number of representatives\n");
+    printf("\tnumPtsPerRep = number of points assigned to each representative\n");
+    printf("\toutFile      = output file (optional); stored in text format\n");
+    printf("\tGPU num      = ID # of the GPU to use (optional) for multi-GPU machines\n");
+    printf("\n\n");
     exit(0);
   }
   
@@ -177,9 +191,9 @@ void parseInput(int argc, char **argv){
 }
 
 
-void readData(char *dataFile, int rows, int cols, real *data){
+void readData(char *dataFile, unint rows, unint cols, real *data){
   FILE *fp;
-  int numRead;
+  unint numRead;
 
   fp = fopen(dataFile,"r");
   if(fp==NULL){
@@ -196,15 +210,37 @@ void readData(char *dataFile, int rows, int cols, real *data){
 }
 
 
+void readDataText(char *dataFile, unint rows, unint cols, real *data){
+  FILE *fp;
+  real t;
+
+  fp = fopen(dataFile,"r");
+  if(fp==NULL){
+    fprintf(stderr,"error opening file.. exiting\n");
+    exit(1);
+  }
+    
+  for(int i=0; i<rows; i++){
+    for(int j=0; j<cols; j++){
+      if(fscanf(fp,"%f ", &t)==EOF){
+       fprintf(stderr,"error reading file.. exiting \n");
+       exit(1);
+      }
+      data[IDX( i, j, cols )]=(real)t;
+    }
+  }
+  fclose(fp);
+}
+
 //This function splits the data into two matrices, x and q, of 
 //their specified dimensions.  The data is split randomly.
 //It is assumed that the number of rows of data (the parameter n)
 //is at least as large as x.r+q.r
-void orgData(real *data, int n, int d, matrix x, matrix q){
+void orgData(real *data, unint n, unint d, matrix x, matrix q){
    
-  int i,fi,j;
-  int *p;
-  p = (int*)calloc(n,sizeof(*p));
+  unint i,fi,j;
+  unint *p;
+  p = (unint*)calloc(n,sizeof(*p));
   
   randPerm(n,p);
 
index 292d941..e5a69ce 100644 (file)
@@ -1,3 +1,7 @@
+/* This file is part of the Random Ball Cover (RBC) library.
+ * (C) Copyright 2010, Lawrence Cayton [lcayton@tuebingen.mpg.de]
+ */
+
 #ifndef KERNELWRAP_CU
 #define KERNELWRAP_CU
 
 #include "kernels.h"
 #include "defs.h"
 
-void dist1Wrap(matrix dq, matrix dx, matrix dD){
+void dist1Wrap(const matrix dq, const matrix dx, matrix dD){
   dim3 block(BLOCK_SIZE,BLOCK_SIZE);
   dim3 grid;
   
-  int todoX, todoY, numDoneX, numDoneY;
+  unint todoX, todoY, numDoneX, numDoneY;
 
   numDoneX = 0;
   while ( numDoneX < dx.pr ){
-    todoX = min( dx.pr - numDoneX, MAX_BS*BLOCK_SIZE );
+    todoX = MIN( dx.pr - numDoneX, MAX_BS*BLOCK_SIZE );
     grid.x = todoX/BLOCK_SIZE;
     numDoneY = 0;
     while( numDoneY < dq.pr ){
-      todoY = min( dq.pr - numDoneY, MAX_BS*BLOCK_SIZE );
+      todoY = MIN( dq.pr - numDoneY, MAX_BS*BLOCK_SIZE );
       grid.y = todoY/BLOCK_SIZE;
       dist1Kernel<<<grid,block>>>(dq, numDoneY, dx, numDoneX, dD);
       numDoneY += todoY;
@@ -30,29 +34,35 @@ void dist1Wrap(matrix dq, matrix dx, matrix dD){
 }
 
 
-void findRangeWrap(matrix dD, real *dranges, int cntWant){
+void findRangeWrap(const matrix dD, real *dranges, unint cntWant){
   dim3 block(4*BLOCK_SIZE,BLOCK_SIZE/4);
   dim3 grid(1,4*(dD.pr/BLOCK_SIZE));
-
-  findRangeKernel<<<grid,block>>>(dD,dranges,cntWant);
-
+  unint numDone, todo;
   
+  numDone=0;
+  while( numDone < dD.pr ){
+    todo = MIN ( dD.pr - numDone, MAX_BS*BLOCK_SIZE/4 );
+    grid.y = 4*(todo/BLOCK_SIZE);
+    findRangeKernel<<<grid,block>>>(dD, numDone, dranges, cntWant);
+    numDone += todo;
+  }
   cudaThreadSynchronize();
 }
 
-void rangeSearchWrap(matrix dD, real *dranges, charMatrix dir){
+
+void rangeSearchWrap(const matrix dD, const real *dranges, charMatrix dir){
   dim3 block(BLOCK_SIZE,BLOCK_SIZE);
-  dim3 grid(dD.pc/BLOCK_SIZE,dD.pr/BLOCK_SIZE);
+  dim3 grid;
 
-  int todoX, todoY, numDoneX, numDoneY;
+  unint todoX, todoY, numDoneX, numDoneY;
   
   numDoneX = 0;
   while ( numDoneX < dD.pc ){
-    todoX = min( dD.pc - numDoneX, MAX_BS*BLOCK_SIZE );
+    todoX = MIN( dD.pc - numDoneX, MAX_BS*BLOCK_SIZE );
     grid.x = todoX/BLOCK_SIZE;
     numDoneY = 0;
     while( numDoneY < dD.pr ){
-      todoY = min( dD.pr - numDoneY, MAX_BS*BLOCK_SIZE );
+      todoY = MIN( dD.pr - numDoneY, MAX_BS*BLOCK_SIZE );
       grid.y = todoY/BLOCK_SIZE;
       rangeSearchKernel<<<grid,block>>>(dD, numDoneX, numDoneY, dranges, dir);
       numDoneY += todoY;
@@ -63,32 +73,57 @@ void rangeSearchWrap(matrix dD, real *dranges, charMatrix dir){
   cudaThreadSynchronize();
 }
 
-void nnWrap(const matrix dx, const matrix dy, real *dMins, int *dMinIDs){
-  dim3 dimBlock(BLOCK_SIZE,BLOCK_SIZE);
-  dim3 dimGrid;
+void nnWrap(const matrix dq, const matrix dx, real *dMins, unint *dMinIDs){
+  dim3 block(BLOCK_SIZE,BLOCK_SIZE);
+  dim3 grid;
+  unint numDone, todo;
   
-  dimGrid.x = 1;
-  dimGrid.y = dx.pr/dimBlock.y + (dx.pr%dimBlock.y==0 ? 0 : 1);
-  nnKernel<<<dimGrid,dimBlock>>>(dx,dy,dMins,dMinIDs);
+  grid.x = 1;
+
+  numDone = 0;
+  while( numDone < dq.pr ){
+    todo = MIN( dq.pr - numDone, MAX_BS*BLOCK_SIZE );
+    grid.y = todo/BLOCK_SIZE;
+    nnKernel<<<grid,block>>>(dq,numDone,dx,dMins,dMinIDs);
+    numDone += todo;
+  }
   cudaThreadSynchronize();
+
 }
 
 
-void rangeCountWrap(const matrix dq, const matrix dx, real *dranges, int *dcounts){
+void planNNWrap(const matrix dq, const unint *dqMap, const matrix dx, const intMatrix dxMap, real *dMins, unint *dMinIDs, compPlan dcP, unint compLength){
   dim3 block(BLOCK_SIZE,BLOCK_SIZE);
-  dim3 grid(1,dq.pr/BLOCK_SIZE);
-
-  rangeCountKernel<<<grid,block>>>(dq,dx,dranges,dcounts);
+  dim3 grid;
+  unint todo;
+
+  grid.x = 1;
+  unint numDone = 0;
+  while( numDone<compLength ){
+    todo = MIN( (compLength-numDone) , MAX_BS*BLOCK_SIZE );
+    grid.y = todo/BLOCK_SIZE;
+    planNNKernel<<<grid,block>>>(dq,dqMap,dx,dxMap,dMins,dMinIDs,dcP,numDone);
+    numDone += todo;
+  }
   cudaThreadSynchronize();
 }
 
 
-/*NOTE: can be deleted */
-void pruneWrap(charMatrix dcM, matrix dD, real *dradiiX, real *dradiiQ){
+void rangeCountWrap(const matrix dq, const matrix dx, real *dranges, unint *dcounts){
   dim3 block(BLOCK_SIZE,BLOCK_SIZE);
-  dim3 grid(dD.pr/BLOCK_SIZE,dD.pc/BLOCK_SIZE);
-  
-  pruneKernel<<<grid,block>>>(dD,dradiiX,dradiiQ,dcM);
+  dim3 grid;
+  unint numDone, todo;
+
+  grid.x=1;
+
+  numDone = 0;
+  while( numDone < dq.pr ){
+    todo = MIN( dq.pr - numDone, MAX_BS*BLOCK_SIZE );
+    grid.y = todo/BLOCK_SIZE;
+    rangeCountKernel<<<grid,block>>>(dq,numDone,dx,dranges,dcounts);
+    numDone += todo;
+  }
   cudaThreadSynchronize();
 }
+
 #endif
index 2d5f80d..b5190fe 100644 (file)
@@ -1,13 +1,17 @@
+/* This file is part of the Random Ball Cover (RBC) library.
+ * (C) Copyright 2010, Lawrence Cayton [lcayton@tuebingen.mpg.de]
+ */
 #ifndef KERNELWRAP_H
 #define KERNELWRAP_H
 
 #include "defs.h"
 
-void dist1Wrap(matrix,matrix,matrix);
+void dist1Wrap(const matrix,const matrix,matrix);
 void kMinsWrap(matrix,matrix,intMatrix);
-void findRangeWrap(matrix,real*,int);
-void rangeSearchWrap(matrix,real*,charMatrix);
-void nnWrap(const matrix,const matrix,real*,int*);
-void rangeCountWrap(const matrix,const matrix,real*,int*);
-void pruneWrap(charMatrix,matrix,real*,real*);
+void findRangeWrap(const matrix,real*,unint);
+void rangeSearchWrap(const matrix,const real*,charMatrix);
+void nnWrap(const matrix,const matrix,real*,unint*);
+void rangeCountWrap(const matrix,const matrix,real*,unint*);
+void planNNWrap(const matrix,const unint*,const matrix,const intMatrix,real*,unint*,compPlan,unint);
+
 #endif
index d73d837..aeccf87 100644 (file)
@@ -1,3 +1,7 @@
+/* This file is part of the Random Ball Cover (RBC) library.
+ * (C) Copyright 2010, Lawrence Cayton [lcayton@tuebingen.mpg.de]
+ */
+
 #ifndef KERNELS_CU
 #define KERNELS_CU
 
 #include<stdio.h>
 
 // This kernel does the same thing as nnKernel, except it only considers pairs as 
-// specified by the compPlan.
-__global__ void planNNKernel(const matrix Q, const matrix X, real *dMins, int *dMinIDs, compPlan cP, intMatrix xmap, int *qIDs, int qStartPos ){
-  int qBlock = qStartPos + blockIdx.y * BLOCK_SIZE;  //indexes Q
-  int xBlock; //indexes X;
-  int colBlock;
-  int offQ = threadIdx.y; //the offset of qPos in this block
-  int offX = threadIdx.x; //ditto for x
-  int i,j,k,l;
-  
+// specified by the compPlan. 
+__global__ void planNNKernel(const matrix Q, const unint *qMap, const matrix X, const intMatrix xMap, real *dMins, unint *dMinIDs, compPlan cP,  unint qStartPos ){
+  unint qB = qStartPos + blockIdx.y * BLOCK_SIZE;  //indexes Q
+  unint xB; //X (DB) Block;
+  unint cB; //column Block
+  unint offQ = threadIdx.y; //the offset of qPos in this block
+  unint offX = threadIdx.x; //ditto for x
+  unint i,j,k;
+  unint groupIts;
   
   __shared__ real min[BLOCK_SIZE][BLOCK_SIZE];
-  __shared__ int minPos[BLOCK_SIZE][BLOCK_SIZE];
+  __shared__ unint minPos[BLOCK_SIZE][BLOCK_SIZE];
 
-  __shared__ real Xb[BLOCK_SIZE][BLOCK_SIZE];
-  __shared__ real Qb[BLOCK_SIZE][BLOCK_SIZE];
+  __shared__ real Xs[BLOCK_SIZE][BLOCK_SIZE];
+  __shared__ real Qs[BLOCK_SIZE][BLOCK_SIZE];
 
-  __shared__ int g; //group of q  
-  __shared__ int numGroups;
-  __shared__ int groupCount;
-  __shared__ int groupOff;
+  unint g; //query group of q
+  unint xG; //DB group currently being examined
+  unint numGroups;
+  unint groupCount;
 
-  if(offX==0 && offQ==0){
-    g = cP.qToGroup[qBlock]; 
-    numGroups = cP.numGroups[g];
-  }
+  g = cP.qToQGroup[qB]; 
+  numGroups = cP.numGroups[g];
+  
   min[offQ][offX]=MAX_REAL;
   __syncthreads();
-      
-  //NOTE: might be better to save numGroup, groupIts in local reg.
-  for(i=0;i<numGroups;i++){ //iterate over groups of X  
-    if(offX==0 && offQ==0){
-      groupCount = cP.groupCountX[IDX(g,i,cP.ld)];
-      groupOff = cP.groupOff[IDX(g,i,cP.ld)];
-    }
-    /* if(qBlock==0 && offQ==0) */
-    /*   printf("g=%d groupCount=%d \n",g,groupCount); */
-    __syncthreads();
-    
-    int groupIts = groupCount/BLOCK_SIZE + (groupCount%BLOCK_SIZE==0? 0 : 1);
-   
-    xBlock=groupOff;
-    for(j=0;j<groupIts;j++){ //iterate over elements of group
-      xBlock=j*BLOCK_SIZE;
-      
-      real ans=0;
-      for(k=0;k<X.pc/BLOCK_SIZE;k++){ // iterate over cols to compute the distances
-       colBlock = k*BLOCK_SIZE;
+  
+
+  for(i=0; i<numGroups; i++){ //iterate over DB groups
+    xG = cP.qGroupToXGroup[IDX( g, i, cP.ld )];
+    groupCount = cP.groupCountX[IDX( g, i, cP.ld )];
+    groupIts = (groupCount+BLOCK_SIZE-1)/BLOCK_SIZE;
 
-       //Each thread loads one element of X and Q into memory.
-       //Note that the indexing is flipped to increase memory
-       //coalescing.
+    for(j=0; j<groupIts; j++){ //iterate over elements of group
+      xB=j*BLOCK_SIZE;
 
-       Xb[offX][offQ] = X.mat[IDX( xmap.mat[IDX( g, xBlock+offQ, xmap.ld)], colBlock+offX, X.ld)];
-       Qb[offX][offQ] = ( (qIDs[qBlock+offQ]==DUMMY_IDX) ? 0 : Q.mat[IDX(qIDs[qBlock+offQ],colBlock+offX,Q.ld)] );
+      real ans=0;
+      for(cB=0; cB<X.pc; cB+=BLOCK_SIZE){ // iterate over cols to compute distances
+
+       Xs[offX][offQ] = X.mat[IDX( xMap.mat[IDX( xG, xB+offQ, xMap.ld )], cB+offX, X.ld )];
+       Qs[offX][offQ] = ( (qMap[qB+offQ]==DUMMY_IDX) ? 0 : Q.mat[IDX( qMap[qB+offQ], cB+offX, Q.ld )] );
        __syncthreads();
        
-       for(l=0;l<BLOCK_SIZE;l++){
-         ans+=abs(Xb[l][offX]-Qb[l][offQ]);
-       }
+       for(k=0; k<BLOCK_SIZE; k++)
+         ans+=DIST( Xs[k][offX], Qs[k][offQ] );
+
        __syncthreads();
       }
      
       //compare to previous min and store into shared mem if needed.
-      if(j*BLOCK_SIZE+offX<groupCount && ans<min[offQ][offX]){
+      if(xB+offX<groupCount && ans<min[offQ][offX]){
        min[offQ][offX]=ans;
-       minPos[offQ][offX]= xmap.mat[IDX( g, xBlock+offX, xmap.ld )];
+       minPos[offQ][offX]= xMap.mat[IDX( xG, xB+offX, xMap.ld )];
       }
       __syncthreads();
-      
     }
   }
   
-  //compare across threads
-  for(k=BLOCK_SIZE/2;k>0;k/=2){
-    if(offX<k){
-      if(min[offQ][offX+k]<min[offQ][offX]){
-       min[offQ][offX] = min[offQ][offX+k];
-       minPos[offQ][offX] = minPos[offQ][offX+k];      
+  //Reduce across threads
+  for(i=BLOCK_SIZE/2; i>0; i/=2){
+    if( offX<i ){
+      if( min[offQ][offX+i] < min[offQ][offX] ){
+       min[offQ][offX] = min[offQ][offX+i];
+       minPos[offQ][offX] = minPos[offQ][offX+i];      
       }
     }
     __syncthreads();
   }
-  
-  if(offX==0 && qIDs[qBlock+offQ]!=DUMMY_IDX){
-    dMins[qIDs[qBlock+offQ]] = min[offQ][0];
-    dMinIDs[qIDs[qBlock+offQ]] = minPos[offQ][0];
-  }
-}
-
 
-
-__global__ void pruneKernel(const matrix D, const real *radiiX, const real *radiiQ, charMatrix cM){
-  int offX = threadIdx.x;
-  int offQ = threadIdx.y;
-
-  int blockX = blockIdx.x * BLOCK_SIZE;
-  int blockQ = blockIdx.y * BLOCK_SIZE;
-  
-  __shared__ real sD[BLOCK_SIZE][BLOCK_SIZE];
-  __shared__ real sRQ[BLOCK_SIZE];
-  __shared__ real sRX[BLOCK_SIZE];
-
-  sD[offQ][offX]=D.mat[IDX(blockQ+offQ,blockX+offX,D.ld)];
-  
-  if(offQ==0)
-    sRX[offX]=radiiX[blockX+offX];
-  if(offX==0)
-    sRQ[offQ]=radiiQ[blockQ+offQ];
-  
-  __syncthreads();
-  
-  if(blockQ+offQ < D.r && blockX+offX < D.c){
-    cM.mat[IDX(blockQ+offQ,blockX+offX,cM.ld)] = (sD[offQ][offX]-sRX[offX]-2*sRQ[offQ] <= 0) ? 1 : 0;
-    //cM.mat[IDX(blockQ+offQ,blockX+offX,cM.ld)] = (sD[offQ][offX]-4*sRQ[offQ] <= 0) ? 1 : 0;
+  if(offX==0 && qMap[qB+offQ]!=DUMMY_IDX){
+    dMins[qMap[qB+offQ]] = min[offQ][0];
+    dMinIDs[qMap[qB+offQ]] = minPos[offQ][0];
   }
 }
 
 
-__global__ void nnKernel(const matrix Q, const matrix X, real *dMins, int *dMinIDs){
 
-  int qBlock = blockIdx.y * BLOCK_SIZE;  //indexes Q
-  int xBlock; //indexes X;
-  int colBlock;
-  int offQ = threadIdx.y; //the offset of qPos in this block
-  int offX = threadIdx.x; //ditto for x
-  int i,j,k;
+__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
+  unint offQ = threadIdx.y; //the offset of qPos in this block
+  unint offX = threadIdx.x; //ditto for x
+  unint i;
+  real ans;
 
   __shared__ real min[BLOCK_SIZE][BLOCK_SIZE];
-  __shared__ int minPos[BLOCK_SIZE][BLOCK_SIZE];
+  __shared__ unint minPos[BLOCK_SIZE][BLOCK_SIZE];
 
-  __shared__ real Xb[BLOCK_SIZE][BLOCK_SIZE];
-  __shared__ real Qb[BLOCK_SIZE][BLOCK_SIZE];
+  __shared__ real Xs[BLOCK_SIZE][BLOCK_SIZE];
+  __shared__ real Qs[BLOCK_SIZE][BLOCK_SIZE];
 
   min[offQ][offX]=MAX_REAL;
   __syncthreads();
 
-  for(i=0;i<X.pr/BLOCK_SIZE;i++){
-    xBlock = i*BLOCK_SIZE;
-    real ans=0;
-    for(j=0;j<X.pc/BLOCK_SIZE;j++){
-      colBlock = j*BLOCK_SIZE;
+  for(xB=0; xB<X.pr; xB+=BLOCK_SIZE){
+    ans=0;
+    for(cB=0; cB<X.pc; cB+=BLOCK_SIZE){
       
       //Each thread loads one element of X and Q into memory.
-      //Note that the indexing is flipped to increase memory
-      //coalescing.
-      Xb[offX][offQ]=X.mat[IDX(xBlock+offQ,colBlock+offX,X.ld)];
-      Qb[offX][offQ]=Q.mat[IDX(qBlock+offQ,colBlock+offX,Q.ld)];
-
+      Xs[offX][offQ] = X.mat[IDX( xB+offQ, cB+offX, X.ld )];
+      Qs[offX][offQ] = Q.mat[IDX( qB+offQ, cB+offX, Q.ld )];
+      
       __syncthreads();
-
-      for(k=0;k<BLOCK_SIZE;k++){
-       ans+=abs(Xb[k][offX]-Qb[k][offQ]);
-      }
+      
+      for(i=0;i<BLOCK_SIZE;i++)
+       ans += DIST( Xs[i][offX], Qs[i][offQ] );
+      
       __syncthreads();
     }
    
-    
-    if( xBlock+offX<X.r && ans<min[offQ][offX] ){
-       minPos[offQ][offX] = xBlock+offX;
-       min[offQ][offX] = ans;
+    if( xB+offX<X.r && ans<min[offQ][offX] ){
+      minPos[offQ][offX] = xB+offX;
+      min[offQ][offX] = ans;
     }
   }
   __syncthreads();
   
   
   //reduce across threads
-  for(j=BLOCK_SIZE/2;j>0;j/=2){
-    if(offX<j){
-      if(min[offQ][offX+j]<min[offQ][offX]){
-       min[offQ][offX] = min[offQ][offX+j];
-       minPos[offQ][offX] = minPos[offQ][offX+j];      
+  for(i=BLOCK_SIZE/2; i>0; i/=2){
+    if(offX<i){
+      if(min[offQ][offX+i]<min[offQ][offX]){
+       min[offQ][offX] = min[offQ][offX+i];
+       minPos[offQ][offX] = minPos[offQ][offX+i];      
       }
     }
     __syncthreads();
   }
   
   if(offX==0){
-    //printf("writing %d, nn= %d, val = %6.4f \n",qBlock+offQ,curMinPos[offQ],curMin[offQ]);
-    dMins[qBlock+offQ] = min[offQ][0];
-    dMinIDs[qBlock+offQ] = minPos[offQ][0];
+    dMins[qB+offQ] = min[offQ][0];
+    dMinIDs[qB+offQ] = minPos[offQ][0];
   }
 }
 
 
-__device__ void dist1Kernel(const matrix Q, int qStart, const matrix X, int xStart, matrix D){
-  int c, i, j;
 
-  int qB = blockIdx.y*BLOCK_SIZE + qStart;
-  int q  = threadIdx.y;
-  int xB = blockIdx.x*BLOCK_SIZE + xStart;
-  int x = threadIdx.x;
+
+__global__ void dist1Kernel(const matrix Q, unint qStart, const matrix X, unint xStart, matrix D){
+  unint c, i, j;
+
+  unint qB = blockIdx.y*BLOCK_SIZE + qStart;
+  unint q  = threadIdx.y;
+  unint xB = blockIdx.x*BLOCK_SIZE + xStart;
+  unint x = threadIdx.x;
 
   real ans=0;
 
@@ -217,8 +176,8 @@ __device__ void dist1Kernel(const matrix Q, int qStart, const matrix X, int xSta
     __syncthreads();
 
     for(j=0 ; j<BLOCK_SIZE ; j++)
-      ans += abs( Qs[j][q] - Xs[j][x] );
-        
+      ans += DIST( Qs[j][q], Xs[j][x] );
+    
     __syncthreads();
   }
   
@@ -228,16 +187,16 @@ __device__ void dist1Kernel(const matrix Q, int qStart, const matrix X, int xSta
 
 
 
-__global__ void findRangeKernel(matrix D, real *ranges, int cntWant){
+__global__ void findRangeKernel(const matrix D, unint numDone, real *ranges, unint cntWant){
   
-  int row = blockIdx.y*(BLOCK_SIZE/4)+threadIdx.y;
-  int ro = threadIdx.y;
-  int co = threadIdx.x;
-  int i, c;
+  unint row = blockIdx.y*(BLOCK_SIZE/4)+threadIdx.y + numDone;
+  unint ro = threadIdx.y;
+  unint co = threadIdx.x;
+  unint i, c;
   real t;
 
-  const int LB = (90*cntWant)/100 ;
-  const int UB = cntWant; 
+  const unint LB = (90*cntWant)/100 ;
+  const unint UB = cntWant; 
 
   __shared__ real smin[BLOCK_SIZE/4][4*BLOCK_SIZE];
   __shared__ real smax[BLOCK_SIZE/4][4*BLOCK_SIZE];
@@ -266,10 +225,10 @@ __global__ void findRangeKernel(matrix D, real *ranges, int cntWant){
 
   //Now start range counting.
 
-  int itcount=0;
-  int cnt;
+  unint itcount=0;
+  unint cnt;
   real rg;
-  __shared__ int scnt[BLOCK_SIZE/4][4*BLOCK_SIZE];
+  __shared__ unint scnt[BLOCK_SIZE/4][4*BLOCK_SIZE];
   __shared__ char cont[BLOCK_SIZE/4];
   
   if(co==0)
@@ -324,30 +283,30 @@ __global__ void findRangeKernel(matrix D, real *ranges, int cntWant){
 }
 
 
-__global__ void rangeSearchKernel(matrix D, int xOff, int yOff, real *ranges, charMatrix ir){
-  int col = blockIdx.x*BLOCK_SIZE + threadIdx.x + xOff;
-  int row = blockIdx.y*BLOCK_SIZE + threadIdx.y + yOff;
+__global__ void rangeSearchKernel(const matrix D, unint xOff, unint yOff, const real *ranges, charMatrix ir){
+  unint col = blockIdx.x*BLOCK_SIZE + threadIdx.x + xOff;
+  unint row = blockIdx.y*BLOCK_SIZE + threadIdx.y + yOff;
 
   ir.mat[IDX( row, col, ir.ld )] = D.mat[IDX( row, col, D.ld )] < ranges[row];
 
 }
 
 
-__global__ void rangeCountKernel(const matrix Q, const matrix X, real *ranges, int *counts){
-  int q = blockIdx.y*BLOCK_SIZE;
-  int qo = threadIdx.y;
-  int xo = threadIdx.x;
+__global__ void rangeCountKernel(const matrix Q, unint numDone, const matrix X, real *ranges, unint *counts){
+  unint q = blockIdx.y*BLOCK_SIZE + numDone;
+  unint qo = threadIdx.y;
+  unint xo = threadIdx.x;
   
   real rg = ranges[q+qo];
   
-  int r,c,i;
+  unint r,c,i;
 
-  __shared__ int scnt[BLOCK_SIZE][BLOCK_SIZE];
+  __shared__ unint scnt[BLOCK_SIZE][BLOCK_SIZE];
 
   __shared__ real xs[BLOCK_SIZE][BLOCK_SIZE];
   __shared__ real qs[BLOCK_SIZE][BLOCK_SIZE];
   
-  int cnt=0;
+  unint cnt=0;
   for( r=0; r<X.pr; r+=BLOCK_SIZE ){
 
     real dist=0;
@@ -357,7 +316,8 @@ __global__ void rangeCountKernel(const matrix Q, const matrix X, real *ranges, i
       __syncthreads();
       
       for( i=0; i<BLOCK_SIZE; i++)
-       dist += abs(xs[i][xo]-qs[i][qo]);
+       dist += DIST( xs[i][xo], qs[i][qo] );
+
       __syncthreads();
 
     }
index cf26734..4d03611 100644 (file)
--- a/kernels.h
+++ b/kernels.h
@@ -1,15 +1,16 @@
+/* This file is part of the Random Ball Cover (RBC) library.
+ * (C) Copyright 2010, Lawrence Cayton [lcayton@tuebingen.mpg.de]
+ */
+
 #ifndef KERNELS_H
 #define KERNELS_H
 
 #include "defs.h"
-__global__ void planNNKernel(const matrix,const matrix,real*,int*,compPlan,intMatrix,int*,int);
-__global__ void pruneKernel(const matrix,const real*,const real*,charMatrix);
-__global__ void dist1Kernel(const matrix,int,const matrix,int,matrix);
-__device__ matrix getSubMat(matrix,int,int);
-__global__ void nnKernel(const matrix,const matrix,real*,int*);
-
-__global__ void findRangeKernel(matrix,real*,int);
-__global__ void rangeSearchKernel(matrix,int,int,real*,charMatrix);
-__global__ void rangeCountKernel(const matrix,const matrix,real*,int*);
+__global__ void planNNKernel(const matrix,const unint*,const matrix,const intMatrix,real*,unint*,compPlan,unint);
+__global__ void dist1Kernel(const matrix,unint,const matrix,unint,matrix);
+__global__ void nnKernel(const matrix,unint,const matrix,real*,unint*);
+__global__ void findRangeKernel(const matrix,unint,real*,unint);
+__global__ void rangeSearchKernel(const matrix,unint,unint,const real*,charMatrix);
+__global__ void rangeCountKernel(const matrix,unint,const matrix,real*,unint*);
 
 #endif
diff --git a/rbc.cu b/rbc.cu
index f6cb34d..51607f9 100644 (file)
--- a/rbc.cu
+++ b/rbc.cu
@@ -1,3 +1,7 @@
+/* This file is part of the Random Ball Cover (RBC) library.
+ * (C) Copyright 2010, Lawrence Cayton [lcayton@tuebingen.mpg.de]
+ */
+
 #ifndef RBC_CU
 #define RBC_CU
 
 #include "rbc.h"
 #include "kernels.h"
 #include "kernelWrap.h"
-#include "sKernel.h"
-
-void rbc(matrix x, matrix q, int numReps, int s, int *NNs){
-  int i;
-  matrix r; // random subset of x
-  matrix dr; //device version of r
-  int *repIDsQ; //assigment of each pt to a rep.
-  real *radiiQ; //radius of each group
-  int *groupCountX, *groupCountQ; //num pts in each group
-  real *distToRepsQ; //distance of each pt to its rep
-  int *qID;
-  charMatrix cM; //compute matrix.  will be treated as a binary matrix
-  intMatrix xmap;
-  compPlan cP; 
-  int compLength;
-  struct timeval tvB, tvE;
-      
-  //convenience variables
-  int n = x.r; //num pts in DB
-  int m = q.r; //num queries
-  int pnr = PAD(numReps);
-  
-  r.r=numReps; r.c=x.c; r.pr=pnr; r.pc=PAD(r.c); r.ld=r.pc;
-  r.mat = (real*)calloc(r.pr*r.pc,sizeof(*(r.mat)));
-
-  repIDsQ = (int*)calloc(m,sizeof(*repIDsQ));
-  radiiQ = (real*)calloc(pnr,sizeof(*radiiQ));
-  groupCountX = (int*)calloc(pnr,sizeof(*groupCountX));
-  groupCountQ = (int*)calloc(pnr,sizeof(*groupCountQ));
-  distToRepsQ = (real*)calloc(m,sizeof(*distToRepsQ));
-
-  qID = (int*)calloc(PAD(m+(BLOCK_SIZE-1)*pnr),sizeof(*qID));
-
-  cM.mat = (char*)calloc(pnr*pnr,sizeof(*cM.mat));
-  cM.r=numReps; cM.c=numReps; cM.pr=pnr; cM.pc=pnr; cM.ld=cM.pc;
-
-  xmap.r=numReps; xmap.pr=PAD(numReps); xmap.c=s; xmap.pc=PAD(s); xmap.ld=xmap.pc;
-  cudaMalloc( (void**)&xmap.mat, xmap.pr*xmap.pc*sizeof(*xmap.mat) );
-
-  //The following lines of code init xmap to 0s.
-  intMatrix tempXmap;
-  tempXmap.r=xmap.r; tempXmap.pr=xmap.pr; tempXmap.c=xmap.c; tempXmap.pc=xmap.pc; tempXmap.ld=xmap.ld;
-  tempXmap.mat = (int*)calloc( tempXmap.pr*tempXmap.pc, sizeof(*tempXmap.mat) );
-  cudaMemcpy(xmap.mat, tempXmap.mat, xmap.pr*xmap.pc*sizeof(*xmap.mat), cudaMemcpyHostToDevice);
-  free(tempXmap.mat);
+#include "sKernelWrap.h"
 
+void queryRBC(const matrix q, const rbcStruct rbcS, unint *NNs){
+  unint m = q.r;
+  unint numReps = rbcS.dr.r;
+  unint compLength;
+  compPlan dcP;
+  unint *qMap, *dqMap;
+  qMap = (unint*)calloc(PAD(m+(BLOCK_SIZE-1)*PAD(numReps)),sizeof(*qMap));
+  matrix dq;
+  copyAndMove(&dq, &q);
   
-  //Choose representatives and move them to device
-  int *randInds;
-  randInds = (int*)calloc(pnr,sizeof(*randInds));
-  subRandPerm(numReps,n,randInds);
-  
-  for(i=0;i<numReps;i++){
-    copyVector(&r.mat[IDX(i,0,r.ld)], &x.mat[IDX(randInds[i],0,x.ld)], x.c);
-  }
-  free(randInds);
-
-  dr.r=r.r; dr.c=r.c; dr.pr=r.pr; dr.pc=r.pc; dr.ld=r.ld;
-  cudaMalloc( (void**)&(dr.mat), dr.pr*dr.pc*sizeof(*(dr.mat)) );
-  cudaMemcpy(dr.mat,r.mat,dr.pr*dr.pc*sizeof(*(dr.mat)),cudaMemcpyHostToDevice);
-
-  
-  matrix dx;
-  copyAndMove(&dx, &x);
-  
-  build(dx, dr, xmap, groupCountX, s); 
+  charMatrix cM;
+  cM.r=cM.c=numReps; cM.pr=cM.pc=cM.ld=PAD(numReps);
+  cM.mat = (char*)calloc( cM.pr*cM.pc, sizeof(*cM.mat) );
   
-  /* tempXmap.r=xmap.r; tempXmap.pr=xmap.pr; tempXmap.c=xmap.c; tempXmap.pc=xmap.pc; tempXmap.ld=xmap.ld; */
-  /* tempXmap.mat = (int*)calloc( tempXmap.pr*tempXmap.pc, sizeof(*tempXmap.mat) ); */
-  /* cudaMemcpy(tempXmap.mat, xmap.mat, xmap.pr*xmap.pc*sizeof(*xmap.mat), cudaMemcpyDeviceToHost); */
-  /* for( i=0; i<16; i++ ) */
-  /*   printf("%d ",tempXmap.mat[i]); */
-  /* printf("\n"); */
-  /* free(tempXmap.mat); */
-
-
-  gettimeofday(&tvB,NULL);  //Start the timer for the queries
+  unint *repIDsQ;
+  repIDsQ = (unint*)calloc( m, sizeof(*repIDsQ) );
+  real *distToRepsQ;
+  distToRepsQ = (real*)calloc( m, sizeof(*distToRepsQ) );
+  unint *groupCountQ;
+  groupCountQ = (unint*)calloc( PAD(numReps), sizeof(*groupCountQ) );
   
-  matrix dqOrig;
-  copyAndMove(&dqOrig, &q);
-
-  computeReps(dqOrig, dr, repIDsQ, distToRepsQ);
+  computeReps(dq, rbcS.dr, repIDsQ, distToRepsQ);
 
   //How many points are assigned to each group?
   computeCounts(repIDsQ, m, groupCountQ);
   
+  //Set up the mapping from groups to queries (qMap).
+  buildQMap(q, qMap, repIDsQ, numReps, &compLength);
   
-  //Set up the mapping from groups to queries (qID).
-  buildQMap(q, qID, repIDsQ, numReps, &compLength);
-
-  // Determine which blocks need to be compared with one another and
-  // store the results in the computation matrix cM.
-  //blockIntersection(cM, dr, radiiX, radiiQ);
+  // 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.  
   idIntersection(cM);
 
-  // Setup the compute plan
-  initCompPlan(&cP, cM, groupCountQ, groupCountX, numReps);
+  initCompPlan(&dcP, cM, groupCountQ, rbcS.groupCount, numReps);
 
-  // Compute the NNs according to the compute plan
-  computeNNs(dx, dqOrig, xmap, cP, qID, NNs, compLength);
+  checkErr( cudaMalloc( (void**)&dqMap, compLength*sizeof(*dqMap) ) );
+  cudaMemcpy( dqMap, qMap, compLength*sizeof(*dqMap), cudaMemcpyHostToDevice );
   
-  gettimeofday(&tvE,NULL);
-  printf("\t.. time elapsed (for queries) = %6.4f \n",timeDiff(tvB,tvE));
-
-  cudaFree(dr.mat); 
-  cudaFree(dqOrig.mat);
-  freeCompPlan(&cP);
+  computeNNs(rbcS.dx, rbcS.dxMap, dq, dqMap, dcP, NNs, compLength);
+  
+  free(qMap);
+  freeCompPlan(&dcP);
+  cudaFree(dq.mat);
   free(cM.mat);
+  free(repIDsQ);
   free(distToRepsQ);
-  free(groupCountX);
   free(groupCountQ);
-  free(qID);
-  free(radiiQ);
-  free(repIDsQ);
-  free(r.mat);
-  cudaFree(xmap.mat);
-  cudaFree(dx.mat);
 }
 
 
-void build(const matrix dx, const matrix dr, intMatrix xmap, int *counts, int s){
-  
-  int n = dx.pr;
-  int p = dr.r;
+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;
 
+  setupReps(x, rbcS, numReps);
+  copyAndMove(&rbcS->dx, &x);
+  
+  xmap.r=numReps; xmap.pr=PAD(numReps); xmap.c=s; xmap.pr=xmap.ld=PAD(s);
+  xmap.mat = (unint*)calloc( xmap.pr*xmap.pc, sizeof(*xmap.mat) );
+  copyAndMoveI(&rbcS->dxMap, &xmap);
+  rbcS->groupCount = (uint*)calloc( PAD(numReps), sizeof(*rbcS->groupCount) );
+  
   //Figure out how much fits into memory
-  unsigned int memFree, memTot;
+  unint memFree, memTot;
   cuMemGetInfo(&memFree, &memTot);
-  memFree = (int)(((float)memFree)*MEM_USABLE);
-     
-  //mem needed per rep = n*sizeof(real)+n*sizeof(char)+n*sizeof(int)+sizeof(real)+sizeof(int)
-  //                   = dist mat      +ir            +dSums        +range       +dCnts
-  int ptsAtOnce = DPAD(memFree/((n+1)*sizeof(real) + n*sizeof(char) +(n+1)*sizeof(int)));
-  if(ptsAtOnce==0){
-    fprintf(stderr,"memfree = %d \n",memFree);
-    fprintf(stderr,"error: not enough memory to build the RBC.. exiting\n");
+  memFree = (unint)(((float)memFree)*MEM_USABLE);
+  /* mem needed per rep:
+   *  n*sizeof(real) - dist mat
+   *  n*sizeof(char) - dir
+   *  n*sizeof(int)  - dSums
+   *  sizeof(real)   - dranges
+   *  sizeof(int)    - dCnts
+   *  MEM_USED_IN_SCAN - memory used internally
+   */
+  unint ptsAtOnce = DPAD(memFree/((n+1)*sizeof(real) + n*sizeof(char) + (n+1)*sizeof(unint) + 2*MEM_USED_IN_SCAN(n)));
+  if(!ptsAtOnce){
+    fprintf(stderr,"error: %d is not enough memory to build the RBC.. exiting\n", memFree);
     exit(1);
   }
 
+  //Now set everything up for the scans
   matrix dD;
-  dD.pr=dD.r=ptsAtOnce; dD.c=dx.r; dD.pc=dx.pr; dD.ld=dD.pc;
-  cudaMalloc( (void**)&dD.mat, dD.pr*dD.pc*sizeof(*dD.mat) );
-
+  dD.pr=dD.r=ptsAtOnce; dD.c=rbcS->dx.r; dD.pc=rbcS->dx.pr; dD.ld=dD.pc;
+  checkErr( cudaMalloc( (void**)&dD.mat, dD.pr*dD.pc*sizeof(*dD.mat) ) );
+  
   real *dranges;
-  cudaMalloc( (void**)&dranges, ptsAtOnce*sizeof(real) );
+  checkErr( cudaMalloc( (void**)&dranges, ptsAtOnce*sizeof(real) ) );
 
   charMatrix ir;
   ir.r=dD.r; ir.pr=dD.pr; ir.c=dD.c; ir.pc=dD.pc; ir.ld=dD.ld;
@@ -164,40 +114,36 @@ void build(const matrix dx, const matrix dr, intMatrix xmap, int *counts, int s)
 
   intMatrix dSums; //used to compute memory addresses.
   dSums.r=dir.r; dSums.pr=dir.pr; dSums.c=dir.c; dSums.pc=dir.pc; dSums.ld=dir.ld;
-  cudaMalloc( (void**)&dSums.mat, dSums.pc*dSums.pr*sizeof(*dSums.mat) );
+  checkErr( cudaMalloc( (void**)&dSums.mat, dSums.pc*dSums.pr*sizeof(*dSums.mat) ) );
 
-  int *dCnts;
-  cudaMalloc( (void**)&dCnts, ptsAtOnce*sizeof(*dCnts) );
-    
-   
-  int numits=0;
-  int numLeft = p; //points left to process
-  int row = 0; //base row for iteration of while loop
-  int pi, pip; //pi=pts per it, pip=pad(pi)
+  unint *dCnts;
+  checkErr( cudaMalloc( (void**)&dCnts, ptsAtOnce*sizeof(*dCnts) ) );
   
+  //Do the scans to build the dxMap
+  unint numLeft = rbcS->dr.r; //points left to process
+  unint row = 0; //base row for iteration of while loop
+  unint pi, pip; //pi=pts per it, pip=pad(pi)
   while( numLeft > 0 ){
-    numits++;
     pi = MIN(ptsAtOnce, numLeft);  //points to do this iteration.
     pip = PAD(pi);
     dD.r = pi; dD.pr = pip; dir.r=pi; dir.pr=pip; dSums.r=pi; dSums.pr=pip;
-  
-    distSubMat(dr, dx, dD, row, pip); //compute the distance matrix
+
+    distSubMat(rbcS->dr, rbcS->dx, dD, row, pip); //compute the distance matrix
     findRangeWrap(dD, dranges, s);  //find an appropriate range
     rangeSearchWrap(dD, dranges, dir); //set binary vector for points in range
-    
     sumWrap(dir, dSums);  //This and the next call perform the parallel compaction.
-
-    buildMapWrap(xmap, dir, dSums, row);
+    buildMapWrap(rbcS->dxMap, dir, dSums, row);
     getCountsWrap(dCnts,dir,dSums);  //How many points are assigned to each rep?  It is not 
                                      //*exactly* s, which is why we need to compute this.
-    cudaMemcpy( &counts[row], dCnts, pi*sizeof(*counts), cudaMemcpyDeviceToHost );
+    cudaMemcpy( &rbcS->groupCount[row], dCnts, pi*sizeof(*rbcS->groupCount), cudaMemcpyDeviceToHost );
     
     numLeft -= pi;
     row += pi;
   }
+  
   cudaFree(dCnts);
   free(ir.mat);
+  free(xmap.mat);
   cudaFree(dranges);
   cudaFree(dir.mat);
   cudaFree(dSums.mat);
@@ -205,13 +151,34 @@ void build(const matrix dx, const matrix dr, intMatrix xmap, int *counts, int s)
 }
 
 
+// Choose representatives and move them to device
+void setupReps(matrix x, rbcStruct *rbcS, int numReps){
+  unint i;
+  unint *randInds;
+  randInds = (unint*)calloc( PAD(numReps), sizeof(*randInds) );
+  subRandPerm(numReps, x.r, randInds);
+  
+  matrix r;
+  r.r=numReps; r.pr=PAD(numReps); r.c=x.c; r.pc=r.ld=PAD(r.c); 
+  r.mat = (real*)calloc( r.pr*r.pc, sizeof(*r.mat) );
+
+  for(i=0;i<numReps;i++)
+    copyVector(&r.mat[IDX(i,0,r.ld)], &x.mat[IDX(randInds[i],0,x.ld)], x.c);
+  
+  copyAndMove(&rbcS->dr, &r);
+
+  free(randInds);
+  free(r.mat);
+}
+
+
 //Assign each point in dq to its nearest point in dr.  
-void computeReps(matrix dq, matrix dr, int *repIDs, real *distToReps){
+void computeReps(matrix dq, matrix dr, unint *repIDs, real *distToReps){
   real *dMins;
-  int *dMinIDs;
+  unint *dMinIDs;
 
-  cudaMalloc((void**)&(dMins), dq.pr*sizeof(*dMins));
-  cudaMalloc((void**)&(dMinIDs), dq.pr*sizeof(*dMinIDs));
+  checkErr( cudaMalloc((void**)&(dMins), dq.pr*sizeof(*dMins)) );
+  checkErr( cudaMalloc((void**)&(dMinIDs), dq.pr*sizeof(*dMinIDs)) );
   
   nnWrap(dq,dr,dMins,dMinIDs);
   
@@ -223,10 +190,9 @@ void computeReps(matrix dq, matrix dr, int *repIDs, real *distToReps){
 }
 
 
-
 //Assumes radii is initialized to 0s
-void computeRadii(int *repIDs, real *distToReps, real *radii, int n, int numReps){
-  int i;
+void computeRadii(unint *repIDs, real *distToReps, real *radii, unint n, unint numReps){
+  unint i;
 
   for(i=0;i<n;i++)
     radii[repIDs[i]] = MAX(distToReps[i],radii[repIDs[i]]);
@@ -234,71 +200,20 @@ void computeRadii(int *repIDs, real *distToReps, real *radii, int n, int numReps
 
 
 //Assumes groupCount is initialized to 0s
-void computeCounts(int *repIDs, int n, int *groupCount){
-  int i;
+void computeCounts(unint *repIDs, unint n, unint *groupCount){
+  unint i;
   
   for(i=0;i<n;i++)
     groupCount[repIDs[i]]++;
 }
 
 
-// This function computes a cumulative sum of the groupCounts.
-// Assumes groupOff is initialized to 0s.  
-void computeOffsets(int *groupCount, int n, int *groupOff){
-  int i;
-
-  for(i=1;i<n;i++)
-    groupOff[i]=groupOff[i-1]+PAD(groupCount[i-1]);
-}
-
-
-// This function sorts points by their repID.  It makes two passes through the 
-// matrix x; one to count the bucket sizes, the next to place points in the 
-// correct bucket.  Note that this function unfortunately allocates a temporary
-// matrix the size of x, then copies the results over to x at the end.  The 
-// sort could be done in place, eg by doing numReps passes through x instead of 2.
-void groupPoints(matrix x, int *xID, int *repIDs, int numReps){
-  matrix y;
-  int n=x.r;
-  int d=x.c;
-  int i;
-  int *gS; //groupSize
-  int *yID;
-
-  yID = (int*)calloc(n,sizeof(*yID));
-  y.mat = (real*)calloc(n*d,sizeof(*y.mat));
-  gS = (int*)calloc(numReps+1,sizeof(*gS));
-
-  y.r=n; y.pr=n; y.c=d; y.pc=d; y.ld=d;
-
-  for(i=0;i<n;i++)
-    gS[repIDs[i]+1]++;
-  for(i=1;i<numReps;i++)
-    gS[i]=gS[i-1]+gS[i];
-  
-  for(i=0;i<n;i++){
-    copyVector(&y.mat[IDX(gS[repIDs[i]],0,y.ld)], &x.mat[IDX(i,0,x.ld)],d);
-    yID[gS[repIDs[i]]]=xID[i];
-    gS[repIDs[i]]++;
-  }
+void buildQMap(matrix q, unint *qMap, unint *repIDs, unint numReps, unint *compLength){
+  unint n=q.r;
+  unint i;
+  unint *gS; //groupSize
   
-  for(i=0;i<n;i++){
-    copyVector(&x.mat[IDX(i,0,x.ld)],&y.mat[IDX(i,0,y.ld)],d);
-    xID[i]=yID[i];
-  }
-  
-  free(yID);
-  free(gS);
-  free(y.mat);
-}
-
-
-void buildQMap(matrix q, int *qID, int *repIDs, int numReps, int *compLength){
-  int n=q.r;
-  int i;
-  int *gS; //groupSize
-  
-  gS = (int*)calloc(numReps+1,sizeof(*gS));
+  gS = (unint*)calloc(numReps+1,sizeof(*gS));
   
   for( i=0; i<n; i++ )
     gS[repIDs[i]+1]++;
@@ -311,10 +226,10 @@ void buildQMap(matrix q, int *qID, int *repIDs, int numReps, int *compLength){
   *compLength = gS[numReps];
   
   for( i=0; i<(*compLength); i++ )
-    qID[i]=DUMMY_IDX;
+    qMap[i]=DUMMY_IDX;
   
   for( i=0; i<n; i++ ){
-    qID[gS[repIDs[i]]]=i;
+    qMap[gS[repIDs[i]]]=i;
     gS[repIDs[i]]++;
   }
 
@@ -322,40 +237,9 @@ void buildQMap(matrix q, int *qID, int *repIDs, int numReps, int *compLength){
 }
 
 
-void blockIntersection(charMatrix cM, matrix dr, real *radiiX, real *radiiQ){
-  matrix dD;
-  real *dradiiX, *dradiiQ;
-  int pnR = dr.pr;
-  charMatrix dcM;
-  
-  dD.r=dD.c=dr.r; dD.pr=dD.pc=dD.ld=dr.pr;
-  dcM.r=cM.r; dcM.c=cM.c; dcM.pr=cM.pr; dcM.pc=cM.pc; dcM.ld=cM.ld;
-  
-  cudaMalloc((void**)&dD.mat, pnR*pnR*sizeof(*dD.mat));
-  cudaMalloc((void**)&dradiiX, pnR*sizeof(*dradiiX));
-  cudaMalloc((void**)&dradiiQ, pnR*sizeof(*dradiiQ));
-  cudaMalloc((void**)&dcM.mat, dcM.pr*dcM.pc*sizeof(*dcM.mat));
-  
-  // Copying over the radii. Note that everything after the first dr.r places 
-  // on the device variables is undefined.
-  cudaMemcpy(dradiiX,radiiX,dr.r*sizeof(*dradiiX),cudaMemcpyHostToDevice);
-  cudaMemcpy(dradiiQ,radiiQ,dr.r*sizeof(*dradiiQ),cudaMemcpyHostToDevice);
-  
-  dist1Wrap(dr, dr, dD);
-  pruneWrap(dcM, dD, dradiiX, dradiiQ);
-
-  cudaMemcpy(cM.mat,dcM.mat,pnR*pnR*sizeof(*dcM.mat),cudaMemcpyDeviceToHost);
-  
-  cudaFree(dcM.mat);
-  cudaFree(dradiiQ);
-  cudaFree(dradiiX);
-  cudaFree(dD.mat);
-}
-
-
 // Sets the computation matrix to the identity.  
 void idIntersection(charMatrix cM){
-  int i;
+  unint i;
   for(i=0;i<cM.r;i++){
     if(i<cM.c)
       cM.mat[IDX(i,i,cM.ld)]=1;
@@ -365,7 +249,7 @@ void idIntersection(charMatrix cM){
 
 
 void fullIntersection(charMatrix cM){
-  int i,j;
+  unint i,j;
   for(i=0;i<cM.r;i++){
     for(j=0;j<cM.c;j++){
       cM.mat[IDX(i,j,cM.ld)]=1;
@@ -374,117 +258,19 @@ void fullIntersection(charMatrix cM){
 }
 
 
-// Danger: this function allocates memory that it does not free.  
-// Use freeCompPlan to clear mem.
-void initCompPlan(compPlan *cP, charMatrix cM, int *groupCountQ, int *groupCountX, int numReps){
-  int i,j,k;
-  int maxNumGroups=0;
-  int *groupOff;
-  
-  groupOff = (int*)calloc(numReps,sizeof(*groupOff));
-
-  cP->sNumGroups=numReps;
-  cP->numGroups = (int*)calloc(cP->sNumGroups,sizeof(*cP->numGroups));
-  
-  for(i=0;i<numReps;i++){
-    cP->numGroups[i] = 0;
-    for(j=0;j<numReps;j++){
-      cP->numGroups[i] += cM.mat[IDX(i,j,cM.ld)];
-    }
-    maxNumGroups = MAX(cP->numGroups[i],maxNumGroups);
-  }
-  
-  cP->ld = maxNumGroups;
-
-  cP->sGroupCountX = maxNumGroups*numReps;
-  cP->groupCountX = (int*)calloc(cP->sGroupCountX,sizeof(*cP->groupCountX));
-  cP->sGroupOff = maxNumGroups*numReps;
-  cP->groupOff = (int*)calloc(cP->sGroupOff,sizeof(*cP->groupOff));
-  
-  computeOffsets(groupCountX,numReps,groupOff);
-  
-  int tempSize=0;
-  for(i=0;i<numReps;i++)
-    tempSize+=PAD(groupCountQ[i]);
-  
-  cP->sQToGroup = tempSize;
-  cP->qToGroup = (int*)calloc(cP->sQToGroup,sizeof(*cP->qToGroup));
-  
-  for(i=0, k=0;i<numReps;i++){
-    for(j=0;j<PAD(groupCountQ[i]);j++){
-      cP->qToGroup[k]=i;
-      k++;
-    }
-  }
-  
-  for(i=0;i<numReps;i++){
-    for(j=0, k=0;j<numReps;j++){
-      if(cM.mat[IDX(i,j,cM.ld)]){
-       cP->groupCountX[IDX(i,k,cP->ld)]=groupCountX[j];
-       cP->groupOff[IDX(i,k,cP->ld)]=groupOff[j];
-       k++;
-      }
-    }
-  }
-  free(groupOff);
-}
-
-
-//Frees memory allocated in initCompPlan.
-void freeCompPlan(compPlan *cP){
-  free(cP->numGroups);
-  free(cP->groupCountX);
-  free(cP->groupOff);
-  free(cP->qToGroup);
-}
-
-
-void computeNNs(matrix dx, matrix dq, intMatrix xmap, compPlan cP, int *qIDs, int *NNs, int compLength){
-  compPlan dcP;
+void computeNNs(matrix dx, intMatrix dxMap, matrix dq, unint *dqMap, compPlan dcP, unint *NNs, unint compLength){
   real *dMins;
-  int *dMinIDs;
-
-  dcP.ld=cP.ld;
-  cudaMalloc((void**)&dcP.numGroups,cP.sNumGroups*sizeof(*dcP.numGroups));
-  cudaMalloc((void**)&dcP.groupCountX,cP.sGroupCountX*sizeof(*dcP.groupCountX));
-  cudaMalloc((void**)&dcP.groupOff,cP.sGroupOff*sizeof(*dcP.groupOff));
-  cudaMalloc((void**)&dcP.qToGroup,cP.sQToGroup*sizeof(*dcP.qToGroup));
-
-  cudaMemcpy(dcP.numGroups,cP.numGroups,cP.sNumGroups*sizeof(*dcP.numGroups),cudaMemcpyHostToDevice);
-  cudaMemcpy(dcP.groupCountX,cP.groupCountX,cP.sGroupCountX*sizeof(*dcP.groupCountX),cudaMemcpyHostToDevice);
-  cudaMemcpy(dcP.groupOff,cP.groupOff,cP.sGroupOff*sizeof(*dcP.groupOff),cudaMemcpyHostToDevice);
-  cudaMemcpy(dcP.qToGroup,cP.qToGroup,cP.sQToGroup*sizeof(*dcP.qToGroup),cudaMemcpyHostToDevice);
-  
-  cudaMalloc((void**)&dMins,compLength*sizeof(*dMins));
-  cudaMalloc((void**)&dMinIDs,compLength*sizeof(*dMinIDs));
-
-  int *dqIDs;
-  cudaMalloc( (void**)&dqIDs, compLength*sizeof(*dqIDs) );
-  cudaMemcpy( dqIDs, qIDs, compLength*sizeof(*dqIDs), cudaMemcpyHostToDevice );
-  
-
-  dim3 dimBlock(BLOCK_SIZE,BLOCK_SIZE);
-  dim3 dimGrid;
-  dimGrid.x = 1;
-  int numDone = 0;
-  while( numDone<compLength ){
-    int todo = MIN( (compLength-numDone) , MAX_BS*BLOCK_SIZE );
-
-    dimGrid.y = todo/BLOCK_SIZE;
-    planNNKernel<<<dimGrid,dimBlock>>>(dq,dx,dMins,dMinIDs,dcP,xmap,dqIDs,numDone);
-    cudaThreadSynchronize();
-    numDone += todo;
-  }
+  unint *dMinIDs;
+  
+  checkErr( cudaMalloc((void**)&dMins,compLength*sizeof(*dMins)) );
+  checkErr( cudaMalloc((void**)&dMinIDs,compLength*sizeof(*dMinIDs)) );
 
+  planNNWrap(dq, dqMap, dx, dxMap, dMins, dMinIDs, dcP, compLength);
   cudaMemcpy( NNs, dMinIDs, dq.r*sizeof(*NNs), cudaMemcpyDeviceToHost);
   
-  cudaFree(dcP.numGroups);
-  cudaFree(dcP.groupCountX);
-  cudaFree(dcP.groupOff);
-  cudaFree(dcP.qToGroup);
+      
   cudaFree(dMins);
   cudaFree(dMinIDs);
-  cudaFree(dqIDs);
 }
 
 
@@ -493,11 +279,85 @@ void computeNNs(matrix dx, matrix dq, intMatrix xmap, compPlan cP, int *qIDs, in
 //only distances from dr[start,:].. dr[start+length-1] to all of x
 //are computed, resulting in a distance matrix of size 
 //length by dx.pr.  It is assumed that length is padded.
-void distSubMat(matrix dr, matrix dx, matrix dD, int start, int length){
+void distSubMat(matrix dr, matrix dx, matrix dD, unint start, unint length){
   dr.r=dr.pr=length;
   dr.mat = &dr.mat[IDX( start, 0, dr.ld )];
   dist1Wrap(dr, dx, dD);
 }
 
 
+void destroyRBC(rbcStruct *rbcS){
+  cudaFree(rbcS->dx.mat);
+  cudaFree(rbcS->dxMap.mat);
+  cudaFree(rbcS->dr.mat);
+  free(rbcS->groupCount);
+}
+
+
+/* Danger: this function allocates memory that it does not free.  
+ * Use freeCompPlan to clear mem.  
+ * See the readme.txt file for a description of why this function is needed.
+ */
+void initCompPlan(compPlan *dcP, charMatrix cM, unint *groupCountQ, unint *groupCountX, unint numReps){
+  unint i,j,k;
+  unint maxNumGroups=0;
+  compPlan cP;
+  
+  unint sNumGroups = numReps;
+  cP.numGroups = (unint*)calloc(sNumGroups, sizeof(*cP.numGroups));
+  
+  for(i=0; i<numReps; i++){
+    cP.numGroups[i] = 0;
+    for(j=0; j<numReps; j++)
+      cP.numGroups[i] += cM.mat[IDX(i,j,cM.ld)];
+    maxNumGroups = MAX(cP.numGroups[i], maxNumGroups);
+  }
+  cP.ld = maxNumGroups;
+  
+  unint sQToQGroup;
+  for(i=0, sQToQGroup=0; i<numReps; i++)
+    sQToQGroup += PAD(groupCountQ[i]);
+  
+  cP.qToQGroup = (unint*)calloc( sQToQGroup, sizeof(*cP.qToQGroup) );
+
+  for(i=0, k=0; i<numReps; i++){
+    for(j=0; j<PAD(groupCountQ[i]); j++)
+      cP.qToQGroup[k++] = i;
+  }
+  
+  unint sQGroupToXGroup = numReps*maxNumGroups;
+  cP.qGroupToXGroup = (unint*)calloc( sQGroupToXGroup, sizeof(*cP.qGroupToXGroup) );
+  unint sGroupCountX = maxNumGroups*numReps;
+  cP.groupCountX = (unint*)calloc( sGroupCountX, sizeof(*cP.groupCountX) );
+  
+  for(i=0; i<numReps; i++){
+    for(j=0, k=0; j<numReps; j++){
+      if( cM.mat[IDX( i, j, cM.ld )] ){
+       cP.qGroupToXGroup[IDX( i, k, cP.ld )] = j;
+       cP.groupCountX[IDX( i, k++, cP.ld )] = groupCountX[j];
+      }
+    }
+  }
+
+  //Move to device
+  checkErr( cudaMalloc( (void**)&dcP->numGroups, sNumGroups*sizeof(*dcP->numGroups) ) );
+  cudaMemcpy( dcP->numGroups, cP.numGroups, sNumGroups*sizeof(*dcP->numGroups), cudaMemcpyHostToDevice );
+  checkErr( cudaMalloc( (void**)&dcP->groupCountX, sGroupCountX*sizeof(*dcP->groupCountX) ) );
+  cudaMemcpy( dcP->groupCountX, cP.groupCountX, sGroupCountX*sizeof(*dcP->groupCountX), cudaMemcpyHostToDevice );
+  checkErr( cudaMalloc( (void**)&dcP->qToQGroup, sQToQGroup*sizeof(*dcP->qToQGroup) ) );
+  cudaMemcpy( dcP->qToQGroup, cP.qToQGroup, sQToQGroup*sizeof(*dcP->qToQGroup), cudaMemcpyHostToDevice );
+  checkErr( cudaMalloc( (void**)&dcP->qGroupToXGroup, sQGroupToXGroup*sizeof(*dcP->qGroupToXGroup) ) );
+  cudaMemcpy( dcP->qGroupToXGroup, cP.qGroupToXGroup, sQGroupToXGroup*sizeof(*dcP->qGroupToXGroup), cudaMemcpyHostToDevice );
+  dcP->ld = cP.ld;
+}
+
+
+//Frees memory allocated in initCompPlan.
+void freeCompPlan(compPlan *dcP){
+  cudaFree(dcP->numGroups);
+  cudaFree(dcP->groupCountX);
+  cudaFree(dcP->qToQGroup);
+  cudaFree(dcP->qGroupToXGroup);
+}
+
 #endif
diff --git a/rbc.h b/rbc.h
index 0d0d517..bee5afb 100644 (file)
--- a/rbc.h
+++ b/rbc.h
@@ -1,26 +1,27 @@
+/* This file is part of the Random Ball Cover (RBC) library.
+ * (C) Copyright 2010, Lawrence Cayton [lcayton@tuebingen.mpg.de]
+ */
+
 #ifndef RBC_H
 #define RBC_H
 
 #include "defs.h"
 
-void rbc(matrix,matrix,int,int,int*);
-void build(const matrix,const matrix,intMatrix,int*,int);
-void distSubMat(matrix,matrix,matrix,int,int);
 
-void computeReps(matrix,matrix,int*,real*);
-void computeRadii(int*,real*,real*,int,int);
-void computeCounts(int*,int,int*);
-void computeOffsets(int*,int,int*);
-void groupPoints(matrix,int*,int*,int);
-void buildQMap(matrix,int*,int*,int,int*);
-void blockIntersection(charMatrix,matrix,real*,real*);
+void buildRBC(const matrix,rbcStruct*,unint, unint);
+void queryRBC(const matrix,const rbcStruct,unint*);
+void destroyRBC(rbcStruct*);
+
+void distSubMat(matrix,matrix,matrix,unint,unint);
+void computeReps(matrix,matrix,unint*,real*);
+void computeRadii(unint*,real*,real*,unint,unint);
+void computeCounts(unint*,unint,unint*);
+void buildQMap(matrix,unint*,unint*,unint,unint*);
 void idIntersection(charMatrix);
 void fullIntersection(charMatrix);
-void initCompPlan(compPlan*,charMatrix,int*,int*,int);
+void initCompPlan(compPlan*,charMatrix,unint*,unint*,unint);
 void freeCompPlan(compPlan*);
-void computeNNs(matrix,matrix,intMatrix, compPlan,int*,int*,int);
-
-
-
+void computeNNs(matrix,intMatrix,matrix,unint*,compPlan,unint*,unint);
+void setupReps(matrix,rbcStruct*,int);
 
 #endif
index bfbb216..0d0244f 100644 (file)
@@ -2,7 +2,7 @@
 Lawrence Cayton
 lcayton@tuebingen.mpg.de
 
-(C) Copyright 2010, Lawrence Cayton
+(C) Copyright 2010, Lawrence Cayton [lcayton@tuebingen.mpg.de]
  
 This program is free software: you can redistribute it and/or modify 
 it under the terms of the GNU General Public License as published by
@@ -18,31 +18,98 @@ You should have received a copy of the GNU General Public License
 along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
 ---------------------------------------------------------------------
+SUMMARY
 
 This is a C and CUDA implementation of the Random Ball Cover data 
-structure for nearest neighbor search.  
+structure for nearest neighbor search described in 
+
+L. Cayton, A nearest neighbor data structure for graphics hardware.
+ADMS, 2010.
+
+
+---------------------------------------------------------------------
+FILES
+
+* brute.{h,cu} -- implementation of brute force search (CPU and GPU
+  versions) 
+* defs.h -- definitions of constants and macros, including the
+  distance metric.
+* driver.cu -- example code for using the RBC data structure.
+* kernels.{h,cu} -- implementation of all the (device) kernel functions,
+  except those related to the scan (see sKernels below)
+* kernelWrap.{h,cu} -- CPU wrapper code around the kernels.
+* rbc.{h,cu} -- the core of the RBC data structure.  Includes the
+  implementation of build and search algorithms.
+* sKernel.{h,cu} -- implementation of the kernel functions related to
+  the parallel scan algorithm (used within the build method).
+* sKernelWrap.{h,cu} -- wrappers for the kernels in sKernel.
+* utils.{h,cu} -- misc utilities used in the code.
+* utilsGPU.{h,cu} -- misc utilities related to the GPU.
+
+
+---------------------------------------------------------------------
+COMPILATION
+
+Type make in a shell.  Requires GCC and NVCC (CUDA).  The code has
+been tested under GCC 4.4 and CUDA 3.1.
+
+
+---------------------------------------------------------------------
+USE
+
+To use the RBC data structure, you will likely need to integrate this
+code into your own.  The driver.cu file provides an example of how to
+use the RBC implementation.  To try it out, type
+>testRBC
+at the prompt and a list of options will be displayed.  Currently, the
+test program assumes that the input is a single binary file, which it
+then splits into queries and a the database randomly.  Clearly, such a
+setup is only useful for testing the performance of the data
+structure.  To use the data structure in a more useful fashion, you
+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.
+
+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.   
+
+
+---------------------------------------------------------------------
+MISC NOTES ON THE CODE
+
+* The code currently computes distance using the L_1 (manhattan)
+  metric.  If you wish to use a different notion of distance, you must
+  modify defs.h.  It is quite simple to switch to any metric that 
+  operates alongs the coordinates independently (eg, any L_p metric),
+  but more complex metrics will require some aditional work.  The L_2
+  metric (standard Euclidean distance) is already defined in defs.h.  
 
-Notes on the code:
 * The code requires that the entire DB and query set fit into the 
-device memory.  
+  device memory.  
 
 * For the most part, device variables (ie arrays residing on the GPU)
-begin with a lowercase d.  For example, the device version of the 
-DB variable x is dx.  
+  begin with a lowercase d.  For example, the device version of the 
+  DB variable x is dx.  
 
 * The computePlan code is a bit more complex than is needed for the 
-version of the RBC search algorithm described in the paper.  The 
-search algorithm described in the paper has two steps: (1) Find the 
-closest representative to the query.  (2) Explore the points owned
-by that representative (ie the s-closest points to the representative
-in the DB).  The computePlan code is more complex to make it easy
-to try out other options.  For example, one could search the points
-owned by the *two* closest representatives to the query instead.  This
-would require only minor changes to the code.
-
-* Currently the software works only in single precision.  If you wish
-to switch to double precision, you must edit the defs.h file.  Simply 
-uncomment the lines
+  version of the RBC search algorithm described in the paper.  The 
+  search algorithm described in the paper has two steps: (1) Find the 
+  closest representative to the query.  (2) Explore the points owned
+  by that representative (ie the s-closest points to the representative
+  in the DB).  The computePlan code is more complex to make it easy
+  to try out other options.  For example, one could search the points
+  owned by the *two* closest representatives to the query instead.  This
+  would require only minor changes to the code, though is currently 
+  untested.
+
+* Currently the software works in single precision.  If you wish to 
+  switch to double precision, you must edit the defs.h file.  Simply 
+  uncomment the lines
 
 typedef double real;
 #define MAX_REAL DBL_MAX
@@ -59,15 +126,23 @@ make clean
 followed by another make.
 
 * This software has been tested on the following graphics cards:
-NVIDIA GTX 285
-NVIDIA Tesla c2050.
+  NVIDIA GTX 285
+  NVIDIA Tesla c2050.
 
 * This sotware has been tested under the following software setup:
-Ubuntu 9.10 (linux)
-gcc 4.4
-cuda 3.1
-
-Please share your experience getting it to work under Windows and
-Mac OSX!
+  Ubuntu 9.10 (linux)
+  gcc 4.4
+  cuda 3.1
+
+  Please share your experience getting it to work under Windows and
+  Mac OSX!
+
+* If you are running this code on a GPU which is also driving your
+  display: A well-known issue with CUDA code in this situation is that 
+  a process within the operating system will automatically kill 
+  kernels that have been running for more than 5-10 seconds or so.
+  You can get around this in Linux by switching out of X-Windows (often 
+  CTRL-ALT-F1 does the trick) and running the code directly from the
+  terminal.
 
 
index a830f5f..f3bdbdc 100644 (file)
@@ -1,3 +1,7 @@
+/* This file is part of the Random Ball Cover (RBC) library.
+ * (C) Copyright 2010, Lawrence Cayton [lcayton@tuebingen.mpg.de]
+ */
+
 #ifndef SKERNEL_CU
 #define SKERNEL_CU
 
 #include "defs.h"
 #include "utils.h"
 
-__global__ void sumKernel(charMatrix in, intMatrix sum, intMatrix sumaux, int n){
-  int id = threadIdx.x;
-  int bo = blockIdx.x*SCAN_WIDTH; //block offset
-  int r = blockIdx.y;
-  int d, t;
+__global__ void sumKernel(charMatrix in, intMatrix sum, intMatrix sumaux, unint n){
+  unint id = threadIdx.x;
+  unint bo = blockIdx.x*SCAN_WIDTH; //block offset
+  unint r = blockIdx.y;
+  unint d, t;
   
-  const int l=SCAN_WIDTH; //length
+  const unint l=SCAN_WIDTH; //length
 
-  int off=1;
+  unint off=1;
 
-  __shared__ int ssum[l];
+  __shared__ unint ssum[l];
 
   ssum[2*id] = (bo+2*id < n) ? in.mat[IDX( r, bo+2*id, in.ld )] : 0;
   ssum[2*id+1] = (bo+2*id+1 < n) ? in.mat[IDX( r, bo+2*id+1, in.ld)] : 0;
@@ -57,122 +61,21 @@ __global__ void sumKernel(charMatrix in, intMatrix sum, intMatrix sumaux, int n)
     sum.mat[IDX( r, bo+2*id, sum.ld )] = ssum[2*id];
   if( bo+2*id+1 < n )
     sum.mat[IDX( r, bo+2*id+1, sum.ld )] = ssum[2*id+1];
-
-}
-
-__global__ void combineSumKernel(intMatrix sum, intMatrix daux, int n){
-  int id = threadIdx.x;
-  int bo = blockIdx.x * SCAN_WIDTH;
-  int r = blockIdx.y;
-  
-  if(bo+2*id < n)
-    sum.mat[IDX( r, bo+2*id, sum.ld )] += daux.mat[IDX( r, blockIdx.x, daux.ld )];
-  if(bo+2*id+1 < n)
-    sum.mat[IDX( r, bo+2*id+1, sum.ld )] += daux.mat[IDX( r, blockIdx.x, daux.ld )];
-  
-}
-
-
-__global__ void getCountsKernel(int *counts, charMatrix ir, intMatrix sums){
-  int r = blockIdx.x*BLOCK_SIZE + threadIdx.x;
-  if ( r < ir.r ){
-    counts[r] = ir.mat[IDX( r, ir.c-1, ir.ld )] ? sums.mat[IDX( r, sums.c-1, sums.ld )]+1 : sums.mat[IDX( r, sums.c-1, sums.ld )];
-  }
-}
-
-
-__global__ void buildMapKernel(intMatrix map, charMatrix ir, intMatrix sums, int offSet){
-  int id = threadIdx.x;
-  int bo = blockIdx.x * SCAN_WIDTH;
-  int r = blockIdx.y;
-
-  if(bo+2*id < ir.c && ir.mat[IDX( r, bo+2*id, ir.ld )])
-    map.mat[IDX( r+offSet, sums.mat[IDX( r, bo+2*id, sums.ld )], map.ld)] = bo+2*id;
-  if(bo+2*id+1 < ir.c && ir.mat[IDX( r, bo+2*id+1, ir.ld )])
-    map.mat[IDX( r+offSet, sums.mat[IDX( r, bo+2*id+1, sums.ld )], map.ld)] = bo+2*id+1;
-}
-
-
-void getCountsWrap(int *counts, charMatrix ir, intMatrix sums){
-  dim3 block(BLOCK_SIZE,1);
-  dim3 grid(ir.pr/BLOCK_SIZE,1);
-  getCountsKernel<<<grid,block>>>(counts, ir, sums);
-}
-
-
-void buildMapWrap(intMatrix map, charMatrix ir, intMatrix sums, int offSet){
-  int numScans = (ir.c+SCAN_WIDTH-1)/SCAN_WIDTH;
-  dim3 block( SCAN_WIDTH/2, 1 );
-  dim3 grid( numScans, ir.r );
-  buildMapKernel<<<grid,block>>>(map, ir, sums, offSet);
-}
-
-
-void sumWrap(charMatrix in, intMatrix sum){
-  int i;
-  int n = in.c;
-  int numScans = (n+SCAN_WIDTH-1)/SCAN_WIDTH;
-  
-
-  int depth = ceil( log(n) / log(SCAN_WIDTH) ) -1 ;
-  int *width = (int*)calloc( depth+1, sizeof(*width) );
-    
-  intMatrix *dAux;
-  dAux = (intMatrix*)calloc( depth+1, sizeof(*dAux) );
-  
-  dAux[0].r=dAux[0].pr=in.r; dAux[0].c=dAux[0].pc=dAux[0].ld=numScans;
-  cudaMalloc( (void**)&dAux[0].mat, dAux[0].pr*dAux[0].pc*sizeof(*dAux[0].mat) );
-   
-  dim3 block( SCAN_WIDTH/2, 1 );
-  dim3 grid( numScans, in.r );
-  
-  sumKernel<<<grid,block>>>(in, sum, dAux[0], n);
-  cudaThreadSynchronize();
-  
-  width[0]=numScans; //Clean up, this is ugly (necc b/c loop might not be entered)
-  for( i=0; i<depth; i++ ){
-    width[i] = numScans;
-    numScans = (numScans+SCAN_WIDTH-1)/SCAN_WIDTH;
-    
-    dAux[i+1].r=dAux[i+1].pr=in.r; dAux[i+1].c=dAux[i+1].pc=dAux[i+1].ld=numScans;
-    cudaMalloc( (void**)&dAux[i+1].mat, dAux[i+1].pr*dAux[i+1].pc*sizeof(*dAux[i+1].mat) );
-        
-    grid.x = numScans;
-    sumKernelI<<<grid,block>>>(dAux[i], dAux[i], dAux[i+1], width[i]);
-    cudaThreadSynchronize();
-  }
-
-    
-  for( i=depth-1; i>0; i-- ){
-    grid.x = width[i];
-    combineSumKernel<<<grid,block>>>(dAux[i-1], dAux[i], width[i-1]);
-    cudaThreadSynchronize();
-  }
-  
-  grid.x = width[0];
-  combineSumKernel<<<grid,block>>>(sum, dAux[0], n);
-  cudaThreadSynchronize();
-  
-  for( i=0; i <= depth; i++)
-    cudaFree(dAux[i].mat);
-  
-  free(dAux);
-  free(width);
 }
 
 
 //This is the same as sumKernel, but takes an int matrix as input.
-__global__ void sumKernelI(intMatrix in, intMatrix sum, intMatrix sumaux, int n){
-  int id = threadIdx.x;
-  int bo = blockIdx.x*SCAN_WIDTH; //block offset
-  int r = blockIdx.y;
-  int d, t;
+__global__ void sumKernelI(intMatrix in, intMatrix sum, intMatrix sumaux, unint n){
+  unint id = threadIdx.x;
+  unint bo = blockIdx.x*SCAN_WIDTH; //block offset
+  unint r = blockIdx.y;
+  unint d, t;
   
-  const int l=SCAN_WIDTH; //length
+  const unint l=SCAN_WIDTH; //length
 
-  int off=1;
+  unint off=1;
 
-  __shared__ int ssum[l];
+  __shared__ unint ssum[l];
 
   ssum[2*id] = (bo+2*id < n) ? in.mat[IDX( r, bo+2*id, in.ld )] : 0;
   ssum[2*id+1] = (bo+2*id+1 < n) ? in.mat[IDX( r, bo+2*id+1, in.ld)] : 0;
@@ -194,7 +97,6 @@ __global__ void sumKernelI(intMatrix in, intMatrix sum, intMatrix sumaux, int n)
     ssum[ l-1 ] = 0;
   }
   
-  
   //down-sweep
   for ( d=1; d<l; d*=2 ){
     off >>= 1;
@@ -205,7 +107,6 @@ __global__ void sumKernelI(intMatrix in, intMatrix sum, intMatrix sumaux, int n)
       ssum[ off*(2*id+1)-1 ] = ssum[ off*(2*id+2)-1 ];
       ssum[ off*(2*id+2)-1 ] += t;
     }
-    
   }
   
   __syncthreads();
@@ -215,8 +116,41 @@ __global__ void sumKernelI(intMatrix in, intMatrix sum, intMatrix sumaux, int n)
   
   if( bo+2*id+1 < n )
     sum.mat[IDX( r, bo+2*id+1, sum.ld )] = ssum[2*id+1];
+}
+
+
+
+__global__ void combineSumKernel(intMatrix sum, unint numDone, intMatrix daux, unint n){
+  unint id = threadIdx.x;
+  unint bo = blockIdx.x * SCAN_WIDTH;
+  unint r = blockIdx.y+numDone;
   
+  if(bo+2*id < n)
+    sum.mat[IDX( r, bo+2*id, sum.ld )] += daux.mat[IDX( r, blockIdx.x, daux.ld )];
+  if(bo+2*id+1 < n)
+    sum.mat[IDX( r, bo+2*id+1, sum.ld )] += daux.mat[IDX( r, blockIdx.x, daux.ld )];
   
 }
 
+
+__global__ void getCountsKernel(unint *counts, unint numDone, charMatrix ir, intMatrix sums){
+  unint r = blockIdx.x*BLOCK_SIZE + threadIdx.x + numDone;
+  if ( r < ir.r ){
+    counts[r] = ir.mat[IDX( r, ir.c-1, ir.ld )] ? sums.mat[IDX( r, sums.c-1, sums.ld )]+1 : sums.mat[IDX( r, sums.c-1, sums.ld )];
+  }
+}
+
+
+__global__ void buildMapKernel(intMatrix map, charMatrix ir, intMatrix sums, unint offSet){
+  unint id = threadIdx.x;
+  unint bo = blockIdx.x * SCAN_WIDTH;
+  unint r = blockIdx.y;
+
+  if(bo+2*id < ir.c && ir.mat[IDX( r, bo+2*id, ir.ld )])
+    map.mat[IDX( r+offSet, sums.mat[IDX( r, bo+2*id, sums.ld )], map.ld)] = bo+2*id;
+  if(bo+2*id+1 < ir.c && ir.mat[IDX( r, bo+2*id+1, ir.ld )])
+    map.mat[IDX( r+offSet, sums.mat[IDX( r, bo+2*id+1, sums.ld )], map.ld)] = bo+2*id+1;
+}
+
+
 #endif
index b38c123..5ede0ab 100644 (file)
--- a/sKernel.h
+++ b/sKernel.h
@@ -1,14 +1,14 @@
+/* This file is part of the Random Ball Cover (RBC) library.
+ * (C) Copyright 2010, Lawrence Cayton [lcayton@tuebingen.mpg.de]
+ */
 #ifndef SKERNEL_H
 #define SKERNEL_H
 
 #include "defs.h"
-__global__ void sumKernelI(intMatrix,intMatrix,intMatrix,int);
-__global__ void sumKernel(charMatrix,intMatrix,intMatrix,int);
-__global__ void combineSumKernel(intMatrix,intMatrix,int);
-__global__ void getCountsKernel(int*,charMatrix,intMatrix);
-void getCountsWrap(int*,charMatrix,intMatrix);
-void buildMapWrap(intMatrix,charMatrix,intMatrix,int);
-__global__ void buildMapKernel(intMatrix,charMatrix,intMatrix,int);
-void sumWrap(charMatrix,intMatrix);
+__global__ void sumKernelI(intMatrix,intMatrix,intMatrix,unint);
+__global__ void sumKernel(charMatrix,intMatrix,intMatrix,unint);
+__global__ void combineSumKernel(intMatrix,unint,intMatrix,unint);
+__global__ void getCountsKernel(unint*,unint,charMatrix,intMatrix);
+__global__ void buildMapKernel(intMatrix,charMatrix,intMatrix,unint);
 
 #endif
diff --git a/sKernelWrap.cu b/sKernelWrap.cu
new file mode 100644 (file)
index 0000000..25d338a
--- /dev/null
@@ -0,0 +1,107 @@
+/* This file is part of the Random Ball Cover (RBC) library.
+ * (C) Copyright 2010, Lawrence Cayton [lcayton@tuebingen.mpg.de]
+ */
+#ifndef SKERNELWRAP_CU
+#define SKERNELWRAP_CU
+
+#include "sKernel.h"
+#include<cuda.h>
+#include "defs.h"
+#include "utilsGPU.h"
+#include<stdio.h>
+
+void getCountsWrap(unint *counts, charMatrix ir, intMatrix sums){
+  dim3 block(BLOCK_SIZE,1);
+  dim3 grid;
+  grid.y=1;
+  unint todo, numDone;
+  
+  numDone = 0;
+  while(numDone < ir.pr){
+    todo = MIN( ir.pr - numDone, MAX_BS*BLOCK_SIZE );
+    grid.x = todo/BLOCK_SIZE;
+    getCountsKernel<<<grid,block>>>(counts, numDone, ir, sums);
+    numDone += todo;
+  }
+}
+
+
+void buildMapWrap(intMatrix map, charMatrix ir, intMatrix sums, unint offSet){
+  unint numScans = (ir.c+SCAN_WIDTH-1)/SCAN_WIDTH;
+  dim3 block( SCAN_WIDTH/2, 1 );
+  dim3 grid;
+  unint todo, numDone;
+
+  grid.x = numScans;
+  numDone = 0;
+  while( numDone < ir.r ){
+    todo = MIN( ir.r-numDone, MAX_BS );
+    grid.y = todo;
+    buildMapKernel<<<grid,block>>>(map, ir, sums, offSet+numDone);
+    numDone += todo;
+  }
+}
+
+
+void sumWrap(charMatrix in, intMatrix sum){
+  int i; 
+  unint todo, numDone, temp;
+  unint n = in.c;
+  unint numScans = (n+SCAN_WIDTH-1)/SCAN_WIDTH;
+  unint depth = ceil( log(n) / log(SCAN_WIDTH) ) -1 ;
+  unint *width = (unint*)calloc( depth+1, sizeof(*width) );
+    
+  intMatrix *dAux;
+  dAux = (intMatrix*)calloc( depth+1, sizeof(*dAux) );
+  
+  for( i=0, temp=n; i<=depth; i++){
+    temp = (temp+SCAN_WIDTH-1)/SCAN_WIDTH;
+    dAux[i].r=dAux[i].pr=in.r; dAux[i].c=dAux[i].pc=dAux[i].ld=temp;
+    checkErr( cudaMalloc( (void**)&dAux[i].mat, dAux[i].pr*dAux[i].pc*sizeof(*dAux[i].mat) ) );
+  }
+
+  dim3 block( SCAN_WIDTH/2, 1 );
+  dim3 grid;
+  
+  numDone=0;
+  while( numDone < in.r ){
+    todo = MIN( in.r - numDone, MAX_BS );
+    numScans = (n+SCAN_WIDTH-1)/SCAN_WIDTH;
+    dAux[0].r=dAux[0].pr=todo;
+    grid.x = numScans;
+    grid.y = todo;
+    sumKernel<<<grid,block>>>(in, sum, dAux[0], n);
+    cudaThreadSynchronize();
+    
+    width[0] = numScans; // Necessary because following loop might not be entered
+    for( i=0; i<depth; i++ ){
+      width[i] = numScans;
+      numScans = (numScans+SCAN_WIDTH-1)/SCAN_WIDTH;
+      dAux[i+1].r=dAux[i+1].pr=todo;
+      
+      grid.x = numScans;
+      sumKernelI<<<grid,block>>>(dAux[i], dAux[i], dAux[i+1], width[i]);
+      cudaThreadSynchronize();
+    }
+  
+    for( i=depth-1; i>0; i-- ){
+      grid.x = width[i];
+      combineSumKernel<<<grid,block>>>(dAux[i-1], numDone, dAux[i], width[i-1]);
+      cudaThreadSynchronize();
+    }
+    
+    grid.x = width[0];
+    combineSumKernel<<<grid,block>>>(sum, numDone, dAux[0], n);
+    cudaThreadSynchronize();
+    
+    numDone += todo;
+  }
+
+  for( i=0; i<=depth; i++)
+   cudaFree(dAux[i].mat);
+  free(dAux);
+  free(width);
+}
+
+
+#endif
diff --git a/sKernelWrap.h b/sKernelWrap.h
new file mode 100644 (file)
index 0000000..275e725
--- /dev/null
@@ -0,0 +1,15 @@
+/* This file is part of the Random Ball Cover (RBC) library.
+ * (C) Copyright 2010, Lawrence Cayton [lcayton@tuebingen.mpg.de]
+ */
+
+#ifndef SKERNELWRAP_H
+#define SKERNELWRAP_H
+
+#include "defs.h"
+
+void getCountsWrap(unint*,charMatrix,intMatrix);
+void buildMapWrap(intMatrix,charMatrix,intMatrix,unint);
+void sumWrap(charMatrix,intMatrix);
+
+
+#endif
index 495d674..2f075dc 100644 (file)
--- a/utils.cu
+++ b/utils.cu
@@ -1,3 +1,7 @@
+/* This file is part of the Random Ball Cover (RBC) library.
+ * (C) Copyright 2010, Lawrence Cayton [lcayton@tuebingen.mpg.de]
+ */
+
 #ifndef UTILS_CU
 #define UTILS_CU
 
@@ -5,13 +9,14 @@
 #include<stdio.h>
 #include "utils.h"
 #include "defs.h"
+#include "kernels.h"
 
 //Returns a length l subset of a random permutation of [0,...,n-1]
 //using the knuth shuffle.
 //Input variable x is assumed to be alloced and of size l.
-void subRandPerm(int l, int n, int *x){
-  int i,ri, *y;
-  y = (int*)calloc(n,sizeof(*y));
+void subRandPerm(unint l, unint n, unint *x){
+  unint i,ri, *y;
+  y = (unint*)calloc(n,sizeof(*y));
     
   struct timeval t3;
   gettimeofday(&t3,NULL);
@@ -33,8 +38,8 @@ void subRandPerm(int l, int n, int *x){
 
 //Generates a random permutation of 0, ... , n-1 using the knuth shuffle.
 //This should probably be merged with subRandPerm. 
-void randPerm(int n, int *x){
-  int i,ri;
+void randPerm(unint n, unint *x){
+  unint i,ri;
   
   struct timeval t3;
   gettimeofday(&t3,NULL);
@@ -50,14 +55,14 @@ void randPerm(int n, int *x){
   }
 }
 
-void swap(int *a, int *b){
-  int t;
+void swap(unint *a, unint *b){
+  unint t;
   t=*a; *a=*b; *b=t;
 }
 
 //generates a rand int in rand [a,b) 
-int randBetween(int a, int b){
-  int val,c;
+unint randBetween(unint a, unint b){
+  unint val,c;
 
   if(b<=a){
     fprintf(stderr,"misuse of randBetween.. exiting\n");
@@ -73,7 +78,7 @@ int randBetween(int a, int b){
 
 
 void printMat(matrix A){
-  int i,j;
+  unint i,j;
   for(i=0;i<A.r;i++){
     for(j=0;j<A.c;j++)
       printf("%6.4f ",(float)A.mat[IDX(i,j,A.ld)]);
@@ -82,8 +87,8 @@ void printMat(matrix A){
 }
 
 
-void printMatWithIDs(matrix A, int *id){
-  int i,j;
+void printMatWithIDs(matrix A, unint *id){
+  unint i,j;
   for(i=0;i<A.r;i++){
     for(j=0;j<A.c;j++)
       printf("%6.4f ",(float)A.mat[IDX(i,j,A.ld)]);
@@ -94,7 +99,7 @@ void printMatWithIDs(matrix A, int *id){
 
 
 void printCharMat(charMatrix A){
-  int i,j;
+  unint i,j;
   for(i=0;i<A.r;i++){
     for(j=0;j<A.c;j++)
       printf("%d ",(char)A.mat[IDX(i,j,A.ld)]);
@@ -103,16 +108,16 @@ void printCharMat(charMatrix A){
 }
 
 void printIntMat(intMatrix A){
-  int i,j;
+  unint i,j;
   for(i=0;i<A.r;i++){
     for(j=0;j<A.c;j++)
-      printf("%d ",(int)A.mat[IDX(i,j,A.ld)]);
+      printf("%d ",(unint)A.mat[IDX(i,j,A.ld)]);
     printf("\n");
   }
 }
 
-void printVector(real *x, int d){
-  int i;
+void printVector(real *x, unint d){
+  unint i;
 
   for(i=0 ; i<d; i++)
     printf("%6.2f ",x[i]);
@@ -120,8 +125,8 @@ void printVector(real *x, int d){
 }
 
 
-void copyVector(real *x, real *y, int d){
-  int i;
+void copyVector(real *x, real *y, unint d){
+  unint i;
   
   for(i=0;i<d;i++)
     x[i]=y[i];
@@ -129,7 +134,7 @@ void copyVector(real *x, real *y, int d){
 
 
 void copyMat(matrix *x, matrix *y){
-  int i,j;
+  unint i,j;
   
   x->r=y->r; x->pr=y->pr; x->c=y->c; x->pc=y->pc; x->ld=y->ld;
   for(i=0; i<y->r; i++){
@@ -140,12 +145,13 @@ void copyMat(matrix *x, matrix *y){
 }
 
 
-real distL1(matrix x, matrix y, int k, int l){
-  int i;
+real distVec(matrix x, matrix y, unint k, unint l){
+  unint i;
   real ans=0;
   
   for(i=0;i<x.c;i++)
-    ans+=abs(x.mat[IDX(k,i,x.ld)]-y.mat[IDX(l,i,x.ld)]);
+    ans += DIST( x.mat[IDX(k,i,x.ld)], y.mat[IDX(l,i,x.ld)] );
+    //ans+=fabs(x.mat[IDX(k,i,x.ld)]-y.mat[IDX(l,i,x.ld)]);
   return ans;
 }
 
diff --git a/utils.h b/utils.h
index 9f9edf1..06f3bda 100644 (file)
--- a/utils.h
+++ b/utils.h
@@ -1,19 +1,23 @@
+/* This file is part of the Random Ball Cover (RBC) library.
+ * (C) Copyright 2010, Lawrence Cayton [lcayton@tuebingen.mpg.de]
+ */
+
 #ifndef UTILS_H
 #define UTILS_H
 
 #include "defs.h"
 
-void swap(int*,int*);
-void randPerm(int,int*);
-void subRandPerm(int,int,int*);
-int randBetween(int,int);
+void swap(unint*,unint*);
+void randPerm(unint,unint*);
+void subRandPerm(unint,unint,unint*);
+unint randBetween(unint,unint);
 void printMat(matrix);
-void printMatWithIDs(matrix,int*);
+void printMatWithIDs(matrix,unint*);
 void printCharMat(charMatrix);
 void printIntMat(intMatrix);
-void printVector(real*,int);
-void copyVector(real*,real*,int);
-real distL1(matrix,matrix,int,int);
+void printVector(real*,unint);
+void copyVector(real*,real*,unint);
+real distVec(matrix,matrix,unint,unint);
 double timeDiff(struct timeval,struct timeval);
 void copyMat(matrix*,matrix*);
 #endif
index a28409e..c45e1d7 100644 (file)
@@ -1,26 +1,14 @@
+/* This file is part of the Random Ball Cover (RBC) library.
+ * (C) Copyright 2010, Lawrence Cayton [lcayton@tuebingen.mpg.de]
+ */
+
 #ifndef UTILSGPU_CU
 #define UTILSGPU_CU
 
 #include<cuda.h>
 #include<stdio.h>
 #include "defs.h"
-
-memPlan createMemPlan(int nPts, int memPerPt){
-  memPlan mp;
-  unsigned int memFree, memTot;
-  int ptsAtOnce;
-
-  cuMemGetInfo(&memFree, &memTot);
-  memFree = (int)(((float)memFree)*MEM_USABLE);
-  printf("memfree = %d \n",memFree);
-  ptsAtOnce = DPAD(memFree/memPerPt); //max number of pts that can be processed at once
-  printf("ptsAtOnce = %d \n",ptsAtOnce);
-  mp.numComputeSegs = nPts/ptsAtOnce + ((nPts%ptsAtOnce==0) ? 0 : 1);
-  mp.normSegSize=PAD(nPts/mp.numComputeSegs); 
-  mp.lastSegSize=PAD(nPts) - mp.normSegSize*(mp.numComputeSegs-1);
-  //Note that lastSegSize is automatically padded if nPts is.
-  return mp;
-}
+#include "utilsGPU.h"
 
 void copyAndMove(matrix *dx, const matrix *x){
   dx->r = x->r; 
@@ -29,7 +17,7 @@ void copyAndMove(matrix *dx, const matrix *x){
   dx->pc = x->pc;
   dx->ld = x->ld;
 
-  cudaMalloc( (void**)&(dx->mat), dx->pr*dx->pc*sizeof(*(dx->mat)) );
+  checkErr( cudaMalloc( (void**)&(dx->mat), dx->pr*dx->pc*sizeof(*(dx->mat)) ) );
   cudaMemcpy( dx->mat, x->mat, dx->pr*dx->pc*sizeof(*(dx->mat)), cudaMemcpyHostToDevice );
   
 }
@@ -42,7 +30,7 @@ void copyAndMoveI(intMatrix *dx, const intMatrix *x){
   dx->pc = x->pc;
   dx->ld = x->ld;
 
-  cudaMalloc( (void**)&(dx->mat), dx->pr*dx->pc*sizeof(*(dx->mat)) );
+  checkErr( cudaMalloc( (void**)&(dx->mat), dx->pr*dx->pc*sizeof(*(dx->mat)) ) );
   cudaMemcpy( dx->mat, x->mat, dx->pr*dx->pc*sizeof(*(dx->mat)), cudaMemcpyHostToDevice );
   
 }
@@ -55,8 +43,24 @@ void copyAndMoveC(charMatrix *dx, const charMatrix *x){
   dx->pc = x->pc;
   dx->ld = x->ld;
 
-  cudaMalloc( (void**)&(dx->mat), dx->pr*dx->pc*sizeof(*(dx->mat)) );
+  checkErr( cudaMalloc( (void**)&(dx->mat), dx->pr*dx->pc*sizeof(*(dx->mat)) ) );
   cudaMemcpy( dx->mat, x->mat, dx->pr*dx->pc*sizeof(*(dx->mat)), cudaMemcpyHostToDevice );
   
 }
+
+
+void checkErr(cudaError_t cError){
+  if(cError != cudaSuccess){
+    fprintf(stderr,"GPU-related error:\n\t%s \n", cudaGetErrorString(cError) );
+    fprintf(stderr,"exiting ..\n");
+    exit(1);
+  }
+  
+}
+
+void checkErr(char* loc, cudaError_t cError){
+  printf("in %s:\n",loc);
+  checkErr(cError);
+}
+
 #endif
index cba437c..9c176a6 100644 (file)
@@ -1,15 +1,18 @@
+/* This file is part of the Random Ball Cover (RBC) library.
+ * (C) Copyright 2010, Lawrence Cayton [lcayton@tuebingen.mpg.de]
+ */
+
 #ifndef UTILSGPU_H
 #define UTILSGPU_H
 
 #include "defs.h"
+#include<cuda.h>
 
-memPlan createMemPlan(int,int);
 
 void copyAndMove(matrix*,const matrix*);
 void copyAndMoveI(intMatrix*,const intMatrix*);
 void copyAndMoveC(charMatrix*,const charMatrix*);
 
-
-
-
+void checkErr(cudaError_t);
+void checkErr(char*,cudaError_t );
 #endif