match.cu

//**********************************************************//
//   Matching test code by Marten Bjorkman aka Celebrandil  //
//                                                          //
//   The code includes an example of gradual optimization   //
//   of a kernel for matching two sets of 16K 128D points.  //
//   You are welcome to the code for educational purposes.  //
//                                                          //
//            Fairlight - When Dreams Come True             //
// https://www.youtube.com/channel/UCdHiji77FlppuNK6xemrsVA //
//**********************************************************//

#include <cuda.h>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <iostream>
#include <vector>
#include <memory>
#include <algorithm>
#include <immintrin.h>
#include "cudautils.h"

#define RUNCPU 1
#define CHECK  1
#define NPTS (2048*8)
#define NDIM 128

#define M1W  128
#define M2W   16
#define M2H   16
#define M5W   16
#define M5H   16
#define M5R    4
#define M7W   32
#define M7H   32
#define M7R    4


/*
Data size:   16 MB
Allocate:    1.01194 ms
Upload:      3.69939 ms  4.32503 MB/ms
MatchCPU1:   34649.6 ms  1.89139 Gflops
MatchCPU2:   3064.36 ms  21.3866 Gflops
MatchCPU3:   184.762 ms  354.706 Gflops
MatchGPU1:   641.828 ms  102.108 Gflops
MatchGPU2:   148.020 ms  442.752 Gflops
MatchGPU3:   31.9609 ms  2050.50 Gflops
MatchGPU4:   29.7891 ms  2200.00 Gflops
MatchGPU5:   17.1484 ms  3821.69 Gflops
MatchGPU6:   16.3516 ms  4007.94 Gflops
MatchGPU7:   14.7995 ms  4428.27 Gflops
MatchGPU8:   10.5291 ms  6224.28 Gflops
Download:    0.16016 ms  0.780488 MB/ms
*/

void MatchC1(float *h_pts1, float *h_pts2, float *h_score, int *h_index)
{
  std::memset(h_score, 0, sizeof(float)*NPTS);
  for (int p1=0;p1<NPTS;p1++) {
    for (int p2=0;p2<NPTS;p2++) {
      float score = 0.0f;
      for (int d=0;d<NDIM;d++)
	score += h_pts1[p1*NDIM + d]*h_pts2[p2*NDIM + d];
      if (score>h_score[p1]) {
	h_score[p1] = score;
	h_index[p1] = p2;
      }
    }
  }
}

void MatchC2(float *h_pts1, float *h_pts2, float *h_score, int *h_index)
{
#define BSIZ  256
  std::memset(h_score, 0, sizeof(float)*NPTS);
  for (int b1=0;b1<NPTS;b1+=BSIZ) {
    for (int b2=0;b2<NPTS;b2+=BSIZ) {
      for (int p1=b1;p1<b1+BSIZ;p1++) {
	float *pt1 = &h_pts1[p1*NDIM];
	for (int p2=b2;p2<b2+BSIZ;p2++) {
	  float *pt2 = &h_pts2[p2*NDIM];
	  __m256 score8 = _mm256_setzero_ps();
	  for (int d=0;d<NDIM;d+=8) {
	    __m256 v1 = _mm256_load_ps(pt1 + d);
	    __m256 v2 = _mm256_load_ps(pt2 + d);
	    score8 = _mm256_fmadd_ps(v1, v2, score8);
	  }
	  score8 = _mm256_add_ps(score8, _mm256_permute2f128_ps(score8, score8, 1));
	  score8 = _mm256_hadd_ps(score8, score8);
	  float score = _mm256_cvtss_f32(_mm256_hadd_ps(score8, score8));
	  if (score>h_score[p1]) {
	    h_score[p1] = score;
	    h_index[p1] = p2;
	  }
	}
      }
    }
  }
}

void MatchC3(float *h_pts1, float *h_pts2, float *h_score, int *h_index)
{
#define BSIZ  256
  std::memset(h_score, 0, sizeof(float)*NPTS);
#pragma omp parallel for
  for (int b1=0;b1<NPTS;b1+=BSIZ) {
    for (int b2=0;b2<NPTS;b2+=BSIZ) {
      for (int p1=b1;p1<b1+BSIZ;p1++) {
	float *pt1 = &h_pts1[p1*NDIM];
	for (int p2=b2;p2<b2+BSIZ;p2++) {
	  float *pt2 = &h_pts2[p2*NDIM];
	  __m256 score8 = _mm256_setzero_ps();
	  for (int d=0;d<NDIM;d+=8) {
	    __m256 v1 = _mm256_load_ps(pt1 + d);
	    __m256 v2 = _mm256_load_ps(pt2 + d);
	    score8 = _mm256_fmadd_ps(v1, v2, score8);
	  }
	  score8 = _mm256_add_ps(score8, _mm256_permute2f128_ps(score8, score8, 1));
	  score8 = _mm256_hadd_ps(score8, score8);
	  float score = _mm256_cvtss_f32(_mm256_hadd_ps(score8, score8));
	  if (score>h_score[p1]) {
	    h_score[p1] = score;
	    h_index[p1] = p2;
	  }
	}
      }
    }
  }
}

