I have an FFT procedure fftconvx
taking two tensors Ttnsr
and S
as input parameters and producing a result into another tensor G
. All tensors are defined as Blitz++ arrays Array<complex<double>, N>
, where N
is the rank of the array. The procedure fftconvx
has to be called multiple times inside of a double loop.
Ideally I would like to pass subarrays Stnsr(ri,rj,rk,0)
or Stnsr(ri,rj,rk,1)
and receive the result into subarrays Gtnsr(t,p,ri,rj,rk,0)
or Gtnsr(t,p,ri,rj,rk,1)
as follows:
fftconvx( Gtnsr(t,p,ri,rj,rk,0), Ttnsr, Stnsr(ri,rj,rk,0) );
Variables ri,rj,rk
are Blitz++ array ranges. Unfortunately this does not work and results in the following compilation error:
error: invalid initialization of non-const reference of type
‘blitz::Array<std::complex<double>, 3>&’ from an rvalue of type
‘blitz::SliceInfo<std::complex<double>, int, int, blitz::Range, blitz::Range,
blitz::Range, int, blitz::nilArraySection, blitz::nilArraySection,
blitz::nilArraySection, blitz::nilArraySection, blitz::nilArraySection>::T_slice
{aka blitz::Array<std::complex<double>, 3>}’
fftconvx(Gtnsr(t,p,ri,rj,rk,0), Ttnsr, Stnsr(ri,rj,rk,0));
Signature of the fftconvx
is:
void fftconvx(Array<complex<double>, 3> &c,
Array<complex<double>, 3> x2,
Array<complex<double>, 3> x1,
...);
There more arrays and variables passed as input parameters, but I omit them for brevity.
So far I have come up with the solution based on temporary arrays S
and G
:
S(ri,rj,rk) = Stnsr(ri,rj,rk,0);
fftconvx(G, Ttnsr, S);
Gtnsr(t,p,ri,rj,rk,0) = G(ri,rj,rk);
I believe there is a more elegant solution.
Without knowing Blitz++ I offer this possible solution.
It looks like Gtnsr is a SliceInfo and not a Array, but that is has a operator Array.
So changing the fftconvx
into
template<class SliceOrArray>
void fftconvx(SliceOrArray &c,
const Array<complex<double>, 3> x2,
const Array<complex<double>, 3> x1,
...);
might work if the operations in fftconvx allows the use of the slice.
If Blitz++ is opdated to C++11 the following might work as well.
G fftconvx( const Array<complex<double>, 3> x2,
const Array<complex<double>, 3> x1,
...) {
G c;
...
return c; // C++11 NRVO
};
and then calling
Gtnsr(t,p,ri,rj,rk,0) = fftconvx( ... );