k-nn brute implemented
authorLawrence Cayton <lcayton@tuebingen.mpg.de>
Wed, 17 Nov 2010 11:49:31 +0000 (12:49 +0100)
committerLawrence Cayton <lcayton@tuebingen.mpg.de>
Wed, 17 Nov 2010 11:49:31 +0000 (12:49 +0100)
Makefile
brute.cu
brute.h
defs.h
driver.cu
kernelWrap.cu
kernelWrap.h
kernels.cu
kernels.h

index f4f1af3..61ccaf8 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -1,9 +1,9 @@
 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 
+LINKFLAGS=-lcuda -lgsl -lgslcblas -lm
 #other linkflags: 
 SOURCES=
 CUSOURCES= driver.cu utils.cu utilsGPU.cu rbc.cu brute.cu kernels.cu kernelWrap.cu sKernel.cu sKernelWrap.cu
index d2377c0..ad106f6 100644 (file)
--- a/brute.cu
+++ b/brute.cu
@@ -14,6 +14,7 @@
 #include "brute.h"
 #include<stdio.h>
 #include<cuda.h>
+#include<gsl/gsl_sort.h>
 
 void bruteRangeCount(matrix x, matrix q, real *ranges, unint *cnts){
   matrix dx, dq;
@@ -67,6 +68,35 @@ void bruteSearch(matrix x, matrix q, unint *NNs){
 }
 
 
+void bruteK(matrix x, matrix q, intMatrix NNs){
+  matrix dMins;
+  intMatrix 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;
+  dMins.r=q.r; dMins.pr=q.pr; dMins.c=K; dMins.pc=K; dMins.ld=dMins.pc;
+  dMinIDs.r=q.r; dMinIDs.pr=q.pr; dMinIDs.c=K; dMinIDs.pc=K; dMinIDs.ld=dMinIDs.pc;
+
+  checkErr( cudaMalloc((void**)&dMins.mat, dMins.pc*dMins.pr*sizeof(*dMins.mat)) );
+  checkErr( cudaMalloc((void**)&dMinIDs.mat, dMinIDs.pc*dMinIDs.pr*sizeof(*dMinIDs.mat)) );
+  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);
+  
+  knnWrap(dq,dx,dMins,dMinIDs);
+
+  cudaMemcpy(NNs.mat,dMinIDs.mat,NNs.pr*NNs.pc*sizeof(*NNs.mat),cudaMemcpyDeviceToHost);
+  
+  cudaFree(dMins.mat);
+  cudaFree(dMinIDs.mat);
+  cudaFree(dx.mat);
+  cudaFree(dq.mat);
+}
+
+
 void bruteCPU(matrix X, matrix Q, unint *NNs){
   real *dtoNNs; 
   real temp;
@@ -89,4 +119,34 @@ void bruteCPU(matrix X, matrix Q, unint *NNs){
   
   free(dtoNNs);  
 }
+
+
+void bruteKCPU(matrix x, matrix q, intMatrix NNs){
+  int i, j;
+
+  float **d;
+  d = (float**)calloc(q.pr, sizeof(*d));
+  size_t **t;
+  t = (size_t**)calloc(q.pr, sizeof(*t));
+  for( i=0; i<q.pr; i++){
+    d[i] = (float*)calloc(x.pr, sizeof(**d));
+    t[i] = (size_t*)calloc(x.pr, sizeof(**t));
+  }
+
+  //#pragma omp parallel for private(j)
+  for( i=0; i<q.r; i++){
+    for( j=0; j<x.r; j++)
+      d[i][j] = distVec( q, x, i, j );
+    gsl_sort_float_index(t[i], d[i], 1, x.r);
+    for ( j=0; j<K; j++)
+      NNs.mat[IDX( i, j, NNs.ld )] = t[i][j];
+  }
+
+  for( i=0; i<q.pr; i++){
+    free(t[i]);
+    free(d[i]);
+  }
+  free(t);
+  free(d);
+}
 #endif
diff --git a/brute.h b/brute.h
index cb25e9f..a0a5bb3 100644 (file)
--- a/brute.h
+++ b/brute.h
@@ -10,5 +10,6 @@
 void bruteRangeCount(matrix,matrix,real*,unint*);
 void bruteSearch(matrix,matrix,unint*);
 void bruteCPU(matrix,matrix,unint*);
-
+void bruteK(matrix,matrix,intMatrix);
+void bruteKCPU(matrix,matrix,intMatrix);
 #endif
diff --git a/defs.h b/defs.h
index 2cbd4d3..c2a6be0 100644 (file)
--- a/defs.h
+++ b/defs.h
@@ -4,6 +4,8 @@
 #include<float.h>
 
 #define FLOAT_TOL 1e-7
+
+#define K 32 //for k-NN
 #define BLOCK_SIZE 16 //must be a power of 2 (current 
 // implementation of findRange requires a power of 4, in fact)
 
@@ -42,8 +44,10 @@ typedef float real;
 #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))
+#define MAXi(i,j,k,l) ((i) > (j) ? (k) : (l)) //indexed version
+#define MINi(i,j,k,l) ((i) <= (j) ? (k) : (l))
 
-#define MIN(i,j) ((i) < (j) ? (i) : (j))
 
 typedef unsigned int unint;
 
index 7bd2711..aa7c06c 100644 (file)
--- a/driver.cu
+++ b/driver.cu
@@ -26,7 +26,8 @@ unint deviceNum=0;
 int main(int argc, char**argv){
   real *data;
   matrix x, q;
-  unint *NNs, *NNsBrute;
+  unint *NNs;
+  intMatrix NNsK, NNsKCPU;
   unint i;
   struct timeval tvB,tvE;
   cudaError_t cE;
@@ -65,75 +66,103 @@ int main(int argc, char**argv){
   NNs = (unint*)calloc( m, sizeof(*NNs) );
   for(i=0; i<m; i++)
     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);
 
   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)); */
