I would add that, in my opinion, if a compiler is going to do implicit initialisation, then zero is just about the worst value to chose. That’s because zero is usually a valid value, which can cause subtle errors to arise when you meant to store something in that variable prior to using it. Trust me - I’ve spent waaaaay too long tracking down errors like this, and I always compile my F90 programs with -stackvar (or its equivalent) for exactly this reason. A good compiler would automatically initialise all floating point variables to signalling NaN, to make sure that any nasty little errors due to assumed implicit initialisation blew the program up immediately. Of course this would be optional behaviour for bug checking purposes, since it would be bad for performance (albeit standard conforming - that’s what ‘undefined’ means).
Edit
I should explain the -stackvar trick for completeness. By default most Fortran compilers will put local variables on the heap and zero initialise, because many a ‘dusty deck’ assumes zero initialisation (not that this is guaranteed by the Fortran standard - like C, the values are undefined). Forcing stack allocation of locals inhibits the zero initialisation. Although not as good as NaN-initialisation, the chances are the the bit pattern ‘left’ on the stack is not a bunch of zeros, and hence obvious errors will appear sooner.
I 100% agree with that. Initializing to invalid values that can be trapped or tested for makes spotting problems much, much easier. I keep homemade 32bit and 64 bit versions of memset around at all times so I can explicitly set floating point values to bit patterns which resolve to Nan,-Inf or +Inf for debugging and exercising error checking code in numerical calculations.
Okay, Okay. It just very surprising to have variables used on a second run with the exact same variable values that were put in on the first run. It just strikes me as a bit odd. Nothing more. I accept that this can happen. I am too familiar with the old programming (FORTAN, PASCAL, LISP) languages where this would be really weird. Then again maybe I do not get around enough in current info science or programmer circle.
This is drifting off topic, but zero initialisation also qualifies as ‘really weird’ in Fortran. I don’t believe that the Fortran standard has ever required it (evidence: Fortran is only just getting rid of punch card compatibility, but does not require zero initialisation), but a fair number of ‘dusty decks’ do.
The problem with the matrix multiplication program is that the row and col equations only hold when the program is within the dimensions of the matrix.
int row = blockIdx.y + blockDim.y + threadIdx.y
int col = blockIdx.x + blockDim.x + threadIdx.x
When in a 10 by 10 product matrix the values of threadIdx.y or threadIdx.x get over 10 (on a c 0 to 9 format) the row and col values from the above equations are no longer accurate. So one must use these equations carefully and never get outside of the matrix dimensions. I think that this is the correct answer.
No. As I have now said three times: the code is only designed to work for matrices which are a round multiple of the block size. If even you can’t understand the code, at least look at the comments:
/ Thread block size
#define BLOCK_SIZE 16
// Forward declaration of the matrix multiplication kernel
__global__ void MatMulKernel(const Matrix, const Matrix, Matrix);
// Matrix Multiplication - Host code
// Matrix dimensions are assumed to be multiples of BLOCK_SIZE
void matrixMulGPU(const Matrix A, const Matrix B, const Matrix C)
ie. this code is designed to only work for matrices whose dimensions are round multiples of 16…
I remember that you said that now look at my output. I ran for two cases a 10 by 10 matrix and a 4 by 4 matrix, the last one fits your stated criteria.
blockIdx.y blockDim.y threadIdx.y blockIdx.x blockDim.x threadIdx.x C.elements Position
0 16 0 0 16 0 40.000000 0
blockIdx.y blockDim.y threadIdx.y blockIdx.x blockDim.x threadIdx.x C.elements Position
0 16 0 0 16 1 40.000000 1
blockIdx.y blockDim.y threadIdx.y blockIdx.x blockDim.x threadIdx.x C.elements Position
0 16 0 0 16 2 40.000000 2
blockIdx.y blockDim.y threadIdx.y blockIdx.x blockDim.x threadIdx.x C.elements Position
0 16 0 0 16 3 40.000000 3
blockIdx.y blockDim.y threadIdx.y blockIdx.x blockDim.x threadIdx.x C.elements Position
0 16 0 0 16 4 40.000000 4
blockIdx.y blockDim.y threadIdx.y blockIdx.x blockDim.x threadIdx.x C.elements Position
0 16 0 0 16 5 40.000000 5
blockIdx.y blockDim.y threadIdx.y blockIdx.x blockDim.x threadIdx.x C.elements Position
0 16 0 0 16 6 40.000000 6
blockIdx.y blockDim.y threadIdx.y blockIdx.x blockDim.x threadIdx.x C.elements Position
0 16 0 0 16 7 40.000000 7
blockIdx.y blockDim.y threadIdx.y blockIdx.x blockDim.x threadIdx.x C.elements Position
0 16 0 0 16 8 40.000000 8
blockIdx.y blockDim.y threadIdx.y blockIdx.x blockDim.x threadIdx.x C.elements Position
0 16 0 0 16 9 40.000000 9
blockIdx.y blockDim.y threadIdx.y blockIdx.x blockDim.x threadIdx.x C.elements Position
0 16 0 0 16 10 36.000000 10
blockIdx.y blockDim.y threadIdx.y blockIdx.x blockDim.x threadIdx.x C.elements Position
0 16 0 0 16 11 36.000000 11
You can see that the correct answer is to not go beyond the dimensions of the matrix. The alignment is a little off, but you can see that once threadIdx.x becomes 10, which is one value above the matrix dimension you start getting incorrect numbers - in this case a 36. So if you only calculate values within the matrix dimensions then you are getting the right answer for a 10 by 10 matrix. Get outside of that and you get wrong answers. Also, the basic row and col formulas do not apply then.
It is hard walking through this program using a gdb debug system. You must ask to print every variable that you are interested in. Thus I just printed out all of the variables that could be the slightest interest to me after each calculation of a c element. This I essentially got everything or most of it. I only copied part if it for you to see. It goes on for several pages.
Now in the case of the 4 by 4 matrix again with A and B matrices having all elements of 2 the output is:
matrix generation
A =
2.000000 2.000000 2.000000 2.000000
2.000000 2.000000 2.000000 2.000000
2.000000 2.000000 2.000000 2.000000
2.000000 2.000000 2.000000 2.000000
B =
2.000000 2.000000 2.000000 2.000000
2.000000 2.000000 2.000000 2.000000
2.000000 2.000000 2.000000 2.000000
2.000000 2.000000 2.000000 2.000000
Wed Dec 16 12:41:43 2009
GPU matrix multiplication
blockIdx.y blockDim.y threadIdx.y blockIdx.x blockDim.x threadIdx.x C.elements Position
0 16 0 0 16 0 16.000000 0
blockIdx.y blockDim.y threadIdx.y blockIdx.x blockDim.x threadIdx.x C.elements Position
0 16 0 0 16 1 16.000000 1
blockIdx.y blockDim.y threadIdx.y blockIdx.x blockDim.x threadIdx.x C.elements Position
0 16 0 0 16 2 16.000000 2
blockIdx.y blockDim.y threadIdx.y blockIdx.x blockDim.x threadIdx.x C.elements Position
0 16 0 0 16 3 16.000000 3
blockIdx.y blockDim.y threadIdx.y blockIdx.x blockDim.x threadIdx.x C.elements Position
0 16 0 0 16 4 12.000000 4
blockIdx.y blockDim.y threadIdx.y blockIdx.x blockDim.x threadIdx.x C.elements Position
0 16 0 0 16 5 12.000000 5
blockIdx.y blockDim.y threadIdx.y blockIdx.x blockDim.x threadIdx.x C.elements Position
0 16 0 0 16 6 12.000000 6
blockIdx.y blockDim.y threadIdx.y blockIdx.x blockDim.x threadIdx.x C.elements Position
0 16 0 0 16 7 12.000000 7
blockIdx.y blockDim.y threadIdx.y blockIdx.x blockDim.x threadIdx.x C.elements Position
0 16 0 0 16 8 8.000000 8
blockIdx.y blockDim.y threadIdx.y blockIdx.x blockDim.x threadIdx.x C.elements Position
0 16 0 0 16 9 8.000000 9
blockIdx.y blockDim.y threadIdx.y blockIdx.x blockDim.x threadIdx.x C.elements Position
0 16 0 0 16 10 8.000000 10
blockIdx.y blockDim.y threadIdx.y blockIdx.x blockDim.x threadIdx.x C.elements Position
0 16 0 0 16 11 8.000000 11
blockIdx.y blockDim.y threadIdx.y blockIdx.x blockDim.x threadIdx.x C.elements Position
0 16 0 0 16 12 4.000000 12
blockIdx.y blockDim.y threadIdx.y blockIdx.x blockDim.x threadIdx.x C.elements Position
0 16 0 0 16 13 4.000000 13
blockIdx.y blockDim.y threadIdx.y blockIdx.x blockDim.x threadIdx.x C.elements Position
0 16 0 0 16 14 4.000000 14
blockIdx.y blockDim.y threadIdx.y blockIdx.x blockDim.x threadIdx.x C.elements Position
0 16 0 0 16 15 4.000000 15
blockIdx.y blockDim.y threadIdx.y blockIdx.x blockDim.x threadIdx.x C.elements Position
0 16 1 0 16 0 16.000000 4
blockIdx.y blockDim.y threadIdx.y blockIdx.x blockDim.x threadIdx.x C.elements Position
0 16 1 0 16 1 16.000000 5
blockIdx.y blockDim.y threadIdx.y blockIdx.x blockDim.x threadIdx.x C.elements Position
0 16 1 0 16 2 16.000000 6
blockIdx.y blockDim.y threadIdx.y blockIdx.x blockDim.x threadIdx.x C.elements Position
0 16 1 0 16 3 16.000000 7
blockIdx.y blockDim.y threadIdx.y blockIdx.x blockDim.x threadIdx.x C.elements Position
0 16 1 0 16 4 12.000000 8
blockIdx.y blockDim.y threadIdx.y blockIdx.x blockDim.x threadIdx.x C.elements Position
0 16 1 0 16 5 12.000000 9
blockIdx.y blockDim.y threadIdx.y blockIdx.x blockDim.x threadIdx.x C.elements Position
0 16 1 0 16 6 12.000000 10
blockIdx.y blockDim.y threadIdx.y blockIdx.x blockDim.x threadIdx.x C.elements Position
0 16 1 0 16 7 12.000000 1
Now once that threadIdx.y or threadIdx.x gets above 3 (0 to 3 are the matrix dimensions) you start getting wrong answers. Thus one must stay within the dimensions. It also make sense: one will only do this to enable faster code execution and if speed is your goal then why do a whole bunch of unnecessary calculations- you are just wasting time. This last output with a square matrix of dimension four fits you criteria exactly. it is a multiple of the block size of 16.
Given that the stated criterion was “matrices whose dimensions are round multiples of 16” perhaps you could enlighten us as to the branch of arithmetic you are using?