Modified the scan kernel so that it can work with > 1024*1024 numbers
authorLawrence Cayton <lcayton@tuebingen.mpg.de>
Mon, 16 Aug 2010 09:52:03 +0000 (11:52 +0200)
committerLawrence Cayton <lcayton@tuebingen.mpg.de>
Mon, 16 Aug 2010 09:52:03 +0000 (11:52 +0200)
Makefile
defs.h
driver.cu
kernelWrap.cu
kernels.cu
rbc.cu
readme.txt
sKernel.cu
todo.txt
utils.cu
utils.h

index 6d58877..30a7799 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -3,7 +3,7 @@ NVCC=nvcc
 CCFLAGS=
 NVCCFLAGS= --ptxas-options=-v
 #other flags: -deviceemu -arch=sm_20 --compiler-bindir=/usr/bin/gcc-4.3
-LINKFLAGS=-lcuda
+LINKFLAGS=-lcuda 
 #other linkflags: 
 SOURCES=
 CUSOURCES=driver.cu utils.cu utilsGPU.cu rbc.cu kernels.cu brute.cu kernelWrap.cu sKernel.cu
diff --git a/defs.h b/defs.h
index 158fccb..ba308f0 100644 (file)
--- a/defs.h
+++ b/defs.h
@@ -6,7 +6,7 @@
 #define FLOAT_TOL 1e-7
 #define BLOCK_SIZE 16 //must be a power of 2
 
-#define MAX_BS 65535 //max block size
+#define MAX_BS 65535 //max block size (specified by CUDA)
 #define SCAN_WIDTH 1024
 
 // Format that the data is manipulated in:
index 3d85631..6ec580f 100644 (file)
--- a/driver.cu
+++ b/driver.cu
@@ -24,6 +24,7 @@ int main(int argc, char**argv){
   int *NNs, *NNsBrute;
   int i;
   struct timeval tvB,tvE;
+  cudaError_t cE;
 
   printf("*****************\n");
   printf("RANDOM BALL COVER\n");
@@ -37,7 +38,8 @@ int main(int argc, char**argv){
     printf("Unable to select device %d.. exiting. \n",deviceNum);
     exit(1);
   }
-
+  
+  
   unsigned int memFree, memTot;
   CUcontext pctx;
   unsigned int flags=0;
@@ -89,7 +91,7 @@ int main(int argc, char**argv){
   printf("\t.. total time elapsed for rbc = %6.4f \n",timeDiff(tvB,tvE));
   printf("finished \n");
   
-  cudaError_t cE = cudaGetLastError();
+  cE = cudaGetLastError();
   if( cE != cudaSuccess ){
     printf("Execution failed; error type: %s \n", cudaGetErrorString(cE) );
   }
@@ -99,9 +101,11 @@ int main(int argc, char**argv){
   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));
-  
+  gettimeofday(&tvB,NULL);
   bruteRangeCount(x,q,ranges,cnts);
+  gettimeofday(&tvE,NULL);
   
   long int nc=0;
   for(i=0;i<m;i++){
@@ -113,7 +117,8 @@ int main(int argc, char**argv){
     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));
-  
+  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) );
index b692f7b..292d941 100644 (file)
@@ -25,7 +25,7 @@ void dist1Wrap(matrix dq, matrix dx, matrix dD){
     }
     numDoneX += todoX;
   }
-  
+
   cudaThreadSynchronize();
 }
 
index d3fc003..d73d837 100644 (file)
@@ -36,13 +36,14 @@ __global__ void planNNKernel(const matrix Q, const matrix X, real *dMins, int *d
   __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);
@@ -68,7 +69,7 @@ __global__ void planNNKernel(const matrix Q, const matrix X, real *dMins, int *d
        }
        __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;
