I'm trying figure out how to get an OpenGL specified object to be displayed correctly according to the device orientation (ie. according to the gravity vector from the accelerometer, and heading from compass).

The GLGravity sample project has an example which is almost like this (despite ignoring heading), but it has some glitches. For example, the teapot jumps 180deg as the device viewing angle crosses the horizon, and it also rotates spuriously if you tilt the device from portrait into landscape. This is fine for the context of this app, as it just shows off an object and it doesn't matter that it does these things. But it means that the code just doesn't work when you attempt to emulate real life viewing of an OpenGL object according to the device's orientation. What happens is that it almost works, but the heading rotation you apply from the compass gets "corrupted" by the spurious additional rotations seen in the GLGravity example project.

Can anyone provide sample code that shows how to adjust correctly for the device orientation (ie. gravity vector), or to fix the GLGravity example so that it doesn't include spurious heading changes?

//Clear matrix to be used to rotate from the current referential to one based on the gravity vector
bzero(matrix, sizeof(matrix));
matrix[3][3] = 1.0;

//Setup first matrix column as gravity vector
matrix[0][0] = accel[0] / length;
matrix[0][1] = accel[1] / length;
matrix[0][2] = accel[2] / length;

//Setup second matrix column as an arbitrary vector in the plane perpendicular to the gravity vector {Gx, Gy, Gz} defined by by the equation "Gx * x + Gy * y + Gz * z = 0" in which we arbitrarily set x=0 and y=1
matrix[1][0] = 0.0;
matrix[1][1] = 1.0;
matrix[1][2] = -accel[1] / accel[2];
length = sqrtf(matrix[1][0] * matrix[1][0] + matrix[1][1] * matrix[1][1] + matrix[1][2] * matrix[1][2]);
matrix[1][0] /= length;
matrix[1][1] /= length;
matrix[1][2] /= length;

//Setup third matrix column as the cross product of the first two
matrix[2][0] = matrix[0][1] * matrix[1][2] - matrix[0][2] * matrix[1][1];
matrix[2][1] = matrix[1][0] * matrix[0][2] - matrix[1][2] * matrix[0][0];
matrix[2][2] = matrix[0][0] * matrix[1][1] - matrix[0][1] * matrix[1][0];

//Finally load matrix
glMultMatrixf((GLfloat*)matrix);