def convolve_spatial3(im, psfs,
                      mode="constant",
                      grid_dim=None,
                      sub_blocks=None,
                      pad_factor=2,
                      plan=None,
                      return_plan=False,
                      verbose=False):
    """
    GPU accelerated spatial varying convolution of an 3d image with a
    (Gz, Gy, Gx) grid of psfs assumed to be equally spaced within the image

    the input image im is subdivided into (Gz, Gy,Gx) blocks, each block is
    convolved with the corresponding psf and linearly interpolated to give the
    final result

    The psfs can be given either in

    A) Stackmode

    psfs.shape =  (Gz, Gy, Gx, Hz, Hy, Hx)
    then psfs[k,j,i] is the psf at the center of each block (i,j,k) in the image

    B) Flatmode

    psfs.shape = im.shape
    then the psfs are assumed to be definied on the gridpoints of the images itself
    in this case grid_dim = (Gz,Gy,Gx) has to be given

    as of now each image dimension has to be divisible by the grid dim, i.e.
    ::
        Nx % Gx == 0
        Ny % Gy == 0
        Nz % Gz == 0


    GPU Memory consumption is of order 8*Nx*Ny*Nz
    If not enough GPU memory is available, consider using
    sub_blocks = (n,m,l)
    then the operation is carried out in a tiled fashion reducing
    memory consumption to 8*Nx*Ny*(1/n+2/Gx)*(1/m+2/Gy)*(1/l+2/Gz)
    (so there is no much use if n>Gx/2...)
    Example
    -------


    im = np.zeros((64,64,64))
    im[::10,::10,::10] = 1.

    # Stackmode
    psfs = np.ones((8,8,8,4,4,4))
    res = convolve_spatial3(im, psfs, mode = "wrap")

    # Flatmode
    _Xs = np.meshgrid(*(np.arange(64),)*2)
    psfs = np.prod([np.clip(np.sin(2*np.pi*_X/8),0,1) for _X in _Xs],axis=0)
    res = convolve_spatial2(im, psfs, grid_dim = (16,16,16))


    Parameters
    ----------
    im: ndarray
        the image to convolve
    psfs: ndarray
        the (Gx,Gy) psf grid, either of shape (Gx,Gy, Hy, Hx) or im.shape

    mode: string, optional
        padding mode, either "constant" or "wrap"
    grid_dim: tuple, optional
        the (Gy,Gx) grid dimension, has to be provided if psfs.shape = im.shape

    sub_blocks: tuple, optional
        tiling mode, give e.g. (2,2) to sequentially operate on quadratnts
    pad_factor: int
        the factor of its size each block get tiled, use pad_factor=2 if the psfs
        are well localized, use pad_factor = 3 if not (e.g. if you experience blocking)_

    plan: fft_plan, optional
        when given use this as the fft plan
    return_plan: bool, optional
        return (res, plan) with plan being the fft plan for further use

    Returns
    -------
    res: ndarray
        the convolved image



    """

    ndim = im.ndim

    if ndim != 3:
        raise ValueError("wrong dimensions of input!")

    if grid_dim:
        if psfs.shape != im.shape:
            raise ValueError("if grid_dim is set, then im.shape = hs.shape !")
    else:
        if not psfs.ndim == 2 * ndim:
            raise ValueError("wrong dimensions of psf grid! should be (Gz,Gy,Gx,Nz,Ny,Nx)")

    if grid_dim:
        Gs = grid_dim
    else:
        Gs = psfs.shape[:ndim]

    if not np.all([n % g == 0 for n, g in zip(im.shape, Gs)]):
        raise NotImplementedError(
            "shape of image has to be divisible by Gx Gy  = %s shape mismatch" % (str(psfs.shape[:2])))

    if sub_blocks == None:
        return _convolve_spatial3(im,
                                  psfs,
                                  mode=mode,
                                  pad_factor=pad_factor,
                                  plan=plan,
                                  return_plan=return_plan,
                                  grid_dim=grid_dim)
    else:
        if not np.all([g % n == 0 for n, g in zip(sub_blocks, Gs)]):
            raise ValueError("psf grid dimension has to be divisible corresponding n_blocks")

        N_sub = [n // s for n, s in zip(im.shape, sub_blocks)]
        Nblocks = [n // g for n, g in zip(im.shape, Gs)]
        Npads = [n * (s > 1) for n, s in zip(Nblocks, sub_blocks)]
        grid_dim_sub = [g // s + 2 * (s > 1) for g, s in zip(Gs, sub_blocks)]

        if grid_dim:
            res = np.empty(im.shape, np.float32)
            plan = None
            for i, ((im_tile, im_s_src, im_s_dest), (hs_tile, hs_s_src, hs_s_dest)) \
                    in enumerate(zip(tile_iterator(im, blocksize=N_sub,
                                                   padsize=Npads,
                                                   mode=mode,
                                                   verbose=verbose), \
                                     tile_iterator(psfs, blocksize=N_sub,
                                                   padsize=Npads,
                                                   mode=mode,
                                                   verbose=verbose
                                                   ))):

                if verbose:
                    print("convolve_spatial3 ... %s\t/ %s" % (i + 1, np.prod(sub_blocks)))

                res_tile, plan = _convolve_spatial3(im_tile.copy(),
                                                    hs_tile.copy(),
                                                    mode=mode,
                                                    pad_factor=pad_factor,
                                                    return_plan=True,
                                                    plan=plan,
                                                    grid_dim=grid_dim_sub)

                res[im_s_src] = res_tile[im_s_dest]
            return res

        else:
            raise NotImplementedError("sub_blocks only implemented for Flatmode")