diff --git a/rbc.cu b/rbc.cu
index 95aedd8..f6cb34d 100644 (file)
--- a/rbc.cu
+++ b/rbc.cu
@@ -77,6 +77,14 @@ void rbc(matrix x, matrix q, int numReps, int s, int *NNs){
   
   build(dx, dr, xmap, groupCountX, s); 
   
+  /* 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
   
@@ -87,7 +95,8 @@ void rbc(matrix x, matrix q, int numReps, int s, int *NNs){
 
   //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);
 
@@ -125,14 +134,14 @@ 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
+  //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);
@@ -156,7 +165,7 @@ 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) );
-  
+
   int *dCnts;
   cudaMalloc( (void**)&dCnts, ptsAtOnce*sizeof(*dCnts) );
     
@@ -175,17 +184,18 @@ void build(const matrix dx, const matrix dr, intMatrix xmap, int *counts, int s)
     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
-    printf("after range search\n");
+    
     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);
@@ -459,6 +469,7 @@ void computeNNs(matrix dx, matrix dq, intMatrix xmap, compPlan cP, int *qIDs, in
   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();
@@ -467,7 +478,6 @@ void computeNNs(matrix dx, matrix dq, intMatrix xmap, compPlan cP, int *qIDs, in
 
   cudaMemcpy( NNs, dMinIDs, dq.r*sizeof(*NNs), cudaMemcpyDeviceToHost);
   
-  
   cudaFree(dcP.numGroups);
   cudaFree(dcP.groupCountX);
   cudaFree(dcP.groupOff);
index 9a792b0..bfbb216 100644 (file)
@@ -37,7 +37,7 @@ 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
+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
index d643927..a830f5f 100644 (file)
@@ -2,8 +2,10 @@
 #define SKERNEL_CU
 
 #include<stdio.h>
+#include<math.h>
 #include "sKernel.h"
 #include "defs.h"
+#include "utils.h"
 
 __global__ void sumKernel(charMatrix in, intMatrix sum, intMatrix sumaux, int n){
   int id = threadIdx.x;
@@ -17,7 +19,6 @@ __global__ void sumKernel(charMatrix in, intMatrix sum, intMatrix sumaux, int n)
 
   __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;
     
@@ -63,7 +64,7 @@ __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)
@@ -74,8 +75,9 @@ __global__ void combineSumKernel(intMatrix sum, intMatrix daux, int n){
 
 __global__ void getCountsKernel(int *counts, charMatrix ir, intMatrix sums){
   int r = blockIdx.x*BLOCK_SIZE + threadIdx.x;
-  if ( r < ir.r )
+  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 )];
+  }
 }
 
 
@@ -95,7 +97,6 @@ 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);
-
 }
 
 
@@ -108,41 +109,55 @@ void buildMapWrap(intMatrix map, charMatrix ir, intMatrix sums, int offSet){
 
 
 void sumWrap(charMatrix in, intMatrix sum){
+  int i;
   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) );
 
+  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, n);
-  cudaThreadSynchronize();
-
-  //While(){}
-  grid.x=1;
-
-  sumKernelI<<<grid,block>>>(dAux, dAux, dDummy, numScans);
+  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();
+  }
 
-  grid.x=numScans;
-  combineSumKernel<<<grid,block>>>(sum,dAux,n);
+    
+  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();
-
-  cudaFree(dAux.mat);
-  cudaFree(dDummy.mat);
-
+  
+  for( i=0; i <= depth; i++)
+    cudaFree(dAux[i].mat);
+  
+  free(dAux);
+  free(width);
 }
 
 
@@ -159,10 +174,9 @@ __global__ void sumKernelI(intMatrix in, intMatrix sum, intMatrix sumaux, int n)
 
   __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();
@@ -180,6 +194,7 @@ __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;
@@ -192,14 +207,16 @@ __global__ void sumKernelI(intMatrix in, intMatrix sum, intMatrix sumaux, int n)
     }
     
   }
-
+  
   __syncthreads();
-  if( bo+2*id < n ) 
+  
+  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
index 8fb52c8..457f219 100644 (file)
--- a/todo.txt
+++ b/todo.txt
@@ -1,3 +1,4 @@
+- make use of void* to clean up code
 - fix problems related to limitations on number of kernel launches
 - Make sure that the code works with r,s that are not powers of 2
 
index 237bce8..495d674 100644 (file)
--- a/utils.cu
+++ b/utils.cu
@@ -102,6 +102,15 @@ void printCharMat(charMatrix A){
   }
 }
 
+void printIntMat(intMatrix A){
+  int 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("\n");
+  }
+}
+
 void printVector(real *x, int d){
   int i;
 
diff --git a/utils.h b/utils.h
index 7882f58..9f9edf1 100644 (file)
--- a/utils.h
+++ b/utils.h
@@ -10,6 +10,7 @@ int randBetween(int,int);
 void printMat(matrix);
 void printMatWithIDs(matrix,int*);
 void printCharMat(charMatrix);
+void printIntMat(intMatrix);
 void printVector(real*,int);
 void copyVector(real*,real*,int);
 real distL1(matrix,matrix,int,int);