Optimization for Real-to-Complex transforms

[Lagrang3]

You probably know that real to complex is very often used in applications of Fourier Methods. This particular use case allows for an optimization in computation and memory of a factor of 2. To exploit the complex conjugate symmetry here, the output is usually an array of floor(N/2)+1 complex numbers when the input consists of N real numbers (see for instance numpy).

The Real-to-Complex transform (and its inverse Complex-to-Real) require for an API that breaks our current design because of the following:

1. this use case applies exclusively to the realm of complex numbers, and not our generalized abstract DFT.
2. input and output have different types, one is Real and the other is Complex.

The first point above makes me think that at this point we should have separate classes to handle two types of transform philosophy: complex dft and ring dft. Moving forward with this proposal will break the nice symmetry that we previously have where ring transforms were handle as if they were complex, with the underlying static type selection making all the work to select the appropriate algorithms. On the bright side having these two api's will give the user more control of what is being done and we no longer need is_complex to select appropriate algorithms. Notice that ring transforms can also apply to complex numbers, with the appropriate selection of the root of unity, this is why I thought they might be unified at first.

If we introduce real-to-complex into the game, we could either expand the complex dft class to accommodate this case or we could introduce a third class specific to it. In any case I think we are likely bound to have template parameters specifying the Real and the Complex types, maybe having by default the following rule Real = Complex::value_type.

I am quite convinced that we should go for the 3 classes solution: ring, complex and real-to-complex dfts. @cosurgi @ckormanyos what do you think?

[ckormanyos]

convinced that we should go for the 3 classes solution: ring, complex and real-to-complex dfts

Hi @Lagrang3 my short answer/feeling is in complete agreement with your conviction. I think the three classes provide clarity of interface that is efficient, understandable and reduces obfuscation.

As we progress, I think we might find some more iinsight into this, but i feel that you are on the right track with that recognition and encourage you to draft out some interface(s) and see how the, let's say, play out.

[cosurgi]

Hi @Lagrang3 I noticed that you have pushed new stuff in real-to-complex branch. I will examine it carefully tomorow and let you know.

[cosurgi]

I notice that you have a little different interface in these tree files:

1. inside bsl_backend.hpp I see:

   • bsl_dft
   • bsl_algebraic_dft
   • bsl_transform
   • bsl_algebraic_transform

2. inside gsl_backend.hpp I see:

   • gsl_dft
   • gsl_rfft
   • gsl_transform

3. inside fftw_backend.hpp I see:

   • fftw_dft
   • fftw_rfft
   • fftw_transform

I understand that "algebraic" are more general, and impossible in other backends. I think that @ckormanyos mentioned somewhere that he has some nice implemented examples of real-to-complex, which are needed fo bsl_rfft, in the wide integer type.

Maybe something similar to this table will end up in documentation, I appreciate if you correct the mistakes. The (…) means that an object (the plan-like API) has to be present. The tr… is a short name for transform to fit in the table.

transform	dft_api.hpp	bsl_backend.hpp	gsl_backend.hpp	fftw_backend.hpp
	transform-like API
ℂ→ℂ fwd	tr…::fwd	bsl_tr…::fwd	gsl_tr…::fwd	fftw_tr…::fwd
ℂ→ℂ bkd	tr…::bkd	bsl_tr…::bkd	gsl_tr…::bkd	fftw_tr…::bkd
ℝ→½ℂ fwd			gsl_rfft::real_to_halfcplx	fftw_rfft::real_to_halfcplx
ℝ→½ℂ bkd			gsl_rfft::halfcplx_to_real	fftw_rfft::halfcplx_to_real
ℝ→ℂ fwd			gsl_rfft::real_to_cplx	fftw_rfft::real_to_cplx
ℝ→ℂ bkd			gsl_rfft::cplx_to_real	fftw_rfft::cplx_to_real
	plan-like API
ℂ→ℂ fwd	cplx_dft<…>(…).fwd	bsl_dft(…).fwd	gsl_dft(…).fwd	fftw_dft(…).fwd
ℂ→ℂ bkd	cplx_dft<…>(…).bkd	bsl_dft(…).bkd	gsl_dft(…).bkd	fftw_dft(…).bkd
ℝ→½ℂ fwd			gsl_rfft(...).real_to_halfcplx	fftw_rfft(…).real_to_halfcplx
ℝ→½ℂ bkd			gsl_rfft(...).halfcplx_to_real	fftw_rfft(…).halfcplx_to_real
ℝ→ℂ fwd			gsl_rfft(...).real_to_cplx	fftw_rfft(…).real_to_cplx
ℝ→ℂ bkd			gsl_rfft(...).cplx_to_real	fftw_rfft(…).cplx_to_real

Also this table is incomplete, because I definitely missed some parts of API when browsing the code.

Maybe this will help to clarify the API.

EDIT: I used following abbrevations so that the table columns would fit:

• forward → fwd
• backward → bkd
• complex → cplx

but in the code I definitely prefer longer names. Just fitting 4 columns was difficult here.

Also please feel free to edit this comment to correct it.

EDIT 2: added real to complex/halfcomplex routines supported by the GSL backend.

[Lagrang3]

Hi @cosurgi and @ckormanyos. While cleaning up the code, I noticed that the real_dft interface is complete on the real to halfcomplex side, but not so much if we want to support real to complex or complex to real, because of the selection of the complex type depending on a real value. These lines imply that the complex type has to be known at the time of instantiation. Yes, we have a complex selector (see issue #16), but I would like to hear your opinion. That approach won't work if the user comes up with a different complex type. Taking into account the deadline ahead I would suggest to do one of the following:

• either use our complex selector from the real type or
• do not provide real to complex nor complex to real or
• deduce the complex type from the iterator traits.

[cosurgi]

I think that using our complex selector boost::multiprecision::complex is presently the best option for a very simple reason: once we decide how to solve #16 we will modify boost::multiprecision::complex into something else. And this place in code which you mention will be affected (e.g. compilation errors) - and hence it will be impossible to forget about it.

In fact we could grep -E "\<std::complex\>" --color . -r --include='*fft*pp' to check if this happens in other places too.

EDIT: and grep -E "\<std::complex\>" --color ./include/boost/math/fft/ -r --include='*pp'

[cosurgi]

Hmm. But yes, if deduction from container is possible here, then perhaps it is an even better solution.