examples/shootout/spectralnorm.pir

   1 #!./parrot
   2 # Copyright (C) 2006-2009, Parrot Foundation.
   3 # $Id$
   4 #
   5 # ./parrot -R jit spectralnorm.pir N         (N = 100 for shootout)
   6 # by Michal Jurosz
   7 # modified by Karl Forner to accept shootout default value of N=100
   8
   9 .sub eval_A
  10         .param int i
  11         .param int j
  12
  13         # return 1.0/((i+j)*(i+j+1)/2+i+1);
  14         $N0 = i + j
  15         $N1 = $N0 + 1
  16         $N0 *= $N1
  17         $N0 /= 2
  18         $N0 += i
  19         $N0 += 1
  20         $N0 = 1 / $N0
  21         .return ($N0)
  22 .end
  23
  24
  25 .sub eval_A_times_u
  26         .param int N
  27         .param pmc u
  28         .param pmc Au
  29
  30         .local int i, j
  31
  32         i = 0
  33 beginfor_i:
  34         unless i < N goto endfor_i
  35                 Au[i] = 0
  36                 j = 0
  37         beginfor_j:
  38                 unless j < N goto endfor_j
  39                         # Au[i]+=eval_A(i,j)*u[j]
  40                         $N0 = eval_A(i,j)
  41                         $N1 = u[j]
  42                         $N0 *= $N1
  43                         $N1 = Au[i]
  44                         $N0 += $N1
  45                         Au[i] = $N0
  46
  47                 inc j
  48                 goto beginfor_j
  49         endfor_j:
  50
  51         inc i
  52         goto beginfor_i
  53 endfor_i:
  54 .end
  55
  56
  57 .sub eval_At_times_u
  58         .param int N
  59         .param pmc u
  60         .param pmc Au
  61
  62         .local int i, j
  63
  64         i = 0
  65 beginfor_i:
  66         unless i < N goto endfor_i
  67                 Au[i] = 0
  68                 j = 0
  69         beginfor_j:
  70                 unless j < N goto endfor_j
  71                         # Au[i]+=eval_A(j,i)*u[j]
  72                         $N0 = eval_A(j,i)
  73                         $N1 = u[j]
  74                         $N0 *= $N1
  75                         $N1 = Au[i]
  76                         $N0 += $N1
  77                         Au[i] = $N0
  78
  79                 inc j
  80                 goto beginfor_j
  81         endfor_j:
  82
  83         inc i
  84         goto beginfor_i
  85 endfor_i:
  86 .end
  87
  88
  89 .sub eval_AtA_times_u
  90         .param int N
  91         .param pmc u
  92         .param pmc AtAu
  93
  94         .local pmc v
  95         v = new 'FixedFloatArray'
  96         v = N
  97
  98         eval_A_times_u(N,u,v)
  99         eval_At_times_u(N,v,AtAu)
 100 .end
 101
 102
 103 .sub main :main
 104         .param pmc argv
 105         .local int argc, N
 106
 107         N = 100
 108         argc = argv
 109         if argc == 1 goto default
 110         $S0 = argv[1]
 111         N = $S0
 112 default:
 113         .local pmc u, v
 114         u = new 'FixedFloatArray'
 115         u = N
 116         v = new 'FixedFloatArray'
 117         v = N
 118
 119         .local int i
 120
 121         i = 0
 122 beginfor_init:
 123         unless i < N goto endfor_init
 124                 u[i] = 1
 125         inc i
 126         goto beginfor_init
 127 endfor_init:
 128
 129
 130         i = 0
 131 beginfor_eval:
 132         unless i < 10 goto endfor_eval
 133             eval_AtA_times_u(N,u,v)
 134             eval_AtA_times_u(N,v,u)
 135         inc i
 136         goto beginfor_eval
 137 endfor_eval:
 138
 139         .local num vBv, vv
 140         vBv = 0.0
 141         vv = 0.0
 142
 143         i = 0
 144 beginfor_calc:
 145         unless i < N goto endfor_calc
 146                 # vBv+=u[i]*v[i]; vv+=v[i]*v[i];
 147                 $N0 = u[i]
 148                 $N1 = v[i]
 149                 $N0 *= $N1
 150                 vBv += $N0
 151                 $N0 = $N1
 152                 $N0 *= $N0
 153                 vv += $N0
 154         inc i
 155         goto beginfor_calc
 156 endfor_calc:
 157
 158         # print "%0.9f" % (sqrt(vBv/vv))
 159         $N0 = vBv / vv
 160         $N0 = sqrt $N0
 161         .local pmc spf
 162         spf = new 'FixedFloatArray'
 163         spf = 1
 164         spf[0] = $N0
 165         $S0 = sprintf "%.9f\n", spf
 166         print $S0
 167         exit 0
 168 .end
 169
 170 # Local Variables:
 171 #   mode: pir
 172 #   fill-column: 100
 173 # End:
 174 # vim: expandtab shiftwidth=4 ft=pir: