Fast gemv()

Accelerated Computing CUDA CUDA Programming and Performance
21

Oh, i found a bug in the code… Change the kernel to this and the results should at least improve by 3-4x:

__global__ static void MY_GEMV(float* A, float* result)

{

	unsigned int t_x = threadIdx.x + blockIdx.x*32;

	

	float sum = 0.0f;

#pragma unroll

	for(int row = 0; row < N; row++)

	{

		sum += A[t_x + row*N]*d_x[row];

	}

	result[t_x] = sum;

}

Seems i was unnecessarily reading and writing to global… But including the transpose i guess your version will still be twice as fast, but what about for larger matrices? I think that is when the payback for transposing might come…

[attachment=16429:MY_GEMV.png]

22

It certainly is faster with the unnecessary global load and stores removed from the sgemv kernel, but it is still twice as slow as what I started with once the transpose is added in. At larger sizes the performance of the gemv kernel is very good, but the overall rate is no better, because the transpose seems to become very costly (more expensive by iteself that my sgemv kernel). These are at 9216x9216:

23

Yes seems it doesn’t scale well for increasing problem sizes.

Out of curiousity, how did you get the x-values to fit into constant memory for such large problem sizes? How much constant memory per SM do you have?

24

Isn’t it 64kb with 8kb cache working set on all current architectures? So a 9216 float vector (about 36kb) should fit fine.

25

My compiler started complaining as soon as i passed 8kb on this card, but those numbers seem familiar. Will check if there is any difference for 1.1 devices.

26

I don’t think there is.

Incidentally, it also just occurred to me that your code is calculating r=Ax, rather than r=Ax-b, like mine is, so there are some additional memory loads that would have to be included somehow as well.

27

Yes Ive been well aware of this, my only interest was in the matrix-vector multiplication… Actually thought maybe you had modified your kernel by removing this code for better comparision… Maybe ‘b’ could be loaded in a similar fashion as ‘x’.

It doesn’t matter though, you already have the fastest implementation!

28

Btw, if you’re running very large problems you will see 1-1.5% improvement if you change your code slightly:

__global__ void JORResidual(const float *A, const float *x, const float *b, float *r)

{

	__shared__ float buff[TPB];

	unsigned int rowIdx = __mul24(LDA_,blockIdx.x);

	// Compute parallel partial sums in shared memory

   float partial_sum = 0.0f;   

	for(int idx = threadIdx.x; idx < N_; idx += blockDim.x) {

		float Aval   = A[rowIdx+idx];

		float xval   = x[idx];

		float partial_sum  += Aval * xval;

		

	}

	

   buff[threadIdx.x] = partial_sum;

	__syncthreads();

shared memory read writes aren’t as fast as using registers instead ( ~20-30 ns). But if you’re running the smaller 1024x1024 this probably won’t be visible.

29

Thanks for looking further - there are a couple of little “1%” things that could be tweaked to wring the last few FLOPs out of it. The results of the original code I came up with really suprised me compared to CUBLAS, which I have always thought was reasonable. The code in the first post is over twice as fast as the transpose CUBLAS sgemv at every size I have tested it at (all in single precision GFLOP/s):

CUBLAS sgemvt_main()   New code

N=1024   10.15				  22.43

N=4096   11.52				  28.51

N=9126   14.29				  30.30
30

Yes, and i think your code should be general enough to be very useful wouldn’t it? Maybe nvidia could include it in their next CUBLAS update!

31

Hey avidday

I think I’ll stick to questions and hints about your jacobi.cu =)

You say that you need to do a dense matrix vector product but what you are actually with is sort of sparse. If you do not plan to change the matrix, you might very well gain alot from taking advantage of that (this relates to JORkernel as well).

You will definiately speed up your algorithm by transposing the matrix throughout all of your kernels: The thing is that every element in x is being accessed as many times as you have rows! By transposing you will be able to load a piece of the vector into shared memory and then treat the corresponding part of the matrix. This will reduce reads from the vector by a factor of 192… And as far as I can see you have the same ‘problem’ in JORkernel as well. Excess global memory access is bad!

