I have a matrix stored in row-major order. I am trying to compute the DCT of a ub-matrix using FFTW, and I get nonsense. In the following paragraphs I will describe the problem and my solution, and I would like your help in understanding why it does not work.

Given some i and l, I would like to compute the DCT of a sub-matrix which consists of all rows k for which k mod l == i. For example, assume that l = 3 and i = 2. In the following matrix, the sub-matrix which I wish to transform is marked in red (2 mod 3 = 2, 5 mod 3 = 2, 8 mod 3 = 2).

The source and destination arrays have the same layout, and the transformed matrix should be stored in the same location in the destination array.

void transform(double* src, double* dest, size_t rows, size_t cols, size_t l, size_t i)
{
    int rank = 2;
    fftw_iodim64 dims[] = {
        { rows / l, l, l },
        { cols, rows, rows } };
    fftw_r2r_kind kind = FFTW_REDFT10;

    fftw_plan plan = fftw_plan_guru64_r2r(rank, dims, 0, NULL, src + l, dest + l, &kind, FFTW_ESTIMATE | FFTW_UNALIGNED | FFTW_PRESERVE_INPUT);
    fftw_execute(plan);
    fftw_destroy_plan(pla);
}

Update

I tested it for the simple case when i=l=1. Even in that case I get nonsense. I tested with a 3x4 matrix which is exactly one of the DCT basis vectors:

A(i,j) = cos((i + 0.5)*2*pi/3) * cos((j + 0.5)*3*pi/4)

I expected to get a result where all elements are (close to) zero except for one. But I get a resulting matrix that looks like this:

0             -2.22045e-016  1.33227e-015  2.22045e-016
2.22045e-016  -2.77556e-016  9.99201e-016  5.55112e-017
-8.88178e-016 -1.62359       7.83938       1.62359

Seems quite weird.

Update 2

I also tested with a simple matrix where the (0,0) element is 1 and the rest are zero. In this case, too, i=l=1 (the sub-matrix is the whole matrix). Here is the result:

      2       2       2       0
1.73205 1.73205 1.73205       0
      1       1       1       0