void CheckMatches(int *h_index, int *h_index2, float *h_score, float *h_score2)
{
  int ndiff = 0;
  for (int i=0;i<NPTS;i++) {
    ndiff += (h_index[i] != h_index2[i]);
    if (h_index[i] != h_index2[i])
      std::cout << "  " << i << " " << h_index[i] << " " << h_index2[i] << " " << h_score[i] << " " << h_score2[i] << std::endl;
  }
  std::cout << "Number of incorrect matches: " << ndiff << std::endl;
}
      

__global__ void Match1(float *d_pts1, float *d_pts2, float *d_score, int *d_index)
{
  int p1 = threadIdx.x + M1W*blockIdx.x;
  float max_score = 0.0f;
  int index = -1;
  
  for (int p2=0;p2<NPTS;p2++) {
    float score = 0.0f;
    for (int d=0;d<NDIM;d++)
      score += d_pts1[p1*NDIM + d]*d_pts2[p2*NDIM + d];
    if (score>max_score) {
      max_score = score;
      index = p2;
    }
  }
  
  d_score[p1] = max_score;
  d_index[p1] = index;
}

__global__ void Match2(float *d_pts1, float *d_pts2, float *d_score, int *d_index)
{
  __shared__ float buffer1[M2W*NDIM];  //%%%%
  __shared__ float buffer2[M2H*NDIM];  //%%%%
  __shared__ float scores[M2W*M2H];    //%%%%
  int tx = threadIdx.x;
  int ty = threadIdx.y;
  int idx = tx + M2W*ty;
  int bp1 = M2W*blockIdx.x;
  if (ty<M2W)
    for (int d=tx;d<NDIM;d+=M2W)
      for (int j=ty;j<M2W;j+=M2H)
	buffer1[j*NDIM + d] = d_pts1[(bp1 + j)*NDIM + d];   //%%%%
  __syncthreads();
  
  float max_score = 0.0f;
  int index = -1;
  for (int bp2=0;bp2<NPTS;bp2+=M2H) {
    for (int d=tx;d<NDIM;d+=M2W)
      buffer2[ty*NDIM + d] = d_pts2[(bp2 + ty)*NDIM + d]; //%%%%
    __syncthreads();

    float score = 0.0f;
    for (int d=0;d<NDIM;d++) 
      score += buffer1[tx*NDIM + d]*buffer2[ty*NDIM + d];   //%%%%
    scores[idx] = score;
    __syncthreads();
    
    if (ty==0) {
      for (int i=0;i<M2H;i++) {
	if (scores[i*M2W + tx]>max_score) {
	  max_score = scores[i*M2W + tx];
	  index = bp2 + i;
	}
      }
    }
    __syncthreads();
  }
  
  if (ty==0) {
    d_score[bp1 + tx] = max_score;
    d_index[bp1 + tx] = index;
  }
}


__global__ void Match3(float *d_pts1, float *d_pts2, float *d_score, int *d_index)
{
  __shared__ float buffer1[M2W*(NDIM + 1)]; //%%%%
  __shared__ float buffer2[M2H*NDIM];
  __shared__ float scores[M2W*M2H];
  int tx = threadIdx.x;
  int ty = threadIdx.y;
  int idx = tx + M2W*ty;
  int bp1 = M2W*blockIdx.x;
  if (ty<M2W)
    for (int d=tx;d<NDIM;d+=M2W)
      for (int j=ty;j<M2W;j+=M2H)
	buffer1[j*(NDIM + 1) + d] = d_pts1[(bp1 + j)*NDIM + d]; //%%%%
  __syncthreads();
  
  float max_score = 0.0f;
  int index = -1;
  for (int bp2=0;bp2<NPTS;bp2+=M2H) {
    for (int d=tx;d<NDIM;d+=M2W)
      buffer2[ty*NDIM + d] = d_pts2[(bp2 + ty)*NDIM + d];
    __syncthreads();

    float score = 0.0f;
    for (int d=0;d<NDIM;d++) 
      score += buffer1[tx*(NDIM + 1) + d]*buffer2[ty*NDIM + d]; //%%%%
    scores[idx] = score;
    __syncthreads();
    
    if (ty==0) {
      for (int i=0;i<M2H;i++) {
	if (scores[i*M2W + tx]>max_score) {
	  max_score = scores[i*M2W + tx];
	  index = bp2 + i;
	}
      }
    }
    __syncthreads();
  }
  
  if (ty==0) {
    d_score[bp1 + tx] = max_score;
    d_index[bp1 + tx] = index;
  }
}


