Homec4science
Diffusion iCAPs public version 2 (master) Browse

iCAPs_v2/MyTemporal_MEX.cpp
e5f5005bb2efmaster

Last Change
Raw File

MyTemporal_MEX.cpp
View Options

/*
Author : E. Orliac, SCITAS, EPFL
Date : 08.11.2017
Purpose: CUDA C MEX implementation of MyTemporal.m loop over voxels.
Remarks:
*/
#include "mex.h"
#include "cuda_runtime.h"
#include "kernel.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <device_launch_parameters.h>
#include <device_functions.h>
#include <math.h>
#define CHECK(call) { \
cudaError_t err; \
if ( (err = (call)) != cudaSuccess) { \
fprintf(stderr, "Got error: %s at %s:%d\n", cudaGetErrorString(err), __FILE__, __LINE__); \
exit(0); \
} \
}
static const double h_6[6] = {
0.33267055295008261599851158914,
0.80689150931109257649449360409,
0.45987750211849157009515194215,
-0.13501102001025458869638990670,
-0.08544127388202666169281916918,
0.03522629188570953660274066472
};
static const double g_6[6] = {
0.03522629188570953660274066472,
0.08544127388202666169281916918,
-0.13501102001025458869638990670,
-0.45987750211849157009515194215,
0.80689150931109257649449360409,
-0.33267055295008261599851158914
};
static const double unity[1] = {1.0};
extern void my_kernel_wrapper(dim3 dimGrid,
dim3 dimBlock,
double *tcIn, const int tcLength,
const int voxelNb,
double *num, const int lnum,
const int lden,
double *den1, const int lden1,
double *den2, const int lden2,
const double lambdaTemp,
double maxeig,
const int cost_save,
const int nit,
const double *noiseFinIn,
double *tcOut,
double *noiseFinOut);
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
double *tcin, *num, *den1, *den2, *tcout, *noiseFinIn, *noiseFinOut;
double *d_tcin, *d_num, *d_den1, *d_den2, *d_tcout, *d_noiseFinIn, *d_noiseFinOut;
double lambdaTemp, /*lambda,*/ maxeig, noise_est;
int lden, cost_save, nit;
const double *g6, *h6, *p_unity;
g6 = g_6;
h6 = h_6;
p_unity = unity;
tcin = mxGetPr(prhs[0]); /* Input data */
num = mxGetPr(prhs[1]); /* Numerator of filter function */
lden = (int) mxGetScalar(prhs[2]);
den1 = mxGetPr(prhs[3]); /* Causal part of denominator */
den2 = mxGetPr(prhs[4]); /* Non-causal part of denominator */
lambdaTemp = mxGetScalar(prhs[5]);
maxeig = mxGetScalar(prhs[6]);
cost_save = (int) mxGetScalar(prhs[7]);
nit = (int) mxGetScalar(prhs[8]);
noiseFinIn = mxGetPr(prhs[9]);
const mwSize *dimsTcin = mxGetDimensions(prhs[0]); /* tc length x number of voxels */
const mwSize *dimsNum = mxGetDimensions(prhs[1]); /* Numerator dimensions*/
const mwSize *dimsDen1 = mxGetDimensions(prhs[3]);
const mwSize *dimsDen2 = mxGetDimensions(prhs[4]);
printf("tcin dimensions = %d x %d\n", dimsTcin[0], dimsTcin[1]);
printf("num dimensions = %d x %d\n", dimsNum[0], dimsNum[1]);
printf("lden = %d\n", lden);
printf("den1 dimensions = %d x %d\n", dimsDen1[0], dimsDen1[1]);
printf("den2 dimensions = %d x %d\n", dimsDen2[0], dimsDen2[1]);
printf("lambdaTemp = %15.10f\n", lambdaTemp);
printf("maxeig = %15.10f\n", maxeig);
printf("cost_save = %d\n", cost_save);
printf("numb. it. nit = %d\n", nit);
// See if that should be flexible or not
if (dimsNum[0] != 1 || dimsNum[1] != 6)
mexErrMsgIdAndTxt("MyTemporal_MEX:num:dims", "Expected dims 1x6");
uint tcinBytes = dimsTcin[0] * dimsTcin[1] * sizeof(double);
uint numBytes = dimsNum[0] * dimsNum[1] * sizeof(double);
uint den1Bytes = dimsDen1[0] * dimsDen1[1] * sizeof(double);
uint den2Bytes = dimsDen2[0] * dimsDen2[1] * sizeof(double);
uint noiseFinBytes = dimsTcin[1] * sizeof(double);
// Set up output for Matlab parallel process
plhs[0] = mxCreateDoubleMatrix(dimsTcin[0], dimsTcin[1], mxREAL);
tcout = mxGetPr(plhs[0]);
plhs[1] = mxCreateDoubleMatrix(dimsTcin[1], 1, mxREAL);
noiseFinOut = mxGetPr(plhs[1]);
// Set device
int dev = 0;
int deviceCount = 0;
CHECK(cudaSetDevice(dev));
CHECK(cudaGetDeviceCount(&deviceCount));
printf("There are %d GPUs available.\n", deviceCount);
size_t size;
CHECK(cudaDeviceGetLimit(&size, cudaLimitMallocHeapSize));
printf("GPU cudaLimitMallocHeapSize = %d\n", size);
CHECK(cudaDeviceSetLimit(cudaLimitMallocHeapSize, size*50));
CHECK(cudaDeviceGetLimit(&size, cudaLimitMallocHeapSize));
printf("GPU cudaLimitMallocHeapSize = %d\n", size);
// Force the creation of the CUDA context so context
// creation overhead is readily distinguishable when profiling
//CHECK(cudaFree(0));
// Allocate device global memory
CHECK(cudaMalloc((double**)&d_tcin, tcinBytes));
CHECK(cudaMalloc((double**)&d_num, numBytes));
CHECK(cudaMalloc((double**)&d_den1, den1Bytes));
CHECK(cudaMalloc((double**)&d_den2, den2Bytes));
CHECK(cudaMalloc((double**)&d_noiseFinIn, noiseFinBytes));
CHECK(cudaMalloc((double**)&d_tcout, tcinBytes));
CHECK(cudaMalloc((double**)&d_noiseFinOut, noiseFinBytes));
// Transfer data from host to device
CHECK(cudaMemcpy(d_tcin, tcin, tcinBytes, cudaMemcpyHostToDevice));
CHECK(cudaMemcpy(d_num, num, numBytes, cudaMemcpyHostToDevice));
CHECK(cudaMemcpy(d_den1, den1, den1Bytes, cudaMemcpyHostToDevice));
CHECK(cudaMemcpy(d_den2, den2, den2Bytes, cudaMemcpyHostToDevice));
CHECK(cudaMemcpy(d_noiseFinIn, noiseFinIn, noiseFinBytes, cudaMemcpyHostToDevice));
//CHECK(cudaDeviceSetSharedMemConfig(cudaSharedMemBankSizeEightByte));
CHECK(cudaDeviceSetCacheConfig(cudaFuncCachePreferL1));
// Transfer common constant information in constant memory:
// - Daubechies D6 coefficients
// - Numerator and denominator filter coefficients
set_constant_memory(g6, h6, p_unity);
// Call CUDA kernel
int blockSize = 5;
int numBlocks = (dimsTcin[1] + blockSize -1) / blockSize;
printf("blocksize is %d and numBlocks is %d\n", blockSize, numBlocks);
dim3 dimGrid(numBlocks, 1);
dim3 dimBlock(blockSize, 1);
my_kernel_wrapper(dimGrid, dimBlock,
d_tcin, dimsTcin[0], dimsTcin[1],
d_num, dimsNum[1],
lden,
d_den1, dimsDen1[1],
d_den2, dimsDen2[1],
lambdaTemp,
maxeig,
cost_save,
nit,
d_noiseFinIn,
d_tcout,
d_noiseFinOut);
// Wait for GPU to finish before accessing on host
CHECK(cudaDeviceSynchronize());
CHECK(cudaMemcpy(tcout, d_tcout, tcinBytes, cudaMemcpyDeviceToHost));
CHECK(cudaMemcpy(noiseFinOut, d_noiseFinOut, noiseFinBytes, cudaMemcpyDeviceToHost));
CHECK(cudaFree(d_tcin));
CHECK(cudaFree(d_num));
CHECK(cudaFree(d_den1));
CHECK(cudaFree(d_den2));
CHECK(cudaFree(d_tcout));
CHECK(cudaFree(d_noiseFinIn));
CHECK(cudaFree(d_noiseFinOut));
// Reset device
CHECK(cudaDeviceReset());
return;
}
c4science · Help