GLSL mat4 * vec4 multiplication different from C++ implementation (Not column vs row-major)

Question

I have defined the corners of a billboard in "parameter" space. To transform from parameter space to world space, I use a matrix called T:

bbCorners[i] = T*pBbCorners[i] + center; // billboard corner i transformed into world space

Now I render these corners as follows:

glBegin(GL_QUADS);
for (unsigned i = 0; i < 4; i++)
  glVertex3f(bbCorners[i].x, bbCorners[i].y, bbCorners[i].z);
glEnd();

And when my vertex shader looks like this, everything renders just fine (I make sure to upload the parameter space corners for the uniforms):

#version 120
#extension GL_EXT_gpu_shader4 : enable

// Parameter space billboard corners
uniform vec4 pBBTopLeft;
uniform vec4 pBBBottomLeft;
uniform vec4 pBBBottomRight;
uniform vec4 pBBTopRight;

// Parameter to world space matrix
uniform mat4 T;

// Center of ellipsoid
uniform vec4 center;

// Parameter space position of fragment on billboard
varying vec4 pPosBB;

void main() {   
  switch (gl_VertexID) {
  case 0:   pPosBB = pBBTopLeft;        break;
  case 1:   pPosBB = pBBBottomLeft;     break;
  case 2:   pPosBB = pBBBottomRight;    break;
  case 3:   pPosBB = pBBTopRight;       break;
  default:  pPosBB = vec4(vec3(0), 1);  break;
  }

  gl_Position = gl_ModelViewProjectionMatrix * gl_Vertex;
}

Now, the problem is the following. When I replace the gl_Vertex with the same calculation for vertex position, but in the shader, I get different results.

vec4 worldPos = transpose(T) * pPosBB + center; // Transpose because column major
gl_Position = gl_ModelViewProjectionMatrix * worldPos;

Does anyone know why the transformation might be different?

Edit: this is how T is calculated:

  // Specify orientation matrix by axes of ellipsoid
  float fO[16] = {
    axis1.x,  axis2.x,  axis3.x,  0.f,
    axis1.y,  axis2.y,  axis3.y,  0.f,
    axis1.z,  axis2.z,  axis3.z,  0.f,
    0.f,      0.f,      0.f,      1.f};
  qfeMatrix4f O;
  qfeMatrix4f::qfeSetMatrixElements(O, fO);

  // Specify radius matrix (lengths of axes)
  float fLambda[16] = {
    lambda1,  0.f,      0.f,      0.f,
    0.f,      lambda2,  0.f,      0.f,
    0.f,      0.f,      lambda3,  0.f,
    0.f,      0.f,      0.f,      1.f};
  qfeMatrix4f Lambda;
  qfeMatrix4f::qfeSetMatrixElements(Lambda, fLambda);

  // Calculate transformation matrix
  T = O*Lambda;

Edit: Solved it!

Turns out that T*pPosBB and center both have a w = 1. Adding them together leads to worldPos = (x, y, z, 2), which messes up our coordinates.