matching.cu

#include "cudaSift.h"
#include "cudautils.h"

//================= Device matching functions =====================//

__global__ void MatchSiftPoints(SiftPoint *sift1, SiftPoint *sift2, float *corrData, int numPts1, int numPts2)
{
  __shared__ float siftPoint[128];
  __shared__ float sums[16*16];
  const int tx = threadIdx.x;
  const int ty = threadIdx.y;
  const int p1 = blockIdx.x;
  const int p2 = blockIdx.y*16 + ty;
  const float *ptr1 = sift1[p1].data;
  const float *ptr2 = sift2[p2].data;
  const int i = 16*ty + tx;
  if (ty<8)
    siftPoint[i] = ptr1[i];
  __syncthreads();
  float sum = 0.0f;
  if (p2<numPts2)
    for (int j=0;j<8;j++)
      sum += siftPoint[16*j+tx] * ptr2[16*j+tx];
  sums[i] = sum;
  __syncthreads();
  if (tx<8)
    sums[i] += sums[i+8];
  __syncthreads();
  if (tx<4)
    sums[i] += sums[i+4];
  __syncthreads();
  if (ty==0) {
    sum = sums[16*tx+0] + sums[16*tx+1] + sums[16*tx+2] + sums[16*tx+3];
    corrData[p1*gridDim.y*16 + blockIdx.y*16 + tx] = sum;
  }
  __syncthreads();
}

__global__ void MatchSiftPoints2(SiftPoint *sift1, SiftPoint *sift2, float *corrData, int numPts1, int numPts2)
{
  __shared__ float siftPoints1[16*128];
  __shared__ float siftPoints2[16*128];
  const int tx = threadIdx.x;
  const int ty = threadIdx.y;
  const float *ptr1 = sift1[min(numPts1-1,blockIdx.x*16 + ty)].data;
  const float *ptr2 = sift2[min(numPts2-1,blockIdx.y*16 + ty)].data;
  for (int i=0;i<8;i++) {
    siftPoints1[128*ty+16*i+tx] = ptr1[16*i+tx];
    siftPoints2[128*ty+16*i+tx] = ptr2[16*i+tx];
  }
  __syncthreads();
  const int p1 = blockIdx.x*16 + ty;
  const int p2 = blockIdx.y*16 + tx;
  const float *pt1 = &siftPoints1[ty*128];
  const float *pt2 = &siftPoints2[tx*128];
  float sum = 0.0f;
  for (int i=0;i<128;i++) {
    int itx = (i + tx)&127; // avoid bank conflicts
    sum += pt1[itx]*pt2[itx];
  }
  if (p1<numPts1)
    corrData[p1*gridDim.y*16 + p2] = (p2<numPts2 ? sum : -1.0f);
}

__global__ void FindMaxCorr(float *corrData, SiftPoint *sift1, SiftPoint *sift2, int numPts1, int corrWidth, int siftSize)
{
  __shared__ float maxScore[16*16];
  __shared__ float maxScor2[16*16];
  __shared__ int maxIndex[16*16];
  const int tx = threadIdx.x;
  const int ty = threadIdx.y;
  const int idx = ty*16 + tx;
  int p1 = blockIdx.x*16 + threadIdx.y;
  p1 = (p1>=numPts1 ? numPts1-1 : p1);
  maxScore[idx] = -1.0f;
  maxScor2[idx] = -1.0f;
  maxIndex[idx] = -1;
  __syncthreads();
  float *corrs = &corrData[p1*corrWidth];
  for (int i=tx;i<corrWidth;i+=16) {
    float val = corrs[i];
    if (val>maxScore[idx]) {
      maxScor2[idx] = maxScore[idx];
      maxScore[idx] = val;
      maxIndex[idx] = i;
    } else if (val>maxScor2[idx])
      maxScor2[idx] = val;
  }
  __syncthreads();
  for (int len=8;len>0;len/=2) {
    if (tx<8) {
      float val = maxScore[idx+len];
      int i = maxIndex[idx+len];
      if (val>maxScore[idx]) {
	maxScor2[idx] = maxScore[idx];
	maxScore[idx] = val;
	maxIndex[idx] = i;
      } else if (val>maxScor2[idx])
	maxScor2[idx] = val;
      float va2 = maxScor2[idx+len];
      if (va2>maxScor2[idx])
	maxScor2[idx] = va2;
    }
    __syncthreads();
  }
  if (tx==0) {
    sift1[p1].score = maxScore[ty*16];
    sift1[p1].ambiguity = maxScor2[ty*16] / (maxScore[ty*16] + 1e-6);
    sift1[p1].match = maxIndex[ty*16];
    sift1[p1].match_xpos = sift2[maxIndex[ty*16]].xpos;
    sift1[p1].match_ypos = sift2[maxIndex[ty*16]].ypos;
  }
}

// Version based on suggestion by Nicholas Lin
__global__ void FindMaxCorr3(float *corrData, SiftPoint *sift1, SiftPoint *sift2, int numPts1, int numPts2)
{
  int block_dim = blockDim.x; // blockDim.x == 16
  const int tx = threadIdx.x;
  const int ty = threadIdx.y;
  const int p1 = blockIdx.x * block_dim + ty;
  const int idx = ty * 16 + tx;
  
  __shared__ int maxIndex[16 * 16];
  maxIndex[idx] = 0;
  __syncthreads();
  
  float *corrs = NULL;
  if (p1 < numPts1) {
    corrs = &corrData[p1 * block_dim * 2];
    corrs[tx] = 0.0f;
    corrs[tx + 16] = 0.0f;
    const float *pt1 = sift1[p1].data;
    for (int p2 = tx; p2 < numPts2; p2 += 16) {
      float *pt2 = sift2[p2].data;
      float sum = 0.0f;
      for (int i = 0; i < 128; i++) 
	sum += pt1[i] * pt2[i];
      if (sum > corrs[tx]) {
	corrs[tx + 16] = corrs[tx];
	corrs[tx] = sum;
	maxIndex[idx] = p2;
      }
      else if (sum > corrs[tx + 16])
	corrs[tx + 16] = sum;
    }
  }
  __syncthreads();
  if (p1 < numPts1) {
    for (int len = 8; len > 0; len /= 2) {
      if (tx < len) {
	float val = corrs[tx + len];
	int i = maxIndex[idx + len];
	if (val > corrs[tx]) {
	  corrs[tx + 16] = corrs[tx];
	  corrs[tx] = val;
	  maxIndex[idx] = i;
	}
	else if (val > corrs[tx + 16])
	  corrs[tx + 16] = val;
	float va2 = corrs[tx + 16 + len];
	if (va2 > corrs[tx + 16])
	  corrs[tx + 16] = va2;
      }
      __syncthreads();
    }
    if (tx==0) {
      sift1[p1].score = corrs[0];
      sift1[p1].ambiguity = corrs[16] / (corrs[0] + 1e-6);
      sift1[p1].match = maxIndex[ty << 4];
      sift1[p1].match_xpos = sift2[maxIndex[ty << 4]].xpos;
      sift1[p1].match_ypos = sift2[maxIndex[ty << 4]].ypos;
    }
  }
}

