Números Aleatórios - Performance

Eu fiz um teste de performance de diferentes implementações de geradores de números aleatórios, comparando a biblioteca GSL e a MKL. O teste consistiu em calcular o tempo necessário para gera 2^30 números aleatórios (float):

1. implementação do Numerical Recipes do ran2                  ->  25 segundos
2. implementação da biblioteca GSL do Mersenne Twister  ->  15 segundos
3. implementação da biblioteca MKL do Mersenne Twister ->  5  segundos
   

 Este é o fonte do teste

Geradores de números pseudo-aleatórios

Geradores de números pseudo-aleatórios GNPA são uma caixa preta, em geral usamos as indicações dos estudantes mais antigos e dos professores. Provavelmente, graças a esse mecanismo, os GNPA mais usados são o ran2 ou o rand, e a implementação destes é encontrada no clássico Numerical Recipes.

O cerne das simulações de Monte-Carlo são justamente os GNPA e dificilmente nossos programas passarão fazendo mais tempo algo que não seja gerar números aleatórios. Por isso é importante conhecer e se atualizar quanto este assunto.

Embora seja importante entender um pouco da minúcias dos GNPA eu não vou nem tocar nesse assunto, para mim eles permaneceram uma caixa preta por muito tempo. Bom..., as atualidades ? Ah, essas são bem mais fáceis. 

O GNPA padrão nos dias de hoje o Mersenne-Twister, introduzido em 1997 http://en.wikipedia.org/wiki/Mersenne_twister (alguem tem que escrever o artigo na wikipedia portuguesa, nem que seja tradução). Este é provavelmente o algoritmo mais popular, ele possui um período de 2^19937 − 1, passou em vários testes e é superior aos ran2,3,4 .

Outro ponto de fundamental importância é a sua performance.

Como implementar o algoritmo:

Eu sugiro fortemente a utilização da biblioteca GSL ( Gnu Scientific Library ) http://www.gnu.org/software/gsl/ e o manual http://www.gnu.org/software/gsl/manual , essa biblioteca já esta instalada no abax e você pode encontra-la nos repositórias da sua distribuição Linux.

E muito fácil usar, segue um exemplo

#include < stdio.h >
#include < gsl/gsl_rng.h >
int main ()
  {
  int i;
  unsigned long a;
/////////////////////////////////////////////////////
  const gsl_rng *r;
  r = gsl_rng_alloc(gsl_rng_mt19937);
////////////////////////////////////////////////////
  gsl_rng_set(r, 123); // semente 123
   for (i = 0; i < 100000000; i++)
    {
    printf ("%u\n",  gsl_rng_get(r));
    }
  printf ("\n");
  printf ("%lu %lu\n", gsl_rng_max(r), gsl_rng_min(r) );
  return 0;
}

Compilar com: gcc rng.c -lgsl -lgslcblas -lm

Comentários:

"const gsl_rng *r;
r = gsl_rng_alloc(gsl_rng_mt19937); "
- Define o nome do gerador, no caso é "r",  e o GNPA usado ( mt19937 = Mersenne Twisted ).

"gsl_rng_set (r, 123);" - Semente (123) do gerador r.

"gsl_rng_get (r);" - Retorna um número inteiro (unsigned long) pseudo-aleatório

"gsl_rng_max(r) e gsl_rng_min(r)" - Valor máximo e mínimo do gerador.


Outras dicas:

"gsl_rng_uniform(r);" - Retorna um real (double) pseudo-aleatório o intervalo [0,1)
"gsl_rng_uniform_int(r,j);" Retorna um número inteiro (unsigned long) pseudo-aleatório entre [0, j-1]

Existe uma lista de extensa de implementações de GNPA nesta biblioteca, inclusive o ran2, basta substituir o "mt19937" pelo algoritmo desejado.

