@@ -31,18 +31,16 @@ int main(int argc, char *argv[]) else cudaSetDevice( cutGetMaxGflopsDeviceId() ); - int result = MLMain(argc, argv); - cutilExit(argc, argv); - return result; + return MLMain(argc, argv); } void cuFourier1D (double *h_A, long n) { - float norm = 1.0/sqrt((float) n); - int mem_size = sizeof(Complex) * n; + double norm = 1.0/sqrt((double) n); + long mem_size = sizeof(Complex) * n; // Allocate host memory for the signal - Complex* h_signal = (Complex*)malloc(sizeof(Complex) * n); + Complex* h_signal = (Complex*)malloc(mem_size); // Initalize the memory for the signal for (long i = 0; i < n; ++i) { @@ -54,35 +52,41 @@ void cuFourier1D (double *h_A, long n) Complex* d_signal; cutilSafeCall(cudaMalloc((void**)&d_signal, mem_size)); // Copy host memory to device - cutilSafeCall(cudaMemcpy(d_signal, h_signal, mem_size, - cudaMemcpyHostToDevice)); + cutilSafeCall(cudaMemcpy(d_signal, h_signal, mem_size, cudaMemcpyHostToDevice)); // CUFFT plan cufftHandle plan; cufftSafeCall(cufftPlan1d(&plan, n, CUFFT_C2C, 1)); // Transform signal - cufftSafeCall(cufftExecC2C(plan, (cufftComplex *)d_signal, (cufftComplex *)d_signal, CUFFT_FORWARD)); + cufftSafeCall(cufftExecC2C(plan, (cufftComplex *)d_signal, (cufftComplex *)d_signal, CUFFT_INVERSE)); // Copy device memory to host Complex* h_convolved_signal = h_signal; - cutilSafeCall(cudaMemcpy(h_convolved_signal, d_signal, mem_size, - cudaMemcpyDeviceToHost)); + cutilSafeCall(cudaMemcpy(h_convolved_signal, d_signal, mem_size, cudaMemcpyDeviceToHost)); + + // Release d_signal + cutilSafeCall(cudaFree(d_signal)); // Destroy CUFFT context cufftSafeCall(cufftDestroy(plan)); // Return transformed signal to Mathematica as a Complex List - MLPutFunction(stdlink,"List",n); - for (long i = 0; i < n; i++) { - MLPutFunction(stdlink,"Complex",2); - MLPutFloat(stdlink,h_convolved_signal[i].x*norm); - MLPutFloat(stdlink,h_convolved_signal[i].y*norm); - } + MLPutFunction(stdlink, "Map", 2); + MLPutFunction(stdlink, "Function", 2); + MLPutFunction(stdlink, "List", 1); + MLPutSymbol(stdlink, "x"); + MLPutFunction(stdlink, "Apply", 2); + MLPutSymbol(stdlink, "Complex"); + MLPutSymbol(stdlink, "x"); + MLPutFunction(stdlink, "Partition", 2); + MLPutFunction(stdlink, "Times", 2); + MLPutReal(stdlink, norm); + MLPutReal32List(stdlink, (float*)h_convolved_signal, 2*n); + MLPutInteger(stdlink, 2); // Cleanup memory free(h_signal); - cutilSafeCall(cudaFree(d_signal)); cudaThreadExit(); -} \ No newline at end of file +} src/cuFourier/cuFourier.cu cc4833f782001f3277ca7121fce531b9897c328b kashif kashif@nomad-labs.com http://github.com/kashif/mathematica_cuda/commit/6b6bda290f010ba087a900466ab43e41e6984815 6b6bda290f010ba087a900466ab43e41e6984815 2009-09-21T13:46:08-07:00 2009-09-21T13:46:08-07:00 Sending Complex numbers as an array. Suggestion from Patrick Scheibe. d2b29b636309bb6a7382fa4a6f577bb69d40a4a7 kashif kashif@nomad-labs.com