#define FMC2W 16
#define FMC2H 4

__global__ void FindMaxCorr2(SiftPoint *sift1, SiftPoint *sift2, int numPts1, int numPts2)
{
  __shared__ float siftPoint[128];
  __shared__ float maxScore[FMC2H]; 
  __shared__ float maxScor2[FMC2H]; 
  __shared__ int maxIndex[FMC2H]; 
  const int p1 = blockIdx.x;
  if (p1>=numPts1)
    return;
  const int tx = threadIdx.x;
  const int ty = threadIdx.y;
  const int idx = ty*FMC2W + tx;
  if (idx<FMC2H) {
    maxScore[idx] = -1.0f;
    maxScor2[idx] = -1.0f;
    maxIndex[idx] = 0;
  }
  __syncthreads();
  const float *pt1 = sift1[p1].data;
  for (int i=idx;i<128;i+=FMC2W*FMC2H)
    siftPoint[i] = pt1[i];
  __syncthreads();
  for (int p2=ty;p2<numPts2;p2+=FMC2H) {
    const float *pt2 = sift2[p2].data;
    float sum = 0.0f;
    for (int j=tx;j<128;j+=FMC2W)
      sum += siftPoint[j] * pt2[j];
    for (int j=FMC2W/2;j>0;j/=2)
      sum += ShiftDown(sum, j);
    if (tx==0) {
      if (sum>maxScore[ty]) {
	maxScor2[ty] = maxScore[ty];
	maxScore[ty] = sum;
	maxIndex[ty] = p2;
      } else if (sum>maxScor2[ty])
	maxScor2[ty] = sum;
    }
  }
  __syncthreads();
  for (int len=FMC2H/2;len>0;len/=2) {
    if (ty==0 && tx<len) {
      float val = maxScore[tx+len];
      int p2 = maxIndex[tx+len];
      if (val>maxScore[tx]) {
	maxScor2[tx] = maxScore[tx];
	maxScore[tx] = val;
	maxIndex[tx] = p2;
      } else if (val>maxScor2[tx])
	maxScor2[tx] = val;
      float va2 = maxScor2[tx+len];
      if (va2>maxScor2[tx])
	maxScor2[tx] = va2;
    }
    __syncthreads();
  }
  if (ty==0 && tx==0) {
    sift1[p1].score = maxScore[0];
    sift1[p1].ambiguity = maxScor2[0] / (maxScore[0] + 1e-6);
    sift1[p1].match = maxIndex[0];
    sift1[p1].match_xpos = sift2[maxIndex[0]].xpos;
    sift1[p1].match_ypos = sift2[maxIndex[0]].ypos;
  }
}

__global__ void FindMaxCorr4(SiftPoint *sift1, SiftPoint *sift2, int numPts1, int numPts2)
{
  __shared__ float siftPoint[128*FMC2H];
  __shared__ float maxScore[FMC2H]; 
  __shared__ float maxScor2[FMC2H]; 
  __shared__ int maxIndex[FMC2H]; 
  const int tx = threadIdx.x;
  const int ty = threadIdx.y;
  if (tx==0) {
    maxScore[ty] = -1.0f;
    maxScor2[ty] = -1.0f;
    maxIndex[ty] = 0;
  }
  const int p1 = blockIdx.x*FMC2H + ty;
  const float *pt1 = sift1[p1].data;
  for (int j=tx;j<128;j+=FMC2W)
    siftPoint[128*ty + j] = pt1[j];
  __syncthreads();
  for (int p2=0;p2<numPts2;p2++) {
    const float *pt2 = sift2[p2].data;
    float sum = 0.0f;
    for (int j=tx;j<128;j+=FMC2W)
      sum += siftPoint[128*ty + j] * pt2[j]; 
    for (int j=FMC2W/2;j>0;j/=2)
      sum += ShiftDown(sum, j);
    if (tx==0) {
      if (sum>maxScore[ty]) {
	maxScor2[ty] = maxScore[ty];
	maxScore[ty] = sum;
	maxIndex[ty] = p2;
      } else if (sum>maxScor2[ty])
	maxScor2[ty] = sum;
    }
  }
  __syncthreads();
  if (tx==0) {
    sift1[p1].score = maxScore[ty];
    sift1[p1].ambiguity = maxScor2[ty] / (maxScore[ty] + 1e-6);
    sift1[p1].match = maxIndex[ty];
    sift1[p1].match_xpos = sift2[maxIndex[ty]].xpos;
    sift1[p1].match_ypos = sift2[maxIndex[ty]].ypos;
  }
}


__global__ void CleanMatches(SiftPoint *sift1, int numPts1)
{
  const int p1 = min(blockIdx.x*64 + threadIdx.x, numPts1-1);
  sift1[p1].score = 0.0f;
}

#define M7W   32
#define M7H   32
#define M7R    4
#define NRX    2
#define NDIM 128