+    NNs[i]=DUMMY_IDX;
   
-  printf("\nrunning rbc..\n");
+  NNsK.r=q.r; NNsK.pr=q.pr; NNsK.pc=NNsK.c=K; NNsK.ld=NNsK.pc;
+  NNsKCPU.r=q.r; NNsKCPU.pr=q.pr; NNsKCPU.pc=NNsKCPU.c=K; NNsKCPU.ld=NNsKCPU.pc;
+  NNsK.mat = (unint*)calloc(NNsK.pr*NNsK.pc, sizeof(*NNsK.mat));
+  NNsKCPU.mat = (unint*)calloc(NNsKCPU.pr*NNsKCPU.pc, sizeof(*NNsKCPU.mat));
+
+  printf("running k-brute force..\n");
   gettimeofday(&tvB,NULL);
-  buildRBC(x, &rbcS, numReps, s);
+  bruteK(x,q,NNsK);
   gettimeofday(&tvE,NULL);
-  printf("\t.. build time for rbc = %6.4f \n",timeDiff(tvB,tvE));
+  printf("\t.. time elapsed = %6.4f \n",timeDiff(tvB,tvE));
 
+  printf("running regular brute force..\n");
   gettimeofday(&tvB,NULL);
-  queryRBC(q, rbcS, NNs);
+  bruteSearch(x,q,NNs);
   gettimeofday(&tvE,NULL);
-  printf("\t.. query time for rbc = %6.4f \n",timeDiff(tvB,tvE));
+  printf("\t.. time elapsed = %6.4f \n",timeDiff(tvB,tvE));
+
+
+  printf("running cpu version..\n");
+  gettimeofday(&tvB,NULL);
+  bruteKCPU(x,q,NNsKCPU);
+  gettimeofday(&tvE,NULL);
+  printf("\t.. time elapsed = %6.4f \n",timeDiff(tvB,tvE));
+
+
+  /* printf("\nrunning rbc..\n"); */
+  /* gettimeofday(&tvB,NULL); */
+  /* buildRBC(x, &rbcS, numReps, s); */
+  /* gettimeofday(&tvE,NULL); */
+  /* printf("\t.. build time for rbc = %6.4f \n",timeDiff(tvB,tvE)); */
 
