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.