Code 3 is a function, written in C++, that performs the first task in the process. As mentioned earlier, elaborating a topological ordering of nodes with a depth-first recursive structure is much like actually computing an FFT with a depth-first recursive program (cf. Listing 3 in Appendix 2). Table 1 lists the nodes contained in the list `ns' after elaborating a size-8 transform by invoking elaborate(8, 0, 0, 0).

  CSplitRadix::elaborate(int N, int ioffset, int offset, int stride) {
    if(N > 4) {
      elaborate(N/2, ioffset, offset, stride+1);
      if(N/4 >= 4) {
        elaborate(N/4, ioffset+(1<<stride), offset+(N/2), stride+2);
        elaborate(N/4, ioffset-(1<<stride), offset+(3*N/4), stride+2);
      }else{
        CNodeLoad *n = new CNodeLoad(this, 4, ioffset, stride, 0);
        ns.push_back(assign_leaf_registers(n));
      }
      for(int k=0;k<N/4;k++) {
        CNodeBfly *n = new CNodeBfly(this, 4, k, stride);
        ns.push_back(assign_body_registers(n,k,N);
      }
    }else if(N==4) {
      CNodeLoad *n = new CNodeLoad(this, 4, ioffset, stride, 1);
      ns.push_back(assign_leaf_registers(n));
    }else if(N==2) {
      CNodeLoad *n = new CNodeLoad(this, 2, ioffset, stride, 1);
      ns.push_back(assign_leaf_registers(n));
    }
  }

A transform is divided into sub-transforms with recursive calls at lines 4, 6 and 7, until the base cases of size 2 or size 4 are reached at the leaves of the elaboration. As well as the size-2 and size-4 base cases, which are handled at lines 20-21 and 17-18 (respectively), there is a special case where two size-2 base cases are handled in parallel at lines 9-10. This special case of handling two size-2 base cases as a larger size-4 node ensures that larger transforms are composed of nodes that are homogeneous in size – this is of little utility when emitting VL=1VL=1 code, but it is exploited in "Other vector lengths" where the topological ordering of nodes is vectorized. The second row of Table 1 is just such a special case, since two size-2 leaf nodes are being computed, and thus the size is listed as 2(x2).

The elaborate function modifies the class member variable `ns' at lines 10, 14, 18 and 21, where it appends a new node to the back of the list. After the function returns, the ns list represents a topological ordering of the computation with CNodeLoad and CNodeBfly nodes. The nodes of type CNodeLoad represent leaf computations: these computations load elements from the input array and perform a small amount of leaf computation, leaving the result in a set of registers. The CNodeBfly nodes represent computations in the body of the transform: these use a twiddle factor to perform a butterfly computation on a vector of registers, leaving the result in the same registers.

The constructor for a CNodeLoad object computes input array addresses for the load operations using the input array offset (ioffset), the input array stride, the size of the node (the nodes instantiated at lines 9 and 17 are size-4, and the node instantiated at line 20 is size-2) and a final parameter that is non-zero if the node is a single node (the nodes instantiated at lines 17 and 20 are single nodes, while the node instantiated at line 9 is composed of two size-2 nodes).

As the newly instantiated CNodeLoad objects are appended to the back of ns at lines 10, 14 and 21, the assign_leaf_registers function assigns registers to the outputs of each instance. Registers are identified with integers beginning at zero, and when each register is created it is assigned an identifier from an auto-incrementing counter (RcounterRcounter). This function also maintains a map of registers to node pointers, referred to as rmap, where the node for a given register is the last node to reference that register.

The constructor for a CNodeBfly object uses k and stride to compute a twiddle factor for the new instance of a butterfly computation node. When the new instance of CNodeBfly is appended to the end of ns at line 14, the assign_body_registers function assigns registers RiRi to a node of size NnodeNnode with the following logic:

