Convert Quaternion rotation to rotation matrix?

One way to do it, which is pretty easy to visualize, is to apply the rotation specified by your quaternion to the basis vectors (1,0,0), (0,1,0), and (0,0,1). The rotated values give the basis vectors in the rotated system relative to the original system. Use these vectors to form the rows of the rotation matrix. The resulting matrix, and its transpose, represent the forward and inverse transformations between the original system and the rotated system.

I'm not familiar with the conventions used by OpenGL, so maybe someone else can answer that part of your question...

The following code is based on a quaternion (qw, qx, qy, qz), where the order is based on the Boost quaternions:

boost::math::quaternion<float> quaternion;
float qw = quaternion.R_component_1();
float qx = quaternion.R_component_2();
float qy = quaternion.R_component_3();
float qz = quaternion.R_component_4();

First you have to normalize the quaternion:

const float n = 1.0f/sqrt(qx*qx+qy*qy+qz*qz+qw*qw);
qx *= n;
qy *= n;
qz *= n;
qw *= n;

Then you can create your matrix:

Matrix<float, 4>(
    1.0f - 2.0f*qy*qy - 2.0f*qz*qz, 2.0f*qx*qy - 2.0f*qz*qw, 2.0f*qx*qz + 2.0f*qy*qw, 0.0f,
    2.0f*qx*qy + 2.0f*qz*qw, 1.0f - 2.0f*qx*qx - 2.0f*qz*qz, 2.0f*qy*qz - 2.0f*qx*qw, 0.0f,
    2.0f*qx*qz - 2.0f*qy*qw, 2.0f*qy*qz + 2.0f*qx*qw, 1.0f - 2.0f*qx*qx - 2.0f*qy*qy, 0.0f,
    0.0f, 0.0f, 0.0f, 1.0f);

Depending on your matrix class, you might have to transpose it before passing it to OpenGL.