__global__ void FindMaxCorr10(SiftPoint *sift1, SiftPoint *sift2, int numPts1, int numPts2)
{
  __shared__ float4 buffer1[M7W*NDIM/4]; 
  __shared__ float4 buffer2[M7H*NDIM/4];       
  int tx = threadIdx.x;
  int ty = threadIdx.y;
  int bp1 = M7W*blockIdx.x;
  for (int j=ty;j<M7W;j+=M7H/M7R) {    
    int p1 = min(bp1 + j, numPts1 - 1);
    for (int d=tx;d<NDIM/4;d+=M7W)
      buffer1[j*NDIM/4 + (d + j)%(NDIM/4)] = ((float4*)&sift1[p1].data)[d];
  }
      
  float max_score[NRX];
  float sec_score[NRX];
  int index[NRX];
  for (int i=0;i<NRX;i++) {
    max_score[i] = 0.0f;
    sec_score[i] = 0.0f;
    index[i] = -1;
  }
  int idx = ty*M7W + tx;
  int ix = idx%(M7W/NRX);
  int iy = idx/(M7W/NRX);
  for (int bp2=0;bp2<numPts2 - M7H + 1;bp2+=M7H) {
    for (int j=ty;j<M7H;j+=M7H/M7R) {      
      int p2 = min(bp2 + j, numPts2 - 1);
      for (int d=tx;d<NDIM/4;d+=M7W)
	buffer2[j*NDIM/4 + d] = ((float4*)&sift2[p2].data)[d];
    }
    __syncthreads();

    if (idx<M7W*M7H/M7R/NRX) {
      float score[M7R][NRX];                                    
      for (int dy=0;dy<M7R;dy++)
	for (int i=0;i<NRX;i++)
	  score[dy][i] = 0.0f;
      for (int d=0;d<NDIM/4;d++) {
	float4 v1[NRX];
	for (int i=0;i<NRX;i++) 
	  v1[i] = buffer1[((M7W/NRX)*i + ix)*NDIM/4 + (d + (M7W/NRX)*i + ix)%(NDIM/4)];
	for (int dy=0;dy<M7R;dy++) {
	  float4 v2 = buffer2[(M7R*iy + dy)*(NDIM/4) + d];    
	  for (int i=0;i<NRX;i++) {
	    score[dy][i] += v1[i].x*v2.x;
	    score[dy][i] += v1[i].y*v2.y;
	    score[dy][i] += v1[i].z*v2.z;
	    score[dy][i] += v1[i].w*v2.w;
	  }
	}
      }
      for (int dy=0;dy<M7R;dy++) {
	for (int i=0;i<NRX;i++) {
	  if (score[dy][i]>max_score[i]) {
	    sec_score[i] = max_score[i];
	    max_score[i] = score[dy][i];     
	    index[i] = min(bp2 + M7R*iy + dy, numPts2-1);
	  } else if (score[dy][i]>sec_score[i])
	    sec_score[i] = score[dy][i]; 
	}
      }
    }
    __syncthreads();
  }

  float *scores1 = (float*)buffer1;
  float *scores2 = &scores1[M7W*M7H/M7R];
  int *indices = (int*)&scores2[M7W*M7H/M7R];
  if (idx<M7W*M7H/M7R/NRX) {
    for (int i=0;i<NRX;i++) {
      scores1[iy*M7W + (M7W/NRX)*i + ix] = max_score[i];  
      scores2[iy*M7W + (M7W/NRX)*i + ix] = sec_score[i];  
      indices[iy*M7W + (M7W/NRX)*i + ix] = index[i];
    }
  }
  __syncthreads();
  
  if (ty==0) {
    float max_score = scores1[tx];
    float sec_score = scores2[tx];
    int index = indices[tx];
    for (int y=0;y<M7H/M7R;y++)
      if (index != indices[y*M7W + tx]) {
	if (scores1[y*M7W + tx]>max_score) {
	  sec_score = max(max_score, sec_score);
	  max_score = scores1[y*M7W + tx]; 
	  index = indices[y*M7W + tx];
	} else if (scores1[y*M7W + tx]>sec_score)
	  sec_score = scores1[y*M7W + tx];
      }
    sift1[bp1 + tx].score = max_score;
    sift1[bp1 + tx].match = index;
    sift1[bp1 + tx].match_xpos = sift2[index].xpos;
    sift1[bp1 + tx].match_ypos = sift2[index].ypos;
    sift1[bp1 + tx].ambiguity = sec_score / (max_score + 1e-6f);
  }
}
  
#define FMC_GH  512
#define FMC_BW   32
#define FMC_BH   32
#define FMC_BD   16
#define FMC_TW    1
#define FMC_TH    4
#define FMC_NW   (FMC_BW/FMC_TW)   //  32
#define FMC_NH   (FMC_BH/FMC_TH)   //   8
#define FMC_NT   (FMC_NW*FMC_NH)   // 256 = 8 warps

__device__ volatile int lock = 0;

__global__ void FindMaxCorr9(SiftPoint *sift1, SiftPoint *sift2, int numPts1, int numPts2)
{
  __shared__ float4 siftParts1[FMC_BW*FMC_BD]; // 4*32*8 = 1024
  __shared__ float4 siftParts2[FMC_BH*FMC_BD]; // 4*32*8 = 1024
  //__shared__ float blksums[FMC_BW*FMC_BH];     // 32*32  = 1024
  const int tx = threadIdx.x;
  const int ty = threadIdx.y;
  const int idx = ty*FMC_NW + tx;
  float4 *pts1 = 0, *pts2 = 0;
  if (idx<FMC_BW) {
    const int p1l = min(blockIdx.x*FMC_BW + idx, numPts1-1);
    pts1 = (float4*)sift1[p1l].data;
  }
  float maxScore = -1.0f;
  float maxScor2 = -1.0f;
  int maxIndex = 0;
  for (int k=0;k<min(FMC_GH, numPts2 - FMC_BH + 1);k+=FMC_BH) {
    if (idx<FMC_BH) {
      const int p2l = min(blockIdx.y*FMC_GH + k + idx, numPts2-1);
      pts2 = (float4*)sift2[p2l].data;
    }
    float sums[FMC_TW*FMC_TH];
    for (int i=0;i<FMC_TW*FMC_TH;i++) 
      sums[i] = 0.0f;

    if (idx<FMC_BW)
      for (int i=0;i<FMC_BD/2;i++) 
	siftParts1[(i + 0)*FMC_BW + idx] = pts1[0 + i];
    if (idx<FMC_BH)
      for (int i=0;i<FMC_BD/2;i++) 
	siftParts2[(i + 0)*FMC_BH + idx] = pts2[0 + i];
    __syncthreads();
    
    int b = FMC_BD/2;
    for (int d=FMC_BD/2;d<32;d+=FMC_BD/2) {
      if (idx<FMC_BW)
	for (int i=0;i<FMC_BD/2;i++) 
	  siftParts1[(i + b)*FMC_BW + idx] = pts1[d + i];
      if (idx<FMC_BH)
	for (int i=0;i<FMC_BD/2;i++) 
	  siftParts2[(i + b)*FMC_BH + idx] = pts2[d + i];

      b ^= FMC_BD/2;
      for (int i=0;i<FMC_BD/2;i++) {
	float4 v1[FMC_TW];
	for (int ix=0;ix<FMC_TW;ix++)
	  v1[ix] = siftParts1[(i + b)*FMC_BW + (tx*FMC_TW + ix)];
	for (int iy=0;iy<FMC_TH;iy++) {
	  float4 v2 = siftParts2[(i + b)*FMC_BH + (ty*FMC_TH + iy)];
	  for (int ix=0;ix<FMC_TW;ix++) {
	    sums[iy*FMC_TW + ix] += v1[ix].x * v2.x;
	    sums[iy*FMC_TW + ix] += v1[ix].y * v2.y;
	    sums[iy*FMC_TW + ix] += v1[ix].z * v2.z;
	    sums[iy*FMC_TW + ix] += v1[ix].w * v2.w;
	  }
	}
      }
      __syncthreads();
    }
    
    b ^= FMC_BD/2;
    for (int i=0;i<FMC_BD/2;i++) {
      float4 v1[FMC_TW];
      for (int ix=0;ix<FMC_TW;ix++)
	v1[ix] = siftParts1[(i + b)*FMC_BW + (tx*FMC_TW + ix)];
      for (int iy=0;iy<FMC_TH;iy++) {
	float4 v2 = siftParts2[(i + b)*FMC_BH + (ty*FMC_TH + iy)];
	for (int ix=0;ix<FMC_TW;ix++) {
	  sums[iy*FMC_TW + ix] += v1[ix].x * v2.x;
	  sums[iy*FMC_TW + ix] += v1[ix].y * v2.y;
	  sums[iy*FMC_TW + ix] += v1[ix].z * v2.z;
	  sums[iy*FMC_TW + ix] += v1[ix].w * v2.w;
	}
      }
    }
    __syncthreads();
    
    float *blksums = (float*)siftParts1;
    for (int iy=0;iy<FMC_TH;iy++) 
      for (int ix=0;ix<FMC_TW;ix++) 
	blksums[(ty*FMC_TH + iy)*FMC_BW + (tx*FMC_TW + ix)] = sums[iy*FMC_TW + ix];
    __syncthreads();
    if (idx<FMC_BW) { 
      for (int j=0;j<FMC_BH;j++) {
	float sum = blksums[j*FMC_BW + idx];
	if (sum>maxScore) { 
	  maxScor2 = maxScore;
	  maxScore = sum;
	  maxIndex = min(blockIdx.y*FMC_GH + k + j, numPts2-1);
	} else if (sum>maxScor2)
	  maxScor2 = sum;
      }
    }
    __syncthreads();
  }
  const int p1 = min(blockIdx.x*FMC_BW + idx, numPts1-1);
  if (idx==0)
    while (atomicCAS((int *)&lock, 0, 1) != 0);
  __syncthreads();
  if (idx<FMC_BW) {
    float maxScor2Old = sift1[p1].ambiguity*(sift1[p1].score + 1e-6f);
    if (maxScore>sift1[p1].score) {
      maxScor2 = max(sift1[p1].score, maxScor2);
      sift1[p1].ambiguity = maxScor2 / (maxScore + 1e-6f);
      sift1[p1].score = maxScore;
      sift1[p1].match = maxIndex;
      sift1[p1].match_xpos = sift2[maxIndex].xpos;
      sift1[p1].match_ypos = sift2[maxIndex].ypos;
    } else if (maxScore>maxScor2Old)
      sift1[p1].ambiguity = maxScore / (sift1[p1].score + 1e-6f);
  }
  __syncthreads();
  if (idx==0)
    atomicExch((int* )&lock, 0);
}

