48076

Generate random number within a function with cuRAND without preallocation

Question:

Is it possible to generate random numbers within a device function without preallocate all the states? I would like to generate and use them in "realtime". I need to use them for Monte Carlo simulations what are the most suitable for this purpose? The number generated below are single precision is it possible to have them in double precision?

#include <iostream> #include "cuda_runtime.h" #include "device_launch_parameters.h" #include <curand_kernel.h> __global__ void cudaRand(float *d_out, unsigned long seed) { int i = blockDim.x * blockIdx.x + threadIdx.x; curandState state; curand_init(seed, i, 0, &state); d_out[i] = curand_uniform(&state); } int main(int argc, char** argv) { size_t N = 1 << 4; float *v = new float[N]; float *d_out; cudaMalloc((void**)&d_out, N * sizeof(float)); // generate random numbers cudaRand << < 1, N >> > (d_out, time(NULL)); cudaMemcpy(v, d_out, N * sizeof(float), cudaMemcpyDeviceToHost); for (size_t i = 0; i < N; i++) { printf("out: %f\n", v[i]); } cudaFree(d_out); delete[] v; return 0; } <hr />

UPDATE

#include <iostream> #include "cuda_runtime.h" #include "device_launch_parameters.h" #include <curand_kernel.h> #include <ctime> __global__ void cudaRand(double *d_out) { int i = blockDim.x * blockIdx.x + threadIdx.x; curandState state; curand_init((unsigned long long)clock() + i, 0, 0, &state); d_out[i] = curand_uniform_double(&state); } int main(int argc, char** argv) { size_t N = 1 << 4; double *h_v = new double[N]; double *d_out; cudaMalloc((void**)&d_out, N * sizeof(double)); // generate random numbers cudaRand << < 1, N >> > (d_out); cudaMemcpy(h_v, d_out, N * sizeof(double), cudaMemcpyDeviceToHost); for (size_t i = 0; i < N; i++) printf("out: %f\n", h_v[i]); cudaFree(d_out); delete[] h_v; return 0; }

Answer1:

How I was dealing with the similar situation in the past, within __device__/__global__ function:

int tId = threadIdx.x + (blockIdx.x * blockDim.x); curandState state; curand_init((unsigned long long)clock() + tId, 0, 0, &state); double rand1 = curand_uniform_double(&state); double rand2 = curand_uniform_double(&state);

So just use curand_uniform_double for generating random doubles and also I believe you don't want the same seed for all of the threads, thats what I am trying to achieve by using clock() + tId instead. This way the odds of having the same rand1/rand2 in any of the two threads are close to nil.

<strong>EDIT:</strong>

However, based on below comments, proposed approach <em>may</em> <em>perhaps</em> lead to biased result:

<ul><li>

JackOLantern pointed me to this part of curand documentation:

<blockquote>

Sequences generated with different seeds usually do not have statistically correlated values, but some choices of seeds <strong>may give statistically correlated sequences</strong>.

</blockquote></li> <li>

Also there is a <a href="https://devtalk.nvidia.com/default/topic/480586/curand-initialization-time/" rel="nofollow">devtalk thread</a> devoted to how to improve performance of curand_init in which the proposed solution to speed up the curand initialization is:

<blockquote>

One thing you can do is use different seeds for each thread and a fixed subsequence of 0 and offset of 0.

</blockquote>

But the same poster is later stating:

<blockquote>

The downside is that you lose some of the nice mathematical properties between threads. It is possible that there is a bad interaction between the hash function that initializes the generator state from the seed and the periodicity of the generators. If that happens, you might get two threads with highly correlated outputs for some seeds. <strong>I don't know of any problems like this, and even if they do exist they will most likely be rare</strong>.

</blockquote></li> </ul>

So it is basically up to you whether you want better performance (as I did) or 1000% unbiased results. If that is what you desire, then solution proposed by JackOLantern is the correct one, i.e. initialize curand as:

