Add cuda intersect and unittest.

src/Function1D.cu

@@ -1253,3 +1253,60 @@ CudaMatrix Aurora::auroraUnion(const CudaMatrix& aMatrix1, const CudaMatrix& aMa
 
     return CudaMatrix::fromRawData(data, endPointer - data);    
 }
+
+CudaMatrix Aurora::intersect(const CudaMatrix& aMatrix1, const CudaMatrix& aMatrix2)
+{
+    if(aMatrix1.isNull() || aMatrix2.isNull() || aMatrix1.isComplex() || aMatrix2.isComplex())
+    {
+        std::cerr<<"intersect not support complex cudamatrix"<<std::endl;
+        return CudaMatrix();
+    }
+
+    size_t size1= aMatrix1.getDataSize();
+    size_t size2= aMatrix2.getDataSize();
+    float* data = nullptr;
+    cudaMalloc((void**)&data, sizeof(float) * (size1 + size2));
+    cudaMemcpy(data, aMatrix1.getData(), sizeof(float) * size1, cudaMemcpyDeviceToDevice);
+    cudaMemcpy(data + size1, aMatrix2.getData(), sizeof(float) * size2, cudaMemcpyDeviceToDevice);
+    thrust::sort(thrust::device, data, data+size1);
+    thrust::sort(thrust::device, data+size1, data+size1+size2);
+    float* end = thrust::set_intersection(thrust::device, data, data+size1,data+size1, data+size1+size2,data);
+
+    return CudaMatrix::fromRawData(data, end - data);
+}
+
+__global__ void intersectKernel(float* aMatrixData, float* aIntersectData, unsigned int aMatrixDataSize, float* aOutputData, unsigned int aOutputDataSize)
+{
+    unsigned int idx = blockIdx.x * blockDim.x + threadIdx.x;
+    if (idx < aOutputDataSize)
+    {
+        for(unsigned int i=0; i<aMatrixDataSize; ++i)
+        {
+            if(aMatrixData[i] == aIntersectData[idx])
+            {
+                aOutputData[idx] = i+1;
+                return;
+            }
+        }
+    }
+}
+
+CudaMatrix Aurora::intersect(const CudaMatrix& aMatrix1, const CudaMatrix& aMatrix2, CudaMatrix& aIa)
+{
+    if(aMatrix1.isNull() || aMatrix2.isNull() || aMatrix1.isComplex() || aMatrix2.isComplex())
+    {
+        std::cerr<<"intersect not support complex cudamatrix"<<std::endl;
+        return CudaMatrix();
+    }
+    CudaMatrix result = intersect(aMatrix1,aMatrix2);
+
+    size_t size = result.getDataSize();
+    float* iaResult = nullptr;
+    cudaMalloc((void**)&iaResult, sizeof(float) * size);
+    int blocksPerGrid = (size + THREADS_PER_BLOCK - 1) / THREADS_PER_BLOCK;
+    intersectKernel<<<blocksPerGrid, THREADS_PER_BLOCK>>>(aMatrix1.getData(), result.getData(), aMatrix1.getDataSize(), iaResult, size);
+    cudaDeviceSynchronize();
+
+    aIa = CudaMatrix::fromRawData(iaResult,size);
+    return result;
+}

src/Function1D.cuh

@@ -75,6 +75,10 @@ namespace Aurora
 
     CudaMatrix auroraUnion(const CudaMatrix& aMatrix1, const CudaMatrix& aMatrix2);
 
+    CudaMatrix intersect(const CudaMatrix& aMatrix1, const CudaMatrix& aMatrix2);
+
+    CudaMatrix intersect(const CudaMatrix& aMatrix1, const CudaMatrix& aMatrix2, CudaMatrix& aIa);
+
     // ------compareSet----------------------------------------------------