__global__ void FindMaxCorr8(SiftPoint *sift1, SiftPoint *sift2, int numPts1, int numPts2)
{
  __shared__ float4 siftParts1[FMC_BW*FMC_BD]; // 4*32*8 = 1024
  __shared__ float4 siftParts2[FMC_BH*FMC_BD]; // 4*32*8 = 1024
  __shared__ float blksums[FMC_BW*FMC_BH];     // 32*32  = 1024
  const int tx = threadIdx.x;
  const int ty = threadIdx.y;
  const int idx = ty*FMC_NW + tx;
  float4 *pts1 = 0, *pts2 = 0;
  if (idx<FMC_BW) {
    const int p1l = min(blockIdx.x*FMC_BW + idx, numPts1-1);
    pts1 = (float4*)sift1[p1l].data;
  }
  float maxScore = -1.0f;
  float maxScor2 = -1.0f;
  int maxIndex = 0;
  for (int k=0;k<min(FMC_GH, numPts2 - FMC_BH + 1);k+=FMC_BH) {
    if (idx<FMC_BH) {
      const int p2l = min(blockIdx.y*FMC_GH + k + idx, numPts2-1);
      pts2 = (float4*)sift2[p2l].data;
    }
    float sums[FMC_TW*FMC_TH];
    for (int i=0;i<FMC_TW*FMC_TH;i++) 
      sums[i] = 0.0f;
    for (int d=0;d<32;d+=FMC_BD) {
      if (idx<FMC_BW)
	for (int i=0;i<FMC_BD;i++) 
	  siftParts1[i*FMC_BW + idx] = pts1[d + i];
      if (idx<FMC_BH)
	for (int i=0;i<FMC_BD;i++) 
	  siftParts2[i*FMC_BH + idx] = pts2[d + i];
      __syncthreads();
      
      for (int i=0;i<FMC_BD;i++) {
	float4 v1[FMC_TW];
	for (int ix=0;ix<FMC_TW;ix++)
	  v1[ix] = siftParts1[i*FMC_BW + (tx*FMC_TW + ix)];
	for (int iy=0;iy<FMC_TH;iy++) {
	  float4 v2 = siftParts2[i*FMC_BH + (ty*FMC_TH + iy)];
	  for (int ix=0;ix<FMC_TW;ix++) {
	    sums[iy*FMC_TW + ix] += v1[ix].x * v2.x;
	    sums[iy*FMC_TW + ix] += v1[ix].y * v2.y;
	    sums[iy*FMC_TW + ix] += v1[ix].z * v2.z;
	    sums[iy*FMC_TW + ix] += v1[ix].w * v2.w;
	  }
	}
      }
      __syncthreads();
    }
    //float *blksums = (float*)siftParts1;
    for (int iy=0;iy<FMC_TH;iy++) 
      for (int ix=0;ix<FMC_TW;ix++) 
	blksums[(ty*FMC_TH + iy)*FMC_BW + (tx*FMC_TW + ix)] = sums[iy*FMC_TW + ix];
    __syncthreads();
    if (idx<FMC_BW) { 
      for (int j=0;j<FMC_BH;j++) {
	float sum = blksums[j*FMC_BW + idx];
	if (sum>maxScore) { 
	  maxScor2 = maxScore;
	  maxScore = sum;
	  maxIndex = min(blockIdx.y*FMC_GH + k + j, numPts2-1);
	} else if (sum>maxScor2)
	  maxScor2 = sum;
      }
    }
    __syncthreads();
  }
  const int p1 = min(blockIdx.x*FMC_BW + idx, numPts1-1);
  if (idx==0)
    while (atomicCAS((int *)&lock, 0, 1) != 0);
  __syncthreads();
  if (idx<FMC_BW) {
    float maxScor2Old = sift1[p1].ambiguity*(sift1[p1].score + 1e-6f);
    if (maxScore>sift1[p1].score) {
      maxScor2 = max(sift1[p1].score, maxScor2);
      sift1[p1].ambiguity = maxScor2 / (maxScore + 1e-6f);
      sift1[p1].score = maxScore;
      sift1[p1].match = maxIndex;
      sift1[p1].match_xpos = sift2[maxIndex].xpos;
      sift1[p1].match_ypos = sift2[maxIndex].ypos;
    } else if (maxScore>maxScor2Old)
      sift1[p1].ambiguity = maxScore / (sift1[p1].score + 1e-6f);
  }
  __syncthreads();
  if (idx==0)
    atomicExch((int* )&lock, 0);
}

