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):
Exemplo simples do uso dos geradores de números aleatórios da biblioteca MKL,
#include <stdio.h> /////////////////////////////// MKL GSL Mersenne Twister #include <mkl_vsl.h> #define BRNG VSL_BRNG_MT19937 #define METHOD 0 //////////////////////////////// main () { int N = 4, sqrtN = 2 , i , j; int seed; seed = 345; float ranfloat[sqrtN]; //cada entrada corresponde a um numero aleatorio VSLStreamStatePtr stream; vslNewStream (&stream, BRNG, seed); //semeando o gerador for (i = 0; i < sqrtN; i++) { vsRngUniform (METHOD, stream, sqrtN, ranfloat, 0.0, 1.0); //carregando sqrtN numeros for (j = 0; j < sqrtN; j++) { printf ("%f\n", ranfloat[sqrtN]); } } return (0); }
Instalando a MKL
A Math Kernel Library MKL pode ser obtida em um pacote com compilador, debugger etc, ou separadamente. Você pode comprar o pacote (licenças acadêmicas têm desconto), testar ou obter uma licença para testes.
Eu optei pelo pacote completo.
Descompacte o arquivo, ex:
tar xvzf l_cproc_p_11.1.064.tgz
Instale com
cd l_cproc_p_11.1.064/
./install
O processo de instalação é muito simples, basta seguir as instruções.
Dentro do diretório de instalação temos a documentação, exemplo, bibliotecas etc. Sugiro um tour para conhecer o que a Intel oferece nesse pacote, provavelmente você vai encontrar muita coisa interessante.
Bom, aqui vão os passos (how-to) para fazer as coisas funcionarem:
- setar as variáveis de ambiente, o pacote traz um script que faz o serviço, basta executar o script iccvars.sh localizado dentro da pasta bin da instalação. No meu caso /home/maicon/opt. Precisamos indicar também a arquitetura desejada (intel64):
source ./iccvars.sh intel64
para evitar ter que fazer este procedimento para cada terminal aberto, vamos adicionar esse comando no bashrc
echo "source /home/maicon/opt/bin/iccvars.sh intel64" >> .bashrc
esse script inclui o compilador icc no PATH.
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
Simple benchmark, just floating point calculations.
This can be use to evaluate the performance on Monte Carlo application and numerical recursions without memory use.
We calculate the time for 10^9 m=tanh(ran2(&seed)*m), m=sin(ran2(&seed)*m) and if(ran2(&seed)). (seed=234)
Used cpus in decreasing performance order :
Intel(R) Core(TM) i5 CPU 750 @ 2.67GHz
AMD Phenom(tm) II X4 945 Processor @ 3GHz
Intel(R) Xeon(R) CPU X5550 @ 2.67GHz
AMD Athlon(tm) 64 X2 Dual Core Processor 4000+ @ 2.1GHz
Intel(R) Core(TM)2 Quad CPU Q6600 @ 2.40GHz - Chipset ICH10
Intel(R) Core(TM)2 Quad CPU Q6600 @ 2.40GHz - Chipset ICH7
Intel(R) Core(TM)2 Quad CPU Q9300 @ 2.50GHz - Chipset ICH9
Mobile AMD Sempron(tm) Processor 3500+ @1.8 GHz
Intel(R) Core 2(TM) CPU T7200@ 2.0Ghz
Below graph (incomplete) shows time versus cpu for the if(ran2(&seed)) routine
Also the complete spreadsheet
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;
}
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;
}