c++multidimensional-arrayprocedureblitz++

How to pass Blitz++ subarray as an input/output parameter of a procedure


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.


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( ... );