__global__ void FindMaxCorr7(SiftPoint *sift1, SiftPoint *sift2, int numPts1, int numPts2)
{
  __shared__ float siftParts1[17*64]; // features in columns
  __shared__ float siftParts2[16*64]; // one extra to avoid shared conflicts
  float4 *pts1 = (float4*)siftParts1;
  float4 *pts2 = (float4*)siftParts2;
  const int tx = threadIdx.x;
  const int ty = threadIdx.y;
  const int p1l = min(blockIdx.x*16 + ty, numPts1-1);
  const float4 *p1l4 = (float4*)sift1[p1l].data;
  float maxScore = -1.0f;
  float maxScor2 = -1.0f;
  int maxIndex = 0;
  for (int k=0;k<512/16;k++) {
    const int p2l = min(blockIdx.y*512 + k*16 + ty, numPts2-1);
    const float4 *p2l4 = (float4*)sift2[p2l].data;
#define NUM 4
    float sum[NUM];
    if (ty<(16/NUM))
      for (int l=0;l<NUM;l++)
	sum[l] = 0.0f;
    __syncthreads();
    for (int i=0;i<2;i++) {
      pts1[17*tx + ty] = p1l4[i*16 + tx];
      pts2[16*ty + tx] = p2l4[i*16 + tx];
      __syncthreads(); 
      if (ty<(16/NUM)) {
#pragma unroll
	for (int j=0;j<16;j++) {
	  float4 p1v = pts1[17* j + tx];
#pragma unroll
	  for (int l=0;l<NUM;l++) {
	    float4 p2v = pts2[16*(ty + l*(16/NUM)) +  j];
	    sum[l] += p1v.x * p2v.x;
	    sum[l] += p1v.y * p2v.y;
	    sum[l] += p1v.z * p2v.z;
	    sum[l] += p1v.w * p2v.w;
	  }
	}
      }
      __syncthreads();
    }
    float *sums = siftParts1;
    if (ty<(16/NUM))
      for (int l=0;l<NUM;l++) 
	sums[16*(ty + l*(16/NUM)) + tx] = sum[l];
    __syncthreads();
    if (ty==0) { 
      for (int j=0;j<16;j++) {
	float sum = sums[16*j + tx];
	if (sum>maxScore) { 
	  maxScor2 = maxScore;
	  maxScore = sum;
	  maxIndex = min(blockIdx.y*512 +  k*16 + j, numPts2-1);
	} else if (sum>maxScor2)
	  maxScor2 = sum;
      }
    }
    __syncthreads();
  }
  const int p1 = min(blockIdx.x*16 + tx, numPts1-1);
  if (tx==0 && ty==0)
    while (atomicCAS((int *)&lock, 0, 1) != 0);
  __syncthreads();
  if (ty==0) {
    float maxScor2Old = sift1[p1].ambiguity*(sift1[p1].score + 1e-6f);
    if (maxScore>sift1[p1].score) {
      maxScor2 = max(sift1[p1].score, maxScor2);
      sift1[p1].ambiguity = maxScor2 / (maxScore + 1e-6f);
      sift1[p1].score = maxScore;
      sift1[p1].match = maxIndex;
      sift1[p1].match_xpos = sift2[maxIndex].xpos;
      sift1[p1].match_ypos = sift2[maxIndex].ypos;
    } else if (maxScore>maxScor2Old)
      sift1[p1].ambiguity = maxScore / (sift1[p1].score + 1e-6f);
  }
  __syncthreads();
  if (tx==0 && ty==0)
    atomicExch((int* )&lock, 0);
}

__global__ void FindMaxCorr6(SiftPoint *sift1, SiftPoint *sift2, int numPts1, int numPts2)
{
  //__shared__ float siftParts1[128*16]; // features in columns
  __shared__ float siftParts2[128*16]; // one extra to avoid shared conflicts
  __shared__ float sums[16*16];
  const int tx = threadIdx.x;
  const int ty = threadIdx.y;
  const int p1l = min(blockIdx.x*16 + ty, numPts1-1);
  float *pt1l = sift1[p1l].data;
  float4 part1 = reinterpret_cast<float4*>(pt1l)[tx];
  float maxScore = -1.0f;
  float maxScor2 = -1.0f;
  int maxIndex = 0;
  for (int k=0;k<512;k+=16) {
    const int p2l = min(blockIdx.y*512 + k + ty, numPts2-1);
    float *pt2l = sift2[p2l].data;
    reinterpret_cast<float4*>(siftParts2)[32*ty + tx] = reinterpret_cast<float4*>(pt2l)[tx];
    __syncthreads();
    for (int i=0;i<16;i++) {
      float4 part2 = reinterpret_cast<float4*>(siftParts2)[32*i  + tx];
      float sum = part1.x*part2.x + part1.y*part2.y + part1.z*part2.z + part1.w*part2.w;
      sum += ShiftDown(sum, 16);
      sum += ShiftDown(sum, 8);
      sum += ShiftDown(sum, 4);
      sum += ShiftDown(sum, 2);
      sum += ShiftDown(sum, 1);
      if (tx==0)
	sums[16*i + ty] = sum;
    }
    __syncthreads();
    if (ty==0 && tx<16) { 
      for (int j=0;j<16;j++) {
	float sum = sums[16*j + tx];
	if (sum>maxScore) { 
	  maxScor2 = maxScore;
	  maxScore = sum;
	  maxIndex = min(blockIdx.y*512 +  k + j, numPts2-1);
	} else if (sum>maxScor2)
	  maxScor2 = sum;
      }
    }
    __syncthreads();
  }
  if (tx==0 && ty==0)
    while (atomicCAS((int *)&lock, 0, 1) != 0);
  __syncthreads();
  if (ty==0 && tx<16) {
    const int p1 = min(blockIdx.x*16 + tx, numPts1-1);
    float maxScor2Old = sift1[p1].ambiguity*(sift1[p1].score + 1e-6f);
    if (maxScore>sift1[p1].score) {
      maxScor2 = max(sift1[p1].score, maxScor2);
      sift1[p1].ambiguity = maxScor2 / (maxScore + 1e-6f);
      sift1[p1].score = maxScore;
      sift1[p1].match = maxIndex;
      sift1[p1].match_xpos = sift2[maxIndex].xpos;
      sift1[p1].match_ypos = sift2[maxIndex].ypos;
    } else if (maxScore>maxScor2Old)
      sift1[p1].ambiguity = maxScore / (sift1[p1].score + 1e-6f);
  }
  __syncthreads();
  if (tx==0 && ty==0)
    atomicExch((int* )&lock, 0);
}
 
