updated text files
[RBC.git] / rbc.cu
1 /* This file is part of the Random Ball Cover (RBC) library.
2  * (C) Copyright 2010, Lawrence Cayton [lcayton@tuebingen.mpg.de]
3  */
4
5 #ifndef RBC_CU
6 #define RBC_CU
7
8 #include<sys/time.h>
9 #include<stdio.h>
10 #include<cuda.h>
11 #include "utils.h"
12 #include "defs.h"
13 #include "utilsGPU.h"
14 #include "rbc.h"
15 #include "kernels.h"
16 #include "kernelWrap.h"
17 #include "sKernelWrap.h"
18
19 void queryRBC(const matrix q, const rbcStruct rbcS, unint *NNs, real* NNdists){
20   unint m = q.r;
21   unint numReps = rbcS.dr.r;
22   unint compLength;
23   compPlan dcP;
24   unint *qMap, *dqMap;
25   qMap = (unint*)calloc(PAD(m+(BLOCK_SIZE-1)*PAD(numReps)),sizeof(*qMap));
26   matrix dq;
27   copyAndMove(&dq, &q);
28   
29   charMatrix cM;
30   cM.r=cM.c=numReps; cM.pr=cM.pc=cM.ld=PAD(numReps);
31   cM.mat = (char*)calloc( cM.pr*cM.pc, sizeof(*cM.mat) );
32   
33   unint *repIDsQ;
34   repIDsQ = (unint*)calloc( m, sizeof(*repIDsQ) );
35   real *distToRepsQ;
36   distToRepsQ = (real*)calloc( m, sizeof(*distToRepsQ) );
37   unint *groupCountQ;
38   groupCountQ = (unint*)calloc( PAD(numReps), sizeof(*groupCountQ) );
39   
40   computeReps(dq, rbcS.dr, repIDsQ, distToRepsQ);
41
42   //How many points are assigned to each group?
43   computeCounts(repIDsQ, m, groupCountQ);
44   
45   //Set up the mapping from groups to queries (qMap).
46   buildQMap(q, qMap, repIDsQ, numReps, &compLength);
47   
48   // Setup the computation matrix.  Currently, the computation matrix is 
49   // just the identity matrix: each query assigned to a particular 
50   // representative is compared only to that representative's points.  
51   idIntersection(cM);
52
53   initCompPlan(&dcP, cM, groupCountQ, rbcS.groupCount, numReps);
54
55   checkErr( cudaMalloc( (void**)&dqMap, compLength*sizeof(*dqMap) ) );
56   cudaMemcpy( dqMap, qMap, compLength*sizeof(*dqMap), cudaMemcpyHostToDevice );
57   
58   computeNNs(rbcS.dx, rbcS.dxMap, dq, dqMap, dcP, NNs, NNdists, compLength);
59   
60   free(qMap);
61   cudaFree(dqMap);
62   freeCompPlan(&dcP);
63   cudaFree(dq.mat);
64   free(cM.mat);
65   free(repIDsQ);
66   free(distToRepsQ);
67   free(groupCountQ);
68 }
69
70 //This function is very similar to queryRBC, with a couple of basic changes to handle
71 //k-nn.  
72 void kqueryRBC(const matrix q, const rbcStruct rbcS, intMatrix NNs, matrix NNdists){
73   unint m = q.r;
74   unint numReps = rbcS.dr.r;
75   unint compLength;
76   compPlan dcP;
77   unint *qMap, *dqMap;
78   qMap = (unint*)calloc(PAD(m+(BLOCK_SIZE-1)*PAD(numReps)),sizeof(*qMap));
79   matrix dq;
80   copyAndMove(&dq, &q);
81   
82   charMatrix cM;
83   cM.r=cM.c=numReps; cM.pr=cM.pc=cM.ld=PAD(numReps);
84   cM.mat = (char*)calloc( cM.pr*cM.pc, sizeof(*cM.mat) );
85   
86   unint *repIDsQ;
87   repIDsQ = (unint*)calloc( m, sizeof(*repIDsQ) );
88   real *distToRepsQ;
89   distToRepsQ = (real*)calloc( m, sizeof(*distToRepsQ) );
90   unint *groupCountQ;
91   groupCountQ = (unint*)calloc( PAD(numReps), sizeof(*groupCountQ) );
92   
93   computeReps(dq, rbcS.dr, repIDsQ, distToRepsQ);
94
95   //How many points are assigned to each group?
96   computeCounts(repIDsQ, m, groupCountQ);
97   
98   //Set up the mapping from groups to queries (qMap).
99   buildQMap(q, qMap, repIDsQ, numReps, &compLength);
100   
101   // Setup the computation matrix.  Currently, the computation matrix is 
102   // just the identity matrix: each query assigned to a particular 
103   // representative is compared only to that representative's points.  
104
105   // NOTE: currently, idIntersection is the *only* computation matrix 
106   // that will work properly with k-nn search (this is not true for 1-nn above).
107   idIntersection(cM);
108
109   initCompPlan(&dcP, cM, groupCountQ, rbcS.groupCount, numReps);
110
111   checkErr( cudaMalloc( (void**)&dqMap, compLength*sizeof(*dqMap) ) );
112   cudaMemcpy( dqMap, qMap, compLength*sizeof(*dqMap), cudaMemcpyHostToDevice );
113   
114   computeKNNs(rbcS.dx, rbcS.dxMap, dq, dqMap, dcP, NNs, NNdists, compLength);
115
116   free(qMap);
117   freeCompPlan(&dcP);
118   cudaFree(dq.mat);
119   free(cM.mat);
120   free(repIDsQ);
121   free(distToRepsQ);
122   free(groupCountQ);
123 }
124
125
126 void buildRBC(const matrix x, rbcStruct *rbcS, unint numReps, unint s){
127   unint n = x.pr;
128   intMatrix xmap;
129
130   setupReps(x, rbcS, numReps);
131   copyAndMove(&rbcS->dx, &x);
132   
133   xmap.r=numReps; xmap.pr=PAD(numReps); xmap.c=s; xmap.pc=xmap.ld=PAD(s);
134   xmap.mat = (unint*)calloc( xmap.pr*xmap.pc, sizeof(*xmap.mat) );
135   copyAndMoveI(&rbcS->dxMap, &xmap);
136   rbcS->groupCount = (unint*)calloc( PAD(numReps), sizeof(*rbcS->groupCount) );
137   
138   //Figure out how much fits into memory
139   size_t memFree, memTot;
140   cudaMemGetInfo(&memFree, &memTot);
141   memFree = (unint)(((float)memFree)*MEM_USABLE);
142   /* mem needed per rep:
143    *  n*sizeof(real) - dist mat
144    *  n*sizeof(char) - dir
145    *  n*sizeof(int)  - dSums
146    *  sizeof(real)   - dranges
147    *  sizeof(int)    - dCnts
148    *  MEM_USED_IN_SCAN - memory used internally
149    */
150   unint ptsAtOnce = DPAD(memFree/((n+1)*sizeof(real) + n*sizeof(char) + (n+1)*sizeof(unint) + 2*MEM_USED_IN_SCAN(n)));
151   if(!ptsAtOnce){
152     fprintf(stderr,"error: %lu is not enough memory to build the RBC.. exiting\n", (unsigned long)memFree);
153     exit(1);
154   }
155
156   //Now set everything up for the scans
157   matrix dD;
158   dD.pr=dD.r=ptsAtOnce; dD.c=rbcS->dx.r; dD.pc=rbcS->dx.pr; dD.ld=dD.pc;
159   checkErr( cudaMalloc( (void**)&dD.mat, dD.pr*dD.pc*sizeof(*dD.mat) ) );
160   
161   real *dranges;
162   checkErr( cudaMalloc( (void**)&dranges, ptsAtOnce*sizeof(real) ) );
163
164   charMatrix ir;
165   ir.r=dD.r; ir.pr=dD.pr; ir.c=dD.c; ir.pc=dD.pc; ir.ld=dD.ld;
166   ir.mat = (char*)calloc( ir.pr*ir.pc, sizeof(*ir.mat) );
167   charMatrix dir;
168   copyAndMoveC(&dir, &ir);
169
170   intMatrix dSums; //used to compute memory addresses.
171   dSums.r=dir.r; dSums.pr=dir.pr; dSums.c=dir.c; dSums.pc=dir.pc; dSums.ld=dir.ld;
172   checkErr( cudaMalloc( (void**)&dSums.mat, dSums.pc*dSums.pr*sizeof(*dSums.mat) ) );
173
174   unint *dCnts;
175   checkErr( cudaMalloc( (void**)&dCnts, ptsAtOnce*sizeof(*dCnts) ) );
176   
177   //Do the scans to build the dxMap
178   unint numLeft = rbcS->dr.r; //points left to process
179   unint row = 0; //base row for iteration of while loop
180   unint pi, pip; //pi=pts per it, pip=pad(pi)
181
182   while( numLeft > 0 ){
183     pi = MIN(ptsAtOnce, numLeft);  //points to do this iteration.
184     pip = PAD(pi);
185     dD.r = pi; dD.pr = pip; dir.r=pi; dir.pr=pip; dSums.r=pi; dSums.pr=pip;
186
187     distSubMat(rbcS->dr, rbcS->dx, dD, row, pip); //compute the distance matrix
188     findRangeWrap(dD, dranges, s);  //find an appropriate range
189     rangeSearchWrap(dD, dranges, dir); //set binary vector for points in range
190     sumWrap(dir, dSums);  //This and the next call perform the parallel compaction.
191     buildMapWrap(rbcS->dxMap, dir, dSums, row);
192     getCountsWrap(dCnts,dir,dSums);  //How many points are assigned to each rep?  It is not 
193                                      //*exactly* s, which is why we need to compute this.
194     cudaMemcpy( &rbcS->groupCount[row], dCnts, pi*sizeof(*rbcS->groupCount), cudaMemcpyDeviceToHost );
195     
196     numLeft -= pi;
197     row += pi;
198   }
199   
200
201   cudaFree(dCnts);
202   free(ir.mat);
203   free(xmap.mat);
204   cudaFree(dranges);
205   cudaFree(dir.mat);
206   cudaFree(dSums.mat);
207   cudaFree(dD.mat);
208 }
209
210
211 // Choose representatives and move them to device
212 void setupReps(matrix x, rbcStruct *rbcS, unint numReps){
213   unint i;
214   unint *randInds;
215   randInds = (unint*)calloc( PAD(numReps), sizeof(*randInds) );
216   subRandPerm(numReps, x.r, randInds);
217   
218   matrix r;
219   r.r=numReps; r.pr=PAD(numReps); r.c=x.c; r.pc=r.ld=PAD(r.c); 
220   r.mat = (real*)calloc( r.pr*r.pc, sizeof(*r.mat) );
221
222   for(i=0;i<numReps;i++)
223     copyVector(&r.mat[IDX(i,0,r.ld)], &x.mat[IDX(randInds[i],0,x.ld)], x.c);
224   
225   copyAndMove(&rbcS->dr, &r);
226
227   free(randInds);
228   free(r.mat);
229 }
230
231
232 //Assign each point in dq to its nearest point in dr.  
233 void computeReps(matrix dq, matrix dr, unint *repIDs, real *distToReps){
234   real *dMins;
235   unint *dMinIDs;
236
237   checkErr( cudaMalloc((void**)&(dMins), dq.pr*sizeof(*dMins)) );
238   checkErr( cudaMalloc((void**)&(dMinIDs), dq.pr*sizeof(*dMinIDs)) );
239   
240   nnWrap(dq,dr,dMins,dMinIDs);
241   
242   cudaMemcpy(distToReps,dMins,dq.r*sizeof(*dMins),cudaMemcpyDeviceToHost);
243   cudaMemcpy(repIDs,dMinIDs,dq.r*sizeof(*dMinIDs),cudaMemcpyDeviceToHost);
244   
245   cudaFree(dMins);
246   cudaFree(dMinIDs);
247 }
248
249
250 //Assumes radii is initialized to 0s
251 void computeRadii(unint *repIDs, real *distToReps, real *radii, unint n, unint numReps){
252   unint i;
253
254   for(i=0;i<n;i++)
255     radii[repIDs[i]] = MAX(distToReps[i],radii[repIDs[i]]);
256 }
257
258
259 //Assumes groupCount is initialized to 0s
260 void computeCounts(unint *repIDs, unint n, unint *groupCount){
261   unint i;
262   
263   for(i=0;i<n;i++)
264     groupCount[repIDs[i]]++;
265 }
266
267
268 void buildQMap(matrix q, unint *qMap, unint *repIDs, unint numReps, unint *compLength){
269   unint n=q.r;
270   unint i;
271   unint *gS; //groupSize
272   
273   gS = (unint*)calloc(numReps+1,sizeof(*gS));
274   
275   for( i=0; i<n; i++ )
276     gS[repIDs[i]+1]++;
277   for( i=0; i<numReps+1; i++ )
278     gS[i]=PAD(gS[i]);
279   
280   for( i=1; i<numReps+1; i++ )
281     gS[i]=gS[i-1]+gS[i];
282   
283   *compLength = gS[numReps];
284   
285   for( i=0; i<(*compLength); i++ )
286     qMap[i]=DUMMY_IDX;
287   
288   for( i=0; i<n; i++ ){
289     qMap[gS[repIDs[i]]]=i;
290     gS[repIDs[i]]++;
291   }
292
293   free(gS);
294 }
295
296
297 // Sets the computation matrix to the identity.  
298 void idIntersection(charMatrix cM){
299   unint i;
300   for(i=0;i<cM.r;i++){
301     if(i<cM.c)
302       cM.mat[IDX(i,i,cM.ld)]=1;
303   }
304 }
305
306
307 void fullIntersection(charMatrix cM){
308   unint i,j;
309   for(i=0;i<cM.r;i++){
310     for(j=0;j<cM.c;j++){
311       cM.mat[IDX(i,j,cM.ld)]=1;
312     }
313   }
314 }
315
316
317 void computeNNs(matrix dx, intMatrix dxMap, matrix dq, unint *dqMap, compPlan dcP, unint *NNs, real *NNdists, unint compLength){
318   real *dNNdists;
319   unint *dMinIDs;
320   
321   checkErr( cudaMalloc((void**)&dNNdists,compLength*sizeof(*dNNdists)) );
322   checkErr( cudaMalloc((void**)&dMinIDs,compLength*sizeof(*dMinIDs)) );
323
324   planNNWrap(dq, dqMap, dx, dxMap, dNNdists, dMinIDs, dcP, compLength );
325   cudaMemcpy( NNs, dMinIDs, dq.r*sizeof(*NNs), cudaMemcpyDeviceToHost );
326   cudaMemcpy( NNdists, dNNdists, dq.r*sizeof(*dNNdists), cudaMemcpyDeviceToHost );
327
328   cudaFree(dNNdists);
329   cudaFree(dMinIDs);
330 }
331
332
333 void computeKNNs(matrix dx, intMatrix dxMap, matrix dq, unint *dqMap, compPlan dcP, intMatrix NNs, matrix NNdists, unint compLength){
334   matrix dNNdists;
335   intMatrix dMinIDs;
336   dNNdists.r=compLength; dNNdists.pr=compLength; dNNdists.c=KMAX; dNNdists.pc=KMAX; dNNdists.ld=dNNdists.pc;
337   dMinIDs.r=compLength; dMinIDs.pr=compLength; dMinIDs.c=KMAX; dMinIDs.pc=KMAX; dMinIDs.ld=dMinIDs.pc;
338
339   checkErr( cudaMalloc((void**)&dNNdists.mat,dNNdists.pr*dNNdists.pc*sizeof(*dNNdists.mat)) );
340   checkErr( cudaMalloc((void**)&dMinIDs.mat,dMinIDs.pr*dMinIDs.pc*sizeof(*dMinIDs.mat)) );
341
342   planKNNWrap(dq, dqMap, dx, dxMap, dNNdists, dMinIDs, dcP, compLength);
343   cudaMemcpy( NNs.mat, dMinIDs.mat, dq.r*KMAX*sizeof(*NNs.mat), cudaMemcpyDeviceToHost );
344   cudaMemcpy( NNdists.mat, dNNdists.mat, dq.r*KMAX*sizeof(*NNdists.mat), cudaMemcpyDeviceToHost );
345
346   cudaFree(dNNdists.mat);
347   cudaFree(dMinIDs.mat);
348 }
349
350
351 //This calls the dist1Kernel wrapper, but has it compute only 
352 //a submatrix of the all-pairs distance matrix.  In particular,
353 //only distances from dr[start,:].. dr[start+length-1] to all of x
354 //are computed, resulting in a distance matrix of size 
355 //length by dx.pr.  It is assumed that length is padded.
356 void distSubMat(matrix dr, matrix dx, matrix dD, unint start, unint length){
357   dr.r=dr.pr=length;
358   dr.mat = &dr.mat[IDX( start, 0, dr.ld )];
359   dist1Wrap(dr, dx, dD);
360 }
361
362
363 void destroyRBC(rbcStruct *rbcS){
364   cudaFree(rbcS->dx.mat);
365   cudaFree(rbcS->dxMap.mat);
366   cudaFree(rbcS->dr.mat);
367   free(rbcS->groupCount);
368 }
369
370
371 /* Danger: this function allocates memory that it does not free.  
372  * Use freeCompPlan to clear mem.  
373  * See the readme.txt file for a description of why this function is needed.
374  */
375 void initCompPlan(compPlan *dcP, charMatrix cM, unint *groupCountQ, unint *groupCountX, unint numReps){
376   unint i,j,k;
377   unint maxNumGroups=0;
378   compPlan cP;
379   
380   unint sNumGroups = numReps;
381   cP.numGroups = (unint*)calloc(sNumGroups, sizeof(*cP.numGroups));
382   
383   for(i=0; i<numReps; i++){
384     cP.numGroups[i] = 0;
385     for(j=0; j<numReps; j++)
386       cP.numGroups[i] += cM.mat[IDX(i,j,cM.ld)];
387     maxNumGroups = MAX(cP.numGroups[i], maxNumGroups);
388   }
389   cP.ld = maxNumGroups;
390   
391   unint sQToQGroup;
392   for(i=0, sQToQGroup=0; i<numReps; i++)
393     sQToQGroup += PAD(groupCountQ[i]);
394   
395   cP.qToQGroup = (unint*)calloc( sQToQGroup, sizeof(*cP.qToQGroup) );
396
397   for(i=0, k=0; i<numReps; i++){
398     for(j=0; j<PAD(groupCountQ[i]); j++)
399       cP.qToQGroup[k++] = i;
400   }
401   
402   unint sQGroupToXGroup = numReps*maxNumGroups;
403   cP.qGroupToXGroup = (unint*)calloc( sQGroupToXGroup, sizeof(*cP.qGroupToXGroup) );
404   unint sGroupCountX = maxNumGroups*numReps;
405   cP.groupCountX = (unint*)calloc( sGroupCountX, sizeof(*cP.groupCountX) );
406   
407   for(i=0; i<numReps; i++){
408     for(j=0, k=0; j<numReps; j++){
409       if( cM.mat[IDX( i, j, cM.ld )] ){
410         cP.qGroupToXGroup[IDX( i, k, cP.ld )] = j;
411         cP.groupCountX[IDX( i, k++, cP.ld )] = groupCountX[j];
412       }
413     }
414   }
415
416   //Move to device
417   checkErr( cudaMalloc( (void**)&dcP->numGroups, sNumGroups*sizeof(*dcP->numGroups) ) );
418   cudaMemcpy( dcP->numGroups, cP.numGroups, sNumGroups*sizeof(*dcP->numGroups), cudaMemcpyHostToDevice );
419   checkErr( cudaMalloc( (void**)&dcP->groupCountX, sGroupCountX*sizeof(*dcP->groupCountX) ) );
420   cudaMemcpy( dcP->groupCountX, cP.groupCountX, sGroupCountX*sizeof(*dcP->groupCountX), cudaMemcpyHostToDevice );
421   checkErr( cudaMalloc( (void**)&dcP->qToQGroup, sQToQGroup*sizeof(*dcP->qToQGroup) ) );
422   cudaMemcpy( dcP->qToQGroup, cP.qToQGroup, sQToQGroup*sizeof(*dcP->qToQGroup), cudaMemcpyHostToDevice );
423   checkErr( cudaMalloc( (void**)&dcP->qGroupToXGroup, sQGroupToXGroup*sizeof(*dcP->qGroupToXGroup) ) );
424   cudaMemcpy( dcP->qGroupToXGroup, cP.qGroupToXGroup, sQGroupToXGroup*sizeof(*dcP->qGroupToXGroup), cudaMemcpyHostToDevice );
425   dcP->ld = cP.ld;
426
427   free(cP.numGroups);
428   free(cP.groupCountX);
429   free(cP.qToQGroup);
430   free(cP.qGroupToXGroup);
431 }
432
433
434 //Frees memory allocated in initCompPlan.
435 void freeCompPlan(compPlan *dcP){
436   cudaFree(dcP->numGroups);
437   cudaFree(dcP->groupCountX);
438   cudaFree(dcP->qToQGroup);
439   cudaFree(dcP->qGroupToXGroup);
440 }
441
442 #endif