nsz Git - libfirm/blob - ir/be/test/langshootout/spectral-norm.c

   1 /* -*- mode: c -*-
   2  *
   3  * The Great Computer Language Shootout
   4  * http://shootout.alioth.debian.org/
   5  *
   6  * Contributed by Sebastien Loisel
   7  */
   8
   9 #include <stdio.h>
  10 #include <stdlib.h>
  11 #include <math.h>
  12 /* Includes for alloca() */
  13 #if defined(__FreeBSD__)
  14 #include <stdlib.h>
  15 #elif defined(_WIN32)
  16 #include <malloc.h>
  17 #else
  18 #include <alloca.h>
  19 #endif
  20
  21 double eval_A(int i, int j) { return 1.0/((i+j)*(i+j+1)/2+i+1); }
  22
  23 void eval_A_times_u(int N, const double u[], double Au[])
  24 {
  25         int i,j;
  26         for(i=0;i<N;i++)
  27     {
  28       Au[i]=0;
  29       for(j=0;j<N;j++) Au[i]+=eval_A(i,j)*u[j];
  30     }
  31 }
  32
  33 void eval_At_times_u(int N, const double u[], double Au[])
  34 {
  35         int i,j;
  36         for(i=0;i<N;i++) {
  37                 Au[i]=0;
  38                 for(j=0;j<N;j++)
  39                         Au[i] += eval_A(j,i)*u[j];
  40         }
  41 }
  42
  43 void eval_AtA_times_u(int N, const double u[], double AtAu[])
  44 {
  45         double v[N];
  46         eval_A_times_u(N,u,v);
  47         eval_At_times_u(N,v,AtAu);
  48 }
  49
  50 int main(int argc, char *argv[])
  51 {
  52         int i;
  53         int N = ((argc == 2) ? atoi(argv[1]) : 2000);
  54
  55         double *u = alloca(sizeof(u[0]) * N);
  56         double *v = alloca(sizeof(v[0]) * N);
  57         double vBv, vv;
  58
  59         for(i=0;i<N;i++)
  60                 u[i]=1;
  61
  62         for(i=0;i<10;i++) {
  63                 eval_AtA_times_u(N,u,v);
  64                 eval_AtA_times_u(N,v,u);
  65         }
  66
  67         vBv = vv = 0;
  68
  69         for(i=0;i<N;i++) {
  70                 vBv+=u[i]*v[i];
  71                 vv+=v[i]*v[i];
  72         }
  73
  74         printf("%0.9f\n", sqrt(vBv/vv));
  75         return 0;
  76 }