where i=0,⋯,Nnode-1i=0,⋯,Nnode-1 and RcounterRcounter is the auto-incrementing register counter. The assign_body_registers functions also updates the map of registers to node pointers by setting rmap[Ri]rmap[Ri] to point to the new instance of CNodeBfly.

  void sfft_dcf8_hc(sfft_plan_t *p, const void *vin, void *vout) {
    const SFFT_D *in = vin;
    SFFT_D *out = vout;
    SFFT_R r0,r1,r2,r3,r4,r5,r6,r7;
 
    L_4(in+0,in+8,in+4,in+12,&r0,&r1,&r2,&r3);
    L_2(in+2,in+10,in+14,in+6,&r4,&r5,&r6,&r7);
    K_0(&r0,&r2,&r4,&r6);
    S_4(r0,r2,r4,r6,out+0,out+4,out+8,out+12);
    K_N(VLIT2(0.7071,0.7071),VLIT2(0.7071,-0.7071),&r1,&r3,&r5,&r7);
    S_4(r1,r3,r5,r7,out+2,out+6,out+10,out+14);
  }

Given a list of nodes, it is a simple process to emit C code that can be compiled to actually compute the transform.

The example in Code 30 would be emitted from the list of four nodes in Table 1. Lines 1–4 are emitted from a function that generates a header, and line 13 is emitted from a function that generates a footer. Lines 6–11 are generated based on the list of nodes.

Code 30 contains references to several types, functions and macros that use upper-case identifiers – these are primitive functions or types that have been predefined as inline functions or macros. A benefit of using primitives in this way is that the details specific to numerical representation and the underlying machine have been abstracted away; thus, the same function can be compiled for a variety of types and machines by simply including a different header file with different primitives. Code 30, for example, could be compiled for double-precision arithmetic on an SSE2 machine by including sse_double.h, or it could be compiled with much slower scalar arithmetic by including scalar.h. The same code can even be used, without modification, to compute forward and backwards transforms, by using C preprocessor directives to conditionally alter the macros.

In order to accommodate mixed numerical representations, the signature of the outermost function references data with void pointers. In the case of the double-precision example in Code 30, SFFT_D would be defined to be double in the appropriate header file, and the void pointers are then cast to SFFT_D pointers.

The size-8 transform in Table 1 uses 8 registers, and thus a declaration of 8 registers of type SFFT_R has been emitted at line 4 in Code 30. In the case of double-precision arithmetic on a SSE2 machine, SFFT_R is defined as __m128d in sse_double.h.

The first two rows of Table 1 correspond to lines 6 and 7 of Code 30, respectively. The L_4 primitive is used to compute the size-4 leaf node in the first row of the table. The second row is a load/leaf node of size 2(x2), indicating two size-2 nodes in parallel, which is computed with the L_2 primitive. The input addresses in the table are the addresses of complex words, while the addresses in the generated code refer to the real and imaginary parts of a complex word, and thus the addresses from Table 1 are multiplied by a factor of 2 to obtain the addresses in Code 30.

The final two CNodeBfly nodes of Table 1 correspond to the K_0 and K_N sub-transform (a.k.a. butterfly) primitives at lines 8 and 10, respectively. Because the node in the third row of Table 1 has a twiddle factor of ω80ω80 (i.e., unity), the computation requires no multiplication, and the K_0 primitive is used for this special case. The K_N primitive at line 10 does require a twiddle factor, which is passed to K_N as two vector literals that represent the twiddle factor in unpacked form.1 Fast interleaved complex multiplication describes how interleaved complex multiplication is faster if one operand is pre-unpacked.

After each node is processed, the registers that have been used by it are checked in a map (rmap) that maps each register to the last node to have used that register. If the current node is the last node to have used a register, the register is stored to memory. In the case of the transform in Code 30, four registers are stored with an instance of the S_4 primitive at lines 9 and 11. In contrast to the load operations at the leaves, which are decimated-in-time and thus effectively pseudo-random memory accesses, the store operations are to linear regions of memory, the addresses of which can be determined from each register's integer identifier. The store address offset for data in register RiRi is simply i×2×VLi×2×VL.

Table 2 lists the nodes that represent a VL-1 size-16 transform. A VL of 2 implies that each vector register contains 2 complex words, and load operations on each of the 4 addresses in the first row of Table 2 will also load the complex words in the adjacent memory locations. Note that the complex words that would be incidentally loaded in the upper half of the VL-2 registers are the complex words that the third CNodeLoad object at row 5 would have loaded. This is exploited to load and compute the first and third CNodeLoad objects in parallel.

