Do we have a vectorized random number generator?

General Discussion
1

For my model problem I need a P0 function with Gaussian random values on each triangle. So I use the code:

fespace Nh(Th, P0);
Nh dW;

for (int j=0; j<Th.nt; ++j) {
   dw[][j] = gslrangaussian(ffrng, 1.0); 
}

This is inside other loops and is the most time consuming part of my code (much longer than assembly and solve). For example, generate 15000x100 numbers takes about 2 seconds, while the randn(15000,100) command in MATLAB takes 0.02. I look through gsl docs and find no way to generate a random vector/matrix without using loops. Importing data from MATLAB will not help because the I/O will become the new bottleneck.

2

No to day but is simple to add this feature.

to speed up you can try

for [j,dwj:dw[]] dwj  = gslrangaussian(ffrng, 1.0);

the code

load "gsl" 
mesh Th = square(100,100); 
fespace Nh(Th, P0);
Nh dW;
real cpu0 = clock(); 
gslrng ffrng=gslrngtype(1);
for (int j=0; j<Th.nt; ++j) {
   dW[][j] = gslrangaussian(ffrng, 1.0); 
}

cout <<  clock() - cpu0 << endl; ; 
 cpu0 = clock(); 
 
 for [i,dWi: dW[]]
   dWi  = gslrangaussian(ffrng, 1.0); 
 cout <<  clock() - cpu0 << endl; ;

The result on my lap top Apple M1

  -- Square mesh : nb vertices  =10201 ,  nb triangles = 20000 ,  nb boundary edges 400
0.017193
0.000999
times: compile 0.010353s, execution 0.019229s,  mpirank:0
3

Thanks a lot! It is quite interesting that the efficiency gain is lost if the statement

dWi  = gslrangaussian(ffrng, 1.0);

is enclosed in a {}. What’s special about the syntax for [i,dWi: dW[]]?

A related question is how to fix the seed for the gsl random number generator? I want the results reproducible at different runs and in different machines by other researchers. The randinit function has no effect.

4

the syntaxe

for [i,dWi: dW[]]  code

Is a faster tool to set array.

5

I am trying to generate normally distributed random numbers using gslrangaussian(ffrng, 1.0); which is using gsl library(load “gsl”). However the numbers appearing on the first run of code are the same as those appearing in the next executions of the code. How to overcome this?

Remark: I have tried using the command srandomdev(); but its showing this error “srandomdev The Identifier srandomdev does not exist”.

Here is the code

load “gsl”
srandomdev(); // set a true ramdom seed …
gslrng ffrng=gslrngtype(1);
for(int i = 0; i < 10; i++){
cout << " ranx = " << gslrangaussian(ffrng, 1.0) << endl;
cout << " \n " << endl;
}

Thank you in advance,
Jitendra

6

the code srandomdev is in plugin ffrandom
see FreeFem-sources/examples/plugin/ffrandom.edp at master · FreeFem/FreeFem-sources · GitHub

7

Dear sir,
Thank you for your response. I have tried using load “ffrandom” but for each run I’m getting the same normally distributed random numbers. Is there any other way or any code in plugin gsl to generate different sets of normally distributed random numbers in each run of the code?

Here is the code I have tried

  load "ffrandom"
  load "gsl"
  srandomdev(); // set a true ramdom seed ..
  gslrng ffrng=gslrngtype(1);
  for(int i = 0; i < 5; i++){
    cout << " ranx = " << gslrangaussian(ffrng, 1.0) << endl;
  }

Thank you.

8

It seems to me that GSL does not depend on ffrandom features, but on its own. You may change the random generator seed by setting the GSL_RNG_SEED environment variable, see GSL documentation for an example : Random Number Distributions — GSL 2.7 documentation

I’m not sure there is an embedded method to change the seed directly in your FreeFem code though - Let me know if you find any !

Best regards,

Caroline.

9

I encountered this issue lately and found a way to set the GSL seed in FreeFem.
gslrng ffrng=gslrngtype(1); // Set up the random number generator
gslrngset(ffrng, myseed); // Seed generator with myseed