openpower/sv/remap_fft_yield.py

   1 # a "yield" version of the REMAP algorithm, for FFT Tukey-Cooley schedules
   2 # original code for the FFT Tukey-Cooley schedul:
   3 # https://www.nayuki.io/res/free-small-fft-in-multiple-languages/fft.py
   4 """
   5         # Radix-2 decimation-in-time FFT
   6         size = 2
   7         while size <= n:
   8                 halfsize = size // 2
   9                 tablestep = n // size
  10                 for i in range(0, n, size):
  11                         k = 0
  12                         for j in range(i, i + halfsize):
  13                 jh = j+halfsize
  14                 jl = j
  15                                 temp1 = vec[jh] * exptable[k]
  16                                 temp2 = vec[jl]
  17                                 vec[jh] = temp2 - temp1
  18                                 vec[jl] = temp2 + temp1
  19                                 k += tablestep
  20                 size *= 2
  21 """
  22
  23 # python "yield" can be iterated. use this to make it clear how
  24 # the indices are generated by using natural-looking nested loops
  25 def iterate_indices(SVSHAPE):
  26     # get indices to iterate over, in the required order
  27     n = SVSHAPE.lims[0]
  28     # createing lists of indices to iterate over in each dimension
  29     # has to be done dynamically, because it depends on the size
  30     # first, the size-based loop (which can be done statically)
  31     x_r = []
  32     size = 2
  33     while size <= n:
  34         x_r.append(size)
  35         size *= 2
  36     # invert order if requested
  37     if SVSHAPE.invxyz[0]: x_r.reverse()
  38
  39     if len(x_r) == 0:
  40         return
  41
  42     # start an infinite (wrapping) loop
  43     skip = 0
  44     while True:
  45         for size in x_r:           # loop over 3rd order dimension (size)
  46             # y_r schedule depends on size
  47             halfsize = size // 2
  48             tablestep = n // size
  49             y_r = []
  50             for i in range(0, n, size):
  51                 y_r.append(i)
  52             # invert if requested
  53             if SVSHAPE.invxyz[1]: y_r.reverse()
  54             for i in y_r:       # loop over 2nd order dimension
  55                 k_r = []
  56                 j_r = []
  57                 k = 0
  58                 for j in range(i, i+halfsize):
  59                     k_r.append(k)
  60                     j_r.append(j)
  61                     k += tablestep
  62                 # invert if requested
  63                 if SVSHAPE.invxyz[2]: k_r.reverse()
  64                 if SVSHAPE.invxyz[2]: j_r.reverse()
  65                 for j, k in zip(j_r, k_r):   # loop over 1st order dimension
  66                     # skip the first entries up to offset
  67                     if skip < SVSHAPE.offset:
  68                         skip += 1
  69                         continue
  70                     # now depending on MODE return the index
  71                     if SVSHAPE.skip == 0b00:
  72                         result = j              # for vec[j]
  73                     elif SVSHAPE.skip == 0b01:
  74                         result = j + halfsize   # for vec[j+halfsize]
  75                     elif SVSHAPE.skip == 0b10:
  76                         result = k              # for exptable[k]
  77
  78                     yield result
  79
  80 def demo():
  81     # set the dimension sizes here
  82     xdim = 8
  83     ydim = 0 # not needed
  84     zdim = 0 # again, not needed
  85
  86     # set total. err don't know how to calculate how many there are...
  87     # do it manually for now
  88     VL = 0
  89     size = 2
  90     n = xdim
  91     while size <= n:
  92         halfsize = size // 2
  93         tablestep = n // size
  94         for i in range(0, n, size):
  95             for j in range(i, i + halfsize):
  96                 VL += 1
  97         size *= 2
  98
  99     # set up an SVSHAPE
 100     class SVSHAPE:
 101         pass
 102     # j schedule
 103     SVSHAPE0 = SVSHAPE()
 104     SVSHAPE0.lims = [xdim, ydim, zdim]
 105     SVSHAPE0.order = [0,1,2]  # experiment with different permutations, here
 106     SVSHAPE0.mode = 0b01
 107     SVSHAPE0.skip = 0b00
 108     SVSHAPE0.offset = 0       # experiment with different offset, here
 109     SVSHAPE0.invxyz = [0,0,0] # inversion if desired
 110     # j+halfstep schedule
 111     SVSHAPE1 = SVSHAPE()
 112     SVSHAPE1.lims = [xdim, ydim, zdim]
 113     SVSHAPE1.order = [0,1,2]  # experiment with different permutations, here
 114     SVSHAPE1.mode = 0b01
 115     SVSHAPE1.skip = 0b01
 116     SVSHAPE1.offset = 0       # experiment with different offset, here
 117     SVSHAPE1.invxyz = [0,0,0] # inversion if desired
 118     # k schedule
 119     SVSHAPE2 = SVSHAPE()
 120     SVSHAPE2.lims = [xdim, ydim, zdim]
 121     SVSHAPE2.order = [0,1,2]  # experiment with different permutations, here
 122     SVSHAPE2.mode = 0b01
 123     SVSHAPE2.skip = 0b10
 124     SVSHAPE2.offset = 0       # experiment with different offset, here
 125     SVSHAPE2.invxyz = [0,0,0] # inversion if desired
 126
 127     # enumerate over the iterator function, getting new indices
 128     schedule = []
 129     for idx, (jl, jh, k) in enumerate(zip(iterate_indices(SVSHAPE0),
 130                                           iterate_indices(SVSHAPE1),
 131                                           iterate_indices(SVSHAPE2))):
 132         if idx >= VL:
 133             break
 134         schedule.append((jl, jh, k))
 135
 136     # ok now pretty-print the results, with some debug output
 137     size = 2
 138     idx = 0
 139     while size <= n:
 140         halfsize = size // 2
 141         tablestep = n // size
 142         print ("size %d halfsize %d tablestep %d" % \
 143                 (size, halfsize, tablestep))
 144         for i in range(0, n, size):
 145             prefix = "i %d\t" % i
 146             k = 0
 147             for j in range(i, i + halfsize):
 148                 jl, jh, ks = schedule[idx]
 149                 print ("  %-3d\t%s j=%-2d jh=%-2d k=%-2d -> "
 150                         "j[jl=%-2d] j[jh=%-2d] exptable[k=%d]" % \
 151                                 (idx, prefix, j, j+halfsize, k,
 152                                       jl, jh, ks))
 153                 k += tablestep
 154                 idx += 1
 155         size *= 2
 156
 157 # run the demo
 158 if __name__ == '__main__':
 159     demo()