The second CNodeLoad object computes two size-2 leaf transforms in parallel, while the last CNodeLoad object computes a size-4 leaf transform. Because the size-4 transform is composed of two size-2 transforms, and memory addresses of the fourth CNodeLoad are adjacent (although permuted), some of the computation can be computed in parallel.

If the CNodeLoad objects at rows 1 and 5 are computed in parallel, the output will be four VL-2 registers: {{0,8}, {1,9}, {2,10}, {3,11}} – i.e., the first register contains what would have been register 0 in the lower half, and what would have been register 8 in the top half etc. Similarly, computing rows 2 and 6 in parallel would yield four VL-2 registers: {{4,14}, {5,15}, {6,12}, {7,13}} – note the permutation of the upper halves in this case. These registers are transposed to {{0,1}, {2,3}, {8,9}, {10,11}} and {{4,5}, {6,7}, {14,15}, {12,13}}, as in row 1 and 2 of Table 3.

With the transposed VL-2 registers, it is now possible to compute CNodeBfly nodes in parallel. For example, rows 2 and 3 of Table 2 can be computed in parallel on four VL-2 registers represented by {{0,1}, {2,3}, {4,5}, {6,7}}, as in row 3 of Table 3.

Code 82 is a C++ implementation of the vectorize_loads function. This function modifies a topological ordering of nodes (the class member variable ns) and uses two other functions: find_parallel_loads, which searches forward from the current node to find another CNodeLoad that shares adjacent memory addresses; and merge_loads(a,b), which adds the addresses, registers and type of b to a. Type introspection is used at lines 7 and 36 (and in other Listings), to differentiate between the two types of object.

Code 81 is a C++ implementation of the vectorize_ks function. For each CNodeBfly node, the function searches forward for another CNodeBfly that does not have a register dependence. Once found, the registers of the latter node are added to the former node, and the latter node erased. Finally, at line 19, the registers of the vectorized CNodeBfly node are merged using a perfect shuffle, which is then recursively applied on each half of the list. The effect is a merge that works for any power of 2 vector length.

  void CSplitRadix::vectorize_ks() {
    vector<CNodeHardCoded *>::iterator i;
    for(i=ns.begin(); i != ns.end();++i) {
      if(!(*i)->type().compare(``blockbfly'')) {
        vector<CNodeHardCoded *>::iterator j = i+1, pj = i;
        int count = 1;
        while(j != ns.end() && count < VL) {
          if(!(*j)->type().compare(``blockbfly'') && !register_dependence(*i, *j)) {
            (*i)->rs.insert( (*i)->rs.end(), (*j)->rs.begin(), (*j)->rs.end());
            ns.erase(j);
            count++;
            j = pj+1;
          }else {
            pj = j; ++j;
          }	
        }
        (*i)->merge_rs();
      }
    }
  }

  CNodeLoad *
  CSplitRadix::find_parallel_load(vector<CNodeHardCoded*>::iterator i){
    CNodeLoad *b = (CNodeLoad *)(*i);	
    for(int k=0;k<((N>2)?4:2);k++) {
      vector<CNodeHardCoded *>::iterator j = i+1;
      while(j != ns.end()) {
        if(!(*j)->type().compare(``blockload'')) {
          CNodeLoad *b2 = (CNodeLoad *)(*j);
          if(b2->iaddrs[k] > b->iaddrs[0] && b2->iaddrs[k] < b->iaddrs[0]+VL) {
            ns.erase(j);
            return b2;
          }
          ++j;
        }
      }
    }
    return NULL;
  }
  void CSplitRadix::merge_loads(CNodeLoad *b1, CNodeLoad *b2) {
    for(int i=0;i<b1->size;i++) {
      for(int j=0;j<b2->iaddrs.size();j++) {
        if(b2->iaddrs[j] > b1->iaddrs[i] && b2->iaddrs[j] < b1->iaddrs[i]+VL) {
          b1->iaddrs.push_back(b2->iaddrs[j]);
          b1->rs.push_back(b2->rs[j]);
          if(rmap[b2->rs[j]] == b2) rmap[b2->rs[j]] = b1;
        }
      }
    }
    b1->types.push_back(b2->types[0]);
  }
  void CSplitRadix::vectorize_loads() {
    vector<CNodeHardCoded *>::iterator i;
    for(i=ns.begin(); i != ns.end();++i) {
      if(!(*i)->type().compare(``blockload'')) {
        while(CNodeLoad *b2 = find_parallel_load(i))
          merge_loads((CNodeLoad *)(*i), b2);
      }
    }
  }