__global__ void Match4(float *d_pts1, float *d_pts2, float *d_score, int *d_index)
{
  __shared__ float4 buffer1[M2W*(NDIM/4 + 1)];  //%%%%
  __shared__ float4 buffer2[M2H*NDIM/4];        //%%%%
  __shared__ float scores[M2W*M2H];
  int tx = threadIdx.x;
  int ty = threadIdx.y;
  int idx = tx + M2W*ty;
  int bp1 = M2W*blockIdx.x;
  if (ty<M2W)
    for (int d=tx;d<NDIM/4;d+=M2W)
      for (int j=ty;j<M2W;j+=M2H)
	buffer1[j*(NDIM/4 + 1) + d] = ((float4*)d_pts1)[(bp1 + j)*(NDIM/4) + d]; //%%%%
  __syncthreads();
  
  float max_score = 0.0f;
  int index = -1;
  for (int bp2=0;bp2<NPTS;bp2+=M2H) {
    for (int d=tx;d<NDIM/4;d+=M2W)
      buffer2[ty*NDIM/4 + d] = ((float4*)d_pts2)[(bp2 + ty)*(NDIM/4) + d]; //%%%%
    __syncthreads();

    float score = 0.0f;
    for (int d=0;d<NDIM/4;d++) {
      float4 v1 = buffer1[tx*(NDIM/4 + 1) + d]; //%%%%
      float4 v2 = buffer2[ty*(NDIM/4) + d];     //%%%%
      score += v1.x*v2.x; score += v1.y*v2.y;
      score += v1.z*v2.z; score += v1.w*v2.w;
    }
    scores[idx] = score;
    __syncthreads();
    
    if (ty==0) {
      for (int i=0;i<M2H;i++) {
	if (scores[i*M2W + tx]>max_score) {
	  max_score = scores[i*M2W + tx];
	  index = bp2 + i;
	}
      }
    }
    __syncthreads();
  }
  
  if (ty==0) {
    d_score[bp1 + tx] = max_score;
    d_index[bp1 + tx] = index;
  }
}

__global__ void Match5(float *d_pts1, float *d_pts2, float *d_score, int *d_index)
{
  __shared__ float4 buffer1[M5W*(NDIM/4 + 1)]; 
  __shared__ float4 buffer2[M5H*NDIM/4];       
  __shared__ float scores[M5W*M5H];
  int tx = threadIdx.x;
  int ty = threadIdx.y;
  int bp1 = M5W*blockIdx.x;
  if (ty<M5W)
    for (int d=tx;d<NDIM/4;d+=M5W)
      for (int j=ty;j<M5W;j+=M5H)
	buffer1[j*(NDIM/4 + 1) + d] = ((float4*)d_pts1)[(bp1 + j)*(NDIM/4) + d];
  __syncthreads();
  
  float max_score = 0.0f;
  int index = -1;
  for (int bp2=0;bp2<NPTS;bp2+=M5H) {
    for (int d=tx;d<NDIM/4;d+=M5W)
      buffer2[ty*NDIM/4 + d] = ((float4*)d_pts2)[(bp2 + ty)*(NDIM/4) + d];
    __syncthreads();

    if (ty<M5H/M5R) {  //%%%%
      float score[M5R];                                    //%%%%
      for (int dy=0;dy<M5R;dy++)
	score[dy] = 0.0f;
      for (int d=0;d<NDIM/4;d++) {
	float4 v1 = buffer1[tx*(NDIM/4 + 1) + d];
	for (int dy=0;dy<M5R;dy++) {
	  float4 v2 = buffer2[(M5R*ty + dy)*(NDIM/4) + d];    //%%%%
	  score[dy] += v1.x*v2.x; score[dy] += v1.y*v2.y;
	  score[dy] += v1.z*v2.z; score[dy] += v1.w*v2.w;
	}
      }
      for (int dy=0;dy<M5R;dy++)
	scores[tx + M5W*(M5R*ty + dy)] = score[dy];
    }
    __syncthreads();
    
    if (ty==0) {
      for (int i=0;i<M5H;i++) {
	if (scores[i*M2W + tx]>max_score) {
	  max_score = scores[i*M5W + tx];
	  index = bp2 + i;
	}
      }
    }
    __syncthreads();
  }

  if (ty==0) {
    d_score[bp1 + tx] = max_score;
    d_index[bp1 + tx] = index;
  }
}


__global__ void Match6(float *d_pts1, float *d_pts2, float *d_score, int *d_index)
{
  __shared__ float4 buffer1[M5W*(NDIM/4 + 1)]; 
  __shared__ float4 buffer2[M5H*NDIM/4];       
  int tx = threadIdx.x;
  int ty = threadIdx.y;
  int bp1 = M5W*blockIdx.x;
  if (ty<M5W)
    for (int d=tx;d<NDIM/4;d+=M5W)
      for (int j=ty;j<M5W;j+=M5H)
	buffer1[j*(NDIM/4 + 1) + d] = ((float4*)d_pts1)[(bp1 + j)*(NDIM/4) + d];
  
  float max_score = 0.0f;
  int index = -1;    
  for (int bp2=0;bp2<NPTS;bp2+=M5H) {
    for (int d=tx;d<NDIM/4;d+=M5W)
      buffer2[ty*NDIM/4 + d] = ((float4*)d_pts2)[(bp2 + ty)*(NDIM/4) + d];
    __syncthreads();

    if (ty<M5H/M5R) {  
      float score[M5R];                                    
      for (int dy=0;dy<M5R;dy++)
	score[dy] = 0.0f;
      for (int d=0;d<NDIM/4;d++) {
	float4 v1 = buffer1[tx*(NDIM/4 + 1) + d];
	for (int dy=0;dy<M5R;dy++) {
	  float4 v2 = buffer2[(M5R*ty + dy)*(NDIM/4) + d];    
	  score[dy] += v1.x*v2.x; score[dy] += v1.y*v2.y;
	  score[dy] += v1.z*v2.z; score[dy] += v1.w*v2.w;
	}
      }
      for (int dy=0;dy<M5R;dy++) {
	if (score[dy]>max_score) {   //%%%%
	  max_score = score[dy];     
	  index = bp2 + M5R*ty + dy;               
	}
      }
    }
    __syncthreads();
  }

  float *scores = (float*)buffer1;
  int *indices = (int*)&scores[M5W*M5H/M5R];
  if (ty<M5H/M5R) {
    scores[ty*M5W + tx] = max_score;  //%%%%
    indices[ty*M5W + tx] = index;     //%%%%
  }
  __syncthreads();
  
  if (ty==0) {
    max_score = scores[tx];
    index = indices[tx];
    for (int y=0;y<M5H/M5R;y++)
      if (scores[y*M5W + tx]>max_score) {
	max_score = scores[y*M5W + tx]; //%%%%
	index = indices[y*M5W + tx];    //%%%%
      }
    d_score[bp1 + tx] = max_score;
    d_index[bp1 + tx] = index;
  }
}