__global__ void FindMaxCorr5(SiftPoint *sift1, SiftPoint *sift2, int numPts1, int numPts2)
{
  __shared__ float siftParts1[17*16]; // features in columns
  __shared__ float siftParts2[17*16]; // one extra to avoid shared conflicts
  const int tx = threadIdx.x;
  const int ty = threadIdx.y;
  const int p1l = min(blockIdx.x*16 + ty, numPts1-1);
  const float *pt1l = sift1[p1l].data;
  float maxScore = -1.0f;
  float maxScor2 = -1.0f;
  int maxIndex = 0;
  for (int k=0;k<512/16;k++) {
    const int p2l = min(blockIdx.y*512 + k*16 + ty, numPts2-1);
    const float *pt2l = sift2[p2l].data;
    float sum = 0.0f;
    for (int i=0;i<8;i++) {
      siftParts1[17*tx + ty] = pt1l[i*16 + tx]; // load and transpose
      siftParts2[17*tx + ty] = pt2l[i*16 + tx];
      __syncthreads();
      for (int j=0;j<16;j++)
	sum += siftParts1[17*j + tx] * siftParts2[17*j + ty];
      __syncthreads();
    }
    float *sums = siftParts1;
    sums[16*ty + tx] = sum;
    __syncthreads();
    if (ty==0) { 
      for (int j=0;j<16;j++) {
	float sum = sums[16*j + tx];
	if (sum>maxScore) { 
	  maxScor2 = maxScore;
	  maxScore = sum;
	  maxIndex = min(blockIdx.y*512 +  k*16 + j, numPts2-1);
	} else if (sum>maxScor2)
	  maxScor2 = sum;
      }
    }
    __syncthreads();
  }
  const int p1 = min(blockIdx.x*16 + tx, numPts1-1);
  if (tx==0 && ty==0)
    while (atomicCAS((int *)&lock, 0, 1) != 0);
  __syncthreads();
  if (ty==0) {
    float maxScor2Old = sift1[p1].ambiguity*(sift1[p1].score + 1e-6f);
    if (maxScore>sift1[p1].score) {
      maxScor2 = max(sift1[p1].score, maxScor2);
      sift1[p1].ambiguity = maxScor2 / (maxScore + 1e-6f);
      sift1[p1].score = maxScore;
      sift1[p1].match = maxIndex;
      sift1[p1].match_xpos = sift2[maxIndex].xpos;
      sift1[p1].match_ypos = sift2[maxIndex].ypos;
    } else if (maxScore>maxScor2Old)
      sift1[p1].ambiguity = maxScore / (sift1[p1].score + 1e-6f);
  }
  __syncthreads();
  if (tx==0 && ty==0)
    atomicExch((int* )&lock, 0);
}
 

template <int size>
__device__ void InvertMatrix(float elem[size][size], float res[size][size]) 
{  
  int indx[size];
  float b[size];
  float vv[size];
  for (int i=0;i<size;i++)
    indx[i] = 0;
  int imax = 0;
  float d = 1.0;
  for (int i=0;i<size;i++) { // find biggest element for each row
    float big = 0.0;
    for (int j=0;j<size;j++) {
      float temp = fabs(elem[i][j]); 
      if (temp>big) 
	big = temp;
    }
    if (big>0.0)
      vv[i] = 1.0/big;
    else
      vv[i] = 1e16;
  }
  for (int j=0;j<size;j++) { 
    for (int i=0;i<j;i++) { // i<j
      float sum = elem[i][j]; // i<j (lower left)
      for (int k=0;k<i;k++) // k<i<j
	sum -= elem[i][k]*elem[k][j]; // i>k (upper right), k<j (lower left)
      elem[i][j] = sum; // i<j (lower left)
    }
    float big = 0.0;
    for (int i=j;i<size;i++) { // i>=j
      float sum = elem[i][j]; // i>=j (upper right)
      for (int k=0;k<j;k++) // k<j<=i
	sum -= elem[i][k]*elem[k][j]; // i>k (upper right), k<j (lower left)
      elem[i][j] = sum; // i>=j (upper right)
      float dum = vv[i]*fabs(sum);
      if (dum>=big) {
	big = dum;
	imax = i;  
      }
    }
    if (j!=imax) { // imax>j
      for (int k=0;k<size;k++) {
	float dum = elem[imax][k]; // upper right and lower left
	elem[imax][k] = elem[j][k];
	elem[j][k] = dum;
      }
      d = -d;
      vv[imax] = vv[j];
    }
    indx[j] = imax;
    if (elem[j][j]==0.0)  // j==j (upper right)
      elem[j][j] = 1e-16;
    if (j!=(size-1)) {
      float dum = 1.0/elem[j][j];
      for (int i=j+1;i<size;i++) // i>j
	elem[i][j] *= dum; // i>j (upper right)
    }
  }
  for (int j=0;j<size;j++) {
    for (int k=0;k<size;k++) 
      b[k] = 0.0;  
    b[j] = 1.0;
    int ii = -1;
    for (int i=0;i<size;i++) {
      int ip = indx[i];
      float sum = b[ip];
      b[ip] = b[i];
      if (ii!=-1)
	for (int j=ii;j<i;j++) 
	  sum -= elem[i][j]*b[j]; // i>j (upper right)
      else if (sum!=0.0)
        ii = i;
      b[i] = sum;
    }
    for (int i=size-1;i>=0;i--) {
      float sum = b[i];
      for (int j=i+1;j<size;j++) 
	sum -= elem[i][j]*b[j]; // i<j (lower left)
      b[i] = sum/elem[i][i]; // i==i (upper right)
    }
    for (int i=0;i<size;i++)
      res[i][j] = b[i];
  }
}

__global__ void ComputeHomographies(float *coord, int *randPts, float *homo, 
  int numPts) 
{
  float a[8][8], ia[8][8];
  float b[8]; 
  const int bx = blockIdx.x;
  const int tx = threadIdx.x;
  const int idx = blockDim.x*bx + tx;
  const int numLoops = blockDim.x*gridDim.x;
  for (int i=0;i<4;i++) {
    int pt = randPts[i*numLoops+idx];
    float x1 = coord[pt+0*numPts];
    float y1 = coord[pt+1*numPts];
    float x2 = coord[pt+2*numPts];
    float y2 = coord[pt+3*numPts];
    float *row1 = a[2*i+0];
    row1[0] = x1;
    row1[1] = y1;
    row1[2] = 1.0;
    row1[3] = row1[4] = row1[5] = 0.0;
    row1[6] = -x2*x1;
    row1[7] = -x2*y1;
    float *row2 = a[2*i+1];
    row2[0] = row2[1] = row2[2] = 0.0;
    row2[3] = x1;
    row2[4] = y1;
    row2[5] = 1.0;
    row2[6] = -y2*x1;
    row2[7] = -y2*y1;
    b[2*i+0] = x2;
    b[2*i+1] = y2;
  }
  InvertMatrix<8>(a, ia);
  __syncthreads();
  for (int j=0;j<8;j++) {
    float sum = 0.0f;
    for (int i=0;i<8;i++) 
      sum += ia[j][i]*b[i];
    homo[j*numLoops+idx] = sum;
  }
  __syncthreads();
}

#define TESTHOMO_TESTS 16 // number of tests per block,  alt. 32, 32
#define TESTHOMO_LOOPS 16 // number of loops per block,  alt.  8, 16 