-  destroyRBC(&rbcS);
-  printf("finished \n");
+  /* 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();
   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++){
-    if(NNs[i]>n) printf("error");
-    ranges[i] = distVec(q,x,i,NNs[i]) - 10e-6;
+    int j;
+    for(j=0;j<K;j++){
+      if(NNsK.mat[IDX(i,j,NNsK.ld)] != NNsKCPU.mat[IDX(i,j,NNsKCPU.ld)])
+       printf("%d %d: (%6.2f %6.2f) ",i,j,distVec(q,x,i,NNsK.mat[IDX(i,j,NNsK.ld)]),distVec(q,x,i,NNsKCPU.mat[IDX(i,j,NNsKCPU.ld)]));
+    }
+ /*    printf("\n"); */
   }
-
-  unint *cnts = (unint*)calloc(q.pr,sizeof(*cnts));
-  gettimeofday(&tvB,NULL);
-  bruteRangeCount(x,q,ranges,cnts);
-  gettimeofday(&tvE,NULL);
+  /* 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++){ */
+  /*   if(NNs[i]>n) printf("error"); */
+  /*   ranges[i] = distVec(q,x,i,NNs[i]) - 10e-6; */
+  /* } */
+
+  /* unint *cnts = (unint*)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++){
-    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));
-  printf("(range count took %6.4f) \n", timeDiff(tvB, tvE));
-
-
-  if(outFile){
-    FILE* fp = fopen(outFile, "a");
-    fprintf( fp, "%d %d %6.5f %6.5f \n", numReps, s, mean, sqrt(var) );
-    fclose(fp);
-  }
-
-  free(ranges);
-  free(cnts);
+  /* 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)); */
+  /* printf("(range count took %6.4f) \n", timeDiff(tvB, tvE)); */
+
+
+  /* if(outFile){ */
+  /*   FILE* fp = fopen(outFile, "a"); */
+  /*   fprintf( fp, "%d %d %6.5f %6.5f \n", numReps, s, mean, sqrt(var) ); */
+  /*   fclose(fp); */
+  /* } */
+
+  /* free(ranges); */
+  /* free(cnts); */
   free(NNs);
-  free(NNsBrute);
+  free(NNsK.mat);
+  free(NNsKCPU.mat);
   free(x.mat);
   free(q.mat);
 }
index e5a69ce..1cee847 100644 (file)
@@ -92,6 +92,25 @@ void nnWrap(const matrix dq, const matrix dx, real *dMins, unint *dMinIDs){
 }
 
 
+void knnWrap(const matrix dq, const matrix dx, matrix dMins, intMatrix dMinIDs){
+  dim3 block(BLOCK_SIZE,BLOCK_SIZE);
+  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;
+    knnKernel<<<grid,block>>>(dq,numDone,dx,dMins,dMinIDs);
+    numDone += todo;
+  }
+  cudaThreadSynchronize();
+
+}
+
+
 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;
index b5190fe..b30bd76 100644 (file)
@@ -11,6 +11,7 @@ void kMinsWrap(matrix,matrix,intMatrix);
 void findRangeWrap(const matrix,real*,unint);
 void rangeSearchWrap(const matrix,const real*,charMatrix);
 void nnWrap(const matrix,const matrix,real*,unint*);
+void knnWrap(const matrix,const matrix,matrix,intMatrix);
 void rangeCountWrap(const matrix,const matrix,real*,unint*);
 void planNNWrap(const matrix,const unint*,const matrix,const intMatrix,real*,unint*,compPlan,unint);
 
index aeccf87..954645f 100644 (file)
@@ -149,6 +149,87 @@ __global__ void nnKernel(const matrix Q, unint numDone, const matrix X, real *dM
 }
 
 
