Shell x Listas x Monte Carlo

Muitas vezes é interessante guardas listas de sítios e sortear somente dentro desta lista para evitar gerar números aleatórios que não serão usados. É fácil usar listas em uma dimensão: quando um sítio sai da lista ele pode ser substituído pelo ultimo da fila e esta perde um elemento.

Tudo muito bom, tudo muito legal, mas há um preço a pagar se estamos interessados no comportamento dinâmico... um não, dois !

Gerando números aleatórios com a biblioteca VSL (MKL)

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);

}

Para compilar basta linkar as bibliotecas adequadas:

ia32,  icc exemplo.c  -lmkl_intel -lmkl_sequential -lmkl_core  -o exemplo.bin

intel64, icc exemplo.c -lmkl_intel_lp64 -lmkl_sequential -lmkl_core  -o exemplo.bin

Pontos importante:
os números são obtidos em stream, a idéia é carrega um vetor a cada chamada da função vsRngUniform, no exemplo carregamos sqrtN números entre 0 e 1 (float) distribuídos uniformente. O desempenho é melhor para sqrtN > 1000.

Outras funções: 
vdRngUniform  -  double
viRngUniform  - int

Outros geradores, funções e mais detalhes podem ser obtidos na documentação da biblioteca Vector Statistical Library VSL

Benchmark - cpus

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

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;
}