Bezier surface without using Evaluators OpenGl

Hi all,
This is my first post that contains some actual code. I wont provide the entire code for any of the posts, as it will be a spoon feeding. I will off-course give substantial part of it. (And to much extent with any optimizations, leaving them to the imagination of the user :))

We will use DeCasteljau Algorithm algorithm to draw the surface.The similar code using evaluators is available in Red Book Chapter 13.

Before going onto the code please learn a bit of opengl, as i will only explain how to calculate the points.

GLfloat ctrlpoints[4][4][3];
//set the control points as u wish, 4 sets, 4 points each, x, y, z
int LOD=20;
//increase or decrease this (may be using +/-) to increase 
number of curves making the surface
/* simple linear interpolation between two points
in three dimension*/
void lerp (GLfloat dest[], GLfloat a[], GLfloat b[], float t)
    { 
        dest[0] = a[0] + (b[0]-a[0])*t;
        dest[1] = a[1] + (b[1]-a[1])*t;
        dest[2] = a[2] + (b[2]-a[2])*t;
    }

// evaluate a point on a bezier-curve. t goes from 0 to 1.0    
void bezier (GLfloat dest[],GLfloat a[],GLfloat b[],
             GLfloat c[],GLfloat d[], float t)
    { 
        GLfloat ab[3],bc[3],cd[3],abbc[3],bccd[3];
        lerp (ab, a,b,t);           // point between a and b 
        lerp (bc, b,c,t);           // point between b and c 
        lerp (cd, c,d,t);           // point between c and d 
        lerp (abbc, ab,bc,t);       // point between ab and bc 
        lerp (bccd, bc,cd,t);       // point between bc and cd 
        lerp (dest, abbc,bccd,t);   // point on the bezier-curve 
    }
/*
Given u and v, compute a point on the surface
*/
void eval(GLfloat dest[],float u,float v,GLfloat pnts[4][4][3])
    {
        GLfloat Q[4][3];
        for(int i=0;i<4;i++)
            {
                bezier(Q[i],pnts[i][0],pnts[i][1],pnts[i][2],pnts[i][3],u);
            }
        bezier(dest,Q[0],Q[1],Q[2],Q[3],v);
    }
/*
For orthogonal lines, change we need column vise points
This function find transpose of the matrix
*/
void get_vertical(GLfloat dest[4][4][3])
    {
        for(int i=0;i<4;i++)//ith column
            {
                for(int j=0;j<4;j++)//row
                    {
                        for(int k=0;k<3;k++)//elem
                            {
                                dest[i][j][k]=ctrlpoints[j][i][k];
                            }
                    }
            }
    }


void display()
{
/*
Intialization etc code here
*/

/*
Compute using algo
Plot them using GL_LINE_STRIP to get a smooth line, instead of discrete points
*/
for(int i=0;i<LOD;i++)
            {
                float u = (float)i/(LOD-1);
                glBegin(GL_LINE_STRIP);
                for(int j=0;j<LOD;j++)
                    {
                        float v = (float)j/(LOD-1);
                        GLfloat p[3];
                        eval(p,u,v,ctrlpoints);
                        
                        glVertex3fv(p);
                        
                    }
                glEnd();
            }
            
        GLfloat dest[4][4][3];
        get_vertical(dest);

        for(int i=0;i<LOD;i++)
            {
                float u = (float)i/(LOD-1);
                glBegin(GL_LINE_STRIP);
                for(int j=0;j<LOD;j++)
                    {
                        float v = (float)j/(LOD-1);
                        GLfloat p[3];
                        eval(p,u,v,dest);
                        
                        glVertex3fv(p);
                        
                    }
                glEnd();
            }
/*
Swapping buffer etc code goes here

*/
}