+__global__ void knnKernel(const matrix Q, unint numDone, const matrix X, matrix dMins, intMatrix dMinIDs){
+  
+  unint qB = blockIdx.y * BLOCK_SIZE + numDone;  //indexes Q
+  unint xB; //indexes X;
+  unint cB; //colBlock
+  unint offQ = threadIdx.y; //the offset of qPos in this block
+  unint offX = threadIdx.x; //ditto for x
+  unint i;
+  real ans;
+
+  __shared__ real Xs[BLOCK_SIZE][BLOCK_SIZE];
+  __shared__ real Qs[BLOCK_SIZE][BLOCK_SIZE];
+  
+  __shared__ real D[BLOCK_SIZE][BLOCK_SIZE];
+  __shared__ unint id[BLOCK_SIZE][BLOCK_SIZE];
+  __shared__ real dNN[BLOCK_SIZE][K+BLOCK_SIZE];
+  __shared__ unint idNN[BLOCK_SIZE][K+BLOCK_SIZE];
+
+  dNN[offQ][offX] = MAX_REAL;
+  dNN[offQ][offX+16] = MAX_REAL;
+  idNN[offQ][offX] = DUMMY_IDX;
+  idNN[offQ][offX+16] = DUMMY_IDX;
+  
+  __syncthreads();
+
+  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.
+      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(i=0;i<BLOCK_SIZE;i++)
+       ans += DIST( Xs[i][offX], Qs[i][offQ] );
+      
+      __syncthreads();
+    }
+    D[offQ][offX] = (xB+offX<X.r)? ans:MAX_REAL;
+    id[offQ][offX] = xB + offX;
+    __syncthreads();
+    
+/*     if(offX==0 && offQ==0 & qB==0){  */
+/*       printf("before sort: \n"); */
+/*       for(i=0;i<BLOCK_SIZE;i++)  */
+/*      printf("%6.2f ",D[0][i]);  */
+/*        printf("\n");  */
+/*      }  */
+
+    sort16( D, id );
+/*      if(offX==0 && offQ==0 & qB==0){  */
+/*        printf("after sort:\n"); */
+/*        for(i=0;i<BLOCK_SIZE;i++)  */
+/*      printf("%6.2f ",D[0][i]);  */
+/*        printf("\n");  */
+/*      }  */
+      
+    __syncthreads();
+    dNN[offQ][offX+32] = D[offQ][offX];
+    idNN[offQ][offX+32] = id[offQ][offX];
+    __syncthreads();
+    merge32x16( dNN, idNN );
+/*     if(offX==0 && offQ==0 & qB==0){  */
+/*       for(i=0;i<K+BLOCK_SIZE;i++) */
+/*     printf("%6.2f ",dNN[0][i]);  */
+/*       printf("\n");  */
+    /*   for(i=0;i<K+BLOCK_SIZE;i++)  */
+/*             printf("%d ",idNN[0][i]);  */
+     /*  printf("\n");  */
+/*     } */
+  }
+  __syncthreads();
+  
+  dMins.mat[IDX(qB+offQ, offX, dMins.ld)] = dNN[offQ][offX];
+  dMins.mat[IDX(qB+offQ, offX+16, dMins.ld)] = dNN[offQ][offX+16];
+  dMinIDs.mat[IDX(qB+offQ, offX, dMins.ld)] = idNN[offQ][offX];
+  dMinIDs.mat[IDX(qB+offQ, offX+16, dMins.ld)] = idNN[offQ][offX+16];
+  
+}
 
 
 __global__ void dist1Kernel(const matrix Q, unint qStart, const matrix X, unint xStart, matrix D){
@@ -340,4 +421,99 @@ __global__ void rangeCountKernel(const matrix Q, unint numDone, const matrix X,
 }
 
 
+__device__ void sort16(real x[][16], unint xi[][16]){
+  int i = threadIdx.x;
+  int j = threadIdx.y;
+
+  if(i%2==0)
+    mmGateI( x[j]+i, x[j]+i+1, xi[j]+i, xi[j]+i+1 );
+  __syncthreads();
+
+  if(i%4<2)
+    mmGateI( x[j]+i, x[j]+i+2, xi[j]+i, xi[j]+i+2 );
+  __syncthreads();
+
+  if(i%4==1)
+    mmGateI( x[j]+i, x[j]+i+1, xi[j]+i, xi[j]+i+1 );
+  __syncthreads();
+  
+  if(i%8<4)
+    mmGateI( x[j]+i, x[j]+i+4, xi[j]+i, xi[j]+i+4 );
+  __syncthreads();
+  
+  if(i%8==2 || i%8==3)
+    mmGateI( x[j]+i, x[j]+i+2, xi[j]+i, xi[j]+i+2 );
+  __syncthreads();
+
+  if( i%2 && i%8 != 7 ) 
+    mmGateI( x[j]+i, x[j]+i+1, xi[j]+i, xi[j]+i+1 );
+  __syncthreads();
+  
+  //0-7; 8-15 now sorted.  merge time.
+  if( i<8)
+    mmGateI( x[j]+i, x[j]+i+8, xi[j]+i, xi[j]+i+8 );
+  __syncthreads();
+  
+  if( i>3 && i<8 )
+    mmGateI( x[j]+i, x[j]+i+4, xi[j]+i, xi[j]+i+4 );
+  __syncthreads();
+  
+  int os = (i/2)*4+2 + i%2;
+  if(i<6)
+    mmGateI( x[j]+os, x[j]+os+2, xi[j]+os, xi[j]+os+2 );
+  __syncthreads();
+  
+  if( i%2 && i<15)
+    mmGateI( x[j]+i, x[j]+i+1, xi[j]+i, xi[j]+i+1 );
+
+}
+
+
+__device__ void merge32x16(real x[][48], unint xi[][48]){
+  int i = threadIdx.x;
+  int j = threadIdx.y;
+
+  mmGateI( x[j]+i, x[j]+i+32, xi[j]+i, xi[j]+i+32 );
+  __syncthreads();
+
+  mmGateI( x[j]+i+16, x[j]+i+32, xi[j]+i+16, xi[j]+i+32 );
+  __syncthreads();
+
+  int os = (i<8)? 24: 0;
+  mmGateI( x[j]+os+i, x[j]+os+i+8, xi[j]+os+i, xi[j]+os+i+8 );
+  __syncthreads();
+  
+  os = (i/4)*8+4 + i%4;
+  mmGateI( x[j]+os, x[j]+os+4, xi[j]+os, xi[j]+os+4 );
+  if(i<4)
+    mmGateI(x[j]+36+i, x[j]+36+i+4, xi[j]+36+i, xi[j]+36+i+4 );
+  __syncthreads();
+
+  os = (i/2)*4+2 + i%2;
+  mmGateI( x[j]+os, x[j]+os+2, xi[j]+os, xi[j]+os+2 );
+  
+  os = (i/2)*4+34 + i%2;
+  if(i<6)
+    mmGateI( x[j]+os, x[j]+os+2, xi[j]+os, xi[j]+os+2 );
+  __syncthreads();
+
+  os = 2*i+1;
+  mmGateI(x[j]+os, x[j]+os+1, xi[j]+os, xi[j]+os+1 );
+
+  os = 2*i+33;
+  if(i<7)
+    mmGateI(x[j]+os, x[j]+os+1, xi[j]+os, xi[j]+os+1 );
+
+}
+
+
+__device__ void mmGateI(real *x, real *y, unint *xi, unint *yi){
+  int ti = MINi( *x, *y, *xi, *yi );
+  *yi = MAXi( *x, *y, *xi, *yi );
+  *xi = ti;
+  real t = MIN( *x, *y );
+  *y = MAX( *x, *y );
+  *x = t;
+}
+
 #endif
index 4d03611..0845ae9 100644 (file)
--- a/kernels.h
+++ b/kernels.h
@@ -9,8 +9,11 @@
 __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 knnKernel(const matrix,unint,const matrix,matrix,intMatrix);
 __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*);
-
+__device__ void sort16(real[][16],unint[][16]);
+__device__ void merge32x16(real[][48],unint[][48]);
+__device__ void mmGateI(real*,real*,unint*,unint*);
 #endif