CMatrix3D Multiplication

[Zoomastigophora]

I'm confused by the matrix multiplication implementation; as it is now, A * B is mathematically B * A, which seems incredibly counterintuitive and prone to errors (at least to me). Was there a reason for it to be implemented as a right multiply?

[goblor]

The common definition of matrix multiplication is not an Abelian group.

=> A*B != B*A in general.

What you see here is a division ring (quaternion). It has all the properties of a field, except the commutation of the multiplication.

More in Wikipedia:

http://en.wikipedia.org/wiki/Division_ring

http://en.wikipedia.org/wiki/Quaternion

http://en.wikipedia.org/wiki/Matrix_multiplication

[Zoomastigophora]

Thanks for the reply, but I wasn't asking about the algebra of matrices. Maybe my phrasing was unclear. From what I've done in linear algebra and what I've seen with most math and computer graphics related work using matrix multiplication, it's defined as C = A x B where C[i,j] = sum(A[i, 1:n] x transpose(B[1:m, j])), the same definition Wikipedia gives. The implementation in CMatrix3D instead gives C[i,j] = sum(B[i, 1:n] x transpose(A[1:m, j])) and I was wondering if that was intentional because it's unconventional for matrix multiplication to be defined like that.

Also, what's a sound guy doing with a theoretical background in abstract mathematics? . I haven't thought of fields, rings, and groups since Category Theory and Number Theory in high school.

Edited June 19, 2010 by Zoomastigophora

[janwas]

hm, the matrix layout doesn't match the common (row, col) indexing - the elements and operator() are indexed by col, row instead. However, even accounting for that, you are right in that operator* seems to have it backwards. HOWEVER, it all seems to be rescued by operator*= again flipping the order (a *= b => a = b*a).

I don't often have contact with this code or this kind of linear algebra, but it does seem unnecessarily contorted.

Unfortunately, we probably won't be able to reconstruct the original author's intent or thinking. Do you think it'd be worthwhile to ensure the previous operator* is only ever called from operator*= and then flip it around?

As to Goblor pulling out his group theory, that is to be expected from (almost) mathematicians - he also happens to be interested in sound

[Zoomastigophora]

I might go ahead and redo the * and *= operators so that they're defined in the conventional manner and see what breaks. I'm probably going to be doing more work in the graphics code later and it'd be nice to not have to worry about multiplication weirdness.

[Ykkrosh]

Looks like we never use operator*= (except in a test where the order doesn't matter), so that should probably just be deleted (or else be implemented in terms of operator* so it's always going to match).

The internal storage of CMatrix3D has to be column-major to match OpenGL, i.e. (_11, _21, _31, _41) is the first column. So operator* is multiplying the first row of 'this' by the first column of 'matrix', to get the first output cell, which I think is the more normal direction of multiplication (and matches what Wikipedia shows). So I'm not seeing the problem here.

[Zoomastigophora]

From my reading, (_11, _12, _13, _14) would be the first column, unless that struct definition is setup so that columns read left to right instead of top to bottom, but this comment would seem to indicate otherwise:

// (Be aware that m(0,2) == _data2d[2][0] == _13, etc. This is to be considered a feature.)

In which case, * is multiplying the first column vector of A with the first row vector of B, when it should be the other way around.

[Ykkrosh]

I think the ordering of the field definitions in the source code is necessarily misleading, given that their memory layout must match what OpenGL expects. _13 is written on the 3rd source line but it corresponds to 3rd column 1st row in the column-major memory layout. _data2d[2] has to be the 3rd column since each column is a contiguous chunk of memory. operator()(col, row) seems to have its names the wrong way round - they're more like (row, col) now, so m(0, 2) is column 2 row 0.

It's definitely confusing . If I remember correctly, the first time I figured out what was going on I decided I didn't want to risk changing anything because I'd probably just cause more bugs and more confusion, so I added that comment claiming the confusingness to be a feature and not a bug to be fixed.

[Zoomastigophora]

Oh god, I think I know where the confusion is coming from. The matrix is actually constructed as row vectors, but since OpenGL uses column-major notation, all the operations are adapted to allow column-major convention. For instance, with column-major matrices, vector-matrix multiplication would have the matrix on the left of the vector, but since CMatrix3D is row-major, doing so would give the wrong result so the operator transposes the calculations (and also why *= flips the order around). I'm pretty sure this is the case since translation should be the 4th column if the matrix was column-major, but it's stored in the 4th row in CMatrix3D, which would be the case in row-major.

Edit: Wait, after reading more usages of CMatrix, I'm even more confused. In some cases, things seem to be stored column-major, others as row major. I don't even...

In any case, OpenGL doesn't seem to actually care which order the matrix is constructed in as long as its located in a contiguous block of memory. If the matrix was constructed using row vectors, you can use glLoadTransposeMatrix().