__global__ void TestHomographies(float *d_coord, float *d_homo, 
  int *d_counts, int numPts, float thresh2)
{
  __shared__ float homo[8*TESTHOMO_LOOPS];
  __shared__ int cnts[TESTHOMO_TESTS*TESTHOMO_LOOPS];
  const int tx = threadIdx.x;
  const int ty = threadIdx.y;
  const int idx = blockIdx.y*blockDim.y + tx;
  const int numLoops = blockDim.y*gridDim.y;
  if (ty<8 && tx<TESTHOMO_LOOPS)
    homo[tx*8+ty] = d_homo[idx+ty*numLoops];
  __syncthreads();
  float a[8];
  for (int i=0;i<8;i++) 
    a[i] = homo[ty*8+i];
  int cnt = 0;
  for (int i=tx;i<numPts;i+=TESTHOMO_TESTS) {
    float x1 = d_coord[i+0*numPts];
    float y1 = d_coord[i+1*numPts];
    float x2 = d_coord[i+2*numPts];
    float y2 = d_coord[i+3*numPts];
    float nomx = __fmul_rz(a[0],x1) + __fmul_rz(a[1],y1) + a[2];
    float nomy = __fmul_rz(a[3],x1) + __fmul_rz(a[4],y1) + a[5];
    float deno = __fmul_rz(a[6],x1) + __fmul_rz(a[7],y1) + 1.0f;
    float errx = __fmul_rz(x2,deno) - nomx;
    float erry = __fmul_rz(y2,deno) - nomy;
    float err2 = __fmul_rz(errx,errx) + __fmul_rz(erry,erry);
    if (err2<__fmul_rz(thresh2,__fmul_rz(deno,deno)))
      cnt ++;
  }
  int kty = TESTHOMO_TESTS*ty;
  cnts[kty + tx] = cnt;
  __syncthreads();
  int len = TESTHOMO_TESTS/2;
  while (len>0) {
    if (tx<len)
      cnts[kty + tx] += cnts[kty + tx + len];
    len /= 2;
    __syncthreads();
  }
  if (tx<TESTHOMO_LOOPS && ty==0)
    d_counts[idx] = cnts[TESTHOMO_TESTS*tx];
  __syncthreads();
}

//================= Host matching functions =====================//

double FindHomography(SiftData &data, float *homography, int *numMatches, int numLoops, float minScore, float maxAmbiguity, float thresh)
{
  *numMatches = 0;
  homography[0] = homography[4] = homography[8] = 1.0f;
  homography[1] = homography[2] = homography[3] = 0.0f;
  homography[5] = homography[6] = homography[7] = 0.0f;
#ifdef MANAGEDMEM
  SiftPoint *d_sift = data.m_data;
#else
  if (data.d_data==NULL)
    return 0.0f;
  SiftPoint *d_sift = data.d_data;
#endif
  TimerGPU timer(0);
  numLoops = iDivUp(numLoops,16)*16;
  int numPts = data.numPts;
  if (numPts<8)
    return 0.0f;
  int numPtsUp = iDivUp(numPts, 16)*16;
  float *d_coord, *d_homo;
  int *d_randPts, *h_randPts;
  int randSize = 4*sizeof(int)*numLoops;
  int szFl = sizeof(float);
  int szPt = sizeof(SiftPoint);
  safeCall(cudaMalloc((void **)&d_coord, 4*sizeof(float)*numPtsUp));
  safeCall(cudaMalloc((void **)&d_randPts, randSize));
  safeCall(cudaMalloc((void **)&d_homo, 8*sizeof(float)*numLoops));
  h_randPts = (int*)malloc(randSize);
  float *h_scores = (float *)malloc(sizeof(float)*numPtsUp);
  float *h_ambiguities = (float *)malloc(sizeof(float)*numPtsUp);
  safeCall(cudaMemcpy2D(h_scores, szFl, &d_sift[0].score, szPt, szFl, numPts, cudaMemcpyDeviceToHost));
  safeCall(cudaMemcpy2D(h_ambiguities, szFl, &d_sift[0].ambiguity, szPt, szFl, numPts, cudaMemcpyDeviceToHost));
  int *validPts = (int *)malloc(sizeof(int)*numPts);
  int numValid = 0;
  for (int i=0;i<numPts;i++) {
    if (h_scores[i]>minScore && h_ambiguities[i]<maxAmbiguity)
      validPts[numValid++] = i;
  }
  free(h_scores);
  free(h_ambiguities);
  if (numValid>=8) {
    for (int i=0;i<numLoops;i++) {
      int p1 = rand() % numValid;
      int p2 = rand() % numValid;
      int p3 = rand() % numValid;
      int p4 = rand() % numValid;
      while (p2==p1) p2 = rand() % numValid;
      while (p3==p1 || p3==p2) p3 = rand() % numValid;
      while (p4==p1 || p4==p2 || p4==p3) p4 = rand() % numValid;
      h_randPts[i+0*numLoops] = validPts[p1];
      h_randPts[i+1*numLoops] = validPts[p2];
      h_randPts[i+2*numLoops] = validPts[p3];
      h_randPts[i+3*numLoops] = validPts[p4];
    }
    safeCall(cudaMemcpy(d_randPts, h_randPts, randSize, cudaMemcpyHostToDevice));
    safeCall(cudaMemcpy2D(&d_coord[0*numPtsUp], szFl, &d_sift[0].xpos, szPt, szFl, numPts, cudaMemcpyDeviceToDevice));
    safeCall(cudaMemcpy2D(&d_coord[1*numPtsUp], szFl, &d_sift[0].ypos, szPt, szFl, numPts, cudaMemcpyDeviceToDevice));
    safeCall(cudaMemcpy2D(&d_coord[2*numPtsUp], szFl, &d_sift[0].match_xpos, szPt, szFl, numPts, cudaMemcpyDeviceToDevice));
    safeCall(cudaMemcpy2D(&d_coord[3*numPtsUp], szFl, &d_sift[0].match_ypos, szPt, szFl, numPts, cudaMemcpyDeviceToDevice));
    ComputeHomographies<<<numLoops/16, 16>>>(d_coord, d_randPts, d_homo, numPtsUp);
    safeCall(cudaDeviceSynchronize());
    checkMsg("ComputeHomographies() execution failed\n");
    dim3 blocks(1, numLoops/TESTHOMO_LOOPS);
    dim3 threads(TESTHOMO_TESTS, TESTHOMO_LOOPS);
    TestHomographies<<<blocks, threads>>>(d_coord, d_homo, d_randPts, numPtsUp, thresh*thresh);
    safeCall(cudaDeviceSynchronize());
    checkMsg("TestHomographies() execution failed\n");
    safeCall(cudaMemcpy(h_randPts, d_randPts, sizeof(int)*numLoops, cudaMemcpyDeviceToHost));
    int maxIndex = -1, maxCount = -1;
    for (int i=0;i<numLoops;i++) 
      if (h_randPts[i]>maxCount) {
	maxCount = h_randPts[i];
	maxIndex = i;
      }
    *numMatches = maxCount;
    safeCall(cudaMemcpy2D(homography, szFl, &d_homo[maxIndex], sizeof(float)*numLoops, szFl, 8, cudaMemcpyDeviceToHost));
  }
  free(validPts);
  free(h_randPts);
  safeCall(cudaFree(d_homo));
  safeCall(cudaFree(d_randPts));
  safeCall(cudaFree(d_coord));
  double gpuTime = timer.read();
#ifdef VERBOSE
  printf("FindHomography time =         %.2f ms\n", gpuTime);
#endif
  return gpuTime;
}