__global__ void Match7(float *d_pts1, float *d_pts2, float *d_score, int *d_index)
{
  __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 d=tx;d<NDIM/4;d+=M7W)
    for (int j=ty;j<M7W;j+=M7H/M7R)      //%%%%
      buffer1[j*NDIM/4 + (d + j)%(NDIM/4)] = ((float4*)d_pts1)[(bp1 + j)*(NDIM/4) + d];
  
  float max_score = 0.0f;
  int index = -1;    
  for (int bp2=0;bp2<NPTS;bp2+=M7H) {
    for (int d=tx;d<NDIM/4;d+=M7W)
      for (int j=ty;j<M7H;j+=M7H/M7R)       //%%%%
	buffer2[j*NDIM/4 + d] = ((float4*)d_pts2)[(bp2 + j)*(NDIM/4) + d];
    __syncthreads();

    float score[M7R];                                    
    for (int dy=0;dy<M7R;dy++)
      score[dy] = 0.0f;
    for (int d=0;d<NDIM/4;d++) {
      float4 v1 = buffer1[tx*NDIM/4 + (d + tx)%(NDIM/4)];
      for (int dy=0;dy<M7R;dy++) {
	float4 v2 = buffer2[(M7R*ty + dy)*(NDIM/4) + d];    
	score[dy] += v1.x*v2.x; score[dy] += v1.y*v2.y;
	score[dy] += v1.z*v2.z; score[dy] += v1.w*v2.w;
      }
    }
    for (int dy=0;dy<M7R;dy++) {
      if (score[dy]>max_score) {   
	max_score = score[dy];     
	index = bp2 + M7R*ty + dy;               
      }
    }
    __syncthreads();
  }

  float *scores = (float*)buffer1;
  int *indices = (int*)&scores[M7W*M7H/M7R];
  scores[ty*M7W + tx] = max_score;  
  indices[ty*M7W + tx] = index;     
  __syncthreads();
  
  if (ty==0) {
    max_score = scores[tx];
    index = indices[tx];
    for (int y=0;y<M7H/M7R;y++)
      if (scores[y*M7W + tx]>max_score) {
	max_score = scores[y*M7W + tx]; 
	index = indices[y*M7W + tx];    
      }
    d_score[bp1 + tx] = max_score;
    d_index[bp1 + tx] = index;
  }
}

__global__ void Match8(float *d_pts1, float *d_pts2, float *d_score, int *d_index)
{
  __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 d=tx;d<NDIM/4;d+=M7W)
    for (int j=ty;j<M7W;j+=M7H/M7R)     
      buffer1[j*NDIM/4 + (d + j)%(NDIM/4)] = ((float4*)d_pts1)[(bp1 + j)*(NDIM/4) + d];

#define NRX 2
  float max_score[NRX];
  int index[NRX];
  for (int i=0;i<NRX;i++) {
    max_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<NPTS;bp2+=M7H) {
    for (int d=tx;d<NDIM/4;d+=M7W)
      for (int j=ty;j<M7H;j+=M7H/M7R)       
	buffer2[j*NDIM/4 + d] = ((float4*)d_pts2)[(bp2 + j)*(NDIM/4) + 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]) {
	    max_score[i] = score[dy][i];     
	    index[i] = bp2 + M7R*iy + dy;
	  }
	}
      }
    }
    __syncthreads();
  }

  float *scores = (float*)buffer1;
  int *indices = (int*)&scores[M7W*M7H/M7R];
  if (idx<M7W*M7H/M7R/NRX) {
    for (int i=0;i<NRX;i++) {
      scores[iy*M7W + (M7W/NRX)*i + ix] = max_score[i];  
      indices[iy*M7W + (M7W/NRX)*i + ix] = index[i];
    }
  }
  __syncthreads();
  
  if (ty==0) {
    float max_score = scores[tx];
    int index = indices[tx];
    for (int y=0;y<M7H/M7R;y++)
      if (scores[y*M7W + tx]>max_score) {
	max_score = scores[y*M7W + tx]; 
	index = indices[y*M7W + tx];    
      }
    d_score[bp1 + tx] = max_score;
    d_index[bp1 + tx] = index;
  }
}

