github.com/coyove/nj@v0.0.0-20221110084952-c7f8db1065c3/tests/bench/spectral-norm.txt (about)

     1  n = 0 
     2  nCPU = 1 -- os.numcpus
     3  
     4  function A(i, j)
     5      ij = i + j
     6      return ((ij)*(ij+1)>>1 + i + 1) 
     7  end
     8  
     9  function TimesVec(k, args)
    10      local v, s, n, u = args
    11      ul = #(u) 
    12      for i = int(s),int(n) do
    13          vi = 0 
    14          for j = 0,ul do
    15              vi = vi + u[j] / A(i, j) 
    16          end
    17          v[i] = vi 
    18      end
    19  end
    20  
    21  function TimesTranspVec(k, args)
    22      local v, s, n, u = args
    23      ul = #(u) 
    24      for i = int(s),int(n) do
    25          vi = 0 
    26          for j = 0,ul do
    27              vi = vi + u[j] / A(j, i) 
    28          end
    29          v[i] = vi 
    30      end
    31  end
    32  
    33  function ATimesTransp(v, u)
    34      x = array.make(#u) 
    35  
    36      payload = []
    37      for i = 0, nCPU do
    38          payload.append([x, i* #(v)/nCPU, (i+1)*#(v)/nCPU, u])
    39      end
    40      TimesVec.map(payload, nCPU)
    41  
    42      payload.clear()
    43      for i = 0, nCPU do
    44          payload.append([v, i*#(v)/nCPU, (i+1)*#(v)/nCPU, x])
    45      end
    46      TimesTranspVec.map(payload, nCPU)
    47  end
    48  
    49  n = 550
    50  u = array.make(n) 
    51  v = array.make(n) 
    52   
    53  for idx = 0, n do
    54      u[idx] = 1 
    55      v[idx] = 0 
    56  end
    57  
    58  for i = 0, 10 do
    59      ATimesTransp(v, u) 
    60      ATimesTransp(u, v) 
    61  end
    62  
    63  vBv = 0
    64  vv = 0 
    65  for i = 0, #(v) do
    66      vBv = vBv + u[i] * v[i]
    67      vv = vv + v[i] * v[i]
    68  end
    69  
    70  println(math.sqrt(vBv/vv)) 
    71