Note that there are some explanatory texts on larger screens.

plurals
  1. PO
    text
    copied!<p>Solution:</p> <pre><code>#include &lt;iostream&gt; #include &lt;curand_kernel.h&gt; using namespace std; /**** ERROR HANDLING ****/ static void HandleError(cudaError_t err,const char *file, int line ) { if (err != cudaSuccess) { printf( "%s in %s at line %d\n", cudaGetErrorString( err ), file, line ); system("pause"); exit( EXIT_FAILURE ); } } #define HANDLE_ERROR( err ) (HandleError( err, __FILE__, __LINE__ )) /**** HOST AND DEVICE CONSTANTS****/ const int hNP=100, hD=31, hN=10; __constant__ int NP, D, N; __constant__ float Cr, F; /*** EVAL FUNCTION******/ __device__ float lennardJones(float a[3], float b[3]) { float distance = sqrt((a[0] - b[0]) * (a[0] - b[0]) + (a[1] - b[1]) * (a[1] - b[1]) + (a[2] - b[2]) * (a[2] - b[2])); float distance6 = distance * distance * distance * distance * distance * distance; float distance12 = distance6 * distance6; return 1/distance12 - 2/distance6; } /**** RANDOM GENERATORS***/ __device__ float rndFloat(curandState* globalState, int id) { curandState localState = globalState[id]; float RANDOM = curand_uniform(&amp;localState); globalState[id] = localState; return RANDOM; } __device__ int rndInt(curandState* globalState, int id, int max) { curandState localState = globalState[id]; float RANDOM = curand_uniform(&amp;localState); globalState[id] = localState; return RANDOM*max; } __device__ float rndFloat(curandState* globalState, int id, int max) { curandState localState = globalState[id]; float RANDOM = curand_uniform(&amp;localState); globalState[id] = localState; return RANDOM*max; } __device__ float rndFloat(curandState* globalState, int id, int min,int max) { curandState localState = globalState[id]; float RANDOM = curand_uniform(&amp;localState); globalState[id] = localState; return min+RANDOM*(max-min); } /*** SEEDS ****/ __global__ void setup_kernel (curandState * state, unsigned long seed) { int id= threadIdx.x+blockIdx.x*blockDim.x; if(id &lt; NP) curand_init(seed, id, 0,&amp;state[id]); } /**** DIFFERENTIAL EVOLUTION: INITIALIZATION ***/ __global__ void kernelE(curandState* globalState, float *population) { int id= threadIdx.x+blockIdx.x*blockDim.x; if(id &lt; NP) { //init, just populating array with some specific numbers population[D*id]=0; population[D*id+N]=0; population[D*id +2*N]=0; population[D*id+1]=rndFloat(globalState,threadIdx.x,4); population[D*id+N+1]=0; population[D*id +2*N+1]=0; for(int i=2; i&lt;N; i++){ float min= -4 - 1/4*abs((int)((i-4)/3)); float max= 4 + 1/4*abs((int)((i-4)/3)); if(i==2) { population[D*id+2]=rndFloat(globalState,threadIdx.x,3.14159265359); population[D*id+N+2]=rndFloat(globalState,threadIdx.x,min,max); population[D*id +2*N+2]=0; } else { population[D*id +i]=rndFloat(globalState,threadIdx.x,min,max); population[D*id+N+i]=rndFloat(globalState,threadIdx.x,min,max); population[D*id +2*N+i]=rndFloat(globalState,threadIdx.x,min,max); } } //eval float e=0; for(int i=0; i&lt;N; i++) { for(int j=0; j&lt;i; j++) { float a[]={population[D*id +i], population[D*id+N+i], population[D*id +2*N+i]}, b[]={population[D*id +j],population[D*id +j+N], population[D*id +2*N+j]}; e += lennardJones(a,b); } } population[D*id + D-1]=e; } } /**** DIFFERENTIAL EVOLUTION: MUTATION INDICES ****/ __global__ void kernelP(curandState* globalState, int *mutation) { int id= threadIdx.x+blockIdx.x*blockDim.x; if(id&lt;NP) { int a = rndInt(globalState, id, NP),b = rndInt(globalState, id, NP),c= rndInt(globalState, id, NP); while(a == id){a = rndInt(globalState, id, NP);} while(b == a &amp;&amp; b==id){b=rndInt(globalState, id, NP);} while(c == a &amp;&amp; c== b &amp;&amp; c ==id){c=rndInt(globalState, id, NP);} mutation[3*id+0]=a; mutation[3*id+1]=b; mutation[3*id+2]=c; } } /**** DIFFERENTIAL EVOLUTION: MUTATION, CROSSOVER, EVALUATION AND SELECTION ***/ __global__ void kernelMCER(curandState* globalState, float *population, int *mutation, float *pop) { int id= threadIdx.x+blockIdx.x*blockDim.x; if(id&lt;NP) { int a=mutation[3*id+0], b=mutation[3*id+1], c=mutation[3*id+2]; //DE mutation and crossover int j=rndInt(globalState, id, NP); for(int i=0; i&lt;D-1; i++) { //DE mutation pop[D*id+i]= population[D*a +i] + F*(population[D*b +i]-population[D*c +i]); //DE crossover if(Cr &gt; rndFloat(globalState, id) &amp;&amp; i!= j) pop[D*id+i]=population[D*id +i]; } // Eval pop[D*id+D-1]=0; for(int i=0; i&lt;N; i++) { for(int j=0; j&lt;i; j++) { float a[]={pop[D*id+i], pop[D*id+N+i], pop[D*id+2*N+i]}, b[]={pop[D*id+j],pop[D*id+N+j], pop[D*id+2*N+j]}; pop[D*id+D-1] += lennardJones(a,b); } } __syncthreads(); //DE selection if(pop[D*id+D-1] &lt; population[D*id +D-1]) { for(int i=0; i&lt;D; i++) population[D*id +i]=pop[D*id+i]; } } } void getBestScore(float *hpopulation) { int max=0; for(int i=1; i&lt;hNP; i++) { if(hpopulation[hD*max+hD-1] &gt; hpopulation[hD*i+hD-1]) max=i; } for(int j=0; j&lt;hN; j++) cout&lt;&lt;"Atom "&lt;&lt;(j+1)&lt;&lt;": ("&lt;&lt;hpopulation[hD*max+j]&lt;&lt;", "&lt;&lt;hpopulation[hD*max+hN+j]&lt;&lt;", "&lt;&lt;hpopulation[hD*max+hN*2+j]&lt;&lt;") "&lt;&lt;endl; cout&lt;&lt;"Result: "&lt;&lt;hpopulation[hD*max+hD-1]&lt;&lt;endl; } int main() { cudaEvent_t start,stop; HANDLE_ERROR(cudaEventCreate(&amp;start)); HANDLE_ERROR(cudaEventCreate(&amp;stop)); HANDLE_ERROR(cudaEventRecord(start,0)); int device, st=100; float hCr=0.6f, hF=0.8f; cudaDeviceProp prop; HANDLE_ERROR(cudaGetDevice(&amp;device)); HANDLE_ERROR(cudaGetDeviceProperties(&amp;prop, device)); // int SN = prop.maxThreadsPerBlock; //512 threads per block //int SB = (hNP+(SN-1))/SN; //constants NP, D, N, Cr, F HANDLE_ERROR(cudaMemcpyToSymbol(N, &amp;hN, sizeof(int))); HANDLE_ERROR(cudaMemcpyToSymbol(NP, &amp;hNP, sizeof(int))); HANDLE_ERROR(cudaMemcpyToSymbol(D, &amp;hD, sizeof(int))); HANDLE_ERROR(cudaMemcpyToSymbol(F, &amp;hF, sizeof(float))); HANDLE_ERROR(cudaMemcpyToSymbol(Cr, &amp;hCr, sizeof(float))); //seeds curandState* devStates; HANDLE_ERROR(cudaMalloc (&amp;devStates, hNP*sizeof(curandState))); setup_kernel &lt;&lt;&lt; 1, hNP&gt;&gt;&gt; (devStates, 50); //population float *population, *pop; float hpopulation[hNP*hD]; HANDLE_ERROR(cudaMalloc((void**)&amp;population, hNP*hD*sizeof(float))); HANDLE_ERROR(cudaMalloc((void**)&amp;pop, hNP*hD*sizeof(float))); //mutation int *mutation, *mutation1; int *hmutation; HANDLE_ERROR(cudaHostAlloc((void**)&amp;hmutation, hNP*3*sizeof(int), cudaHostAllocDefault)); HANDLE_ERROR(cudaMalloc((void**)&amp;mutation, hNP*3*sizeof(int))); HANDLE_ERROR(cudaMalloc((void**)&amp;mutation1, hNP*3*sizeof(int))); //stream cudaStream_t stream_i, stream_j; HANDLE_ERROR(cudaStreamCreate(&amp;stream_i)); HANDLE_ERROR(cudaStreamCreate(&amp;stream_j)); kernelE&lt;&lt;&lt;1,hNP, 0,stream_i&gt;&gt;&gt;(devStates,population); kernelP&lt;&lt;&lt;1,hNP, 0,stream_j&gt;&gt;&gt;(devStates,mutation); while(st != 0) { /*** COPYING MUTATION INDICES***/ HANDLE_ERROR(cudaMemcpyAsync(hmutation, mutation,hNP*3*sizeof(int), cudaMemcpyDeviceToHost, stream_j)); HANDLE_ERROR(cudaMemcpyAsync(mutation1, hmutation,hNP*3*sizeof(int), cudaMemcpyHostToDevice, stream_i)); /**** CALLING KERNELS****/ kernelP&lt;&lt;&lt;1,hNP,0,stream_j&gt;&gt;&gt;(devStates,mutation); kernelMCER&lt;&lt;&lt;1,hNP,0,stream_i&gt;&gt;&gt;(devStates,population,mutation1,pop); st--; //HANDLE_ERROR(cudaStreamSynchronize(stream_i)); //HANDLE_ERROR(cudaMemcpy(hpopulation, population, hNP*hD*sizeof(float), cudaMemcpyDeviceToHost)); //getBestScore(hpopulation); //cin.get(); } HANDLE_ERROR(cudaStreamSynchronize(stream_i)); HANDLE_ERROR(cudaMemcpy(hpopulation, population, hNP*hD*sizeof(float), cudaMemcpyDeviceToHost)); getBestScore(hpopulation); cudaEventRecord(stop,0); cudaEventSynchronize(stop); float time; HANDLE_ERROR(cudaEventElapsedTime(&amp;time, start, stop)); cout&lt;&lt;endl&lt;&lt;"Tme: "&lt;&lt;time/1000&lt;&lt;"s"&lt;&lt;endl; HANDLE_ERROR(cudaEventDestroy(start)); HANDLE_ERROR(cudaEventDestroy(stop)); HANDLE_ERROR(cudaStreamDestroy(stream_i)); HANDLE_ERROR(cudaStreamDestroy(stream_j)); HANDLE_ERROR(cudaFree(population)); HANDLE_ERROR(cudaFree(pop)); HANDLE_ERROR(cudaFreeHost(hmutation)); HANDLE_ERROR(cudaFree(mutation1)); HANDLE_ERROR(cudaFree(devStates)); system("pause"); return 0; } </code></pre>
 

Querying!

 
Guidance

SQuiL has stopped working due to an internal error.

If you are curious you may find further information in the browser console, which is accessible through the devtools (F12).

Reload