Related to your JORkernel, I would highly recommend you to remove the branching in a loop. That’s really not good. Try and see if you can find a way to do it outside the loop instead.

And increase JORkernel executions by a factor 10 - memcpy is rediculously expensive!

32

Dense algebra performs considerably better in this application - the sparse solver goes to great lengths to reorder into dense blocks that can use dense blas for precisely this reason.

I don’t believe that is even close to true. The problem is transposition is very expensive - you replace N^2 LHS vector reads with N^2 LHS matrix reads to do the transpose, which is not any more efficient. Jimmy’s code showed that graphically. Using the accepted optimal transpose code, his version winds up being two to three times slower at large sizes than the non-transposed version.

Branching isn’t necessarily bad. Divergent branching is bad. For a vector length of 9216, only a single warp will diverge exactly once out of the 288 warp executions per block required to retire the loop. The rest don’t, and the overall peformance hit is very small - there is less than 1% difference in overall performance between JORkernel and JORresidual, with one containing a branch in the summation loop, and one not.

Actually it isn’t. The device->host memcpy for a 9216 long single precision vector into pageable memory takes 45 us (that is slightly less that 1GB/s), compared with a kernel run time of about 5500 us. So the copy only costs about 0.8% of additional kernel runs. You can see this in the attached profiler output.

I will grab your column major order kernel from the other thread and benchmark it at some point to see how it performs.

33

Just as a note, I modified my kernel last week to see if it would be faster to do the transpose “in place” in my MY_GEMV kernel since that would cut the load/stores from global memory in half.

The performance was of course better but still around 19% slower than your version…

[attachment=16560:my_mv.png] [attachment=16558:JOR.png]

I haven’t checked for partition camping though…

34

Yeah I was thinking about in-place transformation as well, but I still couldn’t come up with a scheme that looked like it would be better, especially at small problem sizes, where it becomes hard to get enough warps per multitprocessor and total blocks to keep throughput up. Especially on a GT200 with 30 MP.

As for partition camping, I think it should only be an issue in that execution configuration on cards with 8 memory banks (ie.GTX280/285 or C1060/S1070). What card did you run that on?

35

Yes I also think my implementation suffers severly on your GPU since it only has 32 blocks with N = 1024.

GTS250 which i believe has a 256 bit interface… yields 256/64 = 4 partitions ?

36

Yeah, so I did do a fix for partition camping which gave ~22 % performance increase. So it seems my version is actually faster now, at least on my GTS250.

[attachment=16561:my_mv.png]

Will have to check how this scales for larger problems as well… Good thing i checked for the partition camping though!

37

Cool. If you are willing to share the code, I can try running it on some “bigger iron” and see how it goes.

38

Sure! My suspiscion is that a solution with only 32 blocks isn’t going to scale well on a card with 30 SMs though, since it will be harder to hide the off-chip latencies.

It might need a rewrite/modification to work better on such a card. Maybe i will write something for that later today, anyways have a go at it:

[attachment=16566:MY_GEMV.cu]

Oh btw, check the “padding” parameter if it works for your memory interface… This is a bit crude but it was a quick and easy way to solve it :)

39

Here is a new version, it’s about 11 usec slower on my card but might actually scale better on yours leading to overall improvement. There are now 64 blocks for N=1024…

code:
[attachment=16568:MY_GEMV.cu]

profiler:
[attachment=16569:my_mv_2.png]

40

OK, it seems to scale well for the larger problems aswell, this is with the first version:

GFLOP/s

					your_JORResidual	 my_mv

N = 1024			   13.1				   14.3

N = 4096			   14.53				 18.05

N = 9216			   14.66				 20.77

As mentioned before it would be interesting to see how it would work for cards with more SMs.