modifying minTV to minP when calculating p norm of the magnitude gradient images

Hello:
I'm trying to realize p-norm version of the ASD-POCS, and complete it with CUDA. the matlab code works well while the CUDA code comes kinds of errors all the time. There maybe comes NaN voxels, which you said maybe the problem of eps, but it worked nothing when I set eps=0.0001, and maybe came the problem of the input parameter couldn't be the single one, and sometimes, it reminded me that the MinP wasn't a available function although I had compiled and got the mexw64 file in the proper path. could you please tell me where the problem is? the only change of minP.mexw64 to the minTV.mexw64 is an additional parameter p. Here is the related files.

1. OS-ASD-POCS:
`f=minimizeP(f0,dtvg,ng,p);`

2. minimizeP.m:
```
if nargin==1
dtvg=1;ng=30;p=1;
else
if nargin == 4
dtvg=varargin{1};
ng=varargin{2};
p=varargin{3};
else
error('Wrogn amount of inputs');
end
end
img=(permute(img,[3 2 1]));
img=minP(img,dtvg,ng,p);
img=(permute(img,[3 2 1]));
end
```

3. minP.cpp
```
#include "tmwtypes.h"
#include "mex.h"
#include
#include
#include "matrix.h"
#include "POCS_P.hpp"
#include
void mexFunction(int nlhs , mxArray *plhs[],
int nrhs, mxArray const *prhs[])
{
///////// First check if the amount of imputs is rigth.
int maxIter;
float alpha;
float p;
if (nrhs==1){
maxIter=100;
alpha=15.0f;
p=1.0f;
}
if (nrhs==2){
mexErrMsgIdAndTxt("err", "Only 1 POCS hyperparemter inputed");
}
if (nrhs==3){
mexErrMsgIdAndTxt("err", "Only 2 POCS hyperparemter inputed");
}
if (nrhs>4){
mexErrMsgIdAndTxt("err", "Too many imput argumets");
}
if (nrhs==4){
size_t mrows = mxGetM(prhs[1]);
size_t ncols = mxGetN(prhs[1]);
if (mrows!=1 || ncols !=1)
mexErrMsgIdAndTxt("err", "POCS parameters shoudl be 1x1");
mrows = mxGetM(prhs[2]);
ncols = mxGetN(prhs[2]);
if (mrows!=1 || ncols !=1)
mexErrMsgIdAndTxt("err", "POCS parameters shoudl be 1x1");
mrows = mxGetM(prhs[3]);
ncols = mxGetN(prhs[3]);
if (mrows!=1 || ncols !=1)
mexErrMsgIdAndTxt("err", "POCS parameters shoudl be 1x1");
alpha= (float)(mxGetScalar(prhs[1]));
maxIter=(int)floor(mxGetScalar(prhs[2])+0.5);
p=(float)(mxGetScalar(prhs[3]));
}

// First input should be x from (Ax=b), or the image.
mxArray const * const image = prhs[0];
mwSize const numDims = mxGetNumberOfDimensions(image);

// Image should be dim 3
if (numDims!=3){
mexErrMsgIdAndTxt("err", "Image is not 3D");
}
// Now that input is ok, parse it to C data types.
float const * const imgaux = static_cast(mxGetData(image));
const mwSize *size_img= mxGetDimensions(image); //get size of image

float * img = (float*)malloc(size_img[0] *size_img[1] *size_img[2]* sizeof(float));

for(int i=0;i= 0 && z= 0 && y= 0 && x= cols || y >= rows || z >= depth )
return;

float df[3] ={0,0,0};
float dfi[3]={0,0,0}; // dfi== \partial f_{i+1,j,k}
float dfj[3]={0,0,0};
float dfk[3]={0,0,0};
gradient(f,df ,z ,y ,x , depth,rows,cols);
gradient(f,dfi ,z ,y ,x+1, depth,rows,cols);
gradient(f,dfj ,z ,y+1,x , depth,rows,cols);
gradient(f,dfk ,z+1,y ,x , depth,rows,cols);
float eps=0.0001; //% avoid division by zero
dftv[idx]=(df[0]+df[1]+df[2])*((pow(df[0] *df[0] +df[1] *df[1] +df[2] *df[2],(p-2)/2))+eps)
-dfi[2]*((pow(dfi[0]*dfi[0]+dfi[1]*dfi[1]+dfi[2]*dfi[2],(p-2)/2))+eps)
-dfj[1]*((pow(dfj[0]*dfj[0]+dfj[1]*dfj[1]+dfj[2]*dfj[2],(p-2)/2))+eps)
-dfk[0]*((pow(dfk[0]*dfk[0]+dfk[1]*dfk[1]+dfk[2]*dfk[2],(p-2)/2))+eps);

}

__device__ void warpReduce(volatile float *sdata, size_t tid) {
sdata[tid] += sdata[tid + 32];
sdata[tid] += sdata[tid + 16];
sdata[tid] += sdata[tid + 8];
sdata[tid] += sdata[tid + 4];
sdata[tid] += sdata[tid + 2];
sdata[tid] += sdata[tid + 1];
}

__global__ void reduceNorm2(float *g_idata, float *g_odata, size_t n){
extern __shared__ volatile float sdata[];
//http://stackoverflow.com/a/35133396/1485872
size_t tid = threadIdx.x;
size_t i = blockIdx.x*blockDim.x + tid;
size_t gridSize = blockDim.x*gridDim.x;
float mySum = 0;
float value=0;
while (i < n) {
value=g_idata[i]; //avoid reading twice
mySum += value*value;
i += gridSize;
}
sdata[tid] = mySum;
__syncthreads();

if (tid < 512)
sdata[tid] += sdata[tid + 512];
__syncthreads();
if (tid < 256)
sdata[tid] += sdata[tid + 256];
__syncthreads();

if (tid < 128)
sdata[tid] += sdata[tid + 128];
__syncthreads();

if (tid < 64)
sdata[tid] += sdata[tid + 64];
__syncthreads();


#if (__CUDA_ARCH__ >= 300)
if ( tid < 32 )
{
mySum = sdata[tid] + sdata[tid + 32];
for (int offset = warpSize/2; offset > 0; offset /= 2) {
mySum += __shfl_down(mySum, offset);
}
}
#else
if (tid < 32) {
warpReduce(sdata, tid);
mySum = sdata[0];
}
#endif
if (tid == 0) g_odata[blockIdx.x] = mySum;
}
__global__ void reduceSum(float *g_idata, float *g_odata, size_t n){
extern __shared__ volatile float sdata[];
//http://stackoverflow.com/a/35133396/1485872
size_t tid = threadIdx.x;
size_t i = blockIdx.x*blockDim.x + tid;
size_t gridSize = blockDim.x*gridDim.x;
float mySum = 0;
// float value=0;
while (i < n) {
mySum += g_idata[i];
i += gridSize;
}
sdata[tid] = mySum;
__syncthreads();

if (tid < 512)
sdata[tid] += sdata[tid + 512];
__syncthreads();
if (tid < 256)
sdata[tid] += sdata[tid + 256];
__syncthreads();

if (tid < 128)
sdata[tid] += sdata[tid + 128];
__syncthreads();

if (tid < 64)
sdata[tid] += sdata[tid + 64];
__syncthreads();


#if (__CUDA_ARCH__ >= 300)
if ( tid < 32 )
{
mySum = sdata[tid] + sdata[tid + 32];
for (int offset = warpSize/2; offset > 0; offset /= 2) {
mySum += __shfl_down(mySum, offset);
}
}
#else
if (tid < 32) {
warpReduce(sdata, tid);
mySum = sdata[0];
}
#endif
if (tid == 0) g_odata[blockIdx.x] = mySum;
}




// main function
void pocs_p(const float* img,float* dst,float alpha,const long* image_size, int maxIter,float p){

size_t total_pixels = image_size[0] * image_size[1] * image_size[2] ;
size_t mem_size = sizeof(float) * total_pixels;

float *d_image, *d_dimgTV,*d_norm2aux,*d_norm2;
// memory for image
cudaMalloc(&d_image, mem_size);
cudaCheckErrors("Malloc Image error");
cudaMemcpy(d_image, img, mem_size, cudaMemcpyHostToDevice);
cudaCheckErrors("Memory Malloc and Memset: SRC");
// memory for df
cudaMalloc(&d_dimgTV, mem_size);
cudaCheckErrors("Memory Malloc and Memset: TV");

cudaMalloc(&d_norm2, mem_size);
cudaCheckErrors("Memory Malloc and Memset: TV");

// memory for L2norm auxiliar
cudaMalloc(&d_norm2aux, sizeof(float)*(total_pixels + MAXTHREADS - 1) / MAXTHREADS);
cudaCheckErrors("Memory Malloc and Memset: NORMAux");

// For the gradient
dim3 blockGrad(10, 10, 10);
dim3 gridGrad((image_size[0]+blockGrad.x-1)/blockGrad.x, (image_size[1]+blockGrad.y-1)/blockGrad.y, (image_size[2]+blockGrad.z-1)/blockGrad.z);

// For the reduction
float sumnorm2;

for(unsigned int i=0;i>>(d_image,d_dimgTV,image_size[2], image_size[1],image_size[0],p);
cudaCheckErrors("Gradient");

cudaMemcpy(d_norm2, d_dimgTV, mem_size, cudaMemcpyDeviceToDevice);
cudaCheckErrors("Copy from gradient call error");
// Compute the L2 norm of the gradint. For that, reduction is used.
//REDUCE
size_t dimblockRed = MAXTHREADS;
size_t dimgridRed = (total_pixels + MAXTHREADS - 1) / MAXTHREADS;
reduceNorm2 << > >(d_norm2, d_norm2aux, total_pixels);
cudaCheckErrors("reduce1");
if (dimgridRed > 1) {
reduceSum << <1, dimblockRed, MAXTHREADS*sizeof(float) >> >(d_norm2aux, d_norm2, dimgridRed);
cudaCheckErrors("reduce2");
cudaMemcpy(&sumnorm2, d_norm2, sizeof(float), cudaMemcpyDeviceToHost);
cudaCheckErrors("cudaMemcpy");

}
else {
cudaMemcpy(&sumnorm2, d_norm2aux, sizeof(float), cudaMemcpyDeviceToHost);
cudaCheckErrors("cudaMemcpy");
}
//mexPrintf("%f ",sqrt(sumnorm2));
//NOMRALIZE
//in a Tesla, maximum blocks =15 SM * 4 blocks/SM
divideArrayScalar <<<60,MAXTHREADS>>>(d_dimgTV,sqrt(sumnorm2),total_pixels);
cudaCheckErrors("Division error");
//MULTIPLY HYPERPARAMETER
multiplyArrayScalar<<<60,MAXTHREADS>>>(d_dimgTV,alpha, total_pixels);
cudaCheckErrors("Multiplication error");
//SUBSTRACT GRADIENT
substractArrays <<<60,MAXTHREADS>>>(d_image,d_dimgTV, total_pixels);
cudaCheckErrors("Substraction error");
sumnorm2=0;
}

cudaCheckErrors("TV minimization");

cudaMemcpy(dst, d_image, mem_size, cudaMemcpyDeviceToHost);
cudaCheckErrors("Copy result back");

cudaFree(d_image);
cudaFree(d_norm2aux);
cudaFree(d_dimgTV);
cudaFree(d_norm2);

cudaCheckErrors("Memory free");

}

```