__global__ void Match8small(float *d_pts1, float *d_pts2, float *d_score, int *d_index)
{
#define NRX 2
  __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 d=tx;d<NDIM/4;d+=M7W)
    for (int j=ty;j<M7W;j+=M7H/M7R/NRX)     
      buffer1[j*NDIM/4 + (d + j)%(NDIM/4)] = ((float4*)d_pts1)[(bp1 + j)*(NDIM/4) + d];

  float max_score[NRX];
  int index[NRX];
  for (int i=0;i<NRX;i++) {
    max_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<NPTS;bp2+=M7H) {
    for (int d=tx;d<NDIM/4;d+=M7W)
      for (int j=ty;j<M7H;j+=M7H/M7R/NRX)       
	buffer2[j*NDIM/4 + d] = ((float4*)d_pts2)[(bp2 + j)*(NDIM/4) + d];
    __syncthreads();

    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]) {
	  max_score[i] = score[dy][i];     
	  index[i] = bp2 + M7R*iy + dy;
	}
      }
    }
    __syncthreads();
  }

  float *scores = (float*)buffer1;
  int *indices = (int*)&scores[M7W*M7H/M7R];
  if (idx<M7W*M7H/M7R/NRX) {
    for (int i=0;i<NRX;i++) {
      scores[iy*M7W + (M7W/NRX)*i + ix] = max_score[i];  
      indices[iy*M7W + (M7W/NRX)*i + ix] = index[i];
    }
  }
  __syncthreads();
  
  if (ty==0) {
    float max_score = scores[tx];
    int index = indices[tx];
    for (int y=0;y<M7H/M7R;y++)
      if (scores[y*M7W + tx]>max_score) {
	max_score = scores[y*M7W + tx]; 
	index = indices[y*M7W + tx];    
      }
    d_score[bp1 + tx] = max_score;
    d_index[bp1 + tx] = index;
  }
}

__global__ void Match8blocked(float *d_pts1, float *d_pts2, float *d_score, int *d_index)
{
#define NRX 2
#define NUM (NRX*M7R)                       // 32*8 threads
  __shared__ float4 buffer1[M7W*NDIM/4];    // 32*32
  __shared__ float4 buffer2[M7H*NUM];       // 32*8
  int tx = threadIdx.x;
  int ty = threadIdx.y;
  int bp1 = M7W*blockIdx.x;
  for (int d=tx;d<NDIM/4;d+=M7W)
    for (int j=ty;j<M7W;j+=M7H/M7R)     
      buffer1[j*NDIM/4 + (d + j)%(NDIM/4)] = ((float4*)d_pts1)[(bp1 + j)*(NDIM/4) + d];

  float max_score[NRX];
  int index[NRX];
  for (int i=0;i<NRX;i++) {
    max_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<NPTS;bp2+=M7H) {
    float score[M7R][NRX];                                    
    for (int dy=0;dy<M7R;dy++)
      for (int i=0;i<NRX;i++)
	score[dy][i] = 0.0f;

    int d = (idx%NUM);
    int j = (idx/NUM);
    buffer2[j*NUM + d] = ((float4*)d_pts2)[(bp2 + j)*(NDIM/4) + d];
    __syncthreads();
    for (int dp=0;dp<NDIM/4;dp+=NUM) {
      float4 temp;
      if (dp<(NDIM/4-NUM))
	temp = ((float4*)d_pts2)[(bp2 + j)*(NDIM/4) + dp + d + NUM];

      if (idx<M7W*M7H/M7R/NRX) {
	for (int d=0;d<NUM;d++) {
	  float4 v1[NRX];
#pragma unroll
	  for (int i=0;i<NRX;i++) 
	    v1[i] = buffer1[(((M7W/NRX)*i + ix)<<5) + ((dp + d + (M7W/NRX)*i + ix)&31)];
	  //v1[i] = buffer1[((M7W/NRX)*i + ix)*NDIM/4 + (dp + d + (M7W/NRX)*i + ix)%(NDIM/4)];
#pragma unroll
	  for (int dy=0;dy<M7R;dy++) {
	    float4 v2 = buffer2[(M7R*iy + dy)*NUM + d];    
#pragma unroll
	    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;
	    }
	  }
	}
      }
      __syncthreads();

      if (dp<(NDIM/4-NUM)) {
	buffer2[j*NUM + d] = temp;
	__syncthreads();
      }
    }
    for (int dy=0;dy<M7R;dy++) {
      for (int i=0;i<NRX;i++) {
	if (score[dy][i]>max_score[i]) {
	  max_score[i] = score[dy][i];     
	  index[i] = bp2 + M7R*iy + dy;
	}
      }
    }
    __syncthreads();
  }

  float *scores = (float*)buffer1;
  int *indices = (int*)&scores[M7W*M7H/M7R];
  if (idx<M7W*M7H/M7R/NRX) {
    for (int i=0;i<NRX;i++) {
      scores[iy*M7W + (M7W/NRX)*i + ix] = max_score[i];  
      indices[iy*M7W + (M7W/NRX)*i + ix] = index[i];
    }
  }
  __syncthreads();
  
  if (ty==0) {
    float max_score = scores[tx];
    int index = indices[tx];
    for (int y=0;y<M7H/M7R;y++)
      if (scores[y*M7W + tx]>max_score) {
	max_score = scores[y*M7W + tx]; 
	index = indices[y*M7W + tx];    
      }
    d_score[bp1 + tx] = max_score;
    d_index[bp1 + tx] = index;
  }
}

