How can I display information in the parallel.gpu.CUDAKernel?

Dear all,

I am a freshman to GPU computing within Matlab and get one question for help. In the following cudekernel, I want to add some codes to displaydisplay some information, for example, "hello" or "computation start", when debugging this kernel.But the command line

      printf(" hello \n");

does not work within Matlab, while it works in the VS2012 platform. Can somebody help me? Many thanks to your potential suggestions!

Shan Yin

attachments #1 Cudakernel

// My own CUDA code to calculate the bifurcation diagram of sing DOF impact oscillator on a GPU.

__global__ void sub_function( double *Matr_0, double *Matr_1, double *Matr_2, double *d_para, double x0, double v0, 
                              const double w, const double p, const double r, const double theta, 
                              const int date_length, const int para_length)
{
    int total_index_num=blockIdx.x*blockDim.x+threadIdx.x;

    const int NUM=4000;

    if( total_index_num < para_length )
    {
       //printf(" hello \n");  // it does not work!

       double d=d_para[total_index_num];
       
       const double pi=3.141592653589793;
       const double T=2*pi/w;
       
       int    flag0=0;
	double tao0=0;
	double dt0=0;
	double dt=0;
       double t=0;
       double sum=0;
       double a=0;
       double b=0;
       double c1=0;
       double c2=0;       
       double x=0;
       double v=0;

       if( T/2000<=0.005 )
	{
	    dt0=T/2000;
       }
       else
       {
           dt0=0.005;
	}

       for(int ii=0; ii<NUM; ii=ii+1)
       {
           flag0=1;
           tao0=theta;
           dt=dt0;
           t=dt;
           sum=0;

           while(flag0>0)
	   {
              a=(1-w*w)/((1-w*w)*(1-w*w)+(2*p*w)*(2*p*w));
              b=(-2*p*w)/((1-w*w)*(1-w*w)+(2*p*w)*(2*p*w));
              c1=x0-a*sin(tao0)-b*cos(tao0);
              c2=(v0+p*x0+(b*w-a*p)*sin(tao0)-(a*w+b*p)*cos(tao0))/(sqrt(1-p*p));
              x=exp(-p*t)*(c1*cos(sqrt(1-p*p)*t)+c2*sin(sqrt(1-p*p)*t))+a*sin(w*t+tao0)+b*cos(w*t+tao0);
              v=-p*exp(-p*t)*c1*cos(sqrt(1-p*p)*t)-c1*(sqrt(1-p*p))*exp(-p*t)*sin(sqrt(1-p*p)*t)+c2*(sqrt(1-p*p))*cos(sqrt(1-p*p)*t)*exp(-p*t)
                            -p*exp(-p*t)*c2*sin(sqrt(1-p*p)*t)+a*w*cos(w*t+tao0)-b*w*sin(w*t+tao0);
              while(x<=d)
              {
		   if(t<(T-sum))
		   {
			t=t+dt;
                       x=exp(-p*t)*(c1*cos(sqrt(1-p*p)*t)+c2*sin(sqrt(1-p*p)*t))+a*sin(w*t+tao0)+b*cos(w*t+tao0);
                       v=-p*exp(-p*t)*c1*cos(sqrt(1-p*p)*t)-c1*(sqrt(1-p*p))*exp(-p*t)*sin(sqrt(1-p*p)*t)+c2*(sqrt(1-p*p))*cos(sqrt(1-p*p)*t)*exp(-p*t)
                                      -p*exp(-p*t)*c2*sin(sqrt(1-p*p)*t)+a*w*cos(w*t+tao0)-b*w*sin(w*t+tao0);
		    }
	            else
		    {
			  t=T-sum;
		          x=exp(-p*t)*(c1*cos(sqrt(1-p*p)*t)+c2*sin(sqrt(1-p*p)*t))+a*sin(w*t+tao0)+b*cos(w*t+tao0);
                          v=-p*exp(-p*t)*c1*cos(sqrt(1-p*p)*t)-c1*(sqrt(1-p*p))*exp(-p*t)*sin(sqrt(1-p*p)*t)+c2*(sqrt(1-p*p))*cos(sqrt(1-p*p)*t)*exp(-p*t)
                                     -p*exp(-p*t)*c2*sin(sqrt(1-p*p)*t)+a*w*cos(w*t+tao0)-b*w*sin(w*t+tao0);
			   x0=x;
                          v0=v;
                          flag0=-1;
                         if( ii>=NUM-date_length )
                         {
                              Matr_0[total_index_num*date_length+(ii-NUM+date_length)]=d;
                              Matr_1[total_index_num*date_length+(ii-NUM+date_length)]=x;
                              Matr_2[total_index_num*date_length+(ii-NUM+date_length)]=v;
                          }
                          break;
		      }
	         }
		 if(x>d)
		 {
		      while (dt>1e-15 && abs(x-d)>1e-12)
		      {
			   if (x<d)
		           {
			       t=t+dt;
		           }
		           else
		           {
			       t=t-dt;
                              dt=dt/2;
                              t=t+dt;
			    }
			    x=exp(-p*t)*(c1*cos(sqrt(1-p*p)*t)+c2*sin(sqrt(1-p*p)*t))+a*sin(w*t+tao0)+b*cos(w*t+tao0);
		       }
		       sum=sum+t;
		       v=-p*exp(-p*t)*c1*cos(sqrt(1-p*p)*t)-c1*(sqrt(1-p*p))*exp(-p*t)*sin(sqrt(1-p*p)*t)+c2*(sqrt(1-p*p))*cos(sqrt(1-p*p)*t)*exp(-p*t)
                                -p*exp(-p*t)*c2*sin(sqrt(1-p*p)*t)+a*w*cos(w*t+tao0)-b*w*sin(w*t+tao0);
		       tao0=w*t+tao0-floor((w*t+tao0)/(2*pi))*2*pi; //tao0=mod(w*t+tao0,2*pi);
		       dt=dt0;
                       t=dt;
                       x0=x;
		        v0=-r*v;
		  }
	     }
        }
    }
}   

