vec3 U    = normalize ( u );
vec3 dirZ = normalize ( c - p );                 // нормализация обязательна!
vec3 dirX = normalize ( cross ( dirZ, U ) );     // cross - векторное произведение (на выходе - вектор перпендикулярный обоим "множителям")
vec3 dirY = normalize ( cross ( dirX, dirZ ) );

// так как используемые в огл матрицы трансформации столбец-ориентированные,
// то при заполнении такой матрицы необходимо строки менять со столбцами (транспонировать),
// если ты используешь строко-ориентированную матрицу, и оставлять как есть, - если нет:

mat4 camView = mat4 ( dirX.x, dirY.x, -dirZ.x, 0.0,  // на самом деле для огл это - столбец, а не строка (то есть ты на место строки кладешь столбец) !
                      dirX.y, dirY.y, -dirZ.y, 0.0,
                      dirX.z, dirY.z, -dirZ.z, 0.0,
                      P.x   , P.y   ,   P.z   , 1.0 );

// здесь: Р = 0 или P = LinearPart(camView) * (-p), где LinearPart - mat3 верхней левой части.
// таким образом, ты можешь использовать либо:

glMatrixMode  ( GL_MODELVIEW );
glLoadMatrixf ( camView );
glTranslatef  ( -p );

// либо (для второго случая):

glMatrixMode  ( GL_MODELVIEW );
glLoadMatrixf ( camView );

// а например в вершинном шейдере если у тебя вдруг происходит расчет матрицы вида,
// то "LinearPart(camView) * (-p)"  (то есть mat3 * vec3) гммм... затратно, лучше сделать так:

gl_Position = yourProjectionMatrix *  yourModelViewMatrix * ( gl_Vertex - vec4 ( p, 0.0 ) );