First of all: very helpful topic title indeed.

Secondly, if you don't want your program to drag all that GLU junk with it just for glulookat, etc.. you may want to consider this pascal code as alternative, translated from SGI'S opengl implementation sample.

Code:
// this is gluPickMatrix equivalent translated from ogl-sample.

procedure PickMatrix(x, y, deltax, deltay: GLfloat; viewport: array of integer);
begin
  if &#40;deltax <= 0&#41; and &#40;deltay <= 0&#41; then
    exit;

  &#40;* Translate and scale the picked region to the entire window *&#41;
  glTranslatef&#40;
    &#40;viewport&#91;2&#93; - 2 * &#40;x - viewport&#91;0&#93;&#41;&#41; / deltax,
    &#40;viewport&#91;3&#93; - 2 * &#40;y - viewport&#91;1&#93;&#41;&#41; / deltay,
    0&#41;;

  glScalef&#40;viewport&#91;2&#93; / deltax, viewport&#91;3&#93; / deltay, 1.0&#41;;
end;

// this is gluPerspective equivalent translated from ogl-sample.
// matrix should be double but i use single precision

procedure Perspective&#40;fovy, aspect, zNear, zFar&#58; single&#41;;
var
  m&#58; TMatrix4f;
  sine, cotangent, deltaZ, radians&#58; single;

begin
  radians &#58;= fovy / 2 * Pi / 180;

  deltaZ &#58;= zFar - zNear;
  sine   &#58;= sin&#40;radians&#41;;
  if &#40;&#40;deltaZ = 0&#41; and &#40;sine = 0&#41; and &#40;aspect = 0&#41;&#41; then
    exit;

  cotangent &#58;= COS&#40;radians&#41; / sine;

  m &#58;= Identity;
  m&#91;0&#93;&#91;0&#93; &#58;= cotangent / aspect;
  m&#91;1&#93;&#91;1&#93; &#58;= cotangent;
  m&#91;2&#93;&#91;2&#93; &#58;= -&#40;zFar + zNear&#41; / deltaZ;
  m&#91;2&#93;&#91;3&#93; &#58;= -1;
  m&#91;3&#93;&#91;2&#93; &#58;= -2 * zNear * zFar / deltaZ;
  m&#91;3&#93;&#91;3&#93; &#58;= 0;
  glMultMatrixf&#40;@m&#91;0&#93;&#91;0&#93;&#41;;
end;

// this is glulookat equivalent translated from ogl-sample.
// matrix should be single but i use single precision

procedure normalize&#40;var v&#58; array of single&#41;;
var
  r&#58; single;
begin
  r &#58;= sqrt&#40;v&#91;0&#93; * v&#91;0&#93; + v&#91;1&#93; * v&#91;1&#93; + v&#91;2&#93; * v&#91;2&#93;&#41;;

  if &#40;r = 0&#41; then
    exit;

  v&#91;0&#93; &#58;= v&#91;0&#93; / r;
  v&#91;1&#93; &#58;= v&#91;1&#93; / r;
  v&#91;2&#93; &#58;= v&#91;2&#93; / r;
end;

procedure crossp&#40;var v1, v2, Result&#58; array of single&#41;;
begin
  Result&#91;0&#93; &#58;= v1&#91;1&#93; * v2&#91;2&#93; - v1&#91;2&#93; * v2&#91;1&#93;;
  Result&#91;1&#93; &#58;= v1&#91;2&#93; * v2&#91;0&#93; - v1&#91;0&#93; * v2&#91;2&#93;;
  Result&#91;2&#93; &#58;= v1&#91;0&#93; * v2&#91;1&#93; - v1&#91;1&#93; * v2&#91;0&#93;;
end;

// translated from ogl-sample
procedure LookAt&#40;eyex, eyey, eyez, centerx, centery, centerz, upx, upy, upz&#58; single&#41;;
var
  forw, side, up&#58; array&#91;0..2&#93; of single;
  m&#58; TMatrix4f;
begin
  forw&#91;0&#93; &#58;= centerx - eyex;
  forw&#91;1&#93; &#58;= centery - eyey;
  forw&#91;2&#93; &#58;= centerz - eyez;

  up&#91;0&#93; &#58;= upx;
  up&#91;1&#93; &#58;= upy;
  up&#91;2&#93; &#58;= upz;

  normalize&#40;forw&#41;;

  &#40;* Side&#58;= forw x up *&#41;
  crossp&#40;forw, up, side&#41;;
  normalize&#40;side&#41;;

  &#40;* Recompute up as&#58; up&#58;= side x forw *&#41;
  crossp&#40;side, forw, up&#41;;

  m &#58;= Identity;
  m&#91;0&#93;&#91;0&#93; &#58;= side&#91;0&#93;;
  m&#91;0&#93;&#91;1&#93; &#58;= side&#91;1&#93;;
  m&#91;0&#93;&#91;2&#93; &#58;= side&#91;2&#93;;

  m&#91;1&#93;&#91;0&#93; &#58;= up&#91;0&#93;;
  m&#91;1&#93;&#91;1&#93; &#58;= up&#91;1&#93;;
  m&#91;1&#93;&#91;2&#93; &#58;= up&#91;2&#93;;

  m&#91;2&#93;&#91;0&#93; &#58;= -forw&#91;0&#93;;
  m&#91;2&#93;&#91;1&#93; &#58;= -forw&#91;1&#93;;
  m&#91;2&#93;&#91;2&#93; &#58;= -forw&#91;2&#93;;

  glMultMatrixf&#40;@m&#41;;
  glTranslated&#40;-eyex, -eyey, -eyez&#41;;
end;