attachments #2 Main function

clear variables
clc

p=0.01;
r=0.68;
w0=sqrt(1-p*p)*(2.0/3.0); 
w=w0+5e-4;
A=(1-w*w)/((1-w*w)*(1-w*w)+(2*p*w)*(2*p*w));
B=(-2*p*w)/((1-w*w)*(1-w*w)+(2*p*w)*(2*p*w));
tao0=pi+atan(A/B);
d0=sqrt(A*A+B*B);
theta=tao0+pi/2;
xs=[A*sin(theta)+B*cos(theta),A*w*cos(theta)-B*w*sin(theta)];

d_left=d0-2e-5;
d_right=d0+2e-5;

date_length=50;  % dates recorded in each thread.
para_length=100;

d_para=gpuArray(linspace(d_left,d_right,para_length));
Matr_0=zeros(para_length*date_length, 1, 'gpuArray');
Matr_1=Matr_0;
Matr_2=Matr_0;

!nvcc -ptx bifur_kernel.cu

kernel=parallel.gpu.CUDAKernel( 'bifur_kernel.ptx', 'bifur_kernel.cu' );
kernel.ThreadBlockSize=[kernel.MaxThreadsPerBlock,1,1];    
kernel.GridSize=[ceil(para_length/kernel.ThreadBlockSize(1)),1]

[Matr_0, Matr_1, Matr_2]=feval(kernel, Matr_0, Matr_1, Matr_2, d_para, xs(1), xs(2), ...
                                    w, p, r, theta, date_length, para_length);

[Matr_0,Matr_1,Matr_2]=gather(Matr_0,Matr_1,Matr_2); 

figure(1)
plot(Matr_0-d0,Matr_1,'.b')