Optimizing array transposing function

Question

I'm working on a homework assignment, and I've been stuck for hours on my solution. The problem we've been given is to optimize the following code, so that it runs faster, regardless of how messy it becomes. We're supposed to use stuff like exploiting cache blocks and loop unrolling.

Problem:

//transpose a dim x dim matrix into dist by swapping all i,j with j,i
void transpose(int *dst, int *src, int dim) {
    int i, j;

    for(i = 0; i < dim; i++) {
        for(j = 0; j < dim; j++) {
                dst[j*dim + i] = src[i*dim + j];
        }
    }
}

What I have so far:

//attempt 1
void transpose(int *dst, int *src, int dim) {
    int i, j, id, jd;

    id = 0;
    for(i = 0; i < dim; i++, id+=dim) {
        jd = 0;
        for(j = 0; j < dim; j++, jd+=dim) {
                dst[jd + i] = src[id + j];
        }
    }
}

//attempt 2
void transpose(int *dst, int *src, int dim) {
    int i, j, id;
    int *pd, *ps;
    id = 0;
    for(i = 0; i < dim; i++, id+=dim) {
        pd = dst + i;
        ps = src + id;
        for(j = 0; j < dim; j++) {
                *pd = *ps++;
                pd += dim;
        }
    }
}

Some ideas, please correct me if I'm wrong:

I have thought about loop unrolling but I dont think that would help, because we don't know if the NxN matrix has prime dimensions or not. If I checked for that, it would include excess calculations which would just slow down the function.

Cache blocks wouldn't be very useful, because no matter what, we will be accessing one array linearly (1,2,3,4) while the other we will be accessing in jumps of N. While we can get the function to abuse the cache and access the src block faster, it will still take a long time to place those into the dst matrix.

I have also tried using pointers instead of array accessors, but I don't think that actually speeds up the program in any way.

Any help would be greatly appreciated.

Thanks

janneb · Accepted Answer · 2012-05-30 09:14:43Z

7

Cache blocking can be useful. For an example, lets say we have a cache line size of 64 bytes (which is what x86 uses these days). So for a large enough matrix such that it's larger than the cache size, then if we transpose a 16x16 block (since sizeof(int) == 4, thus 16 ints fit in a cache line, assuming the matrix is aligned on a cacheline bounday) we need to load 32 (16 from the source matrix, 16 from the destination matrix before we can dirty them) cache lines from memory and store another 16 lines (even though the stores are not sequential). In contrast, without cache blocking transposing the equivalent 16*16 elements requires us to load 16 cache lines from the source matrix, but 16*16=256 cache lines to be loaded and then stored for the destination matrix.

answered May 30, 2012 at 9:14

janneb

37.4k2 gold badges85 silver badges100 bronze badges

Sign up to request clarification or add additional context in comments.

4 Comments

wildplasser Over a year ago

This is the way to go. "cache oblivious matrix transposition" is the google phrase. Note: by taking 2*2 tiles of 16*16 cache lines you fill 4096 bytes, which is a memory page on (most) x86 machines.

sharptooth Over a year ago

Yes!!! Optimizing memory accesses can yield several times worth improvement from my experience.

Glen Takahashi Over a year ago

so how would I actually go about implementing cache blocking in this code? I was just thinking about if I was given say a 100x100 array, which I didn't think would benefit. Maybe I dont understand what you're saying completely.

janneb Over a year ago

@Glen: You have outer loops where you loop over the blocks, and then the inner loops which loop over every element in a block. So you'd have 4 nested loops. For a 100x100 array the benefit will likely not be that great; try it and measure!

ugoren · Accepted Answer · 2012-05-30 06:34:59Z

Unrolling is useful for large matrixes.
You'll need some code to deal with excess elements if the matrix size isn't a multiple of the times you unroll. But this will be outside the most critical loop, so for a large matrix it's worth it.

Regarding the direction of accesses - it may be better to read linearly and write in jumps of N, rather than vice versa. This is because read operations block the CPU, while write operations don't (up to a limit).

Other suggestions:
1. Can you use parallelization? OpenMP can help (though if you're expected to deliver single CPU performance, it's no good).
2. Disassemble the function and read it, focusing on the innermost loop. You may find things you wouldn't notice in C code.
3. Using decreasing counters (stopping at 0) might be slightly more efficient that increasing counters.
4. The compiler must assume that src and dst may alias (point to the same or overlapping memory), which limits its optimization options. If you could somehow tell the compiler that they can't overlap, it may be great help. However, I'm not sure how to do that (maybe use the restrict qualifier).

Peter G. · Accepted Answer · 2012-05-30 11:38:11Z

1

Messyness is not a problem, so: I would add a transposed flag to each matrix. This flag indicates, whether the stored data array of a matrix is to be interpreted in normal or transposed order.

All matrix operations should receive these new flags in addition to each matrix parameter. Inside each operation implement the code for all possible combinations of flags. Perhaps macros can save redundant writing here.

In this new implementation, the matrix transposition just toggles the flag: The space and time needed for the transpose operation is constant.

edited May 30, 2012 at 11:38

answered May 30, 2012 at 9:58

Peter G.

15.2k7 gold badges60 silver badges76 bronze badges

Comments

SKi · Accepted Answer · 2012-05-30 09:03:37Z

Just an idea how to implement unrolling:

void transpose(int *dst, int *src, int dim) {
    int i, j;
    const int dim1 = (dim / 4) * 4;

    for(i = 0; i < dim; i++) {
        for(j = 0; j < dim1; j+=4) {
                dst[j*dim + i]     = src[i*dim + j];
                dst[(j+1)*dim + i] = src[i*dim + (j+1)];
                dst[(j+2)*dim + i] = src[i*dim + (j+2)];
                dst[(j+3)*dim + i] = src[i*dim + (j+3)];
        }
        for( ; j < dim; j++) {
                dst[j*dim + i] = src[i*dim + j];
        }
        __builtin_prefetch (&src[(i+1)*dim], 0, 1);
    }
}

Of cource you should remove counting ( like i*dim) from the inner loop, as you already did in your attempts.

Cache prefetch could be used for source matrix.

fhtuft · Accepted Answer · 2012-05-30 09:52:22Z

0

you probably know this but register int (you tell the compiler that it would be smart to put this in register). And making the int's unsigned, may make things go little bit faster.

answered May 30, 2012 at 9:52

fhtuft

9765 silver badges8 bronze badges

2 Comments

Joel Falcou Over a year ago

register keyword doesn't really help there. The problem is cache/memory oriented adn micro optimizing register usage won't help.

Ryne Wang Over a year ago

int and unsigned int are the same in instruction level.

Collectives™ on Stack Overflow

Optimizing array transposing function

5 Answers 5

4 Comments

Comments

Comments

Comments

2 Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

5 Answers 5

4 Comments

Comments

Comments

Comments

2 Comments

Your Answer

Sign up or log in

Post as a guest

Related