double MatchSiftData(SiftData &data1, SiftData &data2)
{
  TimerGPU timer(0);
  int numPts1 = data1.numPts;
  int numPts2 = data2.numPts;
  if (!numPts1 || !numPts2) 
    return 0.0;
#ifdef MANAGEDMEM
  SiftPoint *sift1 = data1.m_data;
  SiftPoint *sift2 = data2.m_data;
#else
  if (data1.d_data==NULL || data2.d_data==NULL)
    return 0.0f;
  SiftPoint *sift1 = data1.d_data;
  SiftPoint *sift2 = data2.d_data;
#endif
  
// Original version with correlation and maximization in two different kernels
// Global memory reguirement: O(N^2)
#if 0
  float *d_corrData; 
  int corrWidth = iDivUp(numPts2, 16)*16;
  int corrSize = sizeof(float)*numPts1*corrWidth;
  safeCall(cudaMalloc((void **)&d_corrData, corrSize));
#if 0 // K40c 10.9ms, 1080 Ti 3.8ms
  dim3 blocks1(numPts1, iDivUp(numPts2, 16));
  dim3 threads1(16, 16); // each block: 1 points x 16 points
  MatchSiftPoints<<<blocks1, threads1>>>(sift1, sift2, d_corrData, numPts1, numPts2);
#else // K40c 7.6ms, 1080 Ti 1.4ms
  dim3 blocks(iDivUp(numPts1,16), iDivUp(numPts2, 16));
  dim3 threads(16, 16); // each block: 16 points x 16 points
  MatchSiftPoints2<<<blocks, threads>>>(sift1, sift2, d_corrData, numPts1, numPts2);
#endif
  safeCall(cudaDeviceSynchronize());
  dim3 blocksMax(iDivUp(numPts1, 16));
  dim3 threadsMax(16, 16);
  FindMaxCorr<<<blocksMax, threadsMax>>>(d_corrData, sift1, sift2, numPts1, corrWidth, sizeof(SiftPoint));
  safeCall(cudaDeviceSynchronize());
  checkMsg("FindMaxCorr() execution failed\n");
  safeCall(cudaFree(d_corrData));
#endif

// Version suggested by Nicholas Lin with combined correlation and maximization
// Global memory reguirement: O(N)
#if 0 // K40c 51.2ms, 1080 Ti 9.6ms
  int block_dim = 16;
  float *d_corrData;
  int corrSize = numPts1 * block_dim * 2;
  safeCall(cudaMalloc((void **)&d_corrData, sizeof(float) * corrSize));
  dim3 blocks(iDivUp(numPts1, block_dim));
  dim3 threads(block_dim, block_dim); 
  FindMaxCorr3<<<blocks, threads >>>(d_corrData, sift1, sift2, numPts1, numPts2);
  safeCall(cudaDeviceSynchronize());
  checkMsg("FindMaxCorr3() execution failed\n");
  safeCall(cudaFree(d_corrData));
#endif

// Combined version with no global memory requirement using one 1 point per block
#if 0 // K40c 8.9ms, 1080 Ti 2.1ms, 2080 Ti 1.0ms
  dim3 blocksMax(numPts1);
  dim3 threadsMax(FMC2W, FMC2H);
  FindMaxCorr2<<<blocksMax, threadsMax>>>(sift1, sift2, numPts1, numPts2);
  safeCall(cudaDeviceSynchronize());
  checkMsg("FindMaxCorr2() execution failed\n");
#endif
  
// Combined version with no global memory requirement using one FMC2H points per block
#if 0 // K40c 9.2ms, 1080 Ti 1.3ms, 2080 Ti 1.1ms
  dim3 blocksMax2(iDivUp(numPts1, FMC2H));
  dim3 threadsMax2(FMC2W, FMC2H);
  FindMaxCorr4<<<blocksMax2, threadsMax2>>>(sift1, sift2, numPts1, numPts2);
  safeCall(cudaDeviceSynchronize());
  checkMsg("FindMaxCorr4() execution failed\n");
#endif

// Combined version with no global memory requirement using global locks
#if 1
  dim3 blocksMax3(iDivUp(numPts1, 16), iDivUp(numPts2, 512));
  dim3 threadsMax3(16, 16);
  CleanMatches<<<iDivUp(numPts1, 64), 64>>>(sift1, numPts1);
  int mode = 10;
  if (mode==5)// K40c 5.0ms, 1080 Ti 1.2ms, 2080 Ti 0.83ms
    FindMaxCorr5<<<blocksMax3, threadsMax3>>>(sift1, sift2, numPts1, numPts2);
  else if (mode==6) {                    // 2080 Ti 0.89ms
    threadsMax3 = dim3(32, 16);
    FindMaxCorr6<<<blocksMax3, threadsMax3>>>(sift1, sift2, numPts1, numPts2);
  } else if (mode==7)                    // 2080 Ti 0.50ms  
    FindMaxCorr7<<<blocksMax3, threadsMax3>>>(sift1, sift2, numPts1, numPts2);
  else if (mode==8) {                    // 2080 Ti 0.45ms
    blocksMax3 = dim3(iDivUp(numPts1, FMC_BW), iDivUp(numPts2, FMC_GH));
    threadsMax3 = dim3(FMC_NW, FMC_NH);
    FindMaxCorr8<<<blocksMax3, threadsMax3>>>(sift1, sift2, numPts1, numPts2);
  } else if (mode==9) {                  // 2080 Ti 0.46ms
    blocksMax3 = dim3(iDivUp(numPts1, FMC_BW), iDivUp(numPts2, FMC_GH));
    threadsMax3 = dim3(FMC_NW, FMC_NH);
    FindMaxCorr9<<<blocksMax3, threadsMax3>>>(sift1, sift2, numPts1, numPts2);
  } else if (mode==10) {                 // 2080 Ti 0.24ms
    blocksMax3 = dim3(iDivUp(numPts1, M7W));
    threadsMax3 = dim3(M7W, M7H/M7R);
    FindMaxCorr10<<<blocksMax3, threadsMax3>>>(sift1, sift2, numPts1, numPts2);
  }
  safeCall(cudaDeviceSynchronize());
  checkMsg("FindMaxCorr5() execution failed\n");
#endif

  if (data1.h_data!=NULL) {
    float *h_ptr = &data1.h_data[0].score;
    float *d_ptr = &data1.d_data[0].score;
    safeCall(cudaMemcpy2D(h_ptr, sizeof(SiftPoint), d_ptr, sizeof(SiftPoint), 5*sizeof(float), data1.numPts, cudaMemcpyDeviceToHost));
  }

  double gpuTime = timer.read();
#ifndef VERBOSE
  printf("MatchSiftData time =          %.2f ms\n", gpuTime);
#endif
  return gpuTime;
}