This will make writing that SSE asm matrix multiply even more confusing...

Edited June 20, 2010 by Zoomastigophora

[Ykkrosh]

Hmm, yeah, if the storage is column-major then the constructor arguments are row-major...

Translation is set up with

void CMatrix3D::SetTranslation (float x, float y, float z)
{
	_11=1.0f; _12=0.0f; _13=0.0f; _14=x;
	_21=0.0f; _22=1.0f; _23=0.0f; _24=y;
	_31=0.0f; _32=0.0f; _33=1.0f; _34=z;
	_41=0.0f; _42=0.0f; _43=0.0f; _44=1.0f;
}

which is the usual shape of translation matrices, where _n4 is the fourth column. In glLoadMatrix's column-major convention, that means _14=m[12], _24=m[13], _34=m[14], _44=m[15], which matches the class's storage definition:

struct {
	float _11, _21, _31, _41;
	float _12, _22, _32, _42;
	float _13, _23, _33, _43;
	float _14, _24, _34, _44;
};

where the fourth column (_n4) is stored last. So I think that's all okay.

The constructor uses the left-to-right-then-top-to-bottom matrix order (_11, _12, _13, _14, _21, ...) instead of the storage order. In operator*, that means the first 4 values it computes are the first row of the result matrix, and they're based on the first row of the current matrix (_1n) multiplied by columns of the argument matrix. So I think that's all okay too.

operator*= is wrong because n*=m should be equivalent to n=n*m (regardless of what operator* does), and currently it does n=m*n instead.

(What SSE asm matrix multiply? Why is that needed? I don't remember any profiling data accurately but I think skinning was the main thing that did lots of multiplications, and that maybe doesn't need a general-purpose 4x4 matrix multiply since it's only ever translations and rotations and could probably be done more efficiently in a special-case way.)

[Zoomastigophora]

I was actually looking at:

CVector3D CMatrix3D::GetTranslation() const
{
	CVector3D Temp;

	Temp.X = _14;
	Temp.Y = _24;
	Temp.Z = _34;

	return Temp;
}

and

//Applies a translation to the matrix
void CMatrix3D::Translate(float x, float y, float z)
{
	_14 += x;
	_24 += y;
	_34 += z;
}

and this in the header

// concatenate given translation onto this matrix. Assumes the current
// matrix is an affine transformation (i.e. the bottom row is [0,0,0,1])
// as an optimisation.
void Translate(float x, float y, float z);
void Translate(const CVector3D& vector);

SetTranslation is storing the translation as the 4th row, it just looks like the 4th column because the entries are laid out transposed of the struct. I think that CMatrix should be redone so that it does everything column-major consistently and if that breaks stuff, we'll fix it. I don't think it actually will break much since most usage of CMatrix assumes everything is handled column-major.

(The SSE asm matrix multiply was just something I wanted to do for my own personal amusement. It's how this whole brainhurt of convention inconsistency began. Once the renderer becomes more feature complete, it'll probably be a good idea to go through with a profiler and see what operations would benefit being done in optimized asm.)

Edited June 20, 2010 by Zoomastigophora

[Ykkrosh]

SetTranslation is storing the translation as the 4th row

It's storing the translation in (_14, _24, 34) which I will continue to claim is the fourth column (i.e. the fields are named _{row}{column}), which is consistent with the SetTranslation code and the operator* code and with glLoadMatrix/glUniformMatrix's interpretation of the storage, and with the Translate comment that says the bottom row of the matrix is always (0,0,0,1).

The (almost) only thing that is inconsistent is the "float _11, _21, _31, _41;" etc definition which is visually arranged like a row when actually it's defining a column. That could be made less misleading by putting it all on 16 lines or 1 line instead of the current transposed 4-line layout, or by adding a comment or something. (Also inconsistent is operator()(int col,int row) which should flip its naming of col and row.)

Does this interpretation of rows/columns still not work?

[Zoomastigophora]

Ok. I'll just need to remember that the struct is misleading then.

[janwas]

(The SSE asm matrix multiply was just something I wanted to do for my own personal amusement. It's how this whole brainhurt of convention inconsistency began. Once the renderer becomes more feature complete, it'll probably be a good idea to go through with a profiler and see what operations would benefit being done in optimized asm.)

Just a comment on this: rewriting single operations in asm often turns out not to be worthwhile, partially because VC is really stupid about integrating it into functions. Also, VC x64 doesn't allow inline asm at all :/ We do some external asm where there's no other way, but this probably isn't such a case. If you'd like to try your hand at vector stuff (just for fun, knowing it's probably not necessary nor worthwhile), I'd recommend using intrinsics such as _mm_mul_ps instead.