__global__ void Match8blocked2(float *d_pts1, float *d_pts2, float *d_score, int *d_index)
{
#define NRX 2
#define NUM (NRX*M7R)                       // 32*8 threads
  __shared__ float4 buffer1[M7W*NDIM/4];    // 32*32
  __shared__ float4 buffer2[M7H*NUM];       // 32*8
  int tx = threadIdx.x;
  int ty = threadIdx.y;
  int bp1 = M7W*blockIdx.x;
  for (int d=tx;d<NDIM/4;d+=M7W)
    for (int j=ty;j<M7W;j+=M7H/M7R)     
      buffer1[j*NDIM/4 + (d + j)%(NDIM/4)] = ((float4*)d_pts1)[(bp1 + j)*(NDIM/4) + d];

  float max_score[NRX];
  int index[NRX];
  for (int i=0;i<NRX;i++) {
    max_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<NPTS;bp2+=M7H) {
    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 dp=0;dp<NDIM/4;dp+=NUM) {
      int d = (idx%NUM);
      int j = (idx/NUM);
      buffer2[j*NUM + d] = ((float4*)d_pts2)[(bp2 + j)*(NDIM/4) + dp + d];
      __syncthreads();

      if (idx<M7W*M7H/M7R/NRX) {
	for (int d=0;d<NUM;d++) {
	  float4 v1[NRX];
	  for (int i=0;i<NRX;i++) 
	    v1[i] = buffer1[((M7W/NRX)*i + ix)*NDIM/4 + (dp + d + (M7W/NRX)*i + ix)%(NDIM/4)];
	  for (int dy=0;dy<M7R;dy++) {
	    float4 v2 = buffer2[(M7R*iy + dy)*NUM + 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;
	    }
	  }
	}
      }
      __syncthreads();
    }
    for (int dy=0;dy<M7R;dy++) {
      for (int i=0;i<NRX;i++) {
	if (score[dy][i]>max_score[i]) {
	  max_score[i] = score[dy][i];     
	  index[i] = bp2 + M7R*iy + dy;
	}
      }
    }
    __syncthreads();
  }

  float *scores = (float*)buffer1;
  int *indices = (int*)&scores[M7W*M7H/M7R];
  if (idx<M7W*M7H/M7R/NRX) {
    for (int i=0;i<NRX;i++) {
      scores[iy*M7W + (M7W/NRX)*i + ix] = max_score[i];  
      indices[iy*M7W + (M7W/NRX)*i + ix] = index[i];
    }
  }
  __syncthreads();
  
  if (ty==0) {
    float max_score = scores[tx];
    int index = indices[tx];
    for (int y=0;y<M7H/M7R;y++)
      if (scores[y*M7W + tx]>max_score) {
	max_score = scores[y*M7W + tx]; 
	index = indices[y*M7W + tx];    
      }
    d_score[bp1 + tx] = max_score;
    d_index[bp1 + tx] = index;
  }
}

__global__ void Match9(float *d_pts1, float *d_pts2, float *d_score, int *d_index)
{
#define NRX 2
#define NUM 8
  __shared__ float4 buffer1[M7W*NDIM/4]; 
  __shared__ float4 buffer2[M7H*NUM];       
  int tx = threadIdx.x;
  int ty = threadIdx.y;
  int bp1 = M7W*blockIdx.x;
  for (int d=tx;d<NDIM/4;d+=M7W)
    for (int j=ty;j<M7W;j+=M7H/M7R)     
      buffer1[j*NDIM/4 + (d + j)%(NDIM/4)] = ((float4*)d_pts1)[(bp1 + j)*(NDIM/4) + d];

  float max_score[NRX];
  int index[NRX];
  for (int i=0;i<NRX;i++) {
    max_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<NPTS;bp2+=M7H) {
    
      float score[M7R][NRX];                                    
      if (idx<M7W*M7H/M7R/NRX) {    // 128
	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+=NUM) {
	if (idx<M7H*NUM)            // 256
	  buffer2[idx] = ((float4*)d_pts2)[(bp2 + (idx/NUM))*(NDIM/4) + d + (idx%NUM)];
	__syncthreads();

	if (idx<M7W*M7H/M7R/NRX) {  // 128
	  for (int j=0;j<NUM;j++) {
	    float4 v1[NRX];
	    for (int i=0;i<NRX;i++) 
	      v1[i] = buffer1[((M7W/NRX)*i + ix)*NDIM/4 + (d + j + (M7W/NRX)*i + ix)%(NDIM/4)];
	    for (int dy=0;dy<M7R;dy++) {
	      float4 v2 = buffer2[(M7R*ty + dy)*NUM + j];
	      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;
	      }	      
	    }
	  }
	}
	__syncthreads();
      }
      
      if (idx<M7W*M7H/M7R/NRX) {  // 128
	for (int dy=0;dy<M7R;dy++) {
	  for (int i=0;i<NRX;i++) {
	    if (score[dy][i]>max_score[i]) {
	      max_score[i] = score[dy][i];     
	      index[i] = bp2 + M7R*iy + dy;
	    }
	  }
	}
      }
      __syncthreads();
  }

  float *scores = (float*)buffer1;
  int *indices = (int*)&scores[M7W*M7H/M7R];
  if (idx<M7W*M7H/M7R/NRX) {
    for (int i=0;i<NRX;i++) {
      scores[iy*M7W + (M7W/NRX)*i + ix] = max_score[i];  
      indices[iy*M7W + (M7W/NRX)*i + ix] = index[i];
    }
  }
  __syncthreads();
  
  if (ty==0) {
    float max_score = scores[tx];
    int index = indices[tx];
    for (int y=0;y<M7H/M7R;y++)
      if (scores[y*M7W + tx]>max_score) {
	max_score = scores[y*M7W + tx]; 
	index = indices[y*M7W + tx];    
      }
    d_score[bp1 + tx] = max_score;
    d_index[bp1 + tx] = index;
  }
}


