Okay, I looked around and came up with another way to do it (see end of post if you’re interested), but I’d still like to know what’s going wrong with this method.

Your coordinate system in the first case was non-ortho-normal. I.e. the vector dir was not necessarily orthogonal to up and cross.

What you should have done:

Vector p; // your position
Vector front; // your movement direction
Vector up; // your world's up vector (0, 1, 0)
// construct ortho-normal coordinate system from movement and world-y
Vector z = normalize(front);
Vector x = normalize(cross(z, up));
Vector y = cross(x, z); // already unit length
// construct transformation matrix from base vectors
float tf[16] = {
x.x, x.y, x.z, 0,
y.x, y.y, y.z, 0,
z.x, z.y, z.z, 0,
p.x, p.y, p.z, 1
};
glMultMatrixf(tf);

EDIT: Keep in mind that “front” must not be collinear (parallel or anti-parallel) to “up”.