Veja a lista de performance ( http://www.gnu.org/software/gsl/manual/html_node/Random-Number-Generator-Performance.html ):

     1754 k ints/sec,    870 k doubles/sec, taus
     1613 k ints/sec,    855 k doubles/sec, gfsr4
     1370 k ints/sec,    769 k doubles/sec, mt19937
      565 k ints/sec,    571 k doubles/sec, ranlxs0
      400 k ints/sec,    405 k doubles/sec, ranlxs1
      490 k ints/sec,    389 k doubles/sec, mrg
      407 k ints/sec,    297 k doubles/sec, ranlux
      243 k ints/sec,    254 k doubles/sec, ranlxd1
      251 k ints/sec,    253 k doubles/sec, ranlxs2
      238 k ints/sec,    215 k doubles/sec, cmrg
      247 k ints/sec,    198 k doubles/sec, ranlux389
      141 k ints/sec,    140 k doubles/sec, ranlxd2
    
     1852 k ints/sec,    935 k doubles/sec, ran3
      813 k ints/sec,    575 k doubles/sec, ran0
      787 k ints/sec,    476 k doubles/sec, ran1
      379 k ints/sec,    292 k doubles/sec, ran2


Uma pequena explicação sobre as características de cada algoritmo é encontrada no manual http://www.gnu.org/software/gsl/manual/html_node/Random-number-generator-algorithms.html


Matrices Direct Product C GSL routine

This is my amateur routine of Matrices Direct Product or Kronecker Product using GSL. There is no advantage in doing this with GSL library, this is just a exercise.

#include <stdio.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_math.h>

double *produtodireto (int a,int b, double A[a*a],double B[b*b])
     {
 double *AxB;
 AxB = malloc(a*a*b*b*(sizeof *AxB));
 int i, j, k,l;
 gsl_matrix *M = gsl_matrix_alloc (a,a);
 gsl_matrix *L = gsl_matrix_alloc (b,b);
 for (i = 0; i < a; i++)
  {
  for (j = 0; j < a; j++)
   {
             gsl_matrix_set (M, i, j, A[i+a*j]);
   }
  }
 for (i = 0; i < b; i++)
  {
  for (j = 0; j < b; j++)
   {
      gsl_matrix_set (L, i, j, B[i+b*j]);
   }
  }
 gsl_matrix *R = gsl_matrix_alloc (a*b,a*b);
 for (i = 0; i < a; i++)
  {
  for (j = 0; j < a; j++)
   { 
   for (k = 0; k < b; k++)
    {
    for (l = 0; l < b; l++)
     {
     gsl_matrix_set (R, b*i+k-b+2,b*j+l-b+2, gsl_matrix_get (M, i, j)*gsl_matrix_get (L, k, l));
     }
    }
   }
  }
for (i = 0; i < a*b; i++)  
 {/* OUT OF RANGE ERRORRRRR */
 for (j = 0; j < a*b; j++)
  {
  AxB[i+a*b*j]=gsl_matrix_get (R, i, j);
  }
 }

gsl_matrix_free (M);
gsl_matrix_free (L);
gsl_matrix_free (R);
       return AxB;
     }

Using GSL to calculate Nonsymmetric matrices eigenvalues

Here is a example of using GSL - Gnu Scientific Library - to calculate the eigenvalues of a nonsymmetric matrix.
This is a simple example write in C. GSL uses BLAS and can be optimized with ATLAS, so can be used as a real tool for heavy calculations. To get know more about GSL visit http://www.gnu.org/software/gsl/ or the manual page.

Hint: Beware when using GSL or our own linear algebra routine, remember that greats matrices requires greats responsibilities, so keep track of precision errors. Find precision errors will be a great excuse to know more about multiple precision GMP .

#include <stdio.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_eigen.h>
int main (void)
        {
        int a=2; //matrix a x a
        int i,j;
        double cg[a][a];

        // Defining matrix
        cg[0][0] = 2.;
        cg[0][1] = 3.;
        cg[1][0] = -1.;
        cg[1][1] = 1.;

        gsl_matrix *TOTAL = gsl_matrix_alloc (a,a); //Alloc GSL matrix TOTAL

        //Setting TOTAL
        for (i = 0; i < a; i++)
                {
                for (j = 0; j < a; j++)
                        {
                        gsl_matrix_set (TOTAL, i, j, cg[i][j]);
                        }
                }


        gsl_vector_complex *eval = gsl_vector_complex_alloc (a); //Vector whose entries are the eigenvalues
        gsl_eigen_nonsymm_workspace *K = gsl_eigen_nonsymm_alloc (a); //alloc work space for GSL nonsymmetric calculations
        gsl_eigen_nonsymm (TOTAL,eval,K);
        gsl_eigen_nonsymm_free (K);
        gsl_matrix_free (TOTAL);


        for (i = 0; i < a; i++)  /* OUT OF RANGE ERROR */
                {
                gsl_complex eval_i = gsl_vector_complex_get (eval, i);
                printf ("%g + %g*i\n",GSL_REAL(eval_i),GSL_IMAG(eval_i) );
                }
return 0;
}