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
  13 double eval_A(int i, int j) { return 1.0/((i+j)*(i+j+1)/2+i+1); }
  14
  15 void eval_A_times_u(int N, const double u[], double Au[])
  16 {
  17         int i,j;
  18         for(i=0;i<N;i++)
  19     {
  20       Au[i]=0;
  21       for(j=0;j<N;j++) Au[i]+=eval_A(i,j)*u[j];
  22     }
  23 }
  24
  25 void eval_At_times_u(int N, const double u[], double Au[])
  26 {
  27         int i,j;
  28         for(i=0;i<N;i++) {
  29                 Au[i]=0;
  30                 for(j=0;j<N;j++)
  31                         Au[i] += eval_A(j,i)*u[j];
  32         }
  33 }
  34
  35 void eval_AtA_times_u(int N, const double u[], double AtAu[])
  36 {
  37         double *v[N];
  38         eval_A_times_u(N,u,v);
  39         eval_At_times_u(N,v,AtAu);
  40 }
  41
  42 int main(int argc, char *argv[])
  43 {
  44         int i;
  45         int N = ((argc == 2) ? atoi(argv[1]) : 2000);
  46
  47         double *u = alloca(sizeof(u[0]) * N);
  48         double *v = alloca(sizeof(v[0]) * N);
  49         double vBv, vv;
  50
  51         for(i=0;i<N;i++)
  52                 u[i]=1;
  53
  54         for(i=0;i<10;i++) {
  55                 eval_AtA_times_u(N,u,v);
  56                 eval_AtA_times_u(N,v,u);
  57         }
  58
  59         vBv = vv = 0;
  60
  61         for(i=0;i<N;i++) {
  62                 vBv+=u[i]*v[i];
  63                 vv+=v[i]*v[i];
  64         }
  65
  66         printf("%0.9f\n", sqrt(vBv/vv));
  67         return 0;
  68 }