There is a section in Rabiener and Gold’s 1975 DSP book listed near the index term “twiddle factors” around page 347 that shows how to decompose a DFT in 2 dimensions. I used to use it a lot prior to FFTW to do oddball DFT sizes like 384 without zero padding.
Theory and application of digital signal processing
Prentice-Hall signal processing series
Lawrence R. Rabiner, Bernard Gold, Prentice-Hall, 1975
First you need to factor the desired DFT size into 2 numbers M and N and then allocate a M by N matrix. Write you data to the matrix in row order skipping to the next row as you fill. This is an in place algorithm.
Do M point DFTs on each column. Next Multiply each element in the matrix by its twiddle factor $exp( j 2 \pi (m-1)(n-1))$ where m and m are the index terms given the Fortran convention. Do N point DFTs on each row. Read out the results in column order.
X=reshape(x,N,M).'; % read in as rows
X=fft(X) % fft on each column
X=X.*Twiddle;% element by element product
X=fft(X.').' ; %fft on each row
y=reshape(X,N*M,1); % read out as columns
There is also another order mentioned in Rabiener’s book. Twiddle first, row fft, column fft. He also makes it a point to note that each row or column dft can be implemented, in turn as a composite square as long as N is not a prime number. The row in, column out indexing is also interesting because it generalizes bit reverse addressing.
In summary, you can do a 100k DFT with 4k DFTs if you have enough memory.