Stable release of RBC
authorLawrence Cayton <lcayton@tuebingen.mpg.de>
Wed, 11 Aug 2010 15:27:28 +0000 (17:27 +0200)
committerLawrence Cayton <lcayton@tuebingen.mpg.de>
Wed, 11 Aug 2010 15:27:28 +0000 (17:27 +0200)
19 files changed:
Makefile [new file with mode: 0644]
brute.cu [new file with mode: 0644]
brute.h [new file with mode: 0644]
defs.h [new file with mode: 0644]
deprecated.cu [new file with mode: 0644]
driver.cu [new file with mode: 0644]
kernelWrap.cu [new file with mode: 0644]
kernelWrap.h [new file with mode: 0644]
kernels.cu [new file with mode: 0644]
kernels.h [new file with mode: 0644]
oldKernels.cu [new file with mode: 0644]
rbc.cu [new file with mode: 0644]
rbc.h [new file with mode: 0644]
sKernel.cu [new file with mode: 0644]
sKernel.h [new file with mode: 0644]
utils.cu [new file with mode: 0644]
utils.h [new file with mode: 0644]
utilsGPU.cu [new file with mode: 0644]
utilsGPU.h [new file with mode: 0644]

diff --git a/Makefile b/Makefile
new file mode 100644 (file)
index 0000000..6d58877
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,26 @@
+CC=gcc
+NVCC=nvcc
+CCFLAGS=
+NVCCFLAGS= --ptxas-options=-v
+#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
+OBJECTS=$(SOURCES:.c=.o)
+CUOBJECTS=$(CUSOURCES:.cu=.o)
+EXECUTABLE=testRBC
+all: $(SOURCES) $(CUSOURCES) $(EXECUTABLE)
+
+$(EXECUTABLE): $(OBJECTS) $(CUOBJECTS)
+       $(NVCC) $(NVCCFLAGS) $(OBJECTS) $(CUOBJECTS) -o $@ $(LINKFLAGS)
+
+%.o:%.c
+       $(NVCC) $(NVCCFLAGS) -c $+ 
+
+%.o:%.cu
+       $(NVCC) $(NVCCFLAGS) -c $+
+
+clean:
+       rm -f *.o
+       rm -f $(EXECUTABLE)
diff --git a/brute.cu b/brute.cu
new file mode 100644 (file)
index 0000000..388ee26
--- /dev/null
+++ b/brute.cu
@@ -0,0 +1,88 @@
+#ifndef BRUTE_CU
+#define BRUTE_CU
+
+#include "utilsGPU.h"
+#include "utils.h"
+#include "rbc.h"
+#include "defs.h"
+#include "kernels.h"
+#include "kernelWrap.h"
+#include "brute.h"
+#include<stdio.h>
+#include<cuda.h>
+
+void bruteRangeCount(matrix x, matrix q, real *ranges, int *cnts){
+  matrix dx, dq;
+  real* dranges;
+  int *dcnts;
+  
+  copyAndMove(&dx, &x);
+  copyAndMove(&dq, &q);
+
+  cudaMalloc( (void**)&dranges, q.pr*sizeof(*dranges) );
+  cudaMemcpy( dranges, ranges, q.r*sizeof(*dranges), cudaMemcpyHostToDevice );
+
+  cudaMalloc( (void**)&dcnts, q.pr*sizeof(*dcnts) );
+  
+  rangeCountWrap(dq, dx, dranges, dcnts);
+  
+  cudaMemcpy(cnts, dcnts, q.r*sizeof(*cnts), cudaMemcpyDeviceToHost );
+
+  cudaFree(dcnts);
+  cudaFree(dranges);
+  cudaFree(dx.mat);
+  cudaFree(dq.mat);
+}
+
+
+void bruteSearch(matrix x, matrix q, int *NNs){
+  real *dMins;
+  int *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));
+
+  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);
+  
+  nnWrap(dq,dx,dMins,dMinIDs);
+
+  cudaMemcpy(NNs,dMinIDs,dq.r*sizeof(*NNs),cudaMemcpyDeviceToHost);
+  
+  cudaFree(dMins);
+  cudaFree(dMinIDs);
+  cudaFree(dx.mat);
+  cudaFree(dq.mat);
+
+}
+
+void bruteCPU(matrix X, matrix Q, int *NNs){
+  real *dtoNNs; 
+  real temp;
+
+  int i, j;
+
+  dtoNNs = (real*)calloc(Q.r,sizeof(*dtoNNs));
+  
+  for( i=0; i<Q.r; i++ ){
+    dtoNNs[i] = MAX_REAL;
+    NNs[i] = 0;
+    for(j=0; j<X.r; j++ ){
+      temp = distL1( Q, X, i, j );
+      if( temp < dtoNNs[i]){
+       NNs[i] = j;
+       dtoNNs[i] = temp;
+      }
+    }
+  }
+  
+  free(dtoNNs);  
+}
+#endif
diff --git a/brute.h b/brute.h
new file mode 100644 (file)
index 0000000..534cded
--- /dev/null
+++ b/brute.h
@@ -0,0 +1,10 @@
+#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*);
+
+#endif
diff --git a/defs.h b/defs.h
new file mode 100644 (file)
index 0000000..a517024
--- /dev/null
+++ b/defs.h
@@ -0,0 +1,89 @@
+#ifndef DEFS_H
+#define DEFS_H
+
+#include<float.h>
+
+#define FLOAT_TOL 1e-7
+#define BLOCK_SIZE 16 //must be a power of 2
+
+#define MAX_BS 65535 //max block size
+#define SCAN_WIDTH 1024 
+
+// Format that the data is manipulated in:
+typedef float real;
+#define MAX_REAL FLT_MAX
+
+// To switch to double precision, comment out the above 
+// 2 lines and uncomment the following two lines. 
+
+//typedef double real;
+//#define MAX_REAL DBL_MAX
+
+//Percentage of device mem to use
+#define MEM_USABLE .95
+
+#define DUMMY_IDX INT_MAX
+
+//Row major indexing
+#define IDX(i,j,ld) (((i)*(ld))+(j))
+
+//increase an int to the next multiple of BLOCK_SIZE
+#define PAD(i) ( ((i)%BLOCK_SIZE)==0 ? (i):((i)/BLOCK_SIZE)*BLOCK_SIZE+BLOCK_SIZE ) 
+
+//decrease an int to the next multiple of BLOCK_SIZE
+#define DPAD(i) ( ((i)%BLOCK_SIZE)==0 ? (i):((i)/BLOCK_SIZE)*BLOCK_SIZE ) 
+
+#define MAX(i,j) ((i) > (j) ? (i) : (j))
+
+#define MIN(i,j) ((i) < (j) ? (i) : (j))
+
+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)
+} matrix;
+
+
+typedef struct {
+  char *mat;
+  int r;
+  int c;
+  int pr;
+  int pc;
+  int 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;
+
+
+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)
+} compPlan;
+#endif
diff --git a/deprecated.cu b/deprecated.cu
new file mode 100644 (file)
index 0000000..3dfc315
--- /dev/null
@@ -0,0 +1,316 @@
+/* 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 
+use a simpler version of computeReps. */
+
+//Assumes that dr is a matrix already on the device
+void computeReps(matrix x, matrix dr, int *repIDs, real *distToReps){
+  int memPerX, segSize; //seg==segment
+  int index, tempSize; //temp variables used in the loop
+  int i;
+  matrix dx;
+  real *dMins;
+  int *dMinIDs;
+  memPlan mp;
+
+  int n = x.r; //For convenience
+  
+  // Items that need to go on device: x, repIDs, distToReps.  The "+1" is for the
+  // distance from each point to its nearest rep (distToReps) and the int is for
+  // the ID (repIDs).
+  memPerX = (x.pc+1)*sizeof(real)+sizeof(int);
+  mp = createMemPlan(x.r,memPerX);
+  
+  for(i=0;i<mp.numComputeSegs;i++){
+    if(i==mp.numComputeSegs-1)
+      segSize = mp.lastSegSize;
+    else
+      segSize = mp.normSegSize;
+
+    //Allocate & copy over data
+    index = IDX(mp.normSegSize*i,0,x.ld);
+    tempSize = segSize*x.pc*sizeof(*(dx.mat));
+
+    cudaMalloc((void**)&(dx.mat),tempSize);
+    cudaMemcpy(dx.mat,&(x.mat[index]),tempSize,cudaMemcpyHostToDevice);
+    dx.r=segSize; dx.c=x.c; dx.pr=dx.r; dx.pc=x.pc; dx.ld=x.ld;
+
+    //Allocate matrices to temporarily store mins and IDs (NOTE:MOVE OUT OF LOOP FOR EFFICIENCY)
+    cudaMalloc((void**)&(dMins), PAD(MIN(segSize,n))*sizeof(*dMins));
+    cudaMalloc((void**)&(dMinIDs), PAD(MIN(segSize,n))*sizeof(*dMinIDs));
+    nnWrap(dx,dr,dMins,dMinIDs);
+
+    cudaMemcpy(&distToReps[i*segSize],dMins,MIN(segSize,n)*sizeof(*dMins),cudaMemcpyDeviceToHost);
+    cudaMemcpy(&repIDs[i*segSize],dMinIDs,MIN(segSize,n)*sizeof(*dMinIDs),cudaMemcpyDeviceToHost);
+    
+    cudaFree(dMins);
+    cudaFree(dMinIDs);
+    cudaFree(dx.mat);
+  }
+}
+
+
+__global__ void getMinsKernel(matrix,real*,int*);
+
+// Returns the min of each row of D.  dMins and dMinIDs 
+// are assumed to be (at least) of size D.r.
+__global__ void getMinsKernel(const matrix D, real *dMins, int *dMinIDs){
+  int row, locRow, colOff, i, curCol;
+  real temp;
+
+  row = blockIdx.y*BLOCK_SIZE+threadIdx.y;
+  locRow = threadIdx.y;
+  
+  colOff = threadIdx.x; //column offset of this thread
+  __shared__ float mins[BLOCK_SIZE][BLOCK_SIZE];
+  __shared__ int pos[BLOCK_SIZE][BLOCK_SIZE];
+  
+
+  // This loop finds the minimum of cols 
+  // [colOff, BLOCK_SIZE+colOff, 2*BLOCK_SIZE+colOff,...]
+  // and stores it in mins[locRow][colOff].
+  mins[locRow][colOff]=MAX_REAL;
+  pos[locRow][colOff]=-1;
+  for (i=0;i<D.pc/BLOCK_SIZE;i++){
+    curCol = i*BLOCK_SIZE+colOff;
+    if(curCol < D.c){ //ignore padding
+      temp = D.mat[IDX(row,curCol,D.ld)];
+      if(temp<mins[locRow][colOff]){
+       mins[locRow][colOff]=temp;
+       pos[locRow][colOff]=curCol;
+      }
+    }
+  }
+  __syncthreads();
+    
+  //Now find the min of cols [0, ... , BLOCK_SIZE]
+  for (i=BLOCK_SIZE/2; i>0;i/=2){
+    if(colOff<i){
+      //compare (col) to (col+i)
+      if(mins[locRow][colOff]>mins[locRow][colOff+i]){
+       mins[locRow][colOff]=mins[locRow][colOff+i];
+       pos[locRow][colOff]=pos[locRow][colOff+i];
+      }
+    }
+    __syncthreads();
+  }
+  
+  //arbitrarily use the first thread (along x) to set memory
+  if(threadIdx.x==0){  
+    dMins[row] = mins[locRow][0];
+    dMinIDs[row] = pos[locRow][0];
+  }
+}
+
+
+// Returns the min of each row of D.  dMins and dMinIDs 
+// are assumed to be (at least) of size D.r.
+__global__ void getKMinsKernel(matrix D, matrix dMins, intMatrix NNs, int k){
+  int row, locRow, colOff, i, curCol,j;
+  real temp;
+
+  row = blockIdx.y*BLOCK_SIZE+threadIdx.y;
+  locRow = threadIdx.y;
+
+  //printf("row=%d D.r =%d \n",row,D.r);
+  /* if(row>=D.r) */
+  /*   return; */
+
+  colOff = threadIdx.x; //column offset of this thread
+  __shared__ float mins[BLOCK_SIZE][BLOCK_SIZE];
+  __shared__ int pos[BLOCK_SIZE][BLOCK_SIZE];
+  
+  for(i=0;i<k;i++){
+    // This loop finds the minimum of cols 
+    // [colOff, BLOCK_SIZE+colOff, 2*BLOCK_SIZE+colOff,...]
+    // and stores it in mins[locRow][colOff].
+    mins[locRow][colOff]=MAX_REAL;
+    pos[locRow][colOff]=-1;
+    for (j=0;j<D.pc/BLOCK_SIZE;j++){
+      curCol = j*BLOCK_SIZE+colOff;
+      if(curCol < D.c){ //ignore padding
+       temp = D.mat[IDX(row,curCol,D.ld)];
+       if(temp<mins[locRow][colOff]){
+         mins[locRow][colOff]=temp;
+         pos[locRow][colOff]=curCol;
+       }
+      }
+    }
+    __syncthreads();
+
+
+    //Now find the min of cols [0, ... , BLOCK_SIZE]
+    for (j=BLOCK_SIZE/2; j>0; j/=2){
+      if(colOff<j){    
+       //compare (col) to (col+j)
+       if(mins[locRow][colOff]>mins[locRow][colOff+j]){
+         mins[locRow][colOff]=mins[locRow][colOff+j];
+         pos[locRow][colOff]=pos[locRow][colOff+j];
+       }
+      }
+       __syncthreads();
+    }
+    
+  //arbitrarily use the first thread (along x) to set memory
+    if(threadIdx.x==0 && row<D.r){  
+      dMins.mat[IDX(row,i,dMins.ld)] = mins[locRow][0];
+      NNs.mat[IDX(row,i,NNs.ld)] = pos[locRow][0];
+      D.mat[IDX(row,pos[locRow][0],D.ld)]=MAX_REAL;
+      
+    }
+    __syncthreads();
+  }
+}
+
+size_t countCompute(int*,int*,charMatrix);
+
+//This is used for debugging/research.
+size_t countCompute(int *groupCountQ, int *groupCountX, charMatrix cM){
+  int i,j;
+  size_t ans=0;
+  size_t avgBlocks=0;
+  size_t maxBlocks=0;
+  size_t maxBlocksInd;
+  size_t maxTemp;
+  size_t avgBlockQ=0;
+  size_t avgBlockX=0;
+  size_t maxBlockX=0;
+  size_t maxBlockQ=0;
+
+
+  for(i=0;i<cM.c;i++){
+    maxTemp=0;
+    for(j=0;j<cM.r;j++){
+      //printf("%d ",cM.mat[IDX(i,j,cM.ld)]*PAD(groupCountQ[i])*PAD(groupCountX[j]));
+      ans+=cM.mat[IDX(i,j,cM.ld)]*(groupCountQ[i])*(groupCountX[j]);
+      avgBlocks+=cM.mat[IDX(i,j,cM.ld)];
+      maxTemp+=cM.mat[IDX(i,j,cM.ld)]*PAD(groupCountX[j]);
+    }
+    //    printf("\n");
+    if(maxBlocks < maxTemp){
+      maxBlocks=maxTemp;
+      maxBlocksInd=PAD(groupCountQ[i]);
+    }
+    //maxBlocks=MAX(maxTemp,maxBlocks);
+  }
+  
+  for(i=0;i<cM.c;i++){
+    avgBlockQ+=groupCountQ[i];
+    avgBlockX+=groupCountX[i];
+    maxBlockQ=MAX(maxBlockQ,groupCountQ[i]);
+    maxBlockX=MAX(maxBlockX,groupCountX[i]);
+  }
+  
+  printf("most amt of work for a query: %zu (%zu) ; avg = %6.4f \n",maxBlocks,maxBlocksInd,((double)ans)/((double)cM.c));
+  printf("avg blocks/query block = %6.4f ; \n",((double)avgBlocks)/((double)cM.c));
+  printf("avg blockQ = %6.4f; max = %zu \n",((double)avgBlockQ)/((double)cM.c),maxBlockQ);
+  printf("avg blockX = %6.4f; max = %zu \n",((double)avgBlockX)/((double)cM.c),maxBlockX);
+  
+  return ans;
+}
+
+
+void kMinsWrap(matrix dD, matrix dMins, intMatrix dNNs){
+  dim3 block(BLOCK_SIZE,BLOCK_SIZE);
+  dim3 grid(1,dD.pr/BLOCK_SIZE);
+  
+  kMinsKernel<<<grid,block>>>(dD,dMins,dNNs);
+  cudaThreadSynchronize();
+}
+
+
+__global__ void kMinsKernel(matrix,matrix,intMatrix);
+
+__global__ void kMinsKernel(matrix D, matrix dMins, intMatrix NNs){
+  
+  int row = blockIdx.y*BLOCK_SIZE + threadIdx.y;
+  int ro = threadIdx.y; //row offset
+  int co = threadIdx.x; //col offset
+
+  __shared__ real smin[BLOCK_SIZE][BLOCK_SIZE];
+  __shared__ int smp[BLOCK_SIZE][BLOCK_SIZE];
+
+  real min, t;
+  int mp; //min position
+  int i,j,c;
+  
+   
+  for(i=0 ; i<NNs.c ; i++){
+    
+    min = MAX_REAL;
+    for(j=0 ; j<D.pc/BLOCK_SIZE ; j++){
+      c = j*BLOCK_SIZE;
+      if( c+co < D.c ){
+       t = D.mat[ IDX( row, c+co, D.ld ) ];
+       
+       if( t < min ){
+         min = t;
+         mp = c+co;
+       }
+      }
+    }
+    
+    smin[ro][co] = min;
+    smp[ro][co] = mp;
+    __syncthreads();
+
+
+    for(j=BLOCK_SIZE/2 ; j>0 ; j/=2){
+      if( co < j ){
+       if( smin[ro][co+j] < smin[ro][co] ){
+         smin[ro][co] = smin[ro][co+j];
+         smp[ro][co] = smp[ro][co+j];
+       }
+      }
+      __syncthreads();
+    }
+
+    if(co==0){
+      NNs.mat[ IDX( row, i, NNs.ld ) ] = smp[ro][0];
+      dMins.mat[ IDX( row, i, dMins.ld ) ] = smin[ro][0];
+
+      D.mat[ IDX( row, smp[ro][0], D.ld ) ] = MAX_REAL;
+    }
+
+    __syncthreads();
+  }
+}
+
+void brute2Step(matrix,matrix,intMatrix);
+
+
+void brute2Step(matrix x, matrix q, intMatrix NNs){
+
+  matrix dx, dq, dD, dMins;
+  intMatrix dNNs;
+  real *ranges, *dranges;
+  
+  copyAndMove(&dx, &x);
+  copyAndMove(&dq, &q);
+  copyAndMoveI(&dNNs, &NNs); //NNs.mat is garbage, but no matter.
+  
+  dD.r=q.r; dD.pr=q.pr; dD.c=x.r; dD.pc=x.pr; dD.ld=dD.pc;
+  cudaMalloc( (void**)&dD.mat, dD.pr*dD.pc*sizeof(*dD.mat) );
+  
+  dMins.r=NNs.r; dMins.pr=NNs.pr; dMins.c=NNs.c; dMins.pc=NNs.pc; dMins.ld=NNs.ld;
+  cudaMalloc( (void**)&dMins.mat, dMins.pr*dMins.pc*sizeof(*dMins.mat) );
+  
+  ranges = (real*)calloc( q.pr, sizeof(*ranges) );
+  cudaMalloc( (void**)&dranges, q.pr*sizeof(*dranges) );
+
+  dist1Wrap(dq,dx,dD);
+
+  kMinsWrap(dD, dMins, dNNs);
+  cudaMemcpy(NNs.mat, dNNs.mat, NNs.pr*NNs.pc*sizeof(*NNs.mat), cudaMemcpyDeviceToHost);
+  free(ranges);
+  cudaFree(dranges);
+  cudaFree(dx.mat);
+  cudaFree(dq.mat);
+  cudaFree(dD.mat);
+  cudaFree(dNNs.mat);
+  cudaFree(dMins.mat);
+
+}
diff --git a/driver.cu b/driver.cu
new file mode 100644 (file)
index 0000000..3d85631
--- /dev/null
+++ b/driver.cu
@@ -0,0 +1,220 @@
+#include<stdio.h>
+#include<stdlib.h>
+#include<cuda.h>
+#include<sys/time.h>
+#include<math.h>
+#include "defs.h"
+#include "utils.h"
+#include "utilsGPU.h"
+#include "rbc.h"
+#include "brute.h"
+#include "sKernel.h"
+
+void parseInput(int,char**);
+void readData(char*,int,int,real*);
+void orgData(real*,int,int,matrix,matrix);
+
+
+char *dataFile, *outFile;
+int n=0, m=0, d=0, numReps=0, s=0;
+int deviceNum=0;
+int main(int argc, char**argv){
+  real *data;
+  matrix x, q;
+  int *NNs, *NNsBrute;
+  int i;
+  struct timeval tvB,tvE;
+
+  printf("*****************\n");
+  printf("RANDOM BALL COVER\n");
+  printf("*****************\n");
+  
+  parseInput(argc,argv);
+
+  cuInit(0);
+  printf("Using GPU #%d\n",deviceNum);
+  if(cudaSetDevice(deviceNum) != cudaSuccess){
+    printf("Unable to select device %d.. exiting. \n",deviceNum);
+    exit(1);
+  }
+
+  unsigned int memFree, memTot;
+  CUcontext pctx;
+  unsigned int flags=0;
+  int device;
+  cudaGetDevice(&device);
+  cuCtxCreate(&pctx,flags,device);
+  cuMemGetInfo(&memFree, &memTot);
+  printf("GPU memory free = %u/%u (MB) \n",memFree/(1024*1024),memTot/(1024*1024));
+
+  data  = (real*)calloc( (n+m)*d, sizeof(*data) );
+  x.mat = (real*)calloc( PAD(n)*PAD(d), sizeof(*(x.mat)) );
+
+  //Need to allocate extra space, as each group of q will be padded later.
+  q.mat = (real*)calloc( PAD(m)*PAD(d), sizeof(*(q.mat)) );
+  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) );
+  for(i=0; i<m; i++)
+    NNs[i]=-1;
+  NNsBrute = (int*)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;
+  
+  /* printf("running brute force..\n"); */
+  /* gettimeofday(&tvB,NULL); */
+  /* bruteSearch(x,q,NNsBrute); */
+  /* 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); 
+  gettimeofday(&tvE,NULL);
+  printf("\t.. total time elapsed for rbc = %6.4f \n",timeDiff(tvB,tvE));
+  printf("finished \n");
+  
+  cudaError_t cE = cudaGetLastError();
+  if( cE != cudaSuccess ){
+    printf("Execution failed; error type: %s \n", cudaGetErrorString(cE) );
+  }
+
+  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;
+  
+  int *cnts = (int*)calloc(q.pr,sizeof(*cnts));
+  
+  bruteRangeCount(x,q,ranges,cnts);
+  
+  long int nc=0;
+  for(i=0;i<m;i++){
+    nc += cnts[i];
+  }
+  double mean = ((double)nc)/((double)m);
+  double var = 0.0;
+  for(i=0;i<m;i++) {
+    var += (((double)cnts[i])-mean)*(((double)cnts[i])-mean)/((double)m);
+  }
+  printf("\tavg rank = %6.4f; std dev = %6.4f \n\n", mean, sqrt(var));
+  
+  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(NNsBrute);
+  free(x.mat);
+  free(q.mat);
+}
+
+
+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");
+    exit(0);
+  }
+  
+  while(i<argc){
+    if(!strcmp(argv[i], "-f"))
+      dataFile = argv[++i];
+    else if(!strcmp(argv[i], "-n"))
+      n = atoi(argv[++i]);
+    else if(!strcmp(argv[i], "-m"))
+      m = atoi(argv[++i]);
+    else if(!strcmp(argv[i], "-d"))
+      d = atoi(argv[++i]);
+    else if(!strcmp(argv[i], "-r"))
+      numReps = atoi(argv[++i]);
+    else if(!strcmp(argv[i], "-s"))
+      s = atoi(argv[++i]);
+    else if(!strcmp(argv[i], "-o"))
+      outFile = argv[++i];
+    else if(!strcmp(argv[i], "-g"))
+      deviceNum = atoi(argv[++i]);
+    else{
+      fprintf(stderr,"%s : unrecognized option.. exiting\n",argv[i]);
+      exit(1);
+    }
+    i++;
+  }
+
+  if( !n || !m || !d || !numReps || !s || !dataFile){
+    fprintf(stderr,"more arguments needed.. exiting\n");
+    exit(1);
+  }
+  
+  if(numReps>n){
+    fprintf(stderr,"can't have more representatives than points.. exiting\n");
+    exit(1);
+  }
+}
+
+
+void readData(char *dataFile, int rows, int cols, real *data){
+  FILE *fp;
+  int numRead;
+
+  fp = fopen(dataFile,"r");
+  if(fp==NULL){
+    fprintf(stderr,"error opening file.. exiting\n");
+    exit(1);
+  }
+    
+  numRead = fread(data,sizeof(real),rows*cols,fp);
+  if(numRead != rows*cols){
+    fprintf(stderr,"error reading file.. exiting \n");
+    exit(1);
+  }
+  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){
+   
+  int i,fi,j;
+  int *p;
+  p = (int*)calloc(n,sizeof(*p));
+  
+  randPerm(n,p);
+
+  for(i=0,fi=0 ; i<x.r ; i++,fi++){
+    for(j=0;j<x.c;j++){
+      x.mat[IDX(i,j,x.ld)] = data[IDX(p[fi],j,d)];
+    }
+  }
+
+  for(i=0 ; i<q.r ; i++,fi++){
+    for(j=0;j<q.c;j++){
+      q.mat[IDX(i,j,q.ld)] = data[IDX(p[fi],j,d)];
+    } 
+  }
+
+  free(p);
+}
+
diff --git a/kernelWrap.cu b/kernelWrap.cu
new file mode 100644 (file)
index 0000000..a2f59b6
--- /dev/null
@@ -0,0 +1,62 @@
+#ifndef KERNELWRAP_CU
+#define KERNELWRAP_CU
+
+#include<cuda.h>
+#include<stdio.h>
+#include "kernels.h"
+#include "defs.h"
+
+void dist1Wrap(matrix dq, matrix dx, matrix dD){
+  dim3 block(BLOCK_SIZE,BLOCK_SIZE);
+  dim3 grid(dx.pr/BLOCK_SIZE,dq.pr/BLOCK_SIZE);
+  
+  dist1Kernel<<<grid,block>>>(dq,dx,dD);
+  cudaThreadSynchronize();
+}
+
+
+void findRangeWrap(matrix dD, real *dranges, int cntWant){
+  dim3 block(4*BLOCK_SIZE,BLOCK_SIZE/4);
+  dim3 grid(1,4*(dD.pr/BLOCK_SIZE));
+
+  findRangeKernel<<<grid,block>>>(dD,dranges,cntWant);
+  cudaThreadSynchronize();
+}
+
+void rangeSearchWrap(matrix dD, real *dranges, charMatrix dir){
+  dim3 block(BLOCK_SIZE,BLOCK_SIZE);
+  dim3 grid(dD.pc/BLOCK_SIZE,dD.pr/BLOCK_SIZE);
+  
+  rangeSearchKernel<<<grid,block>>>(dD,dranges,dir);
+  cudaThreadSynchronize();
+}
+
+void nnWrap(const matrix dx, const matrix dy, real *dMins, int *dMinIDs){
+  dim3 dimBlock(BLOCK_SIZE,BLOCK_SIZE);
+  dim3 dimGrid;
+  
+  dimGrid.x = 1;
+  dimGrid.y = dx.pr/dimBlock.y + (dx.pr%dimBlock.y==0 ? 0 : 1);
+  nnKernel<<<dimGrid,dimBlock>>>(dx,dy,dMins,dMinIDs);
+  cudaThreadSynchronize();
+}
+
+
+void rangeCountWrap(const matrix dq, const matrix dx, real *dranges, int *dcounts){
+  dim3 block(BLOCK_SIZE,BLOCK_SIZE);
+  dim3 grid(1,dq.pr/BLOCK_SIZE);
+
+  rangeCountKernel<<<grid,block>>>(dq,dx,dranges,dcounts);
+  cudaThreadSynchronize();
+}
+
+
+/*NOTE: can be deleted */
+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();
+}
+#endif
diff --git a/kernelWrap.h b/kernelWrap.h
new file mode 100644 (file)
index 0000000..2d5f80d
--- /dev/null
@@ -0,0 +1,13 @@
+#ifndef KERNELWRAP_H
+#define KERNELWRAP_H
+
+#include "defs.h"
+
+void dist1Wrap(matrix,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*);
+#endif
diff --git a/kernels.cu b/kernels.cu
new file mode 100644 (file)
index 0000000..f015d3d
--- /dev/null
@@ -0,0 +1,389 @@
+#ifndef KERNELS_CU
+#define KERNELS_CU
+
+#include<cuda.h>
+#include "defs.h"
+#include "kernels.h"
+#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;
+  
+  
+  __shared__ real min[BLOCK_SIZE][BLOCK_SIZE];
+  __shared__ int minPos[BLOCK_SIZE][BLOCK_SIZE];
+
+  __shared__ real Xb[BLOCK_SIZE][BLOCK_SIZE];
+  __shared__ real Qb[BLOCK_SIZE][BLOCK_SIZE];
+
+  __shared__ int g; //group of q  
+  __shared__ int numGroups;
+  __shared__ int groupCount;
+  __shared__ int groupOff;
+
+  if(offX==0 && offQ==0){
+    g = cP.qToGroup[qBlock]; 
+    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)];
+    }
+
+    __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;
+
+       //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( 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)] );
+       __syncthreads();
+       
+       for(l=0;l<BLOCK_SIZE;l++){
+         ans+=abs(Xb[l][offX]-Qb[l][offQ]);
+       }
+       __syncthreads();
+      }
+      
+      //compare to previous min and store into shared mem if needed.
+      if(j*BLOCK_SIZE+offX<groupCount && ans<min[offQ][offX]){
+       min[offQ][offX]=ans;
+       minPos[offQ][offX]= xmap.mat[IDX( g, xBlock+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];      
+      }
+    }
+    __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;
+  }
+}
+
+
+__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;
+
+  __shared__ real min[BLOCK_SIZE][BLOCK_SIZE];
+  __shared__ int minPos[BLOCK_SIZE][BLOCK_SIZE];
+
+  __shared__ real Xb[BLOCK_SIZE][BLOCK_SIZE];
+  __shared__ real Qb[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;
+      
+      //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)];
+
+      __syncthreads();
+
+      for(k=0;k<BLOCK_SIZE;k++){
+       ans+=abs(Xb[k][offX]-Qb[k][offQ]);
+      }
+      __syncthreads();
+    }
+   
+    
+    if( xBlock+offX<X.r && ans<min[offQ][offX] ){
+       minPos[offQ][offX] = xBlock+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];      
+      }
+    }
+    __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];
+  }
+}
+
+
+
+
+
+
+
+
+
+
+__device__ void dist1Kernel(const matrix Q, const matrix X, matrix D){
+  int c, i, j;
+
+  int qB = blockIdx.y*BLOCK_SIZE;
+  int q  = threadIdx.y;
+  int xB = blockIdx.x*BLOCK_SIZE;
+  int x = threadIdx.x;
+
+  real ans=0;
+
+  //This thread is responsible for computing the dist between Q[qB+q] and X[xB+x]
+  
+  __shared__ real Qs[BLOCK_SIZE][BLOCK_SIZE];
+  __shared__ real Xs[BLOCK_SIZE][BLOCK_SIZE];
+
+
+  for(i=0 ; i<Q.pc/BLOCK_SIZE ; i++){
+    c=i*BLOCK_SIZE; //current col block
+
+    Qs[x][q] = Q.mat[ IDX(qB+q, c+x, Q.ld) ];
+    Xs[x][q] = X.mat[ IDX(xB+q, c+x, X.ld) ];
+
+    __syncthreads();
+
+    for(j=0 ; j<BLOCK_SIZE ; j++)
+      ans += abs( Qs[j][q] - Xs[j][x] );
+        
+    __syncthreads();
+  }
+  
+  D.mat[ IDX( qB+q, xB+x, D.ld ) ] = ans;
+
+}
+
+
+__global__ void findRangeKernel(matrix D, real *ranges, int cntWant){
+  
+  int row = blockIdx.y*(BLOCK_SIZE/4)+threadIdx.y;
+  int ro = threadIdx.y;
+  int co = threadIdx.x;
+  int i, c;
+  real t;
+
+  const int LB = (90*cntWant)/100 ;
+  const int UB = cntWant; 
+
+  __shared__ real smin[BLOCK_SIZE/4][4*BLOCK_SIZE];
+  __shared__ real smax[BLOCK_SIZE/4][4*BLOCK_SIZE];
+  
+  real min=MAX_REAL;
+  real max=0;
+  for(c=0 ; c<D.pc ; c+=(4*BLOCK_SIZE)){
+    if( c+co < D.c ){
+      t = D.mat[ IDX( row, c+co, D.ld ) ];
+      min = MIN(t,min);
+      max = MAX(t,max);
+    }
+  }
+  
+  smin[ro][co] = min;
+  smax[ro][co] = max;
+  __syncthreads();
+  
+  for(i=2*BLOCK_SIZE ; i>0 ; i/=2){
+    if( co < i ){
+      smin[ro][co] = MIN( smin[ro][co], smin[ro][co+i] );
+      smax[ro][co] = MAX( smax[ro][co], smax[ro][co+i] );
+    }
+    __syncthreads();
+  }
+
+  //Now start range counting.
+
+  int itcount=0;
+  int cnt;
+  real rg;
+  __shared__ int scnt[BLOCK_SIZE/4][4*BLOCK_SIZE];
+  __shared__ char cont[BLOCK_SIZE/4];
+  
+  if(co==0)
+    cont[ro]=1;
+  
+  do{
+    itcount++;
+    __syncthreads();
+
+    if( cont[ro] )  //if we didn't actually need to cont, leave rg as it was.
+      rg = ( smax[ro][0] + smin[ro][0] ) / ((real)2.0) ;
+
+    cnt=0;
+    for(c=0 ; c<D.pc ; c+=(4*BLOCK_SIZE)){
+      cnt += (c+co < D.c && row < D.r && D.mat[ IDX( row, c+co, D.ld ) ] <= rg);
+    }
+
+    scnt[ro][co] = cnt;
+    __syncthreads();
+    
+    for(i=2*BLOCK_SIZE ; i>0 ; i/=2){
+      if( co < i ){
+       scnt[ro][co] += scnt[ro][co+i];
+      }
+      __syncthreads();
+    }
+    
+    if(co==0){
+      if( scnt[ro][0] < cntWant )
+       smin[ro][0]=rg;
+      else
+       smax[ro][0]=rg;
+    }
+    
+    // cont[ro] == this row needs to continue
+    if(co==0)
+      cont[ro] = row<D.r && ( scnt[ro][0] < LB || scnt[ro][0] > UB ); 
+    __syncthreads();
+
+    // Determine if *any* of the rows need to continue
+    for(i=BLOCK_SIZE/8 ; i>0 ; i/=2){
+      if( ro < i && co==0)
+       cont[ro] |= cont[ro+i];
+      __syncthreads();
+    }
+    
+  } while(cont[0]);
+
+  if(co==0 && row<D.r )
+    ranges[row]=rg;
+  
+}
+
+
+__global__ void rangeSearchKernel(matrix D, real *ranges, charMatrix ir){
+  int col = blockIdx.x*BLOCK_SIZE + threadIdx.x;
+  int row = blockIdx.y*BLOCK_SIZE + threadIdx.y;
+
+  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;
+  
+  real rg = ranges[q+qo];
+  
+  int r,c,i;
+
+  __shared__ int scnt[BLOCK_SIZE][BLOCK_SIZE];
+
+  __shared__ real xs[BLOCK_SIZE][BLOCK_SIZE];
+  __shared__ real qs[BLOCK_SIZE][BLOCK_SIZE];
+  
+  int cnt=0;
+  for( r=0; r<X.pr; r+=BLOCK_SIZE ){
+
+    real dist=0;
+    for( c=0; c<X.pc; c+=BLOCK_SIZE){
+      xs[xo][qo] = X.mat[IDX( r+qo, c+xo, X.ld )];
+      qs[xo][qo] = Q.mat[IDX( q+qo, c+xo, Q.ld )];
+      __syncthreads();
+      
+      for( i=0; i<BLOCK_SIZE; i++)
+       dist += abs(xs[i][xo]-qs[i][qo]);
+      __syncthreads();
+
+    }
+    cnt += r+xo<X.r && dist<rg;
+
+  }
+  
+  scnt[qo][xo]=cnt;
+  __syncthreads();
+  
+  for( i=BLOCK_SIZE/2; i>0; i/=2 ){
+    if( xo<i ){
+      scnt[qo][xo] += scnt[qo][xo+i];
+    }
+    __syncthreads();
+  }
+
+  if( xo==0 && q+qo<Q.r )
+    counts[q+qo] = scnt[qo][0];
+}
+
+
+#endif
diff --git a/kernels.h b/kernels.h
new file mode 100644 (file)
index 0000000..c29b174
--- /dev/null
+++ b/kernels.h
@@ -0,0 +1,15 @@
+#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,const matrix,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,real*,charMatrix);
+__global__ void rangeCountKernel(const matrix,const matrix,real*,int*);
+
+#endif
diff --git a/oldKernels.cu b/oldKernels.cu
new file mode 100644 (file)
index 0000000..b810ed5
--- /dev/null
@@ -0,0 +1,59 @@
+/* ****** Device helper functions ****** */
+
+__device__ matrix getSubMat(matrix A, int row, int col){
+  matrix As;
+  
+  As.c = BLOCK_SIZE;
+  As.r = BLOCK_SIZE;
+  As.ld = A.ld;
+  As.mat = &A.mat[A.ld*row*BLOCK_SIZE + col*BLOCK_SIZE];
+  
+  return As;
+}
+
+
+
+//l_1-norm
+/* __global__ void dist1Kernel(const matrix Q, const matrix X, matrix D){ */
+/*   int qBlockRow = blockIdx.y; */
+/*   int xBlockRow = blockIdx.x; */
+
+/*   //matrix Dsub = getSubMat(D,xBlockRow,qBlockRow); */
+
+/*   int qID = threadIdx.y; */
+/*   int xID = threadIdx.x; */
+
+/*   //  printf("calling (%d,%d) \n",qBlockRow*BLOCK_SIZE+qID,xBlockRow*BLOCK_SIZE+xID); */
+  
+/*   int i,j; */
+/*   real ans=0; */
+  
+/*   //Note: assumes that X is padded. */
+/*   for(i=0;i<X.pc/BLOCK_SIZE;i++){ */
+/*     //    matrix Xs = getSubMat(X,xBlockRow,i); */
+/*     //    matrix Qs = getSubMat(Q,qBlockRow,i); */
+    
+/*     __shared__ real Xb[BLOCK_SIZE][BLOCK_SIZE]; */
+/*     __shared__ real Qb[BLOCK_SIZE][BLOCK_SIZE]; */
+    
+/*     // Each thread loads one element of Xs and Qs into shared mem */
+/*     // Note that the indexing is swapped to increase memory coalescing. */
+/*     //printf("reading x[%d,%d] \n",xBlockRow*BLOCK_SIZE+qID,i*BLOCK_SIZE+xID); */
+/*     Xb[xID][qID]=X.mat[IDX(xBlockRow*BLOCK_SIZE+qID,i*BLOCK_SIZE+xID,X.ld)];// getElement(Xs,qID,xID); */
+/*     //printf("reading q[%d,%d] \n",qBlockRow*BLOCK_SIZE+qID,i*BLOCK_SIZE+xID); */
+/*     Qb[xID][qID]=Q.mat[IDX(qBlockRow*BLOCK_SIZE+qID,i*BLOCK_SIZE+xID,Q.ld)]; */
+    
+/*     __syncthreads(); */
+    
+/*     for(j=0;j<BLOCK_SIZE;j++) */
+/*       ans+=abs(Xb[j][xID]-Qb[j][qID]); */
+        
+/*       __syncthreads(); */
+/*   } */
+  
+/*   //  Dsub.mat[IDX(qID,xID,Dsub.ld)]=ans;//setElement(Dsub,qID,xID,ans); */
+/*   //printf("writing (%d,%d) %6.2f \n",qBlockRow*BLOCK_SIZE+qID,xBlockRow*BLOCK_SIZE+xID,ans); */
+/*   //  Dsub.mat[IDX(qID,xID,Dsub.ld)]=ans; */
+/*   D.mat[IDX(qBlockRow*BLOCK_SIZE+qID,xBlockRow*BLOCK_SIZE+xID,D.ld)]=ans; */
+  
+/* } */
diff --git a/rbc.cu b/rbc.cu
new file mode 100644 (file)
index 0000000..db3f339
--- /dev/null
+++ b/rbc.cu
@@ -0,0 +1,504 @@
+#ifndef RBC_CU
+#define RBC_CU
+
+#include<sys/time.h>
+#include<stdio.h>
+#include<cuda.h>
+#include "utils.h"
+#include "defs.h"
+#include "utilsGPU.h"
+#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));
+
+  for(i=0;i<m;i++)
+    qID[i]=i;
+
+  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);
+
+  
+  //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); 
+  
+
+  gettimeofday(&tvB,NULL);  //Start the timer for the queries
+  
+  matrix dqOrig;
+  copyAndMove(&dqOrig, &q);
+
+  computeReps(dqOrig, dr, repIDsQ, distToRepsQ);
+
+  //How many points are assigned to each group?
+  computeCounts(repIDsQ, m, groupCountQ);
+
+  //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);
+  idIntersection(cM);
+
+  // Setup the compute plan
+  initCompPlan(&cP, cM, groupCountQ, groupCountX, numReps);
+
+  // Compute the NNs according to the compute plan
+  computeNNs(dx, dqOrig, xmap, cP, qID, NNs, compLength);
+  
+  gettimeofday(&tvE,NULL);
+  printf("\t.. time elapsed (for queries) = %6.4f \n",timeDiff(tvB,tvE));
+
+  cudaFree(dr.mat); 
+  cudaFree(dqOrig.mat);
+  freeCompPlan(&cP);
+  free(cM.mat);
+  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;
+  
+  //Figure out how much fits into memory
+  unsigned int memFree, memTot;
+  cuMemGetInfo(&memFree, &memTot);
+  memFree = (int)(((float)memFree)*MEM_USABLE);
+     
+  //memory 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");
+    exit(1);
+  }
+
+  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) );
+
+  real *dranges;
+  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;
+  ir.mat = (char*)calloc( ir.pr*ir.pc, sizeof(*ir.mat) );
+  charMatrix dir;
+  copyAndMoveC(&dir, &ir);
+
+  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) );
+  
+  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)
+  
+  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
+    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);
+    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 );
+
+    numLeft -= pi;
+    row += pi;
+  }
+  
+  cudaFree(dCnts);
+  free(ir.mat);
+  cudaFree(dranges);
+  cudaFree(dir.mat);
+  cudaFree(dSums.mat);
+  cudaFree(dD.mat);
+}
+
+
+//Assign each point in dq to its nearest point in dr.  
+void computeReps(matrix dq, matrix dr, int *repIDs, real *distToReps){
+  real *dMins;
+  int *dMinIDs;
+
+  cudaMalloc((void**)&(dMins), dq.pr*sizeof(*dMins));
+  cudaMalloc((void**)&(dMinIDs), dq.pr*sizeof(*dMinIDs));
+  
+  nnWrap(dq,dr,dMins,dMinIDs);
+  
+  cudaMemcpy(distToReps,dMins,dq.r*sizeof(*dMins),cudaMemcpyDeviceToHost);
+  cudaMemcpy(repIDs,dMinIDs,dq.r*sizeof(*dMinIDs),cudaMemcpyDeviceToHost);
+  
+  cudaFree(dMins);
+  cudaFree(dMinIDs);
+}
+
+
+
+//Assumes radii is initialized to 0s
+void computeRadii(int *repIDs, real *distToReps, real *radii, int n, int numReps){
+  int i;
+
+  for(i=0;i<n;i++)
+    radii[repIDs[i]] = MAX(distToReps[i],radii[repIDs[i]]);
+}
+
+
+//Assumes groupCount is initialized to 0s
+void computeCounts(int *repIDs, int n, int *groupCount){
+  int 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]]++;
+  }
+  
+  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
+  int *yID;
+
+  gS = (int*)calloc(numReps+1,sizeof(*gS));
+  
+  for( i=0; i<n; i++ )
+    gS[repIDs[i]+1]++;
+  for( i=0; i<numReps+1; i++ )
+    gS[i]=PAD(gS[i]);
+  
+  for( i=1; i<numReps+1; i++ )
+    gS[i]=gS[i-1]+gS[i];
+  
+  *compLength = gS[numReps];
+  
+  yID = (int*)calloc((*compLength),sizeof(*yID));
+  for( i=0; i<(*compLength); i++ )
+    yID[i]=DUMMY_IDX;
+  
+
+  for( i=0; i<n; i++ ){
+    yID[gS[repIDs[i]]]=qID[i];
+    gS[repIDs[i]]++;
+  }
+
+  for( i=0; i<(*compLength); i++ ){
+    qID[i]=yID[i];
+  }
+  
+  free(yID);
+  free(gS);
+}
+
+
+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;
+  for(i=0;i<cM.r;i++){
+    if(i<cM.c)
+      cM.mat[IDX(i,i,cM.ld)]=1;
+  }
+}
+
+
+
+void fullIntersection(charMatrix cM){
+  int i,j;
+  for(i=0;i<cM.r;i++){
+    for(j=0;j<cM.c;j++){
+      cM.mat[IDX(i,j,cM.ld)]=1;
+    }
+  }
+}
+
+
+// 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;
+  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;
+  }
+
+  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);
+}
+
+
+//This calls the dist1Kernel wrapper, but has it compute only 
+//a submatrix of the all-pairs distance matrix.  In particular,
+//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){
+  dr.r=dr.pr=length;
+  dr.mat = &dr.mat[IDX( start, 0, dr.ld )];
+  dist1Wrap(dr, dx, dD);
+}
+
+
+#endif
diff --git a/rbc.h b/rbc.h
new file mode 100644 (file)
index 0000000..0d0d517
--- /dev/null
+++ b/rbc.h
@@ -0,0 +1,26 @@
+#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 idIntersection(charMatrix);
+void fullIntersection(charMatrix);
+void initCompPlan(compPlan*,charMatrix,int*,int*,int);
+void freeCompPlan(compPlan*);
+void computeNNs(matrix,matrix,intMatrix, compPlan,int*,int*,int);
+
+
+
+
+#endif
diff --git a/sKernel.cu b/sKernel.cu
new file mode 100644 (file)
index 0000000..12342a8
--- /dev/null
@@ -0,0 +1,204 @@
+#ifndef SKERNEL_CU
+#define SKERNEL_CU
+
+#include<stdio.h>
+#include "sKernel.h"
+#include "defs.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;
+  
+  const int l=SCAN_WIDTH; //length
+
+  int off=1;
+
+  __shared__ int 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;
+    
+  //up-sweep
+  for( d=l>>1; d > 0; d>>=1 ){
+    __syncthreads();
+    
+    if( id < d ){
+      ssum[ off*(2*id+2)-1 ] += ssum[ off*(2*id+1)-1 ];
+    }
+    off *= 2;
+  }
+  
+  __syncthreads();
+  
+  if ( id == 0 ){
+    sumaux.mat[IDX( r, blockIdx.x, sumaux.ld )] = ssum[ l-1 ];
+    ssum[ l-1 ] = 0;
+  }
+  
+  //down-sweep
+  for ( d=1; d<l; d*=2 ){
+    off >>= 1;
+    __syncthreads();
+    
+    if( id < d ){
+      t = ssum[ off*(2*id+1)-1 ];
+      ssum[ off*(2*id+1)-1 ] = ssum[ off*(2*id+2)-1 ];
+      ssum[ off*(2*id+2)-1 ] += t;
+    }
+    
+  }
+
+  __syncthreads();
+  if( bo+2*id < 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 n = in.c;
+  int numScans = (n+SCAN_WIDTH-1)/SCAN_WIDTH;
+
+  if(numScans > SCAN_WIDTH){
+    fprintf(stderr,"scan is too large.. exiting \n");
+    exit(1);
+  }
+  
+  intMatrix dAux;
+  dAux.r=dAux.pr=in.r; dAux.c=dAux.pc=dAux.ld=numScans;
+  intMatrix dDummy;
+  dDummy.r=dDummy.pr=in.r; dDummy.c=dDummy.pc=dDummy.ld=1;
+  
+  cudaMalloc( (void**)&dAux.mat, dAux.pr*dAux.pc*sizeof(*dAux.mat) );
+  cudaMalloc( (void**)&dDummy.mat, dDummy.pr*dDummy.pc*sizeof(*dDummy.mat) );
+
+  dim3 block( SCAN_WIDTH/2, 1 );
+  dim3 grid( numScans, in.r );
+  
+  sumKernel<<<grid,block>>>(in, sum, dAux, n);
+  cudaThreadSynchronize();
+  grid.x=1;
+
+  sumKernelI<<<grid,block>>>(dAux, dAux, dDummy, numScans);
+  cudaThreadSynchronize();
+
+  grid.x=numScans;
+  combineSumKernel<<<grid,block>>>(sum,dAux,n);
+  cudaThreadSynchronize();
+
+  cudaFree(dAux.mat);
+  cudaFree(dDummy.mat);
+
+}
+
+
+//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;
+  
+  const int l=SCAN_WIDTH; //length
+
+  int off=1;
+
+  __shared__ int 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;
+    
+  //up-sweep
+  for( d=l>>1; d > 0; d>>=1 ){
+    __syncthreads();
+    
+    if( id < d ){
+      ssum[ off*(2*id+2)-1 ] += ssum[ off*(2*id+1)-1 ];
+    }
+    off *= 2;
+  }
+  
+  __syncthreads();
+  
+  if ( id == 0 ){
+    sumaux.mat[IDX( r, blockIdx.x, sumaux.ld )] = ssum[ l-1 ];
+    ssum[ l-1 ] = 0;
+  }
+  
+  //down-sweep
+  for ( d=1; d<l; d*=2 ){
+    off >>= 1;
+    __syncthreads();
+    
+    if( id < d ){
+      t = ssum[ off*(2*id+1)-1 ];
+      ssum[ off*(2*id+1)-1 ] = ssum[ off*(2*id+2)-1 ];
+      ssum[ off*(2*id+2)-1 ] += t;
+    }
+    
+  }
+
+  __syncthreads();
+  if( bo+2*id < 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];
+
+}
+
+#endif
diff --git a/sKernel.h b/sKernel.h
new file mode 100644 (file)
index 0000000..b38c123
--- /dev/null
+++ b/sKernel.h
@@ -0,0 +1,14 @@
+#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);
+
+#endif
diff --git a/utils.cu b/utils.cu
new file mode 100644 (file)
index 0000000..237bce8
--- /dev/null
+++ b/utils.cu
@@ -0,0 +1,149 @@
+#ifndef UTILS_CU
+#define UTILS_CU
+
+#include<sys/time.h>
+#include<stdio.h>
+#include "utils.h"
+#include "defs.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));
+    
+  struct timeval t3;
+  gettimeofday(&t3,NULL);
+  srand(t3.tv_usec);
+  
+  for(i=0;i<n;i++)
+    y[i]=i;
+  
+  for(i=0;i<MIN(l,n-1);i++){  //The n-1 bit is necessary because you can't swap the last 
+                              //element with something larger.
+    ri=randBetween(i+1,n);
+    swap(&y[i],&y[ri]);
+  }
+  
+  for(i=0;i<l;i++)
+    x[i]=y[i];
+  free(y);
+}
+
+//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;
+  
+  struct timeval t3;
+  gettimeofday(&t3,NULL);
+  srand(t3.tv_usec);
+  
+  for(i=0;i<n;i++){
+    x[i]=i;
+  }
+  
+  for(i=0;i<n-1;i++){
+    ri=randBetween(i+1,n);
+    swap(&x[i],&x[ri]);
+  }
+}
+
+void swap(int *a, int *b){
+  int t;
+  t=*a; *a=*b; *b=t;
+}
+
+//generates a rand int in rand [a,b) 
+int randBetween(int a, int b){
+  int val,c;
+
+  if(b<=a){
+    fprintf(stderr,"misuse of randBetween.. exiting\n");
+    exit(1);
+  }
+  c= b-a;
+
+  while(c<= (val= rand()/(int)(((unsigned)RAND_MAX + 1) / c)));
+  val=val+a;
+  
+  return val;
+}
+
+
+void printMat(matrix A){
+  int 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)]);
+    printf("\n");
+  }
+}
+
+
+void printMatWithIDs(matrix A, int *id){
+  int 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)]);
+    printf("%d ",id[i]);
+    printf("\n");
+  }
+}
+
+
+void printCharMat(charMatrix A){
+  int 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)]);
+    printf("\n");
+  }
+}
+
+void printVector(real *x, int d){
+  int i;
+
+  for(i=0 ; i<d; i++)
+    printf("%6.2f ",x[i]);
+  printf("\n");
+}
+
+
+void copyVector(real *x, real *y, int d){
+  int i;
+  
+  for(i=0;i<d;i++)
+    x[i]=y[i];
+}
+
+
+void copyMat(matrix *x, matrix *y){
+  int 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++){
+    for(j=0; j<y->c; j++){
+      x->mat[IDX( i, j, x->ld )] = y->mat[IDX( i, j, y->ld )];
+    }
+  }
+}
+
+
+real distL1(matrix x, matrix y, int k, int l){
+  int 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)]);
+  return ans;
+}
+
+double timeDiff(struct timeval start, struct timeval end){
+  return (double)(end.tv_sec+end.tv_usec/1e6 - start.tv_sec - start.tv_usec/1e6); 
+}
+
+
+
+#endif
diff --git a/utils.h b/utils.h
new file mode 100644 (file)
index 0000000..7882f58
--- /dev/null
+++ b/utils.h
@@ -0,0 +1,18 @@
+#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 printMat(matrix);
+void printMatWithIDs(matrix,int*);
+void printCharMat(charMatrix);
+void printVector(real*,int);
+void copyVector(real*,real*,int);
+real distL1(matrix,matrix,int,int);
+double timeDiff(struct timeval,struct timeval);
+void copyMat(matrix*,matrix*);
+#endif
diff --git a/utilsGPU.cu b/utilsGPU.cu
new file mode 100644 (file)
index 0000000..a28409e
--- /dev/null
@@ -0,0 +1,62 @@
+#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;
+}
+
+void copyAndMove(matrix *dx, const matrix *x){
+  dx->r = x->r; 
+  dx->c = x->c;
+  dx->pr = x->pr;
+  dx->pc = x->pc;
+  dx->ld = x->ld;
+
+  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 copyAndMoveI(intMatrix *dx, const intMatrix *x){
+  dx->r = x->r; 
+  dx->c = x->c;
+  dx->pr = x->pr;
+  dx->pc = x->pc;
+  dx->ld = x->ld;
+
+  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 copyAndMoveC(charMatrix *dx, const charMatrix *x){
+  dx->r = x->r; 
+  dx->c = x->c;
+  dx->pr = x->pr;
+  dx->pc = x->pc;
+  dx->ld = x->ld;
+
+  cudaMalloc( (void**)&(dx->mat), dx->pr*dx->pc*sizeof(*(dx->mat)) );
+  cudaMemcpy( dx->mat, x->mat, dx->pr*dx->pc*sizeof(*(dx->mat)), cudaMemcpyHostToDevice );
+  
+}
+#endif
diff --git a/utilsGPU.h b/utilsGPU.h
new file mode 100644 (file)
index 0000000..cba437c
--- /dev/null
@@ -0,0 +1,15 @@
+#ifndef UTILSGPU_H
+#define UTILSGPU_H
+
+#include "defs.h"
+
+memPlan createMemPlan(int,int);
+
+void copyAndMove(matrix*,const matrix*);
+void copyAndMoveI(intMatrix*,const intMatrix*);
+void copyAndMoveC(charMatrix*,const charMatrix*);
+
+
+
+
+#endif