If vectorize_loads and vectorize_ks are invoked with VL=2VL=2 on the topological ordering of nodes in Table 2, the result is the vectorized node list shown in Table 3. As in "Emitting code", emitting code is a fairly simple process, and Code 89 is the code emitted from the node list in (Reference). There are only a few differences to note about the emitted code when VL>1VL>1.

void sfft_fcf16_hc(sfft_plan_t *p, const void *vin, void *vout) {
  const SFFT_D *in = vin;
  SFFT_D *out = vout;
  SFFT_R r0_1,r2_3,r4_5,r6_7,r8_9,r10_11,r12_13,r14_15;
 
  L_4_4(in+0,in+16,in+8,in+24,&r0_1,&r2_3,&r8_9,&r10_11);
  L_2_4(in+4,in+20,in+28,in+12,&r4_5,&r6_7,&r14_15,&r12_13);
  K_N(VLIT4(0.7071,0.7071,1,1),
      VLIT4(0.7071,-0.7071,0,-0),
      &r0_1,&r2_3,&r4_5,&r6_7);
  K_N(VLIT4(0.9239,0.9239,1,1),
      VLIT4(0.3827,-0.3827,0,-0),
      &r0_1,&r4_5,&r8_9,&r12_13);
  S_4(r0_1,r4_5,r8_9,r12_13,out+0,out+8,out+16,out+24);
  K_N(VLIT4(0.3827,0.3827,0.7071,0.7071),
      VLIT4(0.9239,-0.9239,0.7071,-0.7071),
      &r2_3,&r6_7,&r10_11,&r14_15);
  S_4(r2_3,r6_7,r10_11,r14_15,out+4,out+12,out+20,out+28);
}

So far, hard-coded transforms of vector length 1 and 2 have been presented. On Intel machines, VL-1 can be used to compute double-precision transforms with SSE2, while VL-2 can be used to compute double-precision transforms with AVX and single-precision transforms with SSE. The method of vectorization presented in this chapter scales above VL-2, and has been successfully used to compute VL-4 single-precision transforms with AVX.

The leaf primitives were coded by hand in all cases; VL-1 required L_2 and L_4, while VL-2 required L_2_2, L_2_4, L_4_2 and L_4_4. In the case of VL-4, not all permutations of possible leaf primitive were required – only 11 out of 16 were needed for the transforms that were generated.

It is an easy exercise to code the leaf primitives for VL≤4VL≤4 by hand, but for future machines that might feature vector lengths larger than 4, the leaf primitives could be automatically generated (in fact, "Other vector lengths" is concerned with automatic generation of leaf sub-transforms at another level of scale).

For a transform of size NN and leaf node size of SS (S=4S=4 in the examples in this chapter), the following constraint must be satisfied:

If this constraint is not satisfied, the size of either VL or SS must be reduced. In practice, VL and SS are small relative to the size of most transforms, and thus these corner cases typically only occur for very small sized transforms. Such an example is a size-2 transform when VL=2VL=2 and S=4S=4, where in this case the transform is too small to be computed with SIMD operations and should be computed with scalar arithmetic instead.

Code 96 is a VL-1 size-64 hard-coded four-step FFT. Before it can be used, an initialization procedure (not shown) allocates and populates the LUT at line 1 with the twiddle factors that are required for the step 2 multiplications. Line 44 shows the main function that executes the first sub-transform on the first column (line 49), and the second sub-transform on all remaining columns (line 50). Finally, the sub-transforms corresponding to step 4 of the four-step algorithm are executed on all columns in line 51.