int main(int argc, char *argv[])
{
  safeCall(cudaSetDevice(0));

  size_t space = sizeof(float)*NPTS*NDIM*2 + 8;
  std::vector<float> data(NPTS*NDIM*2 + 8);
  void *ptr = (void*)&data[0];
  float *h_pts1 = (float*)std::align(32, sizeof(float)*NPTS*NDIM, ptr, space);
  ptr = (void*)&data[NPTS*NDIM];
  float *h_pts2 = (float*)std::align(32, sizeof(float)*NPTS*NDIM, ptr, space);
  std::vector<int> h_index(NPTS);
  std::vector<float> h_score(NPTS);
  std::vector<int> h_index2(NPTS);
  std::vector<float> h_score2(NPTS);
  
  float *d_pts1, *d_pts2, *d_score;
  int *d_index;
  std::cout << std::endl;
  int psize = sizeof(float)*NPTS;
  std::cout << "Data size:   " << 2.0*psize*NDIM/1024/1024 << " MB" << std::endl;
  TimerGPU time;
  float ltime = time.read();

  safeCall(cudaMalloc((void **)&d_pts1, psize*NDIM));
  safeCall(cudaMalloc((void **)&d_pts2, psize*NDIM));
  safeCall(cudaMalloc((void **)&d_index, psize));
  safeCall(cudaMalloc((void **)&d_score, psize));
  std::cout << "Allocate:    " << time.read() - ltime << " ms" << std::endl;

  for (int i=0;i<NPTS;i++) {
    float sum1 = 0.0f, sum2 = 0.0f;
    for (int d=0;d<NDIM;d++) {
      sum1 += h_pts1[i*NDIM + d] = (float)rand()/RAND_MAX;
      sum2 += h_pts2[i*NDIM + d] = (float)rand()/RAND_MAX;
    }
    sum1 = sqrt(NDIM)/sum1;
    sum2 = sqrt(NDIM)/sum2;
    for (int d=0;d<NDIM;d++) {
      h_pts1[i*NDIM + d] *= sum1;
      h_pts2[i*NDIM + d] *= sum2;
    }
  }
  ltime = time.read();
  safeCall(cudaMemcpy(d_pts1, h_pts1, psize*NDIM, cudaMemcpyHostToDevice));
  safeCall(cudaMemcpy(d_pts2, h_pts2, psize*NDIM, cudaMemcpyHostToDevice));
  float delay = time.read() - ltime;
  std::cout << "Upload:      " << delay << " ms  " << 2*psize*NDIM/delay/1024/1024 << " MB/ms" << std::endl;

  if (RUNCPU) {
#if 0
    ltime = time.read();
    MatchC1(h_pts1, h_pts2, h_score.data(), h_index.data());
    delay = time.read() - ltime;
    std::cout << "MatchCPU1:   " << delay << " ms  " << 2.0*NPTS*NPTS*NDIM/delay/1024/1024 << " Gflops" << std::endl;

    ltime = time.read();
    MatchC2(h_pts1, h_pts2, h_score.data(), h_index.data());
    delay = time.read() - ltime;
    std::cout << "MatchCPU2:   " << delay << " ms  " << 2.0*NPTS*NPTS*NDIM/delay/1024/1024 << " Gflops" << std::endl;
#endif

    ltime = time.read();
    MatchC3(h_pts1, h_pts2, h_score.data(), h_index.data());
    delay = time.read() - ltime;
    std::cout << "MatchCPU3:   " << delay << " ms  " << 2.0*NPTS*NPTS*NDIM/delay/1024/1024 << " Gflops" << std::endl;
  }
  dim3 blocks, threads;
#if 0
  blocks = dim3(NPTS/M1W);
  threads = dim3(M1W);
  ltime = time.read();
  Match1<<<blocks,threads>>>(d_pts1, d_pts2, d_score, d_index);
  delay = time.read() - ltime;
  checkMsg("Match1 error");
  std::cout << "MatchGPU1:   " << delay << " ms  " << 2.0*NPTS*NPTS*NDIM/delay/1024/1024 << " Gflops" << std::endl;

  blocks = dim3(NPTS/M2W);
  threads = dim3(M2W, M2H);
  ltime = time.read();
  Match2<<<blocks,threads>>>(d_pts1, d_pts2, d_score, d_index);
  delay = time.read() - ltime;
  checkMsg("Match2 error");
  std::cout << "MatchGPU2:   " << delay << " ms  " << 2.0*NPTS*NPTS*NDIM/delay/1024/1024 << " Gflops" << std::endl;
#endif  

  blocks = dim3(NPTS/M2W);
  threads = dim3(M2W, M2H);
  ltime = time.read();
  Match3<<<blocks,threads>>>(d_pts1, d_pts2, d_score, d_index);
  delay = time.read() - ltime;
  checkMsg("Match3 error");
  std::cout << "MatchGPU3:   " << delay << " ms  " << 2.0*NPTS*NPTS*NDIM/delay/1024/1024 << " Gflops" << std::endl;
  
  blocks = dim3(NPTS/M2W);
  threads = dim3(M2W, M2H);
  ltime = time.read();
  Match4<<<blocks,threads>>>(d_pts1, d_pts2, d_score, d_index);
  delay = time.read() - ltime;
  checkMsg("Match4 error");
  std::cout << "MatchGPU4:   " << delay << " ms  " << 2.0*NPTS*NPTS*NDIM/delay/1024/1024 << " Gflops" << std::endl;
  
  blocks = dim3(NPTS/M5W);
  threads = dim3(M5W, M5H);
  ltime = time.read();
  Match5<<<blocks,threads>>>(d_pts1, d_pts2, d_score, d_index);
  delay = time.read() - ltime;
  checkMsg("Match5 error");
  std::cout << "MatchGPU5:   " << delay << " ms  " << 2.0*NPTS*NPTS*NDIM/delay/1024/1024 << " Gflops" << std::endl;
  
  blocks = dim3(NPTS/M5W);
  threads = dim3(M5W, M5H);
  ltime = time.read();
  Match6<<<blocks,threads>>>(d_pts1, d_pts2, d_score, d_index);
  delay = time.read() - ltime;
  checkMsg("Match6 error");
  std::cout << "MatchGPU6:   " << delay << " ms  " << 2.0*NPTS*NPTS*NDIM/delay/1024/1024 << " Gflops" << std::endl;

  blocks = dim3(NPTS/M7W);
  threads = dim3(M7W, M7H/M7R);
  ltime = time.read();
  Match7<<<blocks,threads>>>(d_pts1, d_pts2, d_score, d_index);
  delay = time.read() - ltime;
  checkMsg("Match7 error");
  std::cout << "MatchGPU7:   " << delay << " ms  " << 2.0*NPTS*NPTS*NDIM/delay/1024/1024 << " Gflops" << std::endl;

  blocks = dim3(NPTS/M7W);
  threads = dim3(M7W, M7H/M7R);
  ltime = time.read();
  Match8<<<blocks,threads>>>(d_pts1, d_pts2, d_score, d_index);
  delay = time.read() - ltime;
  checkMsg("Match8 error");
  std::cout << "MatchGPU8:   " << delay << " ms  " << 2.0*NPTS*NPTS*NDIM/delay/1024/1024 << " Gflops" << std::endl;
  #if 1
  blocks = dim3(NPTS/M7W);
  threads = dim3(M7W, M7H/M7R/2);
  ltime = time.read();
  Match8small<<<blocks,threads>>>(d_pts1, d_pts2, d_score, d_index);
  delay = time.read() - ltime;
  checkMsg("Match8small error");
  std::cout << "Match8small:   " << delay << " ms  " << 2.0*NPTS*NPTS*NDIM/delay/1024/1024 << " Gflops" << std::endl;
  #endif
  #if 1
  blocks = dim3(NPTS/M7W);
  threads = dim3(M7W, M7H/M7R);
  ltime = time.read();
  Match8blocked<<<blocks,threads>>>(d_pts1, d_pts2, d_score, d_index);
  delay = time.read() - ltime;
  checkMsg("Match8blocked error");
  std::cout << "MatchGPU8blocked:   " << delay << " ms  " << 2.0*NPTS*NPTS*NDIM/delay/1024/1024 << " Gflops" << std::endl;
  #endif
  ltime = time.read();
  safeCall(cudaMemcpy(h_index2.data(), d_index, psize, cudaMemcpyDeviceToHost));
  safeCall(cudaMemcpy(h_score2.data(), d_score, psize, cudaMemcpyDeviceToHost));
  delay = time.read() - ltime;
  std::cout << "Download:    " << delay << " ms  " << 2*psize/delay/1024/1024 << " MB/ms" << std::endl;
  ltime = time.read();

  if (CHECK)
    CheckMatches(h_index.data(), h_index2.data(), h_score.data(), h_score2.data());

  std::cout << std::endl;
  safeCall(cudaFree(d_pts1));
  safeCall(cudaFree(d_pts2));
  safeCall(cudaFree(d_index));
  safeCall(cudaFree(d_score));
}