curand_init((unsigned long long)clock(), tId, 0, &state)

Using not 0 value for offset and subsequence parameters is, however, decreasing performance. For more info on these parameters you may review <a href="https://stackoverflow.com/q/25827825/3242721" rel="nofollow">this SO thread</a> and also <a href="http://docs.nvidia.com/cuda/curand/host-api-overview.html#generator-options" rel="nofollow">curand documentation</a>.

I see that JackOLantern stated in comment that:

<blockquote>

I would say it is not recommandable to call curand_init and curand_uniform_double from within the same kernel from two reasons ........ Second, curand_init initializes the pseudorandom number generator and sets all of its parameters, so I'm afraid your approach will be somewhat slow.

</blockquote>

I was dealing with this in my thesis on several pages, tried various approaches to get different random numbers in each thread and creating curandState in each of the threads turned out to be the most viable solution for me. I needed to generate ~10 random numbers in each thread and among others I tried:

<ul><li>developing my own simple random number generator (Linear Congruential Generator) whose intialization was basically for free, however, the performance suffered greatly when generating numbers, so in the end having curandState in each thread turned out to be superior,</li> <li>pre-allocating curandStates and reusing them - this was memory heavy and when I decreased number of preallocated states then I had to use non zero values for offset/subsequence parameters of curand_uniform_double in order to get rid of bias which led to decreased performance when generating numbers.</li> </ul>

So after making thorough analysis I decided to indeed call curand_init and curand_uniform_double in each thread. The only problem was with the amount of registry that these states were occupying so I had to be careful with the block sizes not to exceed the max number of registry available to each block.

Thats what I have to say about provided solution which I was finally able to test and it is working just fine on my machine/GPU. I run the code from <strong>UPDATE</strong> section in the above question and 16 different random numbers were displayed in the console correctly. Therefore I advise you to properly perform error checking after executing kernel to see what went wrong inside. This topic is very well covered in <a href="https://stackoverflow.com/q/14038589/3242721" rel="nofollow">this SO thread</a>.

Recommend

  • User Input + Random Word and Number printing
  • ruby / sort an array of hashes from a defined index
  • proto2 with Spark cannot run
  • .NET math calculation performances
  • How to restrict template functor return and parameter types
  • Task getting rejected although threads are available
  • Drupal: Hierarchical taxonomical breadcrumb trail
  • Converting a list of strings to ints (or doubles) in Python
  • float vs double comparison [duplicate]
  • PHP curl_multi_exec output to array
  • event not defined in Firefox
  • Post comment to WordPress Blog from iPhone programmatically
  • $this->params returns null in cakephp model
  • Simple regex for domain names
  • Edit assembly language code in Visual Studio while stepping through each statement
  • Serverless Framework Dynamo DB Table Resource Definition with Sort Key
  • Whats the right place for testhelper-classes? (phpunit/best practise)
  • Skip Characters in Oracle TO_DATE function
  • In Akka, is ActorContext thread safe?
  • Efficient & Pythonic way of finding all possible sublists of a list in given range and the minim
  • Cloud Code function running twice
  • Django model inheritance, filtering models
  • Do I need to seed any random number generator before using EVP_PKEY_keygen of OpenSSL?
  • Diff between two dataframes in pandas
  • uniform generation of points on 3D box
  • Do I need to reset a Perl hash index?
  • Use of this Javascript
  • Record samples being played with OpenAL
  • C++ Partial template specialization - design simplification
  • Spray.io: When (not) to use non-blocking route handling?
  • Is possible to count alias result on mysql
  • Symfony2: How to get request parameter
  • How to get next/previous record number?
  • GridView Sorting works once only
  • Transpose CSV data with awk (pivot transformation)
  • Python: how to group similar lists together in a list of lists?
  • WPF Applying a trigger on binding failure
  • Benchmarking RAM performance - UWP and C#
  • How does Linux kernel interrupt the application?
  • git trying to push non-existent file … after clearing cache