The twiddle factor multiplication that corresponds to step 2 of the four-step algorithm takes place in lines 21-23 and lines 26-29. The first register is not multiplied with a twiddle factor because the first row of twiddle factors are ωN0ωN0 (i.e., unity). The other registers are multiplied with two registers loaded from the LUT, which are the unpacked real and imaginary parts (see Fast interleaved complex multiplication for details about unpacked complex multiplication).

Code 97 implements a VL-2 size-64 hard-coded four-step FFT. The main function (line 39) computes 8 FFTs along the columns for step 1 at line 44, and 8 FFTs along the columns for step 4 at line 45. There are only 4 iterations of the loop in each case because two sub-transforms are computed in parallel with each invocation of the sub-transform function.

In the function corresponding to the sub-transforms of step 1 (line 3), two store operations are latched (lines 23 and 24) before emitting code, which includes the preceding transposes (the TX2 operations) and twiddle factor multiplications (lines 13–22).

Code 103 is a size-64 hard-coded leaf transform with size-16 leaves. The first function (lines 1–17) is a size-16 leaf sub-transform, while the second (lines 18–32) consists of two size-8 leaf sub-transforms in parallel. The main function (lines 36–46) invokes four leaf sub-transforms (lines 40, 41, 43 and 44), and two loops of body sub-transforms (lines 42 and 45).

The first parameter to the leaf functions (see lines 1 and 18) is a pointer into an array of precomputed indices for the input data array. At lines 41 and 43–44, the array is incremented before subsequent calls to the leaf functions, and at line 39 the pointer is reset to the base of the array so that the transform can be used repeatedly.

The function used for the body sub-transforms (lines 33–35) is a wrapper for a primitive that computes a radix-2/4 butterfly. The last parameter to this function is a pointer to a precomputed LUT of twiddle factors for a sub-transform of size NN (the second parameter).

Code 104 is the main function for the FFT that corresponds to the leaf node list in Table 5. The first and third loops invoke size-16 sub-transforms at lines 8 and 16, and the second loop invokes 2x size-8 sub-transforms at line 12. Following the leaf sub-transforms, the body sub-transforms are called at lines 19-23.

void sfft_dcf128_shl16_4(sfft_plan_t *p,const void *vin,void *vout){
  const SFFT_D *in = vin;
  SFFT_D *out = vout;
  offset_t *is = p->is_base;
  offset_t *offsets = p->offsets_base;
  int i;
  for(i=3;i>0;--i) {
    sfft_dcf128_shl16_4_e(is, in, out+offsets[0]);
    is += 16; offsets += 1;
  }
  for(i=3;i>0;--i) {
    sfft_dcf128_shl16_4_o(is, in, out+offsets[0]);
    is += 16; offsets += 1;
  }
  for(i=2;i>0;--i) {
    sfft_dcf128_shl16_4_e(is, in, out+offsets[0]);
    is += 16; offsets += 1;
  }
  sfft_dcf128_shl16_4_X_4(out+0, 32, p->ws[0]);
  sfft_dcf128_shl16_4_X_4(out+0, 64, p->ws[1]);
  sfft_dcf128_shl16_4_X_4(out+128, 32, p->ws[0]);
  sfft_dcf128_shl16_4_X_4(out+192, 32, p->ws[0]);
  sfft_dcf128_shl16_4_X_4(out+0, 128, p->ws[2]);
}

In terms of code size, computing the leaf sub-transforms with three loops is economical. As the size of the transform grows, the code size attributed to the leaf sub-transforms remains constant. However, as the size of the transform begins to grow large (e.g., ≥65,536≥65,536), the instructions required for the body sub-transform calls (lines 19-23 in Code 104) begins to dominate the overall program size. "Optimizing the hierarchical structure" describes a method for compressing the code size of the body sub-transform calls while maintaining performance.

Because the input array references between consecutive leaves are now linear, and like types of leaf sub-transforms are grouped together, it is now possible to compute several leaf sub-transforms in parallel, which is fully described in "Other vector lengths".

A size-8192 hard-coded leaf FFT requires 229 calls to radix-2/4 and size-8 body sub-transforms. After optimizing the sequence of calls with Sequitur, the compact set of functions shown in Code 108 replaces a sequence of 229 calls.

Compared to the full list of statically elaborated calls, the optimized set of functions requires less code space while achieving better performance; and compared to a recursive program, the optimized set of function calls is faster (due to lower call and stack overhead) while trading off an acceptably small amount of code space.

The technique presented in this section has been verified for transforms ranging in size from 2626 through to 225225 (32 mega) points. The technique works well up until sizes of about 218218 points, but for larger transforms the elaboration and compile times begin to exceed 1 second or so, and the code size again begins to grow large. For transforms larger than 218218 points, a recursive program can be used until leaves of size 218218 points are reached, at which point the technique presented in this section is used.

void
sfft_fcf128_shl16_8_ee(offset_t *is,const SFFT_D *in,SFFT_D *out){
  SFFT_R r0,r1,r2,r3,r4,r5,r6,r7,r8,r9,r10,r11,r12,r13,r14,r15;
  L_4(in+is[0],in+is[1],in+is[2],in+is[3],&r0,&r1,&r2,&r3);
  L_2(in+is[4],in+is[5],in+is[6],in+is[7],&r4,&r5,&r6,&r7);
  K_0(&r0,&r2,&r4,&r6);
  K_N(VLIT4(0.7071,0.7071,0.7071,0.7071),
      VLIT4(0.7071,-0.7071,0.7071,-0.7071),
      &r1,&r3,&r5,&r7);
  L_4(in+is[8],in+is[9],in+is[10],in+is[11],&r8,&r9,&r10,&r11);
  L_4(in+is[12],in+is[13],in+is[14],in+is[15],&r12,&r13,&r14,&r15);
  K_0(&r0,&r4,&r8,&r12);
  K_N(VLIT4(0.9239,0.9239,0.9239,0.9239),
      VLIT4(0.3827,-0.3827,0.3827,-0.3827),
      &r1,&r5,&r9,&r13);
  TX2(&r0,&r1); TX2(&r4,&r5); TX2(&r8,&r9); TX2(&r12,&r13);
  S_4(r0,r4,r8,r12,out0+0,out0+8,out0+16,out0+24);
  S_4(r1,r5,r9,r13,out1+0,out1+8,out1+16,out1+24);
  K_N(VLIT4(0.7071,0.7071,0.7071,0.7071),
      VLIT4(0.7071,-0.7071,0.7071,-0.7071),
      &r2,&r6,&r10,&r14);
  K_N(VLIT4(0.3827,0.3827,0.3827,0.3827),
      VLIT4(0.9239,-0.9239,0.9239,-0.9239),
      &r3,&r7,&r11,&r15);
  TX2(&r2,&r3); TX2(&r6,&r7); TX2(&r10,&r11); TX2(&r14,&r15);
  S_4(r2,r6,r10,r14,out0+4,out0+12,out0+20,out0+28);
  S_4(r3,r7,r11,r15,out1+4,out1+12,out1+20,out1+28);
}

The vector leaf sub-transforms of a single size/type are handled in the same way as a VL-1 sub-transform, with one difference: the vector registers must be transposed before the data is stored to memory in the output array. In the example shown in Code 109, the transposes take place at lines 16 and 25.

Prior to the store operations, each position of the vector register (each position being a whole complex word) contains an element belonging to each of the leaf sub-transforms composing the vectorized sub-transform. Because each leaf sub-transform is stored sequentially to different locations in memory with aligned vector store operations, sets of registers are transposed such that each vector register contains elements from only one leaf sub-transform.

void
sfft_fcf128_shl16_8_eo(offset_t *is,const SFFT_D *in,SFFT_D *out){
  SFFT_R r0_1,r2_3,r4_5,r6_7,r8_9,r10_11,r12_13,r14_15,
	r16_17,r18_19,r20_21,r22_23,r24_25,r26_27,r28_29,r30_31;
  L_4_4(in+is[0],in+is[1],in+is[2],in+is[3],
      &r0_1,&r2_3,&r16_17,&r18_19);
  L_2_2(in+is[4],in+is[5],in+is[6],in+is[7],
      &r4_5,&r6_7,&r20_21,&r22_23);
  K_N(VLIT4(0.7071,0.7071,1,1),VLIT4(0.7071,-0.7071,0,-0),
      &r0_1,&r2_3,&r4_5,&r6_7);
  L_4_2(in+is[8],in+is[9],in+is[10],in+is[11],
        &r8_9,&r10_11,&r28_29,&r30_31);
  L_4_4(in+is[12],in+is[13],in+is[14],in+is[15],
        &r12_13,&r14_15,&r24_25,&r26_27);
  K_N(VLIT4(0.9239,0.9239,1,1),VLIT4(0.3827,-0.3827,0,-0),
      &r0_1,&r4_5,&r8_9,&r12_13);
  S_4(r0_1,r4_5,r8_9,r12_13,out0+0,out0+8,out0+16,out0+24);
  K_N(VLIT4(0.3827,0.3827,0.7071,0.7071),
      VLIT4(0.9239,-0.9239,0.7071,-0.7071),
      &r2_3,&r6_7,&r10_11,&r14_15);
  S_4(r2_3,r6_7,r10_11,r14_15,out0+4,out0+12,out0+20,out0+28);
  K_N(VLIT4(0.7071,0.7071,1,1),VLIT4(0.7071,-0.7071,0,-0),
      &r16_17,&r18_19,&r20_21,&r22_23);
  S_4(r16_17,r18_19,r20_21,r22_23,out1+0,out1+4,out1+8,out1+12);
  K_N(VLIT4(0.7071,0.7071,1,1),VLIT4(0.7071,-0.7071,0,-0),
      &r24_25,&r26_27,&r28_29,&r30_31);
  S_4(r24_25,r26_27,r28_29,r30_31,out1+16,out1+20,out1+24,out1+28);
}

In the case of a vector comprising heterogeneous leaf sub-transforms, the data is transposed into separate sub-transforms following the primitive leaf operations. The remainder of the computation is carried out separately for each leaf sub-transform in the vector, and no further transposes are required.

When elaborating and generating code for VL-2 transforms, there are only two heterogeneous leaf sub-transforms that might be required, but for other vector lengths the combinations are more complex. During the elaboration process, each unique combination that is encountered in the sorted list of leaf sub-transforms is elaborated into a function with repeated calls to the elaborate function, as was done in "Vector length 1" in order to elaborate a sub-transform composed of two size Nleaf/2Nleaf/2 sub-transforms.

Code 110 is an example of a heterogeneous size-16 VL-2 leaf sub-transform, where one size-16 leaf sub-transform is loaded into the lower halves of the vector registers, and the data from another leaf sub-transform composed of two size-8 sub-transforms is loaded into the upper halves. The primitive leaf operations at lines 5, 7, 11 and 13 transpose each sub-transform's data into separate vector registers, and the remainder of the computation is performed on each sub-transform separately. The size-16 sub-transform is stored to sequential locations in memory at lines 17 and 21, while the sub-transform composed of two size-8 leaf sub-transforms is stored to memory at lines 24 and 27.

Before an application can compute an FFT with SFFT, it must initialize a plan for the specific size, precision and direction of FFT. The library may have several FFTs and configurations that can compute the requested FFT, and it chooses the fastest option by timing each of the candidate configurations, which is at most 8 for any size of transform – a very small space compared to FFTW's exhaustive search of all possible FFT algorithms and configurations. Results and discussion describes an alternative to calibration, where machine learning is used with data collected from benchmarks to build a model that predicts performance.

After determining which implementation and parameters will be used, the initialization code allocates memory and populates any lookup tables that may be required. Before returning the plan to the application, a function pointer in the plan is updated to point to the FFT that has just been initialized.

Applications do not invoke any of the FFTs within SFFT directly. Rather they invoke a dispatch function on an initialized plan, which in turn transfers control to the correct FFT code within SFFT. The use of a dispatch function is purely a matter of convenience, so that users only need to deal with a few simple functions.

SFFT contains calibration code to measure the performance of the possible configurations of FFT on the target machine, which is at most 8 for each size of transform. Following calibration, the timing data is written to a file, which is then used by SFFT to select the